Designing Collective Non-local Responses of Metasurfaces

We propose a numerically efficient `adjoint' inverse design method to optimize a planar structure of dipole scatterers, to manipulate the radiation from an electric dipole emitter. Several examples are presented: modification of the near-field to provide a 3 fold enhancement in power emission; re-structuring the far-field radiation pattern to exhibit chosen directivity; and the design of a discrete `Luneburg lens'. Additionally, we develop a clear physical interpretation of the optimized structure, by extracting `eigen-polarizabilities' of the system. We find that large `eigen-polarizability' corresponds to a large collective response of the scatterers. This framework may find utility in wavefront shaping as well as in the design and characterisation of non-local metasurfaces.


I. INTRODUCTION
Designing the scattering properties of materials is a fundamental challenge in a broad range of disciplines, from metamaterial design [1] to imaging through disordered media [2][3][4][5][6]. Metamaterials are structured at the sub-wavelength scale to control wave propagation, leading imaging beyond the limit set by diffraction [7]. In recent years, there has been increasing interest in how appropriately designed metamaterials can induce virtually any desired wave effect, be that acoustic [8] or electromagnetic [9]. To solve this problem the connection between incident and scattered fields must be established, before the metamaterial design can be engineered. Great progress in mapping an input field to a desired output field has been made for imaging through disordered media. However, this has involved manipulating the incident wave itself, rather than the scattering system [2][3][4][5][6]. In this work, we propose a semi-analytic framework for designing arbitrary scattering properties of systems of discrete scatterers.
Some of the earliest examples of metasurfaces are frequency selective surfaces [10,11].
These are a class of periodically structured two-dimensional metal-dielectric structures designed to have specific reflection and transmission properties that depend upon the frequency of the incident wave. These effects are facilitated by modifying the electromagnetic boundary condition through structuring of the surface to induce near-field and resonant effects.
These boundary conditions are typically both frequency and wave-vector dependent, and described by a complex effective surface impedance. Frequency selective surfaces have since been used for many applications [12], including as electromagnetic filters and as perfect absorbers [13]. Achieving perfect absorption in disordererd media has also attracted recent attention due to its potential use in several applications [14], such as energy harvesting. In this application, rather than designing surface impedance distributions, properties of the incident field are controlled, with gain and loss in the system being tuned. However the aim of the two approaches is identical: arbitrary control of light.
By introducing spatial variation to the periodic basis of frequeny selective surfaces one may achieve inhomogeneous effective properties. Inhomogeneity allows more abrupt phase shifts to be imparted upon the incident wave, providing additional degrees of freedom that can be exploited when manipulating the scattering of light [15]. This has led to the development of 'metalenses', which can have superior bandwidth to traditional refractive lenses [16,17], as well as metasurface antennas [18]. Metasurface antennas have been designed to exhibit many valuable properties, such as the opportunity to engineer bespoke beamshaping, steering, polarization control and improved efficiency [18]. Instead of introducing inhomogeniety to a repeating basis, several classes of metasurface have been designed by 'engineering disorder'. This involves using an algorithm to selectively place scattering elements to form a metasurface with specific properties. This principle has been used to design metasurface holograms [19], and for wavefront shaping [20]. What unites the seemingly disparate applications of holograms, metalenses, wavefront shaping and imaging through disorder is the problem of designing materials to realise a given wave effect.
The materials required for each of these functionalities can be designed using very similar methods. For example, holograms [19], metalenses [16] and beam shaping surfaces [18] have all been designed using the Gercherg-Saxton algorithm [21]. This method has been used extensively to find the maps of phase offset required to convert an incident plane wave into a given output. Generalisations of this method to include control over the amplitude of the wave have also been developed [22]. However, this method neglects multiple scattering interactions and assumes that only a local phase offset is imparted upon the incident field. This is inhibits the application of this design method to problems were non-local interactions are key, for example in achieving perfect anomalous reflection [23,24]. As well as this, the number of design degrees of freedom is reduced.
Due to broad demand for methods to design the scattering properties of materials, the problem of devising design methodologies has attracted recent attention [25]. As well as the Gerchberg-Saxton, there are two other popular inverse design paradigms. Firstly, geometry optimisation based upon the adjoint design method [26] has been used to design many electromagnetic structures [25,27,28]. Typically this procedure involves evaluating a cost function, which is to be extremised, over a given geometry using a full-wave solver.
Changes to the geometry are then made iteratively, so that the figure of merit is improved until a convergence is reached. A key feature of the adjoint method is that the cost function contains both a 'forward' and an 'inverse' contribution [26]. Reciprocity [29] is exploited to allow the 'forward' and 'inverse' parts to be calculated together, reducing the number of numerical simulations required to determine how material parameters should be changed.
In order to further simplify the numerical complexity axis symmetries [30] and the locally periodic approximation [31], which assumes a very sub-wavelength unit cell, are often ex-ploited. Secondly, machine learning [32] and genetic algorithms [33] have become extremely popular for solving the inverse design problem due to their ability to traverse large search spaces. However while machine learning has a role to play in optimisation processes, human intervention can provide more time-efficient design. It is the ambition of our current work to seek a design method that admits a clear physical interpretations of both the optimisation method and of the results, while being numerically efficient.
In this work we present two contributions. Leveraging the benefits of adjoint algorithms, we propose a semi-analytic framework to design the scattering properties of non-periodic arrangements of discrete dipolar scatterers. Due to the efficiency of this method, we do not assume the weak scattering limit [34], which is often employed when it is assumed that multiple-scattering can be neglected [35]. Instead, all interactions are taken into account so that all multiple-scattering effects are considered. By examining these strongly non-local properties of the entire scattering system, we suggest an interpretation of the eigenvalues of the scattering system. This provides explanatory detail on the mechanisms behind the optimisation procedure. This paper is structured as follows. In Section II we briefly review the standard method used to solve Maxwell's equations in the presence of structures of dipolar scatterers. Then, in Section III we discuss how one can use the Local Density of Optical States to characterise the effect of a photonic environment upon an emitter and propose an interesting interpretation of the eigenvalues of the response matrix. In Section IV we propose an iterative technique, based upon perturbation theory, to design the scattering properties of arrangements of scatterers. To demonstrate how this approach may be applied, in Section V we illustrate the versatility of our technique with several numerical examples of manipulating the dipole field. We successfully demonstrate control over the near-field by enhancing power emission, over the far-field by re-shaping the angular distribution of the Poynting vector and we suggest how a discrete Luneburg lens might be realised by multiplexing such designs.

TERERS
We seek to design the scattering properties of an arrangement of magnetodielectric scatterers. Before proposing our solution to the inverse problem, we review the formulation of solutions to Maxwell's equations due to an arrangement of dielectric scatterers.
In general, the scatterering from a polarizible object is described by a summation over all the possible multipolar modes the object may possess [36]. We assume that only the electric and magnetic dipole modes are excited; this is consistent with several experimental observations [37,38]. Formally, the dipole approximation is justified when the scatterers are sufficiently small (up to ∼ 1/k, with k being the wavenumber) and have separation ≥ 3r 0 [39], where r 0 is the radius of the scatterer.
We begin from Maxwell's equations for monochromatic waves of frequency ω in a general linear medium characterised by permittivity ε and permeability µ, where k = ω √ εµ is the wavenumber, E is the electric field, H is the magnetic field, J is the current density of the dipole emitter and P and M are polarization and magnetisation densities respectively. Under the approximation that the scatterers and the source may be treated as points, the currents in Maxwell's equations (1) and (2) may be written explicitly as, wherep is the polarization of the source and is the polarizability tensor, in the absence of bianisotropy. The location of the source is r and of the n th scatterer is r n . It is important that the polarizability tensor obeys the usual requirements for energy conservation; namely that the energy the polarization/magnetization current absorbs from the incident field is equal to or greater than the energy re-radiated by the current. This can be expressed as a constraint upon the elements of the polarizability [40][41][42], where I is the unit tensor. In general, any differential equation of the form may be solved in terms of the dyadic Green function [43,44] ← → G (r, r ) = I + 1 where I = diag(1, 1, 1) is the unit tensor. The two signs in the exponent correspond to advanced and retarded boundary conditions. Integrating the Green function (6) against the source currents (3), it can be shown that the solution to Maxwell's equations (1) and (2) can be written in terms of the outgoing wave solution and its curl [45].
To simplify notation and make it clear that our solutions are length-scale agnostic, we adopt a dimensionless unit system. We choose a characteristic length, a, by which to scale our coordinate system (for example, a natural choice might be a = 1/k). In terms of this, we can define a unitless wavenumber ξ = ka. Making use of the interchangeability of electric and magnetic fields, we work in units where free space impedance is one: Z 0 = 1. In this unit system, after affecting the integration of the Green function (6) over the source currents (3), the solution for the fields may be written as where E s (r) and H s (r) are the source fields. This is not yet a fully specified solution to Maxwell's equations, the total field applied at the location of each scatterer (E(r n ), H(r n )) is not known. These fields include two contributions: the source, and the scattering from all the other scatterers. To determine the fields (E(r n ), H(r n )), we must impose selfconsistency upon the fields (7). To do this, we substitute r = r n into the field solutions (7) and solve for (E(r n ), H(r n )). This yields the following matrix equation where −1 denotes matrix inversion, and R is the response matrix, defined as The response matrix contains information about how each particle interacts with all of the other particles so that imposing this self-consistency condition is equivalent to solving the multiple-scattering problem for all N scatterers at once. It should be noted that this matrix inversion is potentially substantial, as R ∈ C 3×2×N (three spatial dimensions, for the magnetic and electric fields for each of the N scatterers), making it the most numerically demanding part of our design process. As this is a standard Ax = b matrix problem, solving the self-consistency condition (8) is easily facilitated numerically [46]. Once the fields applied to each scatterer have been found in this way, they may be simply substituted into the expression for the full field (7) so that the total fields for any structure may be calculated.

III. THE LOCAL DENSITY OF OPTICAL STATES AND EIGENVALUES OF THE RESPONSE MATRIX
To quantify the effect of the scatterers upon the radiation of the source, we calculate the Partial (or Polarized) Local Density of Optical States (PLDoS). For a fixed source current, this quantity is proportional to the power output [44,47]. In terms of the total electric field, given by (7), the PLDoS is given by [47] ρ(p, r, ω) This quantity gives the number of electromagnetic modes per unit volume, for a given polarizationp, position r and frequency ω. ρ(p, r, ω) characterises both the changes to the radiation properties of the dipole emitter, and the modification of the electromagnetic modes at the emitter location.
To illustrate the physical meaning of the R matrix, we consider the simplest possible case: two scatterers with only an electric polarizability ( ← → α H = 0). In this case, the response matrix connects the field induced by the source to the dipole moment induced in the scatterers as, The eigenvalues of R −1 satisfy the characteristic polynomial which can be solved by making use of the reciprocity of the Green function [43] to find that the eigenvalues of R are One way to interpret the terms of (12) are as multiple scattering events. An interaction between the two electric dipole scatterers is comprised of scattering from r 1 to r 2 then back again. This is what is expressed by the product of the Green functions in the discriminant of the characteristic equation (12). This understanding can be extended to more electric scatterers and scatterers with both electric and magnetic dipoles, whence more complex scattering processes become available. However, this interpretation is still evident in the form of the polynomials. Combining the fact that eigenvalues of R contain information about the collective response of the particles with (13), the value of eigenvalue itself can be interpreted as the collective polarizability of the two scatterers. Therefore, ← → α −1 E · λ gives the enhancement of the single particle polarizability due to the multiple scattering events between the two scatterers. So, if ← → α −1 E · λ = 1 there is no change to the single particle polarizability and multiple scattering events provide no enhancement. On the other hand, a large ← → α −1 E · λ corresponds to a large enhancement to the response of a single scatterer, due to collective behaviour.
Additionally, the eigenmodes of R represent configurations of field that produce a certain collective response of the system. As R has no symmetries these eigenmodes do not form an orthogonal basis [48,49], although the left and right eigenvectors of R do. Using this left and right pair, the source field can be decomposed into the basis of these eigenvectors, w n as where N is the number of scatterers. The expansion coefficient c n indicates which eigenmodes contribute most strongly to the response of the system. Identifying these modes allows the response of the system to be understood and characterised by examining only a few eigenmodes, rather than the whole expansion (14). The expansion coefficient is a useful tool in characterising the response of the system. From this decomposition, one may find which modes are excited and how strongly so that the dominant response of the system may be isolated and examined.
Now that we have outlined how one can calculate, decompose and interpret the fields of a dipole emitter in the vicinity of dipole scatterers, under the assumptions of the point-dipole approximation, we shall proceed to apply a perturbative approach to design the scattering properties of scatterer distributions.

IV. DESIGNING SCATTERING PROPERTIES
We calculate the figure of merit for the optimisation procedure as follows. Firstly, the location of each scatterer and all of the fields are expanded to first order under a small perturbation to the location of each scatterer, r n → r n + ∆r n , δ(r − r n ) → δ(r − r n ) + ∆r n · ∇δ(r − r n ), Upon substituting these expressions into the solutions for the fields (7), retaining only terms to first order, we obtain the following expressions for the corrections to the fields in terms of the unperturbed fields as If one wishes to increase the power emission of the dipole, the figure of merit is ∝ Im[E 1 (r )].
This is equivalent to increasing the PLDoS (10), producing more decay channels for the dipole emitter. To structure the far-field, a point is the far field is chosen so that the figure of merit becomes ∝ Im[E 1 (r farfield )], so that directivity in the direction of r farfield is enhanced.
In this way, the figure of merit may be iteratively increased. It is clear that this is an adjoint calculation: the gradient operators act on the total field applied [50] to the scatterer located at r n , meaning that the effect of the change of the position of r n upon the field at all of the other scatterers is implicitly included in the calculation.
A schematic of how our proposed optimisation procedure works is shown in Figure 1.
It is important to note the checks performed. As the scatterers are modelled as points, it is necessary to ensure they remain properly separated so that the dipole approximation remains valid. To guarantee this, scatterers are not allowed to move within a smoothing distance d 0 of each other. It is also possible that there is no move the scatterer can take which improves the figure of merit, so these moves are also blocked.
To elucidate this process, we now present several examples of the application of the procedure outline in Figure 1.

V. NUMERICAL EXAMPLES
To demonstrate the versatility and strengths of our proposed method, we apply it to solve several inverse design problems. Here, we shall demonstrate the ability of our procedure to design the near field and far field equally, as well as how the results generated can be used to construct devices with more complex purposes. The situation we consider is shown in  Figure 1(c). The optical properties of silicon have been extracted from experimental data [51], then combined with the Mie a 1 and b 1 coefficients [52], to allow the calculation of electric and magnetic polarizabilities [53]. Parameters used for all simulations are shown in Table I. It has been assumed that the scatterers are isotropic, so the electric and magnetic polarizability tensors have the form The complex numbers α E and α H are calculated according to [52] as where a 1 and b 1 are the Mie coefficients for the dipole modes of a spherical scatterer.
Before proceeding to the design of complex structures, we consider the simple case of two electric scatterers. The eigenvalues for this situation are given by (13), and the eigenvectors  can be found from the response matrix in (11). The progression of the optimisation procedure is shown in Figure 2 process [54]). Indeed, we see that the solution presented in Figure 2(a) is probably a local minima, rather than a global one.
To demonstrate the ability of our method to manipulate the near field of a source, we show how structures can be designed to enhance the power emission of a dipole. The results of this are shown in Figure 3. An arrangement of 100 dipolar scatterers has been designed using our proposed framework to provide a factor of ∼ 3 enhancement of the power emission of the dipole emitter. This factor of enhancement is far smaller than can be achieved with 3D bulk structures [28], but is of the order of similar works that have utilised genetic algorithm techniques [33].   Figure 3(f). This large collective response of the scatterer system, as well as being strongly excited by the dipole's field, leads to the enhancement of power emission.
Next, we demonstrate the ability of our method to manipulate the far-field of a dipole emitter. We apply our method to enhance directivity of a dipole emitter along a given direction. The results of this optimisation are shown in Figure 4. We have sought to enhance directivity along the θ = 0 direction, where θ is the polar angle in the x − y plane.
Beginning from a square array, with a far-field radiation pattern shown in Figure 4 in the x − z plane is φ HP BW = 1.6 • . To use the common [55] antenna directivity estimate D = 10 log(41000/(θ HPBW φ HPBW D 0 )), where D 0 is the reference directivity. For an isotropic source D 0 = 1 and for a dipole D 0 = 6. This approximation ignores back-lobes. Using this, we estimate the structure shown in Figure 4(c) has directivity 70 dB above isotropic and 65 dB above a dipole.
In addition to designing a structure with the desired far-field properties, by applying our understanding of the eigenmodes of the system the mechanism for effective performance can be revealed. The change in the eigenvalues and expansion coefficients of the modes is shown in Figure 4(d). We note that the eigenvalues and expansion coefficients do not change magnitude considerably as a result of the optimisation. This is perhaps not very surprising, as the aim of the optimisation was not to enhance power emission but rather to re-structure the field. So, rather than designing modes with a large eigenvalue and therefore a large power emission enhancement due to a collective response, the optimisation procedure has designed modes with a certain shape. This is demonstrated in Figure 4(e-h), where the leading order modes, characterised by expansion coefficient, of the system before and after the optimisation are plotted, with the far-field angular distribution of the associated Poynting vector inset.
Both of the modes have similar expansion coefficients and eigenvalues, but very different spatial distributions. Indeed, for this application the optimisation procedure re-shapes the modes rather than enhancing multiple scattering effects. Now, we shall demonstrate how the results of our optimisation procedure may be utilised to achieve more complex functionality. By rotating the structure designed to enhance directivity, as well as the source, a device can be constructed which has the same functionality as a Luneburg lens [56]. Using this device, a point source may be converted to a plane wave with a specific propagation direction. This is shown in Figure 5. A conventional Luneburg lens is made by grading the refractive index profile, so that a point source placed upon the surface of the lens is converted into a plane wave. Its ray diagram , together with the refractive index profile required to achieve this function, is shown in Figure 5(a).
Exploiting the spherical symmetry of the refractive index profile, by rotating the source around the edge of the lens the plane wave may also be rotated. Multiplexing the result shown in Figure 4 can achieve something qualitatively similar, as demonstrated in Figure  5(b-c). One can use the indicated structure to convert a dipole source into a beam directed along a single angle. By placing the source at different locations within the structure, with the correct orientation, the beam may be rotated. It should be noted that additional side lobes can be seen in the far-field Poynting vector in Figure 5(b) compared to Figure 4(b).
This is due to interaction between the different multiplexed elements that make up the structure and can be reduced if the spacing between elements was increased. The structure proposed in Figure 5(c) can be fabricated without having to grade an index, instead 288 scatterers must be arranged as indicated. The multiplexed device has a radius comparable to conventional Luneburg lenses, at ∼ 6λ. While the fabrication is more straightforward, the angular resolution is not continuous as for the usual Luneburg lens. Instead, an angular resolution of 45 • can be achieved, although by making the structure larger a higher resolution could be achieved. It was found that the relationship between device radius and angular resolution was well approximated by the following power law: resolution (degrees) = 83.6 × (R/λ) −0.66 . Therefore, to achieve a resolution of 2 • , the device would have to be ∼ 50λ.

VI. SUMMARY AND CONCLUSIONS
In this work, we have derived a method of designing metasurfaces comprised of dipolar scatterers. This has been applied to structure both the near-field and the far field of a dipole emitter. In the near-field, power emission has been enhanced by a factor of ∼ 3 and the far-field radiation pattern has been re-structured. We have also demonstrated that structures designed in this way may be multiplexed to achieve more complex functionality.
As an example, we approximate the functionality of a Luneburg lens using an array of dipole scatterers.
As well as an iterative design methodology, we propose an interesting physical interpretation of the eigenvalues and eigenvectors of the matrix defining the electromagnetic response of the scattering system. The eigenvalues of the system correspond to eigen-polarizabilities which we attribute to several scatterers responding collectively. A large collective response corresponds to a large eigenvalue. By analysing how the eigenvalues and eigenvectors change over the optimisation procedure, we have identified that power emission is enhanced by a large collective response of the scatterers, corresponding to a large eigenvalue while directivity is achieved by modifying the spatial distribution of the modes, without significant change to the eigenvalues.
The applicability of both our design technique and theoretical understanding are not limited to engineering dipole radiation. A perturbative approach to designing electromagetic field properties might be applied to engineering mode distributions in optical fibers or metalenses. If more arbitrary field distributions could be successfully designed, then this method might find utility in constructing metasurface holograms or to perform wavefront shaping when imaging through disordered media. As our approach automatically takes non-locality into account, it may be used to develop and provide insight into non-local metasurfaces.