Designing the collective non-local responses of metasurfaces

The ability to design the electromagnetic properties of materials to achieve any given wave scattering effect is key to many technologies, from communications to cloaking and biological imaging. Currently, common design methods either neglect degrees of freedom or are difficult to interpret. Here, we derive a simple and efficient method for designing wave–shaping materials composed of dipole scatterers, taking into account multiple scattering effects and both magnetic and electric polarizabilities. As an application of our theory, we design aperiodic metasurfaces that re-structure the radiation from a dipole emitter: (i) modifying of the near-field to provide a 4-fold enhancement in power emission; (ii) re-shaping the far-field radiation pattern to exhibit chosen directivity; and (iii) the design of a discrete Luneburg–like lens. Additionally, we develop a clear physical interpretation of the optimised structure, by extracting eigen-polarizabilities of the system, finding that a large eigen-polarizability corresponds to a large collective response of the scatterers. Mapping the response required to obtain the desired output field from a metamaterial for a given input remains challenging due to the complexity of scattering processes. Here, a semi-analytical model for the design of aperiodic metasurfaces is presented alongside numerical examples.

D esigning the scattering properties of materials is a fundamental challenge in a broad range of disciplines, from metamaterial design 1,2 to imaging through disordered media [3][4][5][6] . Metamaterials are structured at the sub-wavelength scale to control wave propagation, and their properties can yield imaging beyond the limit set by diffraction 7 , or cloaking from an incident field 1, 8 . In recent years, there has been increasing interest in how appropriately designed metamaterials can induce virtually any desired wave effect, be that acoustic 9 or electromagnetic 10 . To solve this problem the connection between the incident and scattered fields must be established, and then this must be used to design the metamaterial structure. By contrast, great progress in mapping an input field to a desired output field has been made for imaging through disordered media. In this case, the incident wave is modified to give a desired output field, and the properties of the scattering system remain fixed [3][4][5][6]11 . Achieving perfect absorption in disordered media has also attracted recent attention due to its potential use in several applications 12 , such as energy harvesting. However, it is usually completely unclear how the scattering system has to be modified to give a desired output field.
Some of the earliest examples of metasurfaces are frequency selective surfaces 13,14 . 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 the 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 15 , including as electromagnetic filters and as perfect absorbers 16 .
By introducing spatial variation to the periodic basis of frequency 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 17 . This has led to the development of metalenses, which can have superior bandwidth to traditional refractive lenses 18,19 , as well as metasurface antennas 20 . Metasurface antennas have been designed to exhibit many valuable properties, such as the opportunity to engineer bespoke beamshaping, steering, polarisation control and improved efficiency 20 . 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 21 , and for wavefront shaping 22 . 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 21 , metalenses 18 and beam-shaping surfaces 20 have all been designed using the Gercherg-Saxton algorithm 23 . 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 24 . However, this method neglects multiple scattering inhibiting the application of this design method to problems where non-local interactions are key, for example in achieving perfect anomalous reflection 25,26 and refraction 27 . As well as this, the number of design degrees of freedom is reduced. Furthermore, there are several experimental situations where multiple scattering cannot be neglected 28 .
For example, in thin dielectric resonators or resonators with low index contrast, evanescent fields cause non-negligible coupling between the building blocks of the metasurface. Being able to design metalenses using ultra-thin resonators will be key to further flattening optical components.
Due to broad demand for methods to design the scattering properties of materials, the problem of devising such methodologies has attracted recent attention 29 . As well as the Gerchberg-Saxton, there are two other popular inverse design paradigms. Firstly, geometry optimisation based upon the adjoint design method 30 has been used to design many electromagnetic structures 29,31,32 . 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 30 . Reciprocity 33 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 34 and the locally periodic approximation 35 , which assumes a very sub-wavelength unit cell, are often exploited. Secondly, machine learning 36 and genetic algorithms 37 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 arbitrary scattering properties of non-periodic arrangements of discrete dipolar scatterers to give bespoke functionality. To illustrate our framework we shall explore the design of a metasurface composed of a plane of scattering particles arranged in an aperiodic fashion. Due to the efficiency of adjoint methods, we do not assume the weak scattering limit 38 , which is often employed when it is assumed that multiple scattering can be neglected 39 . 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. Together, these techniques provide a paradigm for the inverse design of a wide class of metamaterials. For metamaterials comprised of scatterers that are sub-wavelength but strongly coupled, the methods proposed here provide make the design procedure transparent and efficient as well as providing key physical insight into the function of more exotic designs.

Results and discussion
Solving Maxwell's equations for a system of discrete scatterers. 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, using the discrete dipole approximation 40,41 .
In general, the scatterering from a polarizible object is described by a summation over all the possible multipolar modes the object may possess 42 . We assume that only the electric and magnetic dipole modes are excited; this is consistent with several experimental observations 43,44 . For isotropic spherical scatterers, 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 45 , 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 ¼ ω ffiffiffiffi ffi εμ p 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 polarisation 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 polarisation of the source, E and H are the fields applied to the scatterer and is the polarizability tensor, in the absence of bianisotropy. The location of the source is r 0 and of the nth scatterer is r n . It is important that the polarizability tensor obeys the usual requirements for energy conservation; namely that the energy the polarisation/magnetisation 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 tensor: Im½α
In general, any differential equation of the form may be solved in terms of the dyadic Green function 49,50 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 G $ EH ðr; r 0 Þ ¼ ∇ G $ ðr; r 0 Þ 51 . 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

EðrÞ
HðrÞ 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 self-consistency 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 It should be noted that as α $ is defined as the response to the applied fields, any self-interaction terms of the form G $ ðr i ; r i Þ should be removed from the interaction matrix. The response matrix contains information about how each particle interacts with all of the other particles so that imposing this selfconsistency condition is equivalent to solving the multiplescattering problem for all N scatterers at once. It should be noted that this matrix inversion is potentially substantial, as R 2 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 52 . 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.
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 polarised) local density of optical states (PLDoS). For a fixed source current, this quantity is proportional to the power output 50,53 . In terms of the total electric field, the PLDoS is given by 53 This quantity gives the number of electromagnetic modes per unit volume, for a given polarisationp, position r and frequency ω. The PLDoS ρð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. All multiple scattering effects within the system are encoded by the interaction matrix R. To illustrate the physical meaning of the R matrix (9), we consider the simplest possible case: two scatterers with only an electric polarizability ( α $ H ¼ 0). In this case, the response matrix connects the fields due to the the source, (E s , H s ), to the dipole moment induced in the scatterers p as The eigenvalues of R −1 satisfy the characteristic polynomial which can be solved by making use of the reciprocity of the Green function 49 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 also 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 characteristic polynomials. Combining the fact that eigenvalues of the interaction matrix R contain information about the collective response of the particles with the expression for the eigenvalues (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 the interaction matrix 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 54,55 , although the left and right eigenvectors of R do. Right eigenvectors are defined as satisfying Rw = λ R w, and left eigenvectors satisfy vR = λ L v. 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 scattering 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.
Designing scattering properties. To design the scattering properties of the system of dipole scatterers, we begin by expanding the fields under a small perturbation to the locations 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 0 Þ. 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. The application of this procedure to more complicated figures of merit is given in Supplementary Note 1. For each iteration i, the position of the nth scatter is updated according to 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 56 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 Fig. 1.
To elucidate this process, we now present several numerical examples of the application of the procedure outline in Fig. 1.
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 Fig. 1a. A dipole emitter at 550 nm with polarisationp is located at r 0 . Near this, we place an array of discrete scatterers. Each scatterer is modelled as isotropic silicon sphere of radius 65 nm, and polarizabilities calculated to match this. The polarizabilities of the scatterers are shown in Fig. 1b. A figure of merit is chosen and according to (18) the locations of the scatterers are iteratively updated to improve this figure of merit. The result of this procedure is shown schematically in Fig. 1c. The optical properties of silicon have been extracted from experimental data 57 , then combined with the Mie a 1 and b 1 coefficients 58 , to allow the calculation of electric and magnetic polarizabilities 59 . Parameters used for all simulations are shown in Table 1. The consistency of the discrete dipole approximation with these parameters has been verified by comparison with full-wave solutions, presented in Supplementary Note 2.
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 58 as where a 1 and b 1 are the Mie coefficients for the dipole modes of a spherical scatterer. More complicated scatterers could be used, provided their scattering is dominated by the dipole mode. The extraction of electric and magnetic polarizabilities for complex scatterers can be done numerically 60,61 . Before proceeding to the design of complex structures, we consider the simple case of two electric scatterers, shown in Fig. 2a. 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 Fig. 2b. The initial and final eigenmodes and eigenvalues are shown in Fig. 2c, d. It can be seen that the eigenvalues of the modes change very little as the optimisation process progresses. This is unsurprising in this case, as from Fig. 2a it can be seen that the optimisation has resulted in only a small translation of the scatterers. Instead of the modes themselves being extensively modified, the change occurs in which modes are excited. Examining how the expansion coefficients change, we find that the expansion coefficient of modes 2 and 4, both of which are symmetric, are enhanced. Expansion coefficients of all antisymmetric modes remain small. This is the origin of the power enhancement. The magnitude of the enhancement is small as the size of the eigenvalue of modes 2 and 4 are very close to unity meaning that the collective response of these modes is small, indicating weak coupling between the scatterers.
This simple example demonstrates some important things. Firstly, there are two mechanisms by which the optimisation procedure may increase the figure of merit. The eigenmodes of the system may themselves be changed. This may present as a change of the spatial distribution of the mode, or the increase of eigenvalues. Next, with only two scatterers the enhancements we can achieve are small. Indeed, this is congruent with the interpretation that large eigenvalues, corresponding to large collective responses, represent multiple scattering events.
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 Fig. 3a-c, with the progression of the eigenmodes shown in Fig. 3d. An arrangement of 100 dipolar scatterers has been designed using our proposed framework to provide a factor of~4.5 enhancement of the power emission of the dipole emitter. This factor of enhancement is far smaller than can be achieved with 3D bulk structures 32 , but is of the order of studies with similar assumptions and constraints but which use genetic algorithms 37 . Examining the change in the eigenmode with the Emitter location (0, 0, 0)   largest expansion coefficient, by comparing Fig. 3e with Fig. 3f, two key qualitative features of modes that enhance power emission can be determined. Firstly, the mode has a large eigenvalue corresponding to a large collective response of the scatterers. Secondly, this mode also exhibits a strong localisation at the location of the emitter. Clearly, this mode and the field from the dipole have a large overlap, demonstrated by the large expansion coefficient of the mode shown in Fig. 3f. The 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. The results presented in Fig. 3 show the strength of designing aperiodic structures, rather than being limited to periodic ones. The design of aperiodic structures benefits from an improved exploration of the parameter space due to the removal of the restriction to periodic solutions. While it has been demonstrated that periodic structures can enhance antenna radiation 62 , more compact and better-performing solutions can be obtained with aperiodic structures, as we have demonstrated in Fig. 3. The problem of antenna emission is not periodic, so one should not expect an emission enhancing structure to be periodic. To find an optimal, or even well performing, solution one is forced to consider aperiodic strucutres. Indeed, the bow-tie shape of the structure shown in Fig. 3c is consistent with structures designed to enhance dipole radiation using genetic methods 63 as well as simple phase arguments 32 . The resulting structures have the same symmetry as dipole radiation: they are left-right symmetric, but not periodic. To further accelerate our method, one could impose the symmetry of the problem upon the solution. Such a constraint would be useful when designing metasurfaces with a very large radius~100λ, as only half of the field gradients would need to be calculated 19 . The additional parameters available to optimise when designing aperiodic structures can present a challenge for traditional methods based on gradient or look-up methods. However, in the method we propose here, these additional degrees of freedom can be optimised while keeping the numerical problem efficient.
Next, we demonstrate the ability of our method to manipulate the far-field of a dipole emitter. We have sought to enhance directivity along the θ = 0 direction, where θ is the polar angle in the x-y plane. The results of this optimisation are shown in Fig. 4. Beginning from a square array, with a radiation pattern shown in Fig. 4a, a point in the far-field has been chosen with the power radiated to this location being the figure of merit. The result is a clear enhancement of the directivity along the θ = 0 direction, shown in Fig. 4b, c. The starting and optimised scatterer arrangements are shown in Fig. 4d. The beam-width (FWHM) is~20°.
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 Fig. 4e. The eigenvalues and expansion coefficients do not change magnitude considerably as a result of the optimisation. This is demonstrated in Fig. 4f-i, where the leading order modes, characterised by expansion coefficient, of the system before and after the optimisation are plotted, along with the far-field Poynting vector in each case. Both of the modes have similar expansion coefficients and eigenvalues, but very different spatial distributions. Indeed, for this application the optimisation procedure aims to re-shapes the modes rather than enhancing multiple scattering effects.
The operational bandwidth and sensitivity to disorder of the designs shown in Figs. 3, 4 are given in Supplementary Note 3. In the above results, only radiation properties in the plane have been engineered however out solutions are fully 3D and properties of the system can be engineered in 3D. A demonstration of engineering out-of-plane directivity is given in Supplementary Note 4. 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 similar functionality to a Luneburg lens 64 . Using this device, a point source may be converted to a beam in a specific propagation direction. This is shown in Fig. 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; this is shown in Fig. 5a. Exploiting the circular 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 Fig. 4 can achieve something qualitatively similar, as demonstrated in Fig. 5b, 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. The structure proposed in Fig. 5b 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, in degrees, was well approximated by the following power law: (R/λ) = 83.6 × resolution −0.66 . Therefore, to achieve a resolution of 2°, the device would have to be~50λ. The data used to estimate the power law are given in Supplementary Note 5.

Conclusions
In this work, we have derived a method of designing metamaterials comprised of small scatterers. This has been applied to design aperiodic planar structures that have a predetermined effect on both the near-and the far-field of a dipole emitter. In the near-field, power emission has been enhanced by a factor of 4 and the far-field radiation pattern has been made highly directional. 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 a 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 is not limited to engineering dipole radiation. A perturbative approach to designing electromagnetic field properties might be applied to engineering mode distributions in optical fibres 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 also used to develop and provide insight into non-local metasurfaces.

Data availability
All data created during this research are openly available from the University of Exeter's institutional repository https://doi.org/10.24378/exe.3504.   Fig. 4 to produce a Luneburg-like lens with a discrete angular resolution of 45°. By rotating the source and changing its location inside the array the far-field Poynting vectors indicated in (c) can be observed. For example, placing the emitter at the centre of array 2 in (b), and rotating it by 45°clockwise, results in the far-field directivity labelled as 2 in (c).