Metamaterial applications of T MATSOLVER , an easy-to-use software for simulating multiple wave scattering in two dimensions

rapid prototyping of designs and for validating other solution methods. We develop a general formulation of the multiple scattering problem that facilitates efﬁcient application of the multipole-based method. The shape and morphology of the scatterers is not restricted, provided their T-matrices are available. The multipole method is implemented in the T MATSOLVER software package, which uses our general formulation and the T-matrix methodology to simulate accurately multiple scattering by complex conﬁgurations with a large number of identical or non-identical scatterers that can have complex shapes and/or morphologies. This article provides a mathematical description of the algorithm and demonstrates application of the software to four contemporary metamaterial problems. It concludes with a brief overview of the object-oriented structure of the T MATSOLVER code.

Multiple scattering of waves is eminent in a wide range of applications and extensive research is being undertaken into multiple scattering by ever more complicated structures, with emphasis on the design of metamaterial structures that manipulate waves in a desired fashion.Ongoing research investigates the design of structures and new solution methods for the governing partial differential equations.There is a pressing need for easy-to-use software that empowers

Introduction
We consider the problem of numerically simulating multiple scattering of waves by configurations containing many particles in two spatial dimensions.Such multiple scattering processes have been studied across a broad range of wave phenomena, such as atmospheric science, oceanography, astrophysics (see, for example, Van de Hulst [1]), thermal radiation [2,3], water waves [4,5] and design of phonic and photonic crystals [6].More recently, there has been significant interest in the design of metamaterials whose wave propagation properties stem from multiple scattering effects between their component atoms or particles.Such metamaterials can have interesting and important qualities, such as controlling the wave speed within the material [7][8][9], or having unusual anisotropy [10], large attenuation [11] or large amplification [12].
In practice, design of such metamaterials using advanced mathematical techniques is supported by numerical simulations for validation before experimental/practical realization (see, e.g.[8,9]).Simulations can be used to explore the differences between physical problems and their corresponding model problems, which are often simplified so they can be solved analytically (for example, hypothetical infinite-metamaterials and their finite physical counterparts).Simulations are also important for validating new solution methods (see, e.g.[13][14][15]).
Most numerical simulations are based on either the finite-element method (FEM) (e.g.[9,11,[13][14][15][16]) or multipole expansions (e.g.[8,10]).The FEM is capable of handling complexshaped particles and heterogeneous media, but requires truncation of the infinite domain and approximation of the corresponding radiation condition.Numerical simulations using multipoles incorporate the radiation condition exactly, but have often been restricted to circular-shaped homogeneous particles, primarily because more complicated geometries require sophisticated theory and numerical implementations.However, recent advances in numerical methods for computing the so-called T-matrix [17][18][19] have facilitated a significant expansion in the kinds of metamaterials that can be simulated using multipoles, including the kind of complex-shaped scatterers typically simulated using the FEM.
The T-matrix describes the response of a single scatterer to any incident wave in the sense that it is the solution operator to the general scattering problem of that body, and generalizes to non-spherical scatterers an approach originally developed in 1881 by Rayleigh [20].The T-matrix was introduced for electromagnetics in 1965 by Waterman [21], and subsequently described for the Helmholtz equation in two and three dimensions in 1969 [22].In these seminal papers, the T-matrix was computed using the Null Field Method.Subsequently, many alternative methods for computing the T-matrix have been developed, including the invariant imbedding method [23], discrete sources expansions [24,25], point matching on the scatterer boundary [26][27][28], combined FEM-BEM methods using Green's functions in eigenfunction expansions [29,30], and projection onto the wave function basis in the far field [17][18][19].The extension of the T-matrix to multiple scattering is based on the self-consistent method, which represents the field incident on each scatterer as the sum of the fields from the other scatterers, together with the exciting field.The use of wave function expansions, and changes of expansion origin using the translation-addition theorem, facilitate application of the T-matrix on each particle to give its scattered field.Similar representations on all of the particles lead to a system of coupled equations that is solved for the scattered field.There is an extensive literature on this approach and we refer to [31] for a detailed review.The self-consistent method is formulated in [31,Section 1.1.3]and we refer to [31,Section 2.5] for the translation-addition theorem.
In this work, we develop a general formulation of the multiple scattering problem, with particular focus on efficiently capturing the key properties of the scatterers and their arrangement in space.This formulation facilitates application of a flexible and efficient multipole-based method for solving the multiple scattering problem, which exploits identical scatterers, if they occur, even if they have different orientations.A key feature of this approach is its ability to simulate accurately scattering by circular and non-circular scatterers with complex morphology using the T-matrix.The numerical method is implemented in an accompanying object-oriented MATLAB package TMATSOLVER [32], which has been written to provide an easy-to-use tool for rapid prototyping and validation of metamaterials comprising from a few to a few-hundred individual scatterers.
Convergence of the multipole expansions is super-algebraic and the number of unknowns required for a particular scatterer is independent of its shape [33].A reliable estimate for the number of unknowns required to represent the scattered field for a single scatterer is given by (2.19) in §2, and is approximately ka, where a is the radius of the scatterer and k the wavenumber.Thus the number of unknowns required for simulating scattering by a configuration of N scatterers is proportional to kN, independent of the arrangement of the scatterers.This is a lower-order dependence on k than for the FEM, for which the number of unknowns required is proportional to kA × kB, where A × B are the dimensions of a box containing the scatterers, which is decomposed into triangular elements (and, typically, perfectly matched layer (PML) conditions are applied on the boundary).This condition is equivalent to requiring a fixed number of points per wavelength in each dimension as the wavenumber increases.We have validated the TMATSOLVER software against MieSolver [34] for multiple scattering configurations containing sound-soft circular cylinders at frequencies up to ka = 40π ≈ 125.7, i.e. cylinders with diameter 40 times the wavelength.
Numerous software packages have been developed for multiple scattering simulations and an excellent directory of open-source software is provided at the Scattport [35] website.Packages that provide easy-to-use object-oriented interfaces, similar to TMATSOLVER, include MieSolver [34,36], CELES [37,38] and MultipleScattering.jl[39].These packages are also based on multipole expansions and the self-consistent method, but are restricted to spherical scatterers (although MieSolver also supports layered-scatterers).As outlined above, extension of the self-consistent method to non-spherical scatterers requires the T-matrix.Although there are several packages available that can compute the T-matrix for three-dimensional scatterers, the only package that can compute the T-matrix for scattering by cylinders is TMATROM [40].
To demonstrate our approach, we present results for four canonical metamaterial problems that were discussed at the 'Mathematical Theory and Applications of Multiple Wave Scattering' (MWS) programme at the Isaac Newton Institute for Mathematical Sciences, Cambridge, UK in January-June 2023 (from which the articles in the Special Feature of Proceedings A were derived).The problems demonstrate TMATSOLVER across a spectrum of challenging problems that include large numbers of scatterers, anisotropic scatterers (and so there is strong mode-coupling even for scatterers with circular shape), scatterers with challenging square boundaries, penetrable noncircular scatterers with high contrast, and scatterers in configurations that exhibit resonance.Code demonstrating the last two problems is included in the TMATSOLVER package.
In more detail, the first problem involves simulation of multiple scattering in an array of anisotropic metacylinders.The incident wave induces interior fields in the metacylinders but Downloaded from https://royalsocietypublishing.org/ on 19 June 2024 these are already taken care of within the T-matrix, and interactions between the scatterers are calculated without explicit consideration of the interior fields, substantially reducing the complexity of the simulation.The second problem involves simulation of multiple scattering by an infinite wedge interface, where multipole simulations using TMATSOLVER for a large number of scatterers are used to validate the semi-analytical solution found using the Wiener-Hopf technique.The third problem considers design of metamaterials with unusually large attenuation.The attenuation is shown to be dependent on the structure of the materials, and arises from multiple scattering effects between high-contrast square particles.The fourth problem involves simulation of the wavefield around a finite periodic line array of square scatterers dominated by Rayleigh-Bloch waves.This effect is exploited in metamaterial applications by tuning the period to modulate the group velocity of the waves and control the amplitude at key locations.We are able to tackle the difficult square-shaped scatterers in these problems using a multipole-based approach because we compute the T-matrix of the scatterers using the TMATROM package, which computes the T-matrix using projection of the far field onto the wave function basis [17][18][19].Although the T-matrix is computed from far-field data, the near-field response of the scatterer is accurately captured because of the relation between the expansion of the near field (given in equation (2.6) below) and the corresponding expansion of the far field [18].Convergence of the expansion (2.6) is guaranteed outside the circumscribing sphere of the scatterer.Convergence inside the circumscribing sphere is discussed in detail in [41].The far fields are computed using the coupled FEM-BEM method [42] and MPSPack [43], respectively.
The structure of the paper is as follows.In §2, we formulate the governing equations for the multiple scattering problem in a flexible way that encompasses all of our canonical metamaterial problems, and concisely describe the self-consistent method of solution.In the remainder of the paper, we focus on our canonical metamaterial problems.In §3a, we simulate multiple scattering by arrays of anisotropic metacylinders.In §3b, we validate Wiener-Hopf simulations from infinite arrays.In §3c, we investigate damping and band gaps in metamaterials comprising high-contrast square scatterers.In §3d, we simulate Rayleigh-Bloch waves in finite line-arrays of square cylinders.In §4, we draw some conclusions.Finally, in the appendix, we briefly describe the TMATSOLVER code and describe-with examples-the key functions and methods of the package.

Mathematical description
All of the canonical problems considered in this paper require solving the two-dimensional Helmholtz equation in an unbounded heterogeneous medium containing one or more scattering particles.Here, k = ω/c is the wavenumber, ω is the angular frequency and c the wave speed in the medium.We assume that the scatterer D can be written where D 1 , . . ., D N are the individual scattering particles.We establish a general description of the two-dimensional multiple scattering problem that encompasses all four of the canonical problems.The interaction of an incident wave u inc with the scatterers D 1 , . . ., D N induces a scattered wave u s in the exterior R 2 \ D. Physically important incident waves include the plane wave with incident direction specified by a unit vector d, and the point source with source located at Here, H n denotes the first-kind Hankel function with order n.where the limit holds uniformly with respect to θ and we use polar coordinates (r, θ) for x.
A consequence of (2.4) is that where the modulation function u ∞ is the far-field of u s .
It is convenient to decompose the scattered field as where u s I can be considered the part of the scattered field radiated by D I .Then u s I admits the series expansion where x I ∈ D I is a local origin for the expansions, B I is the smallest closed ball centred at x I that contains D I .We also define |n| (kr) e inθ and Reg where J n is the Bessel function of order n, and Reg ψ n and ψ n are known as regular and radiating (cylindrical) wavefunctions, respectively.The coefficients a n in (2.6) are unknown and to be determined.
We observe that the field incident on D I is which comprises the parts of the scattered field radiated by the other particles, as well as the incident wave u inc itself.In a neighbourhood of D I the field u inc I admits the series expansion (2.9) The coefficients f n are unknown but, crucially, related to the coefficients a n through the T-matrix [40] of the scatterer D I .In particular, there holds the relation where a (I) = (a n ) and T (I) is a matrix known as the T-matrix (of D I with respect to the local origin x I ).
It often happens that a scattering problem contains multiple copies of the same kind of scatterer, with different positions and perhaps with different orientations, but having the same shape and material properties.In this case, the same T-matrix can be used for all copies of the scatterer, after introducing (temporary) translations and rotations of the coordinate system, and it is efficient to exploit this where possible.Subsequently, we say that particles with the same shape and material properties, but perhaps different positions and orientations, have the same design.It follows that any particle in the configuration can be described by its design, its position and its orientation.(The position and orientation are then described with respect to a reference location and origin of a template for the design.) We do not try to describe the design precisely because the range of particles that can be tackled is very broad, including homogeneous particles of various shapes, but also particles comprising heterogeneous media.It is sufficient to note that design is associated with the particle's morphology, including shape and boundary conditions, and the key from the modelling perspective is not to describe the morphology in detail, but that the scatterer's T-matrix can be computed and that the design can be visualized schematically.
It is expedient to assume that the particles D 1 , . . ., D N can be obtained by translations and rotations of M ≤ N designs D 1 , . . ., D M .Then, the Ith particle can be represented by P I = (D K I , x I , φ I ) for 1 ≤ K I ≤ M, where the design D K I represents a template particle with local origin 0, and P I is obtained by translating the template particle with translation vector x I , and rotating it by φ I in the positive direction (figure 1).The T-matrix of the Ith particle is then [29] where T(D) denotes the T-matrix of the design D, and rotation of the coordinate system by φ in the positive direction is effected by the diagonal matrix R(φ) = (r nm (φ)) and r nm (φ) = δ nm e inφ . (2.12) Matrix vector products with the T-matrix (2.11) are computed in factored form using the right-hand side, and the matrix T (I) is not assembled.The T-matrices of discs, with soundsoft, sound-hard and transmission boundary conditions have simple analytical expressions [34].
A numerically stable code for calculating the T-matrix of a wide range of other particles is provided by TMATROM [40].
Returning to the field incident on D I given by equation (2.8), it is helpful to rewrite the parts of the scattered field radiated by each of the other scatterers as The coefficients in (2.13) are given by Graf's addition theorem [31, Theorem 2.12] (see also [44]) (2.15) Substituting (2.14) into (2.8) and comparing with (2.9) gives where b (I) = (b (I) n ) are the regular wave function expansion coefficients of u inc with local origin x I , which are known analytically for plane waves and point sources (see [34]).Then from (2.10), we see that the unknown coefficients a (J) for J = 1, . . ., N satisfy the block linear system a (I) − T (I) (2.17) The matrix on the left-hand side is a perturbation of the identity matrix and the linear system (2.17) is typically well-conditioned provided the circumscribing spheres B 1 , . . ., B N do not intersect, that is, provided Equation (2.18) is effectively a separation condition on the scatterers.However, the radii of B 1 , . . ., B N can be minimized by carefully choosing the local origins of the reference particles when computing their T-matrices (see [40] and [45] for details).The system (2.17) may be solvable when condition (2.18) is violated, and we refer to the discussion and numerical results in [46] for more details.
It is efficient to solve the linear system (2.17) iteratively using GMRES [47], which requires only matrix vector products with the matrix on the left-hand side, which are performed dynamically without assembling the off-diagonal blocks.Because the linear system (2.17) is usually well conditioned, preconditioning is typically not required.In practice, the infinite series expansions and the associated matrices must be truncated, and we replace ∞ n=−∞ • by n max n=−n max • on each particle, where we use Wiscombe's formula for the truncation parameter [48], and a denotes the radius of the particle.The truncation order n max is allowed to differ between different particles in the configuration but for brevity we do not include this explicitly in our notation.
Once the coefficients a (J) in (2.17) have been computed, they can be used in the truncated expansions (2.6) to compute the near field at points x outside the circumscribing spheres Remark 2.1.The expansions (2.6) may diverge inside the circumscribing spheres and we refer to [41,46,49] for detailed discussion on the convergence properties of the series close to the scatterers.
Finally, from the asymptotic behaviour of the radiating wavefunctions (2.7), we observe that the far field of Here, the unit vector x is the observation direction, with polar coordinates (1, θ), and (2.21) normal incidence.The cross-section of the particles is itself occupied by a metamaterial consisting of an array of perfectly conducting plates and interstitial dielectric channels of relative refractive index n r with respect to the exterior bulk.The microstructured plate-array material within the cylindrical inclusions is replaced by an effective medium description [50][51][52] encapsulated by three equations: a one-dimensional wave equation in the interior and a pair of continuity conditions at the surface.

Results
The geometry of the plate-array medium is illustrated in figure 2a.The parallel plate-array is of decreasing length outwards, such that when viewed from above it forms a cylinder of circular cross-section of radius a, and is aligned at cylinder angle δ with respect to the x-axis where δ ∈ [0, 2π ).In the rotated coordinate system (x , y ) (figure 2a), the constitutive equations are governing the longitudinal magnetic field H i interior to the metacylinder, and which comprise the continuity of the field and its gradient through the cylindrical surface, where H e denotes the exterior field.Equation (3.1) applies strictly to the interstitial dielectric channels, and admits the following expansions in Chebyshev polynomials and regular waves: Following §2, in the case of a single metacylinder, the field exterior to the metacylinder is the superposition of an incident plane wave and the outgoing scattered field from the metacylinder, Substitution of (3.3) and (3.4) into conditions (3.2) at the metacylinder boundary yields the following system of equations for the unknown multipole moments a = (a n ) and Chebyshev coefficients c where g = (c, d) encapsulates the coefficients of the interior field.The H and H sub-matrices comprise the exterior field Hankel functions and their derivatives, whereas the C and C submatrices comprise the interior field terms, and their derivatives.The vectors Q and Q consist of the f n moments and Bessel functions of the first kind and their derivatives.In practice, the infinite series expansions must be truncated, and in this application, it is expedient to replace ∞ n=−∞ • by n max n=−n max −1 •, so that each sub-matrix is of size (2n max + 2) × (2n max + 2).To compute the T-matrix of the metacylinder, we multiply (3.5) by the inverse of the scattering matrix to give where again each sub-matrix is of size (2n max + 2) × (2n max + 2).We write the incident wave vector (Q, Q ) as where the matrices J, J are diagonal matrices of elements J n (ka) and J n (ka), respectively.Combining (3.6) and (3.7), we have The T-matrix is then, using the top row, We use the submatrix corresponding to n = −n max , . . ., n max in our subsequent calculations.
The T-matrix approach represents a rapid technique for the investigation of the eigenmodes of large, but finite, arrangements of scatterers.A prime example of this type of system is the quasi-periodic array, in which a long chain of periodically arranged inclusions undergoes some perturbation to the positions of the scatterers.In figure 2c and index.The provided solution converged exactly to that computed with the multipole method detailed in [10] (not shown), but required less than half of the computational time, due to the large number of expansion terms necessary to the multipole technique.Thus, we envisage the use of T-matrices as a faster method for the exploration of the parameter space of systems consisting of many hundreds, or even thousands, of inclusions.

(b) Application: infinite wedge interface of point scatterers
In this section, we consider scattering by line arrays with large numbers of scatterers, which provide finite approximations to infinite line arrays.Of particular, interest is the configuration with two semi-infinite arrays positioned to form a wedge interface.The scatterers are illuminated by a plane wave with direction d = (cos(π + θ I ), sin(π + θ I )), where θ I denotes the angle from which the plane wave is incident.The scatterers are identical sound-soft circular particles with small radius a and positions To prevent the scatterers from overlapping, we require a < s/2 and a < s| sin(α 1 + α 2 )/2|.When α 1 = 0 and α 2 = π the configuration becomes an infinite line array, for which an analytical solution is available in the low-frequency case.
The analytical Wiener-Hopf-based techniques in [13,16] can be applied to the semi-infinite arrays under the assumption that the scatterers are small in comparison to the wavelength of the incident wave, that is ka 1.On the other hand, truncating the arrays leads to a scattering problem with a finite number of scatterers that can be solved directly using TMATSOLVER, which provides a useful comparison with the analytical techniques for validation.
We briefly review the method in [13,16].For sound-soft boundary conditions, we can model the scatterers as isotropic point scatterers in the spirit of [53] (known as Foldy's approximation).Using a similar expansion to the one in (2.6), but using only monopole terms (i.e.only the ψ 0 (r, θ) terms) the approximation to the total scattered field is where A (1) J and A (2) J are the unknown scattering coefficients to be determined.The procedure of [13,16] is to find two coupled systems of equations, one for each semi-infinite array, and solve them using the discrete analogue of the Wiener-Hopf technique, in a manner similar to the semiinfinite array problem [54,55].The Wiener-Hopf solution is written in the form of a block matrix equation, which is then inverted to determine the scattering coefficients contained in the column vectors A (l) .Here, A (l) 0 is a column vector containing the scattering coefficients that solve the equivalent independent semi-infinite array problems, the entries of which are given by eqn (25) in [16].The entries of the infinite block matrices M (l ,l) are given by eqn (29) in [16].
We will directly compare the scattered wave fields u s computed by the Wiener-Hopf method and TMATSOLVER for two test cases.For the Wiener-Hopf method, we truncate the matrix in (3.12) so that each block is an M × M matrix.For the TMATSOLVER method, we restrict the number of multipole terms to just the monopole and dipole terms (i.e.truncation parameter n max = 1), which is appropriate for the small radius a = 0.001 of the scatterers considered.The number of scatterers in the finite arrays is chosen to match the number of coefficients evaluated through the Wiener-Hopf method.The remaining parameters are as given in fig.6 of [13].In figure 3, we demonstrate the excellent visual agreement in the scattered fields computed using both methods.In figure 4, we examine in detail the discrepancy between the two methods by plotting the difference in the monopole scattering coefficients at the different scatterers for a wedge and an infinite line array (with α 1 = 0 and α 2 = π ), respectively.We order the coefficients such that A −101 , . . .
100 respectively.We note the similarities between the errors in figure 4 for the two configurations.
The advantage of considering the infinite line array is that there is an analytical solution for the scattering coefficients that is independent of the truncation, K(e iks cos(θ I ) ) , J ≥ 0 and A (2) e iks(J+1) cos(θ I ) K(e iks cos(θ I ) ) , J ≥ 0.
Because this exact solution is available, we include in figure 4c,d comparisons with (3.13).This allows us to observe the errors in each method for simulating the infinite arrays.The similarity of these plots with fig.6 in [16] implies similar conclusions.Firstly, that the Wiener-Hopf method is weaker in the middle of the infinite line array (J ≈ 0), where the ends of the two semi-infinite arrays meet, and stronger where the arrays are truncated (J ≈ ±100).Conversely, TMATSOLVER is stronger in the middle and weaker where the arrays are truncated.The reason for this is how each method models the geometry of the problem.The Wiener-Hopf method considers each semi-infinite array individually and adds the interaction between them when the big matrix is inverted.The TMATSOLVER method considers each scatterer individually and solves for the individual interactions but does not model the infinite number of scatterers.
(c) Application: high-contrast heterogeneous media In this section, we consider wave propagation through high-contrast metamaterials comprising a large number of small, heterogeneous, penetrable scatterers.The scatterers are square particles centred at positions x ∈ εZ 2 and with side length ε/2 for the small parameter 0 < ε < 1.We chose ε = 1/8 in the experiments below.Metamaterials can induce astonishing effective wave phenomena such as exponential damping and, related, the occurrence of all-angle frequency band gaps.The latter means that, independent of the angle of an incoming plane wave, the metamaterial dampens the waves at certain frequencies so that for such frequencies no waves can propagate through the material.For this, a high material contrast is crucial to incite special resonance effects, see [56,57].The high material contrast means that the refractive index is chosen as ε 2 inside the small particles.
We consider a plane incident wave with wavenumber k = 28, which is known to induce resonance effects in this configuration [11].The T-matrix of a single small square particle is computed using TMATROM [40] from far-field data computed using a coupled FEM-BEM solver [42].Having fixed k and ε and having computed the T-matrix, TMATSOLVER allows us to calculate the total field for various placements of the small particles quickly.We use this computational efficiency of TMATSOLVER to investigate the following pre-requisites for exponential damping of the incident plane wave: (i) the necessary thickness of the metamaterial slab and (ii) the importance of a periodic arrangement of the small particles.Code for this application is included in the TMATSOLVER package.
We position our particles at x lm := ((2.5 + l)/8, (0.5 + m)/8) for m = 0, . . ., 7 and l = 0, . . ., n and vary the number of layers n + 1.The computed total field for one, two, three and four layers is visualized in figure 5.A certain number of layers is required to dampen the wave significantly, as we can observe, for instance, by comparing the results for one and three layers.Furthermore, the effect of adding more layers decreases rapidly with the number of layers already present.For our setup, two to three layers seem to be sufficient to get close to maximum damping.
Similar to the previous experiment, we observe the desired wave damping for the particles positioned at x lm := ((2.5 + l)/8, (2.5 + m)/8) for l = 0, . . ., 3 and m = 0, . . ., 3, see figure 6a.Since the incident wave can propagate around the metamaterial in this example, the damping expresses itself in a 'shadow region' directly behind (i.e. to the right of) the material.We compare this periodic configuration with three different disordered configurations for the particles.For the latter, we randomly choose 16 centres x I out of the 32 possibilities x lm := ((2.5 + l)/8, (0.5 + m)/8), m = 0, . . ., 7 and l = 0, . . ., 3. The periodic configuration shown in figure 6 is one of 32 16 = 601080390 possible particle geometries.As figure 6 illustrates, the periodic configuration is clearly distinguished with respect to its damping potential in comparison to the disordered configurations.The latter show a wide variation of wave behaviour due to the different scatterer geometries, but no clear shadow region can be observed for any of them.
(d) Application: Rayleigh-Bloch waves In this section, we consider simulation of wavefields dominated by Rayleigh-Bloch waves in finite periodic line arrays of scatterers.In particular, we consider arrays comprising N identical, equally spaced sound-hard obstacles, D 1 , . . ., D N , which are aligned along the x-axis with uniform spacing.Rayleigh-Bloch waves propagate along the corresponding infinite array with wavenumber β(k) > k, and decay exponentially away from it [58].Nevertheless, Rayleigh-Bloch waves can dominate the response along the finite array [59][60][61].They can also be used to engineer desired responses along the array [62,63].
Rayleigh-Bloch waves are unforced solutions of the corresponding infinite-array problem (i.e. with an infinite number of scatterers).They are a class of Bloch wave, familiar in the analysis of doubly periodic structures such as photonic/phononic crystals, and they exist quite generally along infinite array-like structures [64].They are the localized solutions obtained in the limit that the unit cell for a doubly periodic structure tends to infinity in one dimension [65], and, thus, they are also a class of trapped mode [58].Plane waves cannot excite Rayleigh-Bloch waves along an infinite array as β(k) > k, but they can excite Rayleigh-Bloch waves along a semi-infinite array that propagate away from the end [66,67], or in both directions along a finite array [60].
In the acoustic/water-wave setting, Rayleigh-Bloch waves have been computed predominantly for circular scatterers [59][60][61][65][66][67][68][69], with notable exceptions being elliptical scatterers [58], rectangular scatterers [70] and C-shaped resonators [62,71].For circular cylinders of radius a and centre-to-centre spacing R, Rayleigh-Bloch waves that are symmetric about the x-axis exist for all k ≤ k c < π/R, where k c (a) is known as the cut-off frequency.They have been extended to frequencies above the cut-off, for which the wavenumber β becomes complex and, hence, the Rayleigh-Bloch waves attenuate along the array [61].They have also been extended to multi-line arrays [71].
As mentioned above, calculation of Rayleigh-Bloch waves has mostly been restricted to circular obstacles, primarily because more complicated geometries require sophisticated numerical implementations.Moreover, on finite (or semi-infinite) arrays, they are usually forced by simple plane waves.By contrast, we use TMATSOLVER to observe Rayleigh  scatterers are computed using TMATROM [40] using far-field data computed using MPSPack [43].
Code for this application is included in the TMATSOLVER package.The presence of Rayleigh-Bloch waves is detected by performing a frequency sweep of the response of the array at its centre point x = x * = ((N + 1) R/2, 0).There is a near-resonant peak in the response for k < Rπ (figure 7a), which is caused by constructive interference between Rayleigh-Bloch waves after reflections and re-reflections by the array ends [60].The resonance occurs at a frequency just below the cut-off, which is indicated by the sudden drop in the response at higher frequencies, and can be confirmed by calculations of the spectrum for the infinite array, similar to [68] (not part of TMATSOLVER).The profile at the near-resonant frequency (figure 7b) is symmetrical with a peak at the centre, which is associated with primary resonances caused by Rayleigh-Bloch waves [59][60][61]72].The wavefield at the near-resonant frequency is dominated by the Rayleigh-Bloch waves, particularly around the centre of the array (figure 7c).At this frequency just below the cut-off, the Rayleigh-Bloch waves are approximately anti-symmetric with respect to the vertical axes passing through the scatterers, with blue on one side and yellow on the other side of a given scatterer (compare with fig.3c in [61]).Although the near-resonant field contains both rightward and leftward propagating Rayleigh-Bloch waves, they are close to standing waves (near-zero group velocity) and have similar shapes.

Conclusion
We have addressed the challenge of rapidly simulating multiple scattering processes in metamaterials by developing a general formulation of the multiple scattering problem that can simply and easily describe typical metamaterial structures, while facilitating efficient implementation of the multipole-based self-consistent method.Complex shapes and morphologies of the scatterers are encapsulated in their T-matrices, which allows our formulation to model metamaterials comprising complex particles.Our formulation is implemented in Downloaded from https://royalsocietypublishing.org/ on 19 June 2024 the TMATSOLVER software, which provides a tool for researchers working on metamaterials to prototype their metamaterial designs or validate numerical methods quickly and easily.Numerical results have demonstrated the application of the software for challenging problems arising in recent metamaterials research, which include large numbers of scatterers, anisotropy, complex-shaped boundaries, high-contrast materials and resonances.In two dimensions, the particle's position vector (x, y) can be represented by the complex number z = x + yi and its orientation can be represented by a real number.The particle class represents a particle by associating a design with a position and orientation.For example, the following code creates a particle with position (2, 3) and orientation π/4: z = 2+3i; % position phi = pi/4; % orientation p = particle(d, z, phi); The particle can be visualized using p.plot( ) This uses the design.plotmethod.The user may define their own child classes (of the design class) to implement their own scatterers.This has the particular advantage over the extern_design class in allowing overloading of the plot method to specify how their scatterer should be plotted.

(b) The solver
Solving the scattering problem and visualization of the solution is performed using the tmatsolver class.As outlined in §2, the scattering problem involves an incident wave u inc and a collection of particles.In TMATSOLVER, the incident wave is represented by the TMATROM incident class and its subclasses.For example, the plane wave u inc (x) = e ikx•(cos θ,sin θ) with direction θ = π/4 and wavenumber k = 1 is represented as follows: theta = pi/4; % incident wave direction kwave = 1; % wavenumber uinc = plane_wave(theta, kwave); We refer to the TMATROM documentation [45] for more information about incident waves that are linear combinations of plane waves and point sources.
The solver can be set up using the incident wave and the particles, for example, obj = tmatsolver(uinc, p1, p2, p3, p4) ; or using just the incident wave, with the particles added subsequently The scattering problem is solved using obj.solve( ) The linear system (2.17) is solved iteratively using GMRES [47].The default solver tolerance is 10 −8 and the default maximum number of iterations is the dimension of (2.17).The solver Downloaded from https://royalsocietypublishing.org/ on 19 June 2024 tolerance and maximum number of iterations can be set to other values using, for example obj.setSolverTol(1e-4); obj.setSolverIterations(100); We refer to the online documentation (type help tmatsolver in the MATLAB command window) for methods to check the convergence of the GMRES iteration.
The total field, computed from the solution of the scattering problem, can be plotted using, for example We refer to the online documentation (type help tmatsolver in the MATLAB command window) for the full range of plotting and visualization options.

Figure 1 .
Figure 1.Schematic showing a template particle D at 0 and the associated scatterer P I with origin x I and rotation φ I .

. 6 )Figure 2 .
Figure 2. (a) Coordinate systems of the metacylinder.Cylinder coordinates are rotated through δ to primed plate array coordinates.The plate array is merely illustrative and is not representative of the number of plates, the plate thickness or the inter-plate spacing.(b) Normalized longitudinal magnetic fields (0, 0, H z ) due to the scattering by a plane-wave of incident angle ψ inc = π/2 by a single metacylinder of radius a = 1 and relative index n r = 1.(c) Scattering of a plane-wave of incident angle ψ inc = π/2 by a quasi-periodic chain of metacylinders of uniform angle δ = −π/4.(a) Coordinate systems of the metacylinder, (b) Scattering for δ = π/4, n r ka = 1 and (c) Scattering by a quasi-periodic chain.

Figure 4 .
Figure 4. Absolute difference between the scattering coefficients for a wedge array (a,b) and an infinite line array (c,d).(a) Absolute error of A J for point scatterer wedge (θ I = π/4), (b) absolute error of A J for point scatterer wedge (θ I = π/12), (c) absolute error of A J for infinite array (θ I = π/4), (d) absolute error of A J for infinite array (θ I = π/12).

Figure 5 .Figure 6 .
Figure 5. Visualization of the total field Re{u(x)} from a periodic high-contrast material with varying thickness.The field is excited by a plane wave incoming from the left.The field inside the circumscribing circles is not shown (see remark 2.1).

Figure 7 .
Figure 7. Detection of Rayleigh-Bloch waves along a line array of N = 21 sound-hard square scatterers forced by a point source.(a) Response at the mid-point of the array versus frequency, where the near-resonant response is indicated by the red circle.(b) Profile at the near-resonant frequency, at the mid-points between the scatterers.(c) The wave field, Re(u), at the nearresonant frequency.

Figure 8 .Figure 9 .
Figure 8. UML class diagram showing the structure relating the solver to the particles and incident wave.Key class attributes and methods are shown but others may be omitted.