Designing the collective non-local responses of metasurfaces

Abstract

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.

Introduction

Designing 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⁷, 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⁹ or electromagnetic¹⁰. 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¹², 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¹⁵, including as electromagnetic filters and as perfect absorbers¹⁶.

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¹⁷. This has led to the development of metalenses, which can have superior bandwidth to traditional refractive lenses^18,19, as well as metasurface antennas²⁰. Metasurface antennas have been designed to exhibit many valuable properties, such as the opportunity to engineer bespoke beam-shaping, steering, polarisation control and improved efficiency²⁰. 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²¹, and for wavefront shaping²². 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²¹, metalenses¹⁸ and beam-shaping surfaces²⁰ have all been designed using the Gercherg–Saxton algorithm²³. 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²⁴. 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²⁷. 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²⁸. 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²⁹. As well as the Gerchberg–Saxton, there are two other popular inverse design paradigms. Firstly, geometry optimisation based upon the adjoint design method³⁰ 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³⁰. Reciprocity³³ 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³⁴ and the locally periodic approximation³⁵, which assumes a very sub-wavelength unit cell, are often exploited. Secondly, machine learning³⁶ and genetic algorithms³⁷ 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³⁸, which is often employed when it is assumed that multiple scattering can be neglected³⁹. 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⁴². 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₀⁴⁵, where r₀ 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 μ,

$${{{{{{{\boldsymbol{\nabla }}}}}}}}\times {{{{{{{\boldsymbol{\nabla }}}}}}}}\times {{{{{{{\boldsymbol{E}}}}}}}}-{k}^{2}{{{{{{{\boldsymbol{E}}}}}}}}={{{{{\rm{i}}}}}}\omega \mu {{{{{{{\boldsymbol{J}}}}}}}}+{\omega }^{2}\mu {{{{{{{\boldsymbol{P}}}}}}}}+{{{{{\rm{i}}}}}}\omega \mu {{{{{{{\boldsymbol{\nabla }}}}}}}}\times {{{{{{{\boldsymbol{M}}}}}}}},$$

(1)

$${{{{{{{\boldsymbol{\nabla }}}}}}}}\times {{{{{{{\boldsymbol{\nabla }}}}}}}}\times {{{{{{{\boldsymbol{H}}}}}}}}-{k}^{2}{{{{{{{\boldsymbol{H}}}}}}}}={{{{{{{\boldsymbol{\nabla }}}}}}}}\times {{{{{{{\boldsymbol{J}}}}}}}}+{k}^{2}{{{{{{{\boldsymbol{M}}}}}}}}-{{{{{\rm{i}}}}}}\omega {{{{{{{\boldsymbol{\nabla }}}}}}}}\times {{{{{{{\boldsymbol{P}}}}}}}},$$

(2)

where $k=\omega \sqrt{\varepsilon \mu }$ 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

$${{{{{{{\boldsymbol{J}}}}}}}}=-{{{{{\rm{i}}}}}}\omega \hat{{{{{{{{\boldsymbol{p}}}}}}}}}\delta ({{{{{{{\boldsymbol{r}}}}}}}}-{{{{{{{\boldsymbol{r}}}}}}}}^{\prime} ),\left(\begin{array}{l}{{{{{{{\boldsymbol{P}}}}}}}}\\ {{{{{{{\boldsymbol{M}}}}}}}}\end{array}\right)={\mathop{\alpha}\limits^\leftrightarrow}\left(\begin{array}{l}{{{{{{{\boldsymbol{E}}}}}}}}\\ {{{{{{{\boldsymbol{H}}}}}}}}\end{array}\right)\delta ({{{{{{{\boldsymbol{r}}}}}}}}-{{{{{{{{\boldsymbol{r}}}}}}}}}_{n}),$$

(3)

where $\hat{{{{{{{{\boldsymbol{p}}}}}}}}}$ is the polarisation of the source, E and H are the fields applied to the scatterer and

$${\mathop{\alpha}\limits^\leftrightarrow}=\left(\begin{array}{ll}{\mathop{\alpha}\limits^\leftrightarrow}_{E}&0\\ 0&{\mathop{\alpha}\limits^\leftrightarrow}_{H}\end{array}\right),$$

(4)

is the polarizability tensor, in the absence of bianisotropy. The location of the source is ${{{{{{{\boldsymbol{r}}}}}}}}^{\prime}$ 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: ${{{{{{{\rm{Im}}}}}}}}[{\mathop{\alpha}\limits^\leftrightarrow}_{E,H}^{-1}]\ge -{\mathbb{I}}{k}^{3}/(6\pi \varepsilon )$^46,47,48, where ${\mathbb{I}}$ is the unit tensor.

In general, any differential equation of the form

$${{{{{{{\boldsymbol{\nabla }}}}}}}}\times {{{{{{{\boldsymbol{\nabla }}}}}}}}\times {\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}({{{{{{{\boldsymbol{r}}}}}}}},{{{{{{{\boldsymbol{r}}}}}}}}^{\prime} )-{k}^{2}{\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}({{{{{{{\boldsymbol{r}}}}}}}},{{{{{{{\boldsymbol{r}}}}}}}}^{\prime} )=\delta ({{{{{{{\boldsymbol{r}}}}}}}}-{{{{{{{\boldsymbol{r}}}}}}}}^{\prime} ){\mathbb{I}},$$

(5)

may be solved in terms of the dyadic Green function^49,50

$${\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}({{{{{{{\boldsymbol{r}}}}}}}},{{{{{{{\boldsymbol{r}}}}}}}}^{\prime} )=\left({\mathbb{I}}+\frac{1}{{k}^{2}}{{{{{{{\boldsymbol{\nabla }}}}}}}}\otimes {{{{{{{\boldsymbol{\nabla }}}}}}}}\right)\frac{{{{{{{\rm{e}}}}}}}^{\pm {{{{{\rm{i}}}}}}k| {{{{{{{\boldsymbol{r}}}}}}}}-{{{{{{{\boldsymbol{r}}}}}}}}^{\prime} | }}{4\pi | {{{{{{{\boldsymbol{r}}}}}}}}-{{{{{{{\boldsymbol{r}}}}}}}}^{\prime} | },$$

(6)

where ${\mathbb{I}}={{{{{{{\rm{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 ${{\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}}_{EH}({{{{{{{\boldsymbol{r}}}}}}}},{{{{{{{\boldsymbol{r}}}}}}}}^{\prime} )={{{{{{{\boldsymbol{\nabla }}}}}}}}\times {\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}({{{{{{{\boldsymbol{r}}}}}}}},{{{{{{{\boldsymbol{r}}}}}}}}^{\prime} )$⁵¹.

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₀ = 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

$$\left(\begin{array}{l}{{{{{{{\boldsymbol{E}}}}}}}}({{{{{{{\boldsymbol{r}}}}}}}})\\ {{{{{{{\boldsymbol{H}}}}}}}}({{{{{{{\boldsymbol{r}}}}}}}})\end{array}\right)=\left(\begin{array}{l}{{{{{{{{\boldsymbol{E}}}}}}}}}_{s}({{{{{{{\boldsymbol{r}}}}}}}})\\ {{{{{{{{\boldsymbol{H}}}}}}}}}_{s}({{{{{{{\boldsymbol{r}}}}}}}})\end{array}\right)+\mathop{\sum}\limits_{n}\left(\begin{array}{ll}{\xi }^{2}{\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}({{{{{{{\boldsymbol{r}}}}}}}},{{{{{{{{\boldsymbol{r}}}}}}}}}_{n}){\mathop{\alpha}\limits^\leftrightarrow}_{E}&{{{{{\rm{i}}}}}}\xi {{\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}}_{EH}({{{{{{{\boldsymbol{r}}}}}}}},{{{{{{{{\boldsymbol{r}}}}}}}}}_{n}){\mathop{\alpha}\limits^\leftrightarrow}_{H}\\ -{{{{{\rm{i}}}}}}\xi {{\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}}_{EH}({{{{{{{\boldsymbol{r}}}}}}}},{{{{{{{{\boldsymbol{r}}}}}}}}}_{n}){\mathop{\alpha}\limits^\leftrightarrow}_{E}&{\xi }^{2}{\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}({{{{{{{\boldsymbol{r}}}}}}}},{{{{{{{{\boldsymbol{r}}}}}}}}}_{n}){\mathop{\alpha}\limits^\leftrightarrow}_{H}\end{array}\right)\left(\begin{array}{l}{{{{{{{\boldsymbol{E}}}}}}}}({{{{{{{{\boldsymbol{r}}}}}}}}}_{n})\\ {{{{{{{\boldsymbol{H}}}}}}}}({{{{{{{{\boldsymbol{r}}}}}}}}}_{n})\end{array}\right),$$

(7)

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:

$$\left(\begin{array}{l}{{{{{{{\boldsymbol{E}}}}}}}}({{{{{{{{\boldsymbol{r}}}}}}}}}_{1})\\ {{{{{{{\boldsymbol{H}}}}}}}}({{{{{{{{\boldsymbol{r}}}}}}}}}_{1})\\ {{{{{{{\boldsymbol{E}}}}}}}}({{{{{{{{\boldsymbol{r}}}}}}}}}_{2})\\ {{{{{{{\boldsymbol{H}}}}}}}}({{{{{{{{\boldsymbol{r}}}}}}}}}_{2})\\ \vdots \end{array}\right)={{{{{{{{\boldsymbol{R}}}}}}}}}^{-1}\left(\begin{array}{l}{{{{{{{{\boldsymbol{E}}}}}}}}}_{s}({{{{{{{{\boldsymbol{r}}}}}}}}}_{1})\\ {{{{{{{{\boldsymbol{H}}}}}}}}}_{s}({{{{{{{{\boldsymbol{r}}}}}}}}}_{1})\\ {{{{{{{{\boldsymbol{E}}}}}}}}}_{s}({{{{{{{{\boldsymbol{r}}}}}}}}}_{2})\\ {{{{{{{{\boldsymbol{H}}}}}}}}}_{s}({{{{{{{{\boldsymbol{r}}}}}}}}}_{2})\\ \vdots \end{array}\right),$$

(8)

where ⁻¹ denotes matrix inversion, and R is the response matrix, defined as

$${{{{{{{\boldsymbol{R}}}}}}}}=\left(\begin{array}{lll}{{{{{{{{\boldsymbol{R}}}}}}}}}_{11}&{{{{{{{{\boldsymbol{R}}}}}}}}}_{12}&\cdots \\ {{{{{{{{\boldsymbol{R}}}}}}}}}_{21}&{{{{{{{{\boldsymbol{R}}}}}}}}}_{22}&\cdots \\ \vdots &\vdots &\ddots \end{array}\right),{{{{{{{{\boldsymbol{R}}}}}}}}}_{ij}=\left(\begin{array}{ll}{\mathbb{I}}{\delta }_{ij}-{\xi }^{2}{\mathop{\alpha}\limits^\leftrightarrow}_{E}{\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}({{{{{{{{\boldsymbol{r}}}}}}}}}_{i},{{{{{{{{\boldsymbol{r}}}}}}}}}_{j})&-{{{{{\rm{i}}}}}}\xi {\mathop{\alpha}\limits^\leftrightarrow}_{H}{\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}_{EH}({{{{{{{{\boldsymbol{r}}}}}}}}}_{i},{{{{{{{{\boldsymbol{r}}}}}}}}}_{j})\\ {{{{{\rm{i}}}}}}\xi {\mathop{\alpha}\limits^\leftrightarrow}_{E}{\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}_{EH}({{{{{{{{\boldsymbol{r}}}}}}}}}_{i},{{{{{{{{\boldsymbol{r}}}}}}}}}_{j})&{\mathbb{I}}{\delta }_{ij}-{\xi }^{2}{\mathop{\alpha}\limits^\leftrightarrow}_{H}{\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}({{{{{{{{\boldsymbol{r}}}}}}}}}_{i},{{{{{{{{\boldsymbol{r}}}}}}}}}_{j})\end{array}\right).$$

(9)

It should be noted that as ${\mathop{\alpha}\limits^\leftrightarrow}$ is defined as the response to the applied fields, any self–interaction terms of the form ${\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}({{{{{{{{\boldsymbol{r}}}}}}}}}_{i},{{{{{{{{\boldsymbol{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 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 ${{{{{{{\boldsymbol{R}}}}}}}}\in {{\mathbb{C}}}^{3\times 2\times 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⁵². 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⁵³

$$\rho (\hat{{{{{{{{\bf{p}}}}}}}}},{{{{{{{\bf{r}}}}}}}},\omega )=\frac{2{\epsilon }_{0}{n}^{2}}{\pi \omega }{{{{{{{\rm{Im}}}}}}}}[\hat{{{{{{{{\bf{p}}}}}}}}}\cdot {{{{{{{\bf{E}}}}}}}}({{{{{{{\bf{r}}}}}}}})].$$

(10)

This quantity gives the number of electromagnetic modes per unit volume, for a given polarisation $\hat{{{{{{{{\bf{p}}}}}}}}}$, position r and frequency ω. The PLDoS $\rho (\hat{{{{{{{{\bf{p}}}}}}}}},{{{{{{{\bf{r}}}}}}}},\omega )$ 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 (${\mathop{\alpha}\limits^\leftrightarrow}_{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

$$\left(\begin{array}{l}{{{{{\bf{p}}}}}}({{{{{{\bf{r}}}}}}}_{1})\\ {{{{{\bf{p}}}}}}({{{{{{\bf{r}}}}}}}_{2})\end{array}\right)={\mathop{\alpha}\limits^\leftrightarrow}_{E}\cdot {\left(\begin{array}{rc}{\mathbb{I}}&-{\xi }^{2}{\mathop{{{{{{\bf{G}}}}}}}\limits^\leftrightarrow}({{{{{{\bf{r}}}}}}}_{1},{{{{{{\bf{r}}}}}}}_{2})\cdot {\mathop{\alpha}\limits^\leftrightarrow}_{E}\\ -{\xi }^{2}{\mathop{{{{{{\bf{G}}}}}}}\limits^\leftrightarrow}({{{{{{\bf{r}}}}}}}_{2},{{{{{{\bf{r}}}}}}}_{1})\cdot {\mathop{\alpha}\limits^\leftrightarrow}_{E}&{\mathbb{I}}\end{array}\right)}^{-1}\cdot \left(\begin{array}{l}{{{{{{\bf{E}}}}}}}_{s}({{{{{{\bf{r}}}}}}}_{1})\\ {{{{{{\bf{E}}}}}}}_{s}({{{{{{\bf{r}}}}}}}_{2})\end{array}\right).$$

(11)

The eigenvalues of R⁻¹ satisfy the characteristic polynomial

$${{{{{{{{\boldsymbol{\lambda }}}}}}}}}^{-1}={\mathbb{I}}\pm \frac{1}{2}\sqrt{4\left[{\xi }^{2}{\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}({{{{{{{{\bf{r}}}}}}}}}_{1},{{{{{{{{\bf{r}}}}}}}}}_{2})\cdot {\mathop{\alpha}\limits^\leftrightarrow}_{E}\right]\cdot \left[{\xi }^{2}{\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}({{{{{{{{\bf{r}}}}}}}}}_{2},{{{{{{{{\bf{r}}}}}}}}}_{1})\cdot {\mathop{\alpha}\limits^\leftrightarrow}_{E}\right]},$$

(12)

which can be solved by making use of the reciprocity of the Green function⁴⁹ to find that the eigenvalues of R are

$${{{{{{{\boldsymbol{\lambda }}}}}}}}={{{{{{{\rm{diag}}}}}}}}\left\{\frac{{\mathop{\alpha}\limits^\leftrightarrow}_{E}}{{\mathbb{I}}\pm {\xi }^{2}{\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}({{{{{{{{\bf{r}}}}}}}}}_{1},{{{{{{{{\bf{r}}}}}}}}}_{2})\cdot {\mathop{\alpha}\limits^\leftrightarrow}_{E}}\right\}.$$

(13)

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₁ to r₂ 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, ${\mathop{\alpha}\limits^\leftrightarrow}_{E}^{-1}\cdot {{{{{{{\boldsymbol{\lambda }}}}}}}}$ gives the enhancement of the single-particle polarizability due to the multiple scattering events between the two scatterers. So, if ${\mathop{\alpha}\limits^\leftrightarrow}_{E}^{-1}\cdot {{{{{{{\boldsymbol{\lambda }}}}}}}}=1$ there is no change to the single-particle polarizability and multiple scattering events provide no enhancement. On the other hand, a large ${\mathop{\alpha}\limits^\leftrightarrow}_{E}^{-1}\cdot {{{{{{{\boldsymbol{\lambda }}}}}}}}$ 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 = λ_Rw, and left eigenvectors satisfy vR = λ_Lv. Using this left and right pair, the source field can be decomposed into the basis of these eigenvectors, w_n as

$${\left({{{{{{{{\bf{E}}}}}}}}}_{s},{{{{{{{{\bf{H}}}}}}}}}_{s}\right)}^{{{{{{\rm{T}}}}}}}=\mathop{\sum }\limits_{n=1}^{N}{c}_{n}{{{{{{{{\bf{w}}}}}}}}}_{n},$$

(14)

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,

$$ {{{{{{\bf{r}}}}}}}_{n}\to {{{{{{\bf{r}}}}}}}_{n}+{{\Delta }}{{{{{{\bf{r}}}}}}}_{n},\quad\delta ({{{{{\bf{r}}}}}}-{{{{{{\bf{r}}}}}}}_{n})\to \delta ({{{{{\bf{r}}}}}}-{{{{{{\bf{r}}}}}}}_{n})+{{\Delta }}{{{{{{\bf{r}}}}}}}_{n}\cdot \nabla \delta ({{{{{\bf{r}}}}}}-{{{{{{\bf{r}}}}}}}_{n}),\\ {{{{{\bf{E}}}}}}\to {{{{{{\bf{E}}}}}}}_{0}+{{{{{{\bf{E}}}}}}}_{1},\qquad\quad{{{{{\bf{H}}}}}}\to {{{{{{\bf{H}}}}}}}_{0}+{{{{{{\bf{H}}}}}}}_{1}.$$

(15)

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

$$\begin{array}{lll}{{{{{{{{\bf{E}}}}}}}}}_{1}({{{{{{{\bf{r}}}}}}}})&=&{\xi }^{2}\left[{\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}({{{{{{{\bf{r}}}}}}}},{{{{{{{{\bf{r}}}}}}}}}_{n}){\mathop{\alpha}\limits^\leftrightarrow}_{E}{{{{{{{{\bf{E}}}}}}}}}_{1}({{{{{{{{\bf{r}}}}}}}}}_{n})-{\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}({{{{{{{\bf{r}}}}}}}},{{{{{{{{\bf{r}}}}}}}}}_{n}){\mathop{\alpha}\limits^\leftrightarrow}_{E}{{\Delta }}{{{{{{{{\bf{r}}}}}}}}}_{n}\cdot \nabla {{{{{{{{\bf{E}}}}}}}}}_{0}({{{{{{{{\bf{r}}}}}}}}}_{n})\right]\hfill\\ &+&{{{{{\rm{i}}}}}}\xi \left[{{\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}}_{EH}({{{{{{{\bf{r}}}}}}}},{{{{{{{{\bf{r}}}}}}}}}_{n}){\mathop{\alpha}\limits^\leftrightarrow}_{H}{{{{{{{{\bf{H}}}}}}}}}_{1}({{{{{{{{\bf{r}}}}}}}}}_{n})-{{\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}}_{EH}({{{{{{{\bf{r}}}}}}}},{{{{{{{{\bf{r}}}}}}}}}_{n}){\mathop{\alpha}\limits^\leftrightarrow}_{H}{{\Delta }}{{{{{{{{\bf{r}}}}}}}}}_{n}\cdot \nabla {{{{{{{{\bf{H}}}}}}}}}_{0}({{{{{{{{\bf{r}}}}}}}}}_{n})\right],\hfill\end{array}$$

(16)

$$\begin{array}{lll}{{{{{{{{\bf{H}}}}}}}}}_{1}({{{{{{{\bf{r}}}}}}}})&=&{\xi }^{2}\left[{\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}({{{{{{{\bf{r}}}}}}}},{{{{{{{{\bf{r}}}}}}}}}_{n}){\mathop{\alpha}\limits^\leftrightarrow}_{H}{{{{{{{{\bf{H}}}}}}}}}_{1}({{{{{{{{\bf{r}}}}}}}}}_{n})-{\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}({{{{{{{\bf{r}}}}}}}},{{{{{{{{\bf{r}}}}}}}}}_{n}){\mathop{\alpha}\limits^\leftrightarrow}_{H}{{\Delta }}{{{{{{{{\bf{r}}}}}}}}}_{n}\cdot \nabla {{{{{{{{\bf{H}}}}}}}}}_{0}({{{{{{{{\bf{r}}}}}}}}}_{n})\right]\hfill\\ &-&{{{{{\rm{i}}}}}}\xi \left[{{\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}}_{EH}({{{{{{{\bf{r}}}}}}}},{{{{{{{{\bf{r}}}}}}}}}_{n}){\mathop{\alpha}\limits^\leftrightarrow}_{E}{{{{{{{{\bf{E}}}}}}}}}_{1}({{{{{{{{\bf{r}}}}}}}}}_{n})-{{\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}}_{EH}({{{{{{{\bf{r}}}}}}}},{{{{{{{{\bf{r}}}}}}}}}_{n}){\mathop{\alpha}\limits^\leftrightarrow}_{E}{{\Delta }}{{{{{{{{\bf{r}}}}}}}}}_{n}\cdot \nabla {{{{{{{{\bf{E}}}}}}}}}_{0}({{{{{{{{\bf{r}}}}}}}}}_{n})\right].\hfill\end{array}$$

(17)

If one wishes to increase the power emission of the dipole, the figure of merit is $\propto {{{{{{{\rm{Im}}}}}}}}[{{{{{{{{\bf{E}}}}}}}}}_{1}({{{{{{{\bf{r}}}}}}}}^{\prime} )]$. 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 $\propto {{{{{{{\rm{Im}}}}}}}}[{{{{{{{{\bf{E}}}}}}}}}_{1}({{{{{{{{\bf{r}}}}}}}}}_{{{{{{{{\rm{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

$${{{{{{{{\bf{r}}}}}}}}}_{n}^{i+1}={{{{{{{{\bf{r}}}}}}}}}_{n}^{i}+{{\Delta }}{{{{{{{{\bf{r}}}}}}}}}_{n}\times {{{{{{{\rm{sign}}}}}}}}\left[{{{{{{{\rm{Im}}}}}}}}\left\{{\xi }^{2}{\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}({{{{{{{\bf{r}}}}}}}},{{{{{{{{\bf{r}}}}}}}}}_{n})\cdot {\mathop{\alpha}\limits^\leftrightarrow}_{E}\cdot \nabla {{{{{{{{\bf{E}}}}}}}}}_{0}({{{{{{{{\bf{r}}}}}}}}}_{n})+{{{{{\rm{i}}}}}}\xi {{\mathop{{{{{{{{\bf{G}}}}}}}}}\limits^\leftrightarrow}}_{EH}({{{{{{{\bf{r}}}}}}}},{{{{{{{{\bf{r}}}}}}}}}_{n})\cdot {\mathop{\alpha}\limits^\leftrightarrow}_{H}\cdot \nabla {{{{{{{{\bf{H}}}}}}}}}_{0}({{{{{{{{\bf{r}}}}}}}}}_{n})\right\}\right].$$

(18)

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⁵⁶ 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.

Fig. 1: The system considered in our numerical examples.

A schematic of the problem is shown in (a). A dipole emitter with polarisation $\hat{{{{{{{{\bf{p}}}}}}}}}$ is located at r′. This generates source fields (E_s, H_s), which are scattered by N silicon spheres located at r_n. The scatters have optical properties shown in (b). The scatterers are assumed to be isotropic silicon spheres of radius 65 nm, with the frequency dispersion of the polarizability extracted from experimental data⁵⁷. A simple example of the design process is shown in (c). The power emission of the emitter (red arrow) is enhanced by iteratively moving passive scatterers (black dots). Colour plots show ${{{{{{{\rm{Im}}}}}}}}[{{{{{{{\bf{E}}}}}}}}({{{{{{{\bf{r}}}}}}}})\cdot \hat{{{{{{{{\bf{p}}}}}}}}}]$, which evaluated at the emitter’s location is proportional to power emission. The enhancement of this quantity at the location of the emitter is evident.

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 polarisation $\hat{{{{{{{{\bf{p}}}}}}}}}$ is located at ${{{{{{{\bf{r}}}}}}}}^{\prime}$. 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⁵⁷, then combined with the Mie a₁ and b₁ coefficients⁵⁸, to allow the calculation of electric and magnetic polarizabilities⁵⁹. 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.

Table 1 Parameters used for all numerical calculations.

It has been assumed that the scatterers are isotropic, so the electric and magnetic polarizability tensors have the form

$${\mathop{\alpha}\limits^\leftrightarrow}_{E}=\left(\begin{array}{lll}1&0&0\\ 0&1&0\\ 0&0&1\end{array}\right){\alpha }_{E},{\mathop{\alpha}\limits^\leftrightarrow}_{H}=\left(\begin{array}{lll}1&0&0\\ 0&1&0\\ 0&0&1\end{array}\right){\alpha }_{H}.$$

(19)

The complex numbers α_E and α_H are calculated according to⁵⁸ as

$${\alpha }_{E}={{{{{{{\rm{i}}}}}}}}\frac{6\pi }{{k}^{3}}{a}_{1},{\alpha }_{H}={{{{{{{\rm{i}}}}}}}}\frac{6\pi }{{k}^{3}}{b}_{1},$$

(20)

where a₁ and b₁ 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 anti-symmetric 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.

Fig. 2: A more detailed look at a simple example optimisation, to introduce the tools that will be used to characterise more complex structures.

The aim is to enhance the power emission of the dipole emitter at the origin (black arrow). We consider two scatterers, labelled A and B, that support only an electric dipole (α_H = 0). (a) shows the configuration change during the optimisation. (b) shows the evolution of the power emission. (c) Shows how the optimisation process changes the relative dominance of the modes in the response of the system, as characterised by their expansion coefficient. (d) Shows how the eigenmodes of the system, with their associated eigenvalues, are modified by the design procedure. Note the two classes of mode: anti-symmetric (i.e. initial mode 1) shown in blue and symmetric (i.e. initial mode 2) shown in red. Eigenmodes are indicated as pairs of 3D vectors, representing the components of the electric field in each scatterer.

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³², but is of the order of studies with similar assumptions and constraints but which use genetic algorithms³⁷. Examining the change in the eigenmode with the 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.

Fig. 3: The result of applying our design methodology to enhance power emission of a dipole using 100 scatterers.

In all plots, scatterers are shown as black circles and the dipole emitter as a black arrow. (a) Shows the $\hat{{{{{{{{\bf{y}}}}}}}}}$ component of the electric field in the initial configuration, (b) shows the progress of the power enhancement as the optimisation progresses and (c) shows E_y in the optimised configuration. (d) shows how the eigenvalues and expansion coefficients change due to the optimisation. (e) and (f) show ${{{{{{{\rm{Re}}}}}}}}[{{{{{{{\bf{E}}}}}}}}\cdot \hat{{{{{{{{\bf{y}}}}}}}}}]$ of the mode with the largest expansion coefficient in the initial and optimised structure respectively. The final mode, plotted in (f), with eigenvalue ~8 and expansion coefficient ~0.5 is responsible for the power enhancement.

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⁶², 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⁶³ as well as simple phase arguments³². 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¹⁹. 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°.

Fig. 4: Re-structuring the far-field of a dipole emitter to radiated power in a particular direction.

(a) shows ∣E∣ for the initial configuration (regular array), viewed from the far-field. (b) shows a comparison of the far-field distribution of the Poynting vector for the initial configuration (black dashed line) and the optimised configuration (red line). All far-field distributions are normalised to have unit amplitude. The width of the beam in the optimised structure is ~20°. (c) shows ∣E∣ of the optimised structure, viewed from the far-field. The initial and final scatterer locations are shown in (d), and the change in eigenmodes and their expansion coefficients over the optimisation are shown in (e). The evolution of the properties of the eigenmode with the largest expansion coefficient is shown in (f–i). Both the distribution of E_y in the plane of the scatterers and the normalised far-field Poynting vector are shown. (f–g) show the mode with the largest expansion coefficient in the initial strucutre, with (h–i) showing the mode with the largest expansion coefficient in the final structure. It is clear that the mode shown in (h–i) is responsible for the strong directivity along θ = 0° in the optimised structure.

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⁶⁴. 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.

Fig. 5: Multiplexing the design shown in Fig. 4 to construct a Luneburg-like lens.

(a) demonstrates the function of a normal Luneburg lens of radius R, with a refractive index graded according to the inset equation. A point source is converted into a beam in a single direction. b, c The result of multiplexing the structure proposed in 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).

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.