Designing Multi-functional Metamaterials

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 structure 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 design a device that operates at optical wavelengths and beams light into different directions depending on the source polarisation and a device that works at microwave wavelengths and sorts plane waves by their angle of incidence.


I. INTRODUCTION
Much of 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 permittivty ε(r). This wave obeys the vector Helmholtz equation, ∇ × ∇ × φ(r) + k 2 0 ε(r) φ(r) = 0.
From this statement, we can find a condition upon the two wave-fields for this to be possible.
Integrating this over all space, then using Green's vector identity [13], we find that ∂V φ 2 (r) × ∇ × φ 1 (r) − φ 1 (r) × ∇ × φ 2 (r) · dS = 0. (4) This places a stringent condition 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]. Eq. (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 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 where η 0 is the impedance of free space, we find which is the usual expression for energy conservation expressed in terms of the Poynting vector S = 1 2 Re [E × H * ]. A direct application of Eq. (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 fullwave solver to find the fields then 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 electromagnetics problems. The discrete dipole approximation is briefly described in Section II then a summary of how metamaterials can be designed within this approximation is provided in in Section III. 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 IV. Several examples of utilising this framework are then shown in Section V. 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 polarisaion 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.

II. 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 then 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 distribution can be written 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 G(r, r ) = 1 + 1 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] R mn φ m = φ i,n , 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.

III. 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, with respect to the scatterer locations ∂F ∂rn 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 diving 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 explore in the next section.

IV. 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 emission of the emitter, The second is the overlap integral between the angular distribution of the Poynting vector in the far-field, |S(θ)| and the 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.
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.

V. 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 figure of merits are the power emission (21) and the overlap with the desired radiation pattern (22) and we choose the weights according to (24). Analytic expression for the gradients of these figures are merit can be found analytically. 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 2a, with the path of the optimisation in solution space shown in Figure 2b and the two resulting structures shown in Figure 2c 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 2b 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 3a shows the radiation patterns 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] that may be converted to 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 complicated scatterers, provided they can be approximated as dipolar, 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 C. 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 4a and for a plane wave at −20 • in Figure 4b. The two different focus points are clearly visible in the fields. The path in solution space is shown in Figure 4c, where again the choice of weighting has ensured roughly equal performance of F 1 and F 2 .
Slices of the fields from Figure 4a,b are shown in Figure 4d, indicating the large main peaks at the desired focus locations. The validity of the multi-functionality condition is shown in Figure 4e. Here, the numerical error is larger due to the highly oscillatory nature of the integrands. 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.

VI. 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 semianalytic 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 IV, 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.
supervised the project. All authors commented on the manuscript. where 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, A2), the gradient of the overlap integral can be written analytically.

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, A2) gives analytic expressions for ∇ rn F i .

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 , The positions are different for each of the incidence directions. This figure of merit can be easily expanded to find its gradient analytically,

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.

Appendix C: Validity of the Discrete Dipole Approximation
To verify the validity of the discrete dipole approximation, we compare our results with full-wave solutions using a finite element method numerical solver (COMSOL Multiphysics) [28].
The beaming device shown in Figure 3 has been validated by considering the nearest 20 scatterers to the source, due to memory considerations. For this reduced system, the comparison between the discrete dipole approximation and COMSOL is shown in Figure 6.
The scatterers here are silicon spheres of radius 65 nm, for which the electric and magnetic polarisabilities can be found analytically.
Validation of the 'lensing' device shown in Figure 4 of the main paper is shown in Figure   7. We consider a TE plane wave incident upon a small number of metacubes. Comparing the near and far fields in Figure 7 the main difference is in the field at the location of the