Metaclusters for the Full Control of Mechanical Waves

We present a new method for the control of waves based on inverse multiple scattering theory. Conceived as a generalization of the concept of metagrating, we call metaclusters to a finite set of scatterers whose position and properties are obtained by inverse design once we have defined their response to some external incident field. The method is applied to the propagation of flexural waves in thin plates, and to the design of far field patterns, although its generalization to acoustic or electromagnetic waves is straightforward. Numerical examples are presented to the design of uni and bidirectional ``anomalous scatterers'', which will bend the scattering energy along a specific direction, ``odd pole'' scatterers, whose radiation pattern presents an odd number of poles and to the generation of vortical patterns. Finally, some considerations about the optimal design of these metaclusters are discussed.


I. INTRODUCTION
Active and passive control of the energy transfer in electromagnetic and mechanical waves is a challenging problem with a large number of applications, such as focusing, imaging, beam forming, cloaking and energy harvesting, among others. 1 The advent of so-called metamaterials 2,3 provided a new perspective since these artificial structures allow the design of materials with extraordinary properties capable of manipulating the flow of energy in ways that would be impossible with common materials, enlarging in this manner the number of devices for the control of electromagnetic and mechanical waves.
More recently, the concept of "metasurface", conceived as artificial planar metamaterials, has attracted an increasing interest. Being thinner and less dissipative than bulk metamaterials, these structures allow for more efficient ways of manipulating the wave energy, with the additional simplification in fabrication that planar structures present in comparison with bulk structures. [4][5][6] However, the major drawback of both metamaterials and metasurfaces is that their functionality is based on the extraordinary refractive/reflective properties they present, and most of the devices designed in this framework require a large number of scattering elements in order to form an "effective" material whose effective physical properties provide metamaterials of their extraordinary properties. In the case of metasurfaces, the surface has to be gradually structured so that the effective gradient in the surface impedance allows for the manipulation of the energy flow. This large number of scattering elements is an important limitation in the efficiency of metamaterials and metasurfaces, since in practice the number of different scattering elements will be limited, especially in the micro or nano-scale.
To overcome these difficulties, several approaches have been explored recently to simplify the design of metasurfaces by means of diffraction gratings, [7][8][9][10][11][12] in which it has been possible to find a complex scatterer or unit cell performing the same functionality as some metasurfaces. However, the design process is still complex and functionality is limited to the control of the propagation direction of waves. [13][14][15] In this work, we present a generalization of the concept of a metagrating but for finite structures. The objective is to show how, for a given incident field, we can obtain a cluster of scatterers and their physical properties such that the scattered field presents a pre-selected shape. If a particular diffraction pattern is desired for a specific type of incident wave, we provide a method to design a cluster of scatterers capable of transferring the energy along the desired directions. The inverse design method presented is based on multiple scattering theory 16 and the general principle is applicable to any kind of classical wave, including acoustic and electromagnetic waves. This work therefore provides a general principle for the full control of mechanical and electromagnetic waves based on scattering elements.
The paper is organized as follows: After this introduction, section II develops the idea of the direct and inverse multiple scattering problem. Section III explains how the method can be applied to the design of far field patterns and section IV shows numerical examples of specific patterns. Finally, section V summarizes the work and some mathematical results are given in Appendices A and B.

II. DIRECT AND INVERSE MULTIPLE SCATTERING PROBLEM
When some incident field ψ 0 impinges on a cluster of N point-like scatterers the total field ψ(r) can be expressed as the sum of the incident plus the scattered fields, ψ(r) = ψ 0 (r) + ψ s (r).
(1) arXiv:2009.13376v1 [physics.app-ph] 28 Sep 2020 The scattered field is where G(r) = G(|r|) is the Green's function and the coefficients B β are obtained from the multiple scattering This provides a system of N equations with N unknowns. The quantity t α is the strength of each point-like scatterer and it is the only quantity that contains information about its physical properties. This describes the direct multiple scattering problem, in which the number of scatterers, N , their strengths t α and locations R β are known, from which we compute the B α coefficients to finally determine the field in all of space. The inverse problem is as follows: assume that the scattered field can be expressed as a linear combination of basis functions φ n such that then we specify the inverse problem as determining a finite number N p of A n coefficients, for n = 1 . . . N p , so that the scattered field will have a specified radiation pattern in the far-field. In general there will be a matrix S such that therefore, if we select the number N of particles in the cluster equal to the number N p of modes to design, equation (5) constitutes a determinate system of N equations with N unknowns from which we can solve for the B β coefficients. Once these are known, we can obtain the t α elements from equation (4) as Thus we can obtain the physical properties of each particle. The main challenge is to find a cluster configuration giving physically acceptable particles. For the case of flexural waves on thin elastic plates, ψ is the plate deflection, G is the solution for a point force per unit area applied in the positive ψ-direction, and is the point force per unit area of scatterer α, see Appendix A. The parameter t α is an effective point impedance which can be interpreted in terms of a single degree of freedom system with mass, stiffness and damping. Physically acceptable particles cannot supply energy, i.e. they must be passive. Assuming time dependence e −iωt , the passivity constraints require that one or other of the following is met Equation (8a) requires that the cluster be globally passive, meaning that some of the scatterers can provide energy but there should be a negative energy balance adding all the contributions of the scatterers. Equation (8b), or equivalently Im t α ≥ 0, is a more restrictive condition, since it requires that all scatterers be passive systems (see Appendix A for details). The equality holds for zero dissipation in both equations. The goal of the inverse multiple scattering problem is to obtain a set of particles all simultaneously satisfying the first or both constraints. For the first, global passivity, constraint we assume that although some scatterers may require energy supply, this energy can be transferred from other, locally passive, ones (see Appendix A). The specific problem addressed below is to engineer the cluster of point scatterers to provide a close approximation to a desired far field scattering response. In the next section we outline the steps necessary to achieve this in an optimal sense.

A. Direct far field solution
The functions φ n of (4) are chosen as the infinite set φ n (r) = G(r)e inθ , n ∈ Z, where the position is expressed in polar coordinates r = (r, θ) with respect to an origin at r = 0. This allows to uniquely identify the coefficients A n of (4) as far field amplitudes of the scattered wave. In order to see this, first note that the far-field for a source at R β = (R β , θ β ) is This approximation holds whether the Green's function is for the Helmholtz equation or for the Kirchhoff plate equation. In both cases, the far-field response depends only on the large argument approximation of H 0 (x). The scattered far-field of the cluster follows from (2) and (10) as with the far-field radiation function Alternatively, where the infinite set of coefficients {A n } is related to the N coefficients {B β } by (5) with For a unit amplitude incident plane wave propagating in the direction θ = 0 the radiation pattern function satisfies the optical theorem 18 where the scattering cross-section σ sca and absorption cross-section σ abs are defined in Eq. (A6). Further details can be found in Appendix A. The cross-sections can also be expressed directly in terms of the coefficients {A n } and {B β }, see Eq. (A8), leading to the explicit form of the optical theorem

B. Inverse problem
In the inverse source problem we are given f (θ) and seek the cluster that optimally reproduces this scattering pattern. The radiation pattern, defined by the coefficients {A n } in the form (13), is infinite dimensional, whereas the cluster comprises a finite set of N sources. We define the error where ||X|| 2 = X † X with X † the Hermitian transpose of vector X. Minimizing E for given A and S yields the solution where S † S −1 S † may be identified as the Moore-Penrose inverse of S. The approximated radiation pattern is f (N ) (θ) where A (N ) n , n ∈ Z, are the elements of and the non-negative definite Hermitian matrix P is It is shown in Appendix B that the matrix P is infinite dimensional but finite rank with N non-zero eigenvalues equal to +1, see Eq. (B5). It therefore acts as a projection from the infinite dimensional space of far-field pattern functions to the N −dimensional set of approximate pattern functions: The optimal solution (20) yields an error where In practice we will be interested in the relative error C. Inverse design algorithm Based on the above findings, the inverse scattering design can be formulated as follows.
2. The desired far field pattern f (θ) is specified, or equivalently the set of far field modal amplitudes {A n , n ∈ Z} are given (see (13)).
5. The source strengths B α , the optimal approximation to the far field pattern n , n ∈ Z}, and the relative error E rel are calculated (see (18), (14), (24)).
The first two items are geometrical, independent of frequency and the incident wave. Once the frequency is defined, the approximation f (N ) (θ) to the scattered far field is optimal in the sense of an N -dimensional solution according to the setup, and it is independent of the incident field. The form of the incident wave, combined with the source amplitudes B α , defines the required particle impedances t α in Eq. (6).
The question that must be addressed is whether or not all of the scatterer impedances satisfy the passivity constraints (8a) or (8b).

IV. APPLICATIONS
A. Far field patterns and the matrix P Two groups of cluster patterns are considered, namely regular polygons, where scatterers are uniformly distributed over a circle, and finite lattices, where scatterers are regularly distributed in a 2D finite grid. We describe how different arrangements of the scatterers influence the matrix P of Eq. (21) which defines the optimal approximation to the desired scattering pattern.

Scatterers on a regular polygon
Let us assume that the N scatterers lie on the circle of radius R at θ β = 2πβ/N . We consider A corresponding to each of the modes e imθ , m ∈ Z, so that ||A|| 2 = 1 and E ≤ 1 with E 1 indicating the desired scattering mode is well approximated. The results of numerical experimentation are as follows. For small kR relative to N , E is small for modes m = 0, ±1, . . . , N −1 The accuracy diminishes as kR increases. In other words, for small kR the N unit eigenvalues of P correspond to modes m = 0, ±1, . . . , N −1 2 if N > 1 is odd, with analogous association for N even. Since the modes are multiply degenerate (all of eigenvalue unity) it follows that any linear combination of these modes is an eigenvector.

Scatterers on a finite square lattice
We now assume that the scatterers are distributed in a square but finite lattice. The lattice is M × M ≡ N with lattice spacing a. For instance, with M = 3 and ka = 1 we find that the 9 eigenvectors of P of eigenvalue +1 span the space Ω 4 ≡ {e imθ , m = 0, ±1, . . . , ±4}. This result is arrived at by inspecting the error E for each mode, and noting that it is small, on the order of 1e − 4 typically, while higher modes have error of approximately unity.

General properties of the P matrix
Numerical experiments on matrix P for different spatial configurations of the clusters show that for large and moderate kR there are exactly N eigenvalues of P with values close to 1. For large kR, the corresponding eigenvectors (i.e. patterns of scattering modes of the cluster) are highly irregular and sensitive to both kR and scatterers positions (while the number of eigenvalues of value 1 equals N ). For kR c ≈ 0.5 and smaller, where R c is the characteristic size of the cluster, the eigenvalues of P begin to differ and assume values other than 1. For kR c 1 the number of nonzero eigenvalues reduces and low order scattering patterns are preferred.
Some general remarks on the number of scatterers (N ) and scattering properties of the cluster can be formulated as follows. The larger N , the larger number of cluster modes, thus more complex scattering patterns can be reproduced accurately. Large number of scatterers in the cluster, on the other hand, may result in overconstraining the minimization problem and lack of locally or globally passive solutions. For moderate kR c typically the number of regular patterns (eigenvectors of P) is similar to N , while more degenerate patterns and/or smaller number of similar eigenvalues (≈ 1) are observed for large or small kR c .

B. Scattering patterns
The inverse design of metaclusters is illustrated with the scatterers arranged on regular polygons or square lattices, as outlined in Sec. IV A. Here we present the target scattering patterns that will be later reproduced by proper selection of passive impedances.

Odd-pole patterns
Odd-pole patterns havep scattering lobes directing energy towards that preferential directions. The odd-pole scattering pattern and the corresponding A n coefficients are given as where θ = θ mod 2π is used to ensure that f (θ) is a 2π−periodic function. a. A tripole For a tripolar pattern,p = 3, there are three main lobes spaced every 2/3π. An example of a tripolar pattern is shown in Fig. 1c. b. A pentapole Similarly, forp = 5, a pentapole scattering pattern is obtained. Figure 1d illustrates this type of pattern.

A vortex
A vortex generates uniform constant amplitude pattern with angle-dependent linearly changing phase behavior. The corresponding formulas for the vortex of orderp ∈ Z are Directional characteristics of amplitudes and phases for the vortex pattern are shown in Fig. 1e.

C. Full metacluster designs
Designing a metacluster requires finding all t α for a given cluster topology and the desired scattering pattern. The procedure outlined in Sec. III C is employed here to find t α . We first present metacluster scattering patterns corresponding to the desired patterns from Sec. IV B, obtained for different clusters configurations. Since the inverse procedure frequently leads to active particles, we next impose the condition (8b) to find locally passive optimal metaclusters and present their scattering responses. For all presented examples we introduce the incident wave -without loss of generality -assumed to be a plane wave in the −x direction (θ = π).

Scattering patterns for optimal metaclusters
Scattering patterns obtained for selected cluster topologies are shown in Fig. 2. Very good agreement between the desired patterns of Fig. 1 can be seen, proving the effectiveness of the design procedure. However, some of the corresponding impedances -computed using the inverse approach of Sec. III C -are active, hence require energy supply. We next analyze and adopt the inverse procedure for seeking only locally passive solutions.

An optimization problem for passive metaclusters
Our design objective is the set of point impedances {t α , α = 1, . . . , N }. We aim at fulfilling the local passivity condition, Eq. (8b). Define then Eq. (6) becomes u α = p α + s α , α = 1, . . . , N . Consider plane wave incidence ψ 0 (r) = p 0 e ik·r , for some wavenumber k. There is a further degree of freedom that  ) and (20). Incidence angle θ = π (the −x direction). For the vortex, Fig. 2e, circular shapes are amplitude profiles while the spirals in the center are phases of the scattering pattern. All patterns are normalized.
has not been used. This could be considered as the amplitude and phase of the incident wave, i.e. the complex number p 0 . Alternatively, if we fix p 0 = 1, then there is a similar degree of freedom in how we normalize the far field pattern function f (θ). This has the effect of scaling A and hence B by a complex number. This scaling redefines p α but has no effect on s α of (29). Therefore, with no loss in generality we assume the incident wave has unit amplitude, and rewrite Eq. (6) as where the complex number z defines the scaling of the pattern function, which goes as z −1 . The important point is that z can be chosen arbitrarily; in particular, we use it as an optimization parameter. The fact that the pattern function amplitude is inversely proportional to |z| for unit amplitude incident wave suggests that smaller |z| is preferred. The optimization problem is as follows: given the N complex numbers p α associated with the incident wave and the N complex numbers s α associated with the point sources, find z that ensures Im u α ≤ 0 for all α. If this can be achieved then the optimal solution is the one with minimum value of |z|, ensuring maximum amplitude for the pattern function. It might not be possible using the single complex number z to obtain all of the complex numbers u α in the negative imaginary half-plane. If this is not achievable in practical examples then the constraint may be relaxed, for instance, to minimize the maximum instance of positive Im u α . Then the "nearest" passive configuration can be identified by setting Im u α to zero for those particles with positive Im u α . Another alternative could be based on condition (8a), i.e. when the metacluster is globally passive -meaning that the net energy supplied to the cluster non-positive.
In cases where the search procedure for t α failed to find locally passive metaclusters, a rigid rotation was applied to the cluster (equivalent to changing the incidence angle) and the search was repeated.
Numerical experimentation shows there are metacluster configurations for which the inverse impedance solutions are all passive. Examples of the uni-and bidirectional scattering patterns for a square lattice metacluster are shown in Figure 2. More detailed investigations show that -for instance -a square array with lattice parameter a designed to direct a wave incident from the θ = π direction into a scattered wave preferentially directed toward θ = 3/4π has totally passive solutions for 1.9 ≤ ka ≤ 2.8. The optimal passive admittances t −1 α are frequency dependent, with values at the end of the passive interval shown in Table I. The associated optimal scattering patterns are shown in Figure 3. In all examples we take a = 1 and D = 1.
The examples in Figure 3 and Table I are based on the value of z in (31) for which the largest value of Im t −1 α is zero. This optimizes the passive array in terms of its efficiency in converting the incident energy into a directed far field pattern. The metacluster dissipates wave energy but in a way that is most efficient among all passive options.
Similarly, Fig. 4 shows a passive optimal 2 × 2 metacluster that converts the incident plane wave into a symmetric bi-directional pattern. Although the metacluster scattering pattern roughly approximates the desired far field function, all admittances are purely real indicating no dissipation in the system.  Table I.  α for the 2 × 2 array of passive particles sending the wave incident at θ = π into the θ = 3/4π direction at two frequencies, see Figure 3.
Further experimentation shows that obtained optimal solutions are very sensitive to the scatterers' positions and impedances. Also, requirements of symmetric clusters are overconstrained, most often resulting in at least one active particle, especially for large number of particles N .  The admittances t −1 α for the 2 × 2 array of passive particles transforming the incident wave into a bi-directional symmetric pattern at θ = 3/4π, see Figure 4. Admittances of the cluster are give in Table II.

Example: A passive optimal metacluster for odd-polar patterns
Analogously to the previous search, we look for optimal passive clusters capable of generating a scattering tripole. Figure 5 shows the target and the actual scattering patterns for the tripole obtained for a square 2 × 2 cluster of scatterers. The optimal positions and admittances of the scatterers are shown in Table III. The admittances have nearly the same passive damping properties. The corresponding displacement field pattern generated by the metacluster is shown in Fig. 5.  Table III. Figure 6 shows a metacluster designed for generating a pentapole pattern. The cluster consists of a circular arrangement of five scatterers with optimal positions and impedances listed in Table IV. Clearly, the cluster is locally passive. It is important to note that this metacluster setup, resulting in nearly perfect pentapole, has been obtained accidentally when looking for other scattering patterns (different that the pentapole pattern). The lat-  The admittances t −1 α for the 2 × 2 array of passive particles sending the wave incident at θ = π for ka = 5.4 into the tripole pattern, see Figure 5.  α for the five-element circular array of passive particles sending the wave incident at θ = π for ka = 5.3 into the pentapole pattern, see Figure 6.
ter is a consequence of relaxing the requirement of enforcing the target phase of the scattered field and indicates that much more complex scattering patterns that are still locally passive may be obtained for desired amplitudeonly rather than amplitude-and-phase target fields.  Table IV.

Example: A passive optimal metacluster for a vortex pattern
Finally, we present a locally passive metacluster capable of transforming the incident wavefield into the first-  α for the2 × 2 square array of passive particles sending the wave incident at θ = π for ka = 5.5 into the vortex pattern, see Figure 7.
order vortex,p = 1, as shown in Fig. 7. It can be seen from Fig. 7a that despite the amplitude pattern is not perfectly preserved, the phase behavior (Fig. 7b) recovers the linearly-dependent angular characteristic of the vortex. Figures 7c and 7d show displacements and phases of the wavefields generated by the metacluster.
Admittances of the cluster are give in Table V.

V. SUMMARY
In summary, we have shown that an inverse multiple scattering method can be applied for the design of the radiation patterns of clusters of scatterers. While the design process is complex and passive solutions are not easy to find, the approach has still more degrees of freedom to explore. The specific case of flexural waves in thin elastic plates has been considered, with the target problem of designing the far field patterns, although it is easy to show that near field patterns can also be considered. Similarly, the multiple scattering formulation is not unique of flexural waves, for which this approach can be easily extended to other classical waves, like optical or acoustical.
The proposed metaclusters can be considered a generalization of the notion of metagrating, where the inverse design is performed in the amplitude of the diffraction orders and the structures are infinite periodic gratings. With the alternative presented here we could design not only finite gratings but also flat lenses, beam splitters and even cloaking devices. We consider therefore that this work opens a new direction towards the design of passive devices for the control of mechanic and electromagnetic energy. The plate has thickness h, bending stiffness D (= EI/(1−ν 2 )) and density ρ. Time harmonic motion e −iωt is assumed, so that the flexural wavenumber k is defined by k 4 = ω 2 ρh/D. We assume there are N point scatterers located at R α with impedances t α , α = 1, 2, . . . , N . The total displacement ψ satisfies A generic attachment impedance t may be modeled as single degree of freedom with mass M , spring stiffness κ and damping coefficient ν. Two possible configurations are In case (a) the mass is attached to the plate by a spring and damper in parallel. Model (b) assumes the mass is rigidly attached to the plate, and both are attached to a rigid foundation by the spring and damper in parallel. 19 An important limit is a plate pinned at R α , ψ(R α ) = 0, which corresponds to t → ∞. The (a) and (b) oscillators could also be attached in parallel. e.g. on either side of the plate, to give t = t a + t b . The Green's function is the particular solution ψ = G for a single source of the form δ(r) on the right hand side of (A1), The following identity may be found starting from the plate equation (A1) using the procedure of Norris and Vemula 18 for the analogous case without source terms, Taking the limit as the bounding surface ∂A tends to infinity, and using Eqs. (1), (11) and (30) yields Im t α )|ψ(R α )| 2 , then the energy balance becomes Note that where the infinite vector A (P ) and N −vector B (P ) are the solutions for the passive set of impedances. A sort of equivalent reasoning can be derived from Eq. (6) by rewriting it in the matrix form where G = {G (R α − R β )} and ψ 0 = {ψ 0 (R α )}. From (A9) we are interested only in the imaginary part, as it defines the passive or active character of the cluster. Note that the imaginary part of the quadratic form in (A9) is B † Im (G) B and Im (G) ∝ J 0 , hence Im (G) is real-valued and symmetric. Finally, for a globally passive cluster we require Satisfying Im t −1 α ≤ 0 for all individual particles α corresponds to a locally passive metacluster. From Eq. (A1) and (7) it may be noted that B α is a complex force amplitude that acts on the plate. The passivity of a single scatterer can be seen through the Poynting vector -characteristic of the direction of energy flow. Time-averaged energy flow through the point at which a scatterer is placed is where for Im t α ≥ 0 we have Φ ≥ 0, so energy flows from the plate towards the scatterer, i.e. the scatterer is passive. Φ can be seen as the power absorbed by the scatterer.

Appendix B: Some matrix properties
It follows from the definition of S in (14) and Graf's addition theorem for Bessel functions, Eq. where R αβ = |R α − R β |. Note that J 0 (kR αβ ) ≈ 1 at low frequency, indicating that S † S becomes singular in this limit. Numerical examples shows this in terms of the matrix condition number which becomes large at low frequency.
The N × N matrix S † S is therefore real, symmetric and non-negative definite, and can be expressed with positive eigenvalues λ α > 0 and normalized eigenvectors of length N , u † α u β = δ αβ . Using (B2) in the definition of P, Eq. (21), yields where the infinite dimensional vectors V α are V α = Su α , α = 1, . . . , N.
These are orthogonal, V † α V β = λ α δ αβ , but not orthonormal. We define the orthonormal set U α = λ −1/2 α Su α , α = 1, . . . , N, so that P is in canonical form, Hence P is finite rank with N non-zero eigenvalues equal to +1. Alternatively, P is a projection onto the N −dimensional subspace span{U α , α = 1, . . . , N }, and satisfies the projector property We note some other properties of P and related matrices. Multiplying (21) on the right by S and on the left by S † gives PS = S, S † P = S † . (B8) The fundamental matrix S of (14) has an interesting form in terms of the finite and infinite dimensional normalized eigenvectors: The Moore-Penrose inverse of S is Similarly, the matrix Q of (23) is It follows from this, or from its definition in (23), that the matrix Q satisfies where I N is the identity on span{U α , α = 1, . . . , N } Finally, the physical vectors for the far field pattern function and source strengths of Eqs. (20) and (18) are, respectively, where