Designing disordered multi-functional metamaterials using the discrete dipole approximation

The ability to design passive structures that perform different operations on different electromagnetic fields is key to many technologies, from beam-steering to optical computing. While many techniques have been developed to optimise structures to achieve specific functionality through inverse design, designing multi-function materials remains challenging. We present a semi-analytic method, based on the discrete dipole approximation, to design multi-functional metamaterials. To demonstrate the generality of our method, we present two key examples. Firstly, we work at optical wavelengths to design a disordered 2D arrangement of silicon spheres that beams light into different directions depending on the source polarisation. Secondly, we design a 3D device that works at microwave wavelengths and sorts plane waves by their angle of incidence. In this case, the scatterers are more complicated meta-atoms, with a strong dipole resonance at microwave frequencies.


Introduction
Much modern technology is enabled by the control of electromagnetic waves. One way to manipulate electromagnetic radiation is to structure materials in space [1][2][3], and more recently time [4,5], to control electromagnetic fields. Metamaterials, man-made materials that have electromagnetic properties associated with their structure rather than chemistry, have demonstrated extraordinary control over all types of wave [6]. Typically metamaterials are built from 'meta-atoms', sub-wavelength resonant scattering elements with scattering properties that depend upon their structure. By tuning the geometry of the meta-atoms, almost any wave scattering effect can be realised. Of particular interest is the ability to design multi-functional metamaterials: materials that scatter differently in response to different input fields. Being able to passively perform different operations upon different waves is key to several problems across physics; from multiplexing output beams [7] to mode sorting [8,9] and beam steering [10]. Furthermore, metamaterials comprising several discrete meta-atoms that can be arranged anywhere in space, and that are explicitly aperiodic, are typically designed using genetic algorithms [11]. These are extremely effective at exploring large search spaces, but are numerically expensive and produce results that can be difficult to interpret [12].
To understand why the design of multi-functional materials is challenging, we consider an electric field φ(r) of a fixed frequency ω = ck 0 , with wave number k 0 , in a material with a spatially varying permittivity ε(r). This wave obeys the vector Helmholtz equation, The difficulty in the design of multi-functional materials is the problem of finding a single material distribution (here the permittivity, ε) that performs two (or more) desired wave transformations. This means that both of the desired wave behaviours, φ 1 (r) and φ 2 (r) must be solutions to the same Helmholtz equation, From this statement, we can find a condition upon the two wave-fields for this to be possible. Multiplying the first of these by φ 2 (r) and the second by φ 1 (r), then taking the difference eliminates the material properties such that φ 2 (r) · ∇ × ∇ × φ 1 (r) − φ 1 (r) · ∇ × ∇ × φ 2 (r) = 0.
Integrating this over all space, then using Green's vector identity [13], we find that This condition, equivalent to reciprocity, places a stringent constraint on the two wave-fields if they are to be supported by the same material, which can be used to derive fundamental bounds on the performance of multi-functional devices [14]. Equation (4) is a generalization of Poynting's theorem, representing the conservation of the norm of the system modes; ensuring for example, their orthogonality. To better understand the connection between the above equation and energy conservation, consider the special case where we demand the same permittivity distribution supports the solution φ 1 (r) = E and its complex conjugate (time reverse) φ 2 (r) = E * . The surface integral (4) can then be re-written using the divergence theorem, Applying Maxwell's equations to convert the curls into magnetic fields which is the usual expression for energy conservation expressed in terms of the Poynting vector A direct application of equation (4) to design multi-functional materials is generally difficult. Instead several other design methodologies have emerged recently [15]. Where the function of the metasurface, a 2D metamaterial, is to impart a known phase and amplitude offset, the Gerchberg-Saxton algorithm [16] is commonly used [17]. While simple and efficient this method can struggle to capture the coupling between elements of the metamaterial, requiring either strong field confinement within, or large spacing between the elements. As well as this, several full-wave simulations are required to build up a library of phase and amplitude changing meta-atoms. For problems where a continuous permittivity distribution is required, 'topology optimisation' [18] can be used to design structures that sort waveguide modes [8] or perform wavelength-dependant behaviour [9]. This method typically uses a full-wave solver to find the fields, then employs gradient descent optimisation to design a structure that optimises a figure of merit. For large structures this can be computationally demanding, however the adjoint method can be used to optimise a figure of merit by converting a shape derivative for the entire design into two field calculations [19]. Even with this greatly improved numerical efficiency, many full-wave simulations are still required over the course of the optimisation.
We present a general framework to design multi-functional metamaterials in electromagnetism in order to overcome some of the aforementioned difficulties. Based on the discrete dipole approximation [20], the framework proposed in this paper is numerically efficient while still being easily applicable to a wide range of electromagnetic problems. The discrete dipole approximation is briefly described in section 2, then a summary of how metamaterials can be designed within this approximation is provided in section 3. Considering the problem of shaping the far-field of a dipole emitter while also increasing its efficiency, this is then extended to design multi-functional materials in section 4. Several examples of utilising this framework are then shown in section 5. The first device we present operates at optical wavelengths λ = 550 nm, with silicon nanospheres used as the scattering elements, and beams light into different directions based upon the polarisation of the source. The second device we show works at microwave wavelengths ∼20 mm (15.5 GHz), using 'metacubes' [21] as the scatterers, and sorts signals by their incidence direction.

Modelling metamaterials
Before trying to design the scattering properties of a metamaterial, it is necessary to characterise the effect of the material upon an electromagnetic wave. Maxwell's equations for a fixed frequency ω = ck 0 , where k 0 is the wave-number, can be written as where the material is characterised by the polarisation density P and magnetisation density M and the incident field is (E inc , H inc ) T . The incident field could be due to an emitter near the structure, or might be a plane wave. Assuming that the scatterers are sub-wavelength in size means that the field can be treated as constant over the scatterer, and so the scatterer's spatial distribution can be treated as a delta-like point at its center r n . The scatterer then acquires an electric and magnetic dipole moment in response to the applied fields where α E is the electric polarisability tensor and α H is the magnetic polarisability tensor. The wave-equation with delta-like source terms can be solved with the appropriate Green's function [22,23] as We have adopted the compact notation where is the dyadic Green's function and ξ is a dimensionless wave-number. To completely specify the field solution (9), the fields applied to each scatterer φ(r n ) must be determined. This includes the source field, as well as all orders of multiple scattering interactions between the scatterers themselves. The applied fields φ(r n ) can be found through the self-consistency condition [24] with R nm = 1δ nm − G(r n , r m ), φ m = φ(r m ) and φ i,n = φ i (r n ). This forms a linear system that can be solved for the fields applied to the scatterers φ m using standard matrix methods. Once these are found, the fields (9) are fully specified.

Designing metamaterials
The scattering properties of a metamaterial made of several sub-wavelength scatterers can be designed using the methodology presented in [25]. We briefly review this here, before extending the method to multi-functional devices. By Taylor expanding the delta function sources in the wave equation (8) under small changes in the position of a scatterer yields, and by keeping only terms linear in δr n an expression for how the field φ(r) changes due to a small change in the location of scatterer n can be deduced, A figure of merit F can be written as a functional of the field F = F [φ]. Under small changes in the field as a result of a small change in the position of a scatterer, the figure of merit is changed by a small amount. Given a particular figure of merit, an analytic expression for this change δF [φ, δφ] can be derived. Since this will be linear in the change in position of a scatterer δr n , the resulting expression allows for the derivative of the figure of merit with respect to the scatterer locations ∂F ∂r n to be calculated analytically. As an example, say the figure of merit is the amplitude squared of the field at a particular location |φ(r )| 2 . Expanding under small changes in the field we can find how the figure of merit changes, Substituting into this in the variation of the field (13) and dividing by δr n we find that This approach gives an analytic expression for the derivative of a figure of merit that can be evaluated for all of the scatterers at the same time then used in a gradient descent optimisation [26] r i+1 where γ is the learning rate and i is the iteration number. While the derivatives of the fields still need to be evaluated, the derivative of the figure of merit, which is more numerically expensive, does not. Extending this to apply to multi-functional metamaterials, requires one to seek to increase some set of figures of merit {F 1 , F 2 , F 3 . . .}. A composite figure of merit can be constructed that is a weighted sum of these, where w i are the weights for each figure of merit. This composite figure of merit can be optimised in the same way as a single figure of merit, where the gradient in equation (18) is now Key to the success of this method is a sensible choice of the weights, w i . An appropriate choice can be informed by considering the desired properties of the resulting device, as we explore in the next section.

Multi-objective optimisation considerations
Consider a simple example problem, shown in figure 1. The goal is to distribute scatterers around a point emitter at location r with polarisation p =ẑ such that two figures of merit are simultaneously maximised. More specifically, the goal is to re-shape the radiation pattern of an emitter, while simultaneously increasing the efficiency. The first figure of merit is the power emitted from the dipole, The second is the overlap integral between the angular distribution of the Poynting vector in the far-field, |S(θ)| and a desired angular distribution ψ T (θ) where the angle θ is in the same plane as the metasurface. In the following examples, the target distribution is Both of these can be expanded to first order to find the gradient of the figure of merit with respect to the scatterer locations [25]. These expansions are given in appendix A. For convenience, both will be normalised by their free-space values, P 0 and Choosing the weights to be proportional to 1/F i means that when the figure of merit is small the contribution of the gradient associated with that figure of merit to the sum (20) is large, but when the figure of merit is large the contribution is suppressed. Note that the figures of merit must be normalised so that their magnitudes can be meaningfully compared, for example by division by a free space value. Choosing the weights in this way allows for the design of multi-functional metamaterials built from discrete scatterers, for a variety of applications. A few examples are offered in the following section.

Multi-functional devices
Continuing with the example developed in the previous section and shown in figure 1, the multi-objective problem of designing a structure that re-shapes the radiation pattern of an emitter in the plane of the metasurface, while also increasing efficiency, is addressed. We work at a wavelength of λ = 550 nm and the scatterers are small silicon spheres of radius 65 nm. For this system, the polarisability tensor can be found analytically from the Mie a 1 and b 1 coefficients. Our figures of merit are the power emission (21) and the overlap with the desired radiation pattern (22) and we choose the weights according to (24). Expressions for the gradients of these figures are merit can be found analytically and are given in appendix A. Using these gradients, weights (24) and the gradient descent method (18), we design the structures shown in figure 2.
For comparison, we also consider the single-objective case where only the far-field radiation pattern is shaped. The resulting far-field radiation patterns are shown in figure 2(a), with the path of the optimisation in solution space shown in figure 2(b) and the two resulting structures shown in figures 2(c) and (d).
Examining first the radiation pattern, we note that the multi-objective optimisation produces a slightly worse match to the target distribution than the single-objective case. This is due to the trade-off between the two figures of merit we seek to optimise. In solution space, figure 2(b), we see that in the case where only the radiation pattern is shaped (blue line) only a very small change in emitted power is seen. Conversely, when both power and radiation pattern are optimised, the emitted power approximately doubles while the overlap integral also increases. Unlike the case for a single scatterer shown in figure 1, it is impossible to plot the whole search space and visualise the location of the Pareto front. The scatterers have a diameter ∼ λ/4 and our solution box has size 8λ, meaning that there are 32 2 = 1024 possible 'pixels' a scatterer could occupy. This means that for N scatterers, the number of possible solutions is 1024 ! (N ! (1024 − N)!). For N = 64, this is ≈10 21 . Despite the exceedingly large search space, which even a genetic algorithm would explore only a very small portion of, our method has found a solution that performs well. A comparison between solving this problem using the method we present and a genetic algorithm is given in appendix B. The second example we consider is manipulating the radiation pattern based on source polarisation. Again, we work at λ = 550 nm and use 65 nm silicon spheres as the scatterers. We aim to create beams at angles θ i , associated with source polarisation p i . Our figures of merit are therefore The expansion of this to find analytically the gradient is given in appendix A. We consider the source polarisation being either left or right circularly polarised, i.e.
The Poynting vector can then be expanded to first order to find the derivatives of the figures of merit for the optimisation procedure. Figure 3(a) shows the radiation patterns of the designed structure excited by each of the two different sources we consider. For a right-handed source, the target angle is θ = 0 • and for a left-handed source, θ = 100 • . The far-field Poynting vector in the plane of the metamaterial, figure 3(a), also shows clear peaks at the desired locations, which are also evident in the near-fields shown in figures 3(b) and (c). The full far-field sphere is given in appendix C. The path in solution space, figure 3(d), shows that the choice of weighting has ensured that the performance of both figures of merit remain similar over the optimisation and in the final result. The multi-functionality condition (4) is considered in figure 3(e). Forming the vector fields φ i × ∇ × φ j on the surface of a sphere enclosing the structure, integrating over the surface and taking the difference yields a result of the order 10 −8 , within expected numerical error. The third and final example we consider is designing a device for beam sorting. Working at 15.5 GHz and using 'metacubes' [21] as the scattering element. The metacubes, formed of six metal faces joined by three connecting spokes, exhibit a strong dipole resonance at 15.5 GHz. Due to their complexity the polarisability tensor cannot be found analytically. Instead, one can model a single scatterer under plane-wave incidence using a full-wave solver such as COMSOL [28] and integrate over the currents to find the electric and magnetic dipole moments [29,30], from which one can find the polarisability tensors. Optimisation of a structure of many complex scatterers using such full-wave methods quickly becomes intractable. Our method presents the key benefit of being able to model large systems of potentially . Design of a multi-functional device, with its operation depending upon the direction of incidence of a plane wave. If the wave is incident from ±20 • , the wave is focused to different locations. The performance of the device under −20 • incident is shown in (a) and under 20 • incidence in (b), with direction of incidence shown as a white arrow. We work at 15.5 GHz using metacubes, shown inset in (a), as the scatterers. Panel (c), shows cuts of the fields under the two different incidence angles, demonstrating peaks at the target positions. The path in solution space, (d), shows that the two figures of merit progress at the same rate over the optimisation, leading to a structure with roughly equal performance for each figure of merit. The initial and final structure are shown in (e) and (f). The multi-functionality condition (4) is verified in (g) by evaluating the surface integrals. The numerical error here is larger than the results in figure 3 due to the strongly oscillatory nature of the integrand, making evaluation of the surface integral more sensitive. complicated scatterers, provided they can be approximated as dipoles, although it is possible to include higher order multipoles into the formalism [31][32][33]. The validity of using the discrete dipole approximation to describe these systems is verified by comparison with full-wave solution in appendix D. We seek a structure of metacubes that takes plane waves from different directions and focuses them to distinct points. A device of this sort could be used, for example, to detect from which direction a signal is coming. The figures of merit for this problem are: where r i denotes the location to focus the wave at for incident direction i and E i is the electric field produced by the structure under incidence from direction i. The gradients of these figures of merit are given in appendix A. The structure resulting from this optimisation is shown in figure 4. Operation of the device when driven by a TE plane wave incident at 20 • is shown in figure 4(a) and for a plane wave at −20 • in figure 4(b). The two different focus points (b) are clearly visible in the fields. The path in solution space is shown in figure 4(d), where again the choice of weighting has ensured roughly equal performance of F 1 and F 2 . Slices of the fields from figures 4(a) and (b) are shown in figure 4(c), indicating the large main peaks at the desired focus locations. The initial and final structures are shown in figures 4(e) and (f). We begin from an ordered 3D arrangement of metacubes and the optimisation procedure introduces disorder to achieve the desired functionality. The validity of the multi-functionality condition is shown in figure 4(g). Here, the numerical error is larger due to the highly oscillatory nature of the integrands.

Conclusions & outlook
Beginning from general considerations of vector fields, we have derived a condition that the fields must satisfy if they are to be supported by the same material distribution. Interestingly, this can be expressed as a surface integral so only the fields on a boundary are needed to determine whether a particular multi-functional device is feasible. While interesting, this condition is difficult to solve in general so instead we present an efficient and versatile semi-analytic method for designing multi-functional metamaterials made from discrete scattering elements. Demonstrating the generality of our formulation, we apply our method to design multi-functional devices at both optical and microwave wavelengths. We show structures that: (i) enhance the efficiency of an emitter while shaping its radiation pattern; (ii) beam in different directions based on the source polarisation; and (iii) sort waves by their incidence direction. In addition to the design method being simple, the structures we design are easy to fabricate as permittivity does not need to be graded. Instead many identical resonators must be distributed in space. Our approach could be utilised to design very wide classes of multi-functional devices. For example, the inclusion of higher-order multipoles could allow more complex resonators to be used, proving more degrees of freedom. Due to the generality of the two central ideas behind our framework: expanding figures of merit analytically to avoid expensive numerical derivatives and the formulation of the multi-objective problem in section 4, a wide range of devices could be designed using this approach. Graded index structures as well as propagation in waveguides and fibre optical cables could all be engineered to depend upon polarisation or wavelength for communications, cloaking or sensing applications using the methodology we propose here.

Conflict of interest
The authors declare no competing interests.

Author contributions
JRC and SARH conceived the idea. JRC derived the analytic expressions, performed the numerical simulations and wrote the manuscript. SARH, APH and SJB supervised the project. All authors commented on the manuscript.

Acknowledgments
JRC would like to thank Dean Patient for many useful discussions and for his careful proof-reading of the manuscript, as well as Alex Powell for providing the COMSOL model for the metacubes. We acknowledge financial support from the Engineering and Physical Sciences Research Council (EPSRC) of the United Kingdom, via the EPSRC Centre for Doctoral Training in Metamaterials (Grant No. EP/L015331/1). JRC also wishes to acknowledge financial support from Defence Science Technology Laboratory (DSTL). SARH acknowledges financial support from the Royal Society (URF\R\211033).

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.

Appendix A. Figure of merit expansions
In this section we give expansions of the figures of merit used in the main text and the gradients that result from these expansions. To write these, we will need to use the expressions for the changes in the fields that are given by (14) in the main text. Unpacking the notation, we have, δE(r) = ξ 2 G(r, r n )α E ∇E(r n ) + iξ∇ × G(r, r n )α H ∇H(r n ) δr n , ( A 1 )

A.1. Emitted power
Beginning from the usual expression for emitted power, expanding under small changes in the fields and then substituting (A1) gives the gradient of the power emission with respect to scatterer positions analytically,

A.2. Overlap integral
The figure of merit used to manipulate the shape of the Poynting vector in the far-field is the normalised overlap integral between the current angular distribution of the Poynting vector |S(θ)| and the target angular distribution ψ T (θ) .
The modulus of the Poynting vector can be expanded as Using these expressions to expand both the numerator and denominator in (A6), we find an expression for how the overlap integral changes when the fields change by a small amount, Substituting into this the expressions of the variations of the fields (A1) and (A2), the gradient of the overlap integral can be written analytically.

A.3. Modulus of Poynting vector in the far-field
When designing the beaming device, shown in figure 3, the figure of merit is the modulus of the Poynting vector in the far-field when the structure is excited by a source of polarisation p i . In this example, we have chosen the two circular polarisations. For each polarisation, we write the figure of merit as is, This means that for different source polarisations, power will be beamed into different directions. The expansion of this follows the same procedure as the expansion of the modulus of the Poynting vector did in the overlap integral, giving the result, Substituting into this the expressions for the field variations (A1) and (A2) gives analytic expressions for ∇ r n F i .

A.4. Modulus of electric field
For the 'lensing' problem, demonstrated in figure 4, we seek to increase the modulus of the electric field at particular positions r , F = |E(r )| 2 . (A13) Figure 5. A comparison of the results of our optimisation and a genetic algorithm seeking to shape a far-field radiation pattern while also improving efficiency. The far-field radiation patterns are compared in (a), and the solution space paths are shown in (b). The progress of our method is shown in red, and the progress of the genetic algorithm as green dots. Each dot represents a single population member. The final result of the genetic algorithm is shown as a blue star. The resulting structures are shown in (c) and (d). The positions are different for each of the incidence directions. This figure of merit can be easily expanded to find its gradient analytically, δF = 2 Re E * (r ) · δE(r ) , ∇ r n F = 2 Re E * (r ) · ξ 2 G(r , r n )α E ∇E(r n ) + iξ∇ × G(r , r n )α H ∇H(r n ) .

Appendix B. Comparison with a genetic algorithm
We compare the results of our optimisation for both power emission and directivity with the results of a genetic algorithm solving the same problem. Using the differential evolution algorithm [34], with a population size of 20, and a maximum allowed iterations of 5000. The differential weight parameter is F = 0.5 and the crossover probability is CR = 0.7. This genetic algorithm was run several times and the best solution selected. The comparison between this result and the result of our local optimisation is shown in figure 5. The genetic algorithm produces a slightly higher power emission but a slightly lower value of overlap integral. From the scatter of the solutions generated by the genetic algorithm in solution space, shown in figure 5, it is evident that the genetic algorithm explores more of the search space than our local optimisation. However, due to the size of the search space for multi-functional problems, this does not provide much advantage.