Inverse design of light–matter interactions in macroscopic QED

Inverse design represents a paradigm shift in the development of nanophotonic devices, where optimal geometries and materials are discovered by an algorithm rather than symmetry considerations or intuition. Here we present a very general formulation of inverse design that is applicable to atomic interactions in external environments, and derive from this some explicit formulae for optimisation of spontaneous decay rates, Casimir–Polder forces and resonant energy transfer. Using Purcell enhancement of the latter as a simple example, we employ finite-difference time-domain techniques in a proof-of-principle demonstration of our formula, finding enhancement of the rate many orders of magnitude larger than a selection of traditional designs.


Introduction
Traditional design methods work by specifying a device, then investigating its properties. By contrast, in inverse design the desired property is specified, and an algorithm is left to find a device which fulfils the desired criteria. A naive approach to this would be simply trying all devices that fulfil some set of design constraints. The large space of possible designs renders this numerically unrealistic, meaning that a pre-determined set of designs must be optimised over, at least in the earliest applications of inverse methods to electromagnetic problems [1,2]. The development of adjoint methods [3] originally used in aerodynamics have made unconstrained inverse design computationally feasible, with the first application in photonics being to low-loss waveguide bends [4]. Adjoint methods were subsequently applied to band gaps [5], solar cells [6], on-chip demultiplexers [7] and many more diverse systems-see the recent review articles [8,9] and references therein.
An area in which inverse design has not yet been applied is virtual-photon mediated interactions, such as Casimir-Polder [10] forces and resonant energy transfer [11]. These processes can be described within a very general formalism known as macroscopic quantum electrodynamics (QED) [12], where they can be reduced to various functionals of the classical dyadic Green's tensor G for a source at r , observation point at r and frequency ω defined to satisfy ∇ × ∇ × G r, r , ω − ω 2 c 2 ε (r, ω) G r, r , ω = Iδ r − r (1) subject to given boundary conditions. Though G is a completely classical quantity, it is a vital ingredient for the description of quantised electromagnetic fields in media [12,13], and thereby quantum light-matter interactions in their vicinity. It appears, for example, in the coupling constants found in master equations governing weak light-matter interactions (e.g. those for atomic waveguide QED [14]), as well as in matrix elements in the strong coupling regime (e.g. references [15][16][17]). The quantised electric field in a region with permittivity ε (r, ω) and unit permeability is given via G as [12], whereb † ,b are a set of bosonic creation and annihilation operators for the medium-assisted quantised electromagnetic field 4 . The Green's tensor G takes into account both the geometry and material response of an arbitrarily-shaped medium, meaning that an optimal geometry for particular r, r and ω is represented by a particular functional form of G. It follows that G is the fundamental object which is to be worked with in inverse design of macroscopic QED. In this article we begin in section 2 by introducing the underlying formulae for inverse design of light-matter interactions. In section 3 use the specific example of resonant energy transfer combined with finite-difference time domain (FDTD) techniques to demonstrate that the efficiencies achievable in this method are far beyond those found from 'by-hand' constructions, opening up a new direction in the design and optimisation of a broad class of light-matter interactions. We conclude with some discussion and suggestions for further work in section 4.

General formulation
In order to carry out any optimisation, we need to define a merit function F which we intend to maximise. In traditional presentations of adjoint optimisation, this function is taken to depend on the electromagnetic fields E, D, B and H, but all of these are of course deducible from the dyadic Green's tensor so we consider F as being dependent on only G r, r , ω . The merit function should be an observable quantity, so we take it to be a real-valued functional of G r, r , ω , integrated over all its arguments: The integrals allow us to take into account a delocalised source and extended observation volume, as well as multi-mode effects. The entries of the tensor G are in general complex-valued, so in principle one could consider variations in the real and imaginary parts separately. However, it is more convenient to consider the complex tensors G and G * as independent, in which case the variation of the merit function with G is; where represents the Frobenius product (A B = i j A ij B ij ) and δG is a change in the Green's function brought about by an infinitesimal change in the environment. If this change can be considered as being confined to a small volume V containing a number density n(r ) of atoms with polarisabilities α(r ), as shown in appendix A we can write G in terms of the following Born series; where μ 0 is the vacuum permeability. The change in the merit function due to an additional small piece of dielectric material is then given by; where Lorentz reciprocity G(r, r , ω) = G T (r , r, ω) has been used. Merit functions for observables that depend on ∇ r G r, r , ω can be obtained via the replacements G r, r , ω → ∇ r G r, r , ω and G T r , r, ω → G T r , r, ω ← − ∇ r , with the most general form being a sum over all such derivative contributions.
Equation (6) is our main result, representing a formula of remarkable power and utility, ready for direct application a variety of processes. Before discussing some particular cases, there are several features of equation (6) worth commenting on. In traditional presentations of adjoint optimisation, the equivalent of equation (6) is expressed as the product of two electric fields. The first is the 'direct' field, which is simply the electric field induced by the sources present in the system. The second is the adjoint field, which is that generated by a dipole oscillator at the observation point with an amplitude given by the electric-field derivative of the merit function. The advantage of adjoint methods is that the optimal value of the merit function can be found with only two simulations (rather than a brute force method entailing placement of a dielectric inclusion at each possible point in the optimisation region and repeatedly simulating for each). This is reflected our version of the merit function change shown in (6); once the two independent Green's Table 1. Non-exhaustive list of observables f expressible in terms of G and their associated merit function gradients δF. In each case r A and r D are the positions of any atoms involved. The expression of the merit function change for spontaneous decay is equivalent to that used in reference [35], as should be expected.

Merit function Merit function integrand f change δF
Spontaneous tensors for a source at r and a source at r in a given environment (e.g. vacuum) have been calculated, δF is known at all points. The link with the adjoint electric field is simply that one of the Green's tensors in (6) has been transposed by taking advantage of Lorentz reciprocity. Equation (6) can be directly applied to any quantity that can be expressed in terms of the Green's dyadic G. This includes Casimir [18,19] and Casimir-Polder forces [10,20], spontaneous decay (Purcell factor) [21,22], quantum friction [23,24], interatomic Coulombic decay [25,26], radiative heat transfer [27,28], van der Waals forces [29], non-linear optical processes [30,31], environment-induced coherence [32,33], resonant energy transfer [11,34] and many more (the latter reference for each of these is where the formula in terms of G can be found). The merit functions for a selection of these are shown in table 1.
Given a merit function gradient δF, one has (at least) two choices for practical implementation of an optimisation-the simplest is an additive scheme illustrated in figure 1. Here a small block of material is added at the point of maximal δF, then the two Green's tensors in the new geometry are recalculated and combined to find a the next optimal point, an so on as indicated in figure 1. It is worth noting that the size and shape of the added blocks can be adjusted to take into account constraints from a given additive manufacturing process.
The second way to implement the optimisation consists of gradually optimising the shape of an initial object by changing its boundary, known as the level-set method [36]. Here, the initial boundary shape (as well as its subsequent evolution) is encoded by a chosen function φ. This is defined as negative inside the boundary, zero on it and positive outside, as indicated in figure 1. Introducing a 'time' parameter t representing iteration, one is led to the following equation of motion governing the shape of the boundary [36]:φ(r(t), t) + v n |∇φ(r(t), t)| = 0 where v n is the velocity of motion normal to the surface. Formally, this is an advection equation which can be solved using techniques from fluid dynamics. Taking the volume V in (6) to be that defined by the function φ, we can let; V d 3 r → ∂V dA δx(r ) = ∂V dA v n δt, where the shape deformation has been assumed to be small, as can be ensured by a sufficiently small time step δt in the evolution process. If the integrand of the r integral in (6) is positive at each iteration, the value of the merit function will continually increase. Positivity of (6) can then be ensured by using a velocity such that ∂F = ∂V dA v 2 n δt, which means identifying; This velocity can be directly calculated for a given G, then inserted into the advection equation, after which φ is evolved for δt. This delivers a new φ, which defines a new geometry, for which the new G can be calculated and the process iterates. Complex manufacturing constraints such elimination of highly curved regions, gaps, or 'bridges' can be enforced within in the level set method, as discussed in detail in reference [37].

Example implementation
In order to succinctly demonstrate the use of (6), we make some simplifying assumptions. We assume that the dielectric additions are homogenous and sufficiently small that the integral over r can be approximated by the value at its centre s: In practice, quantities which depend on the field at a single frequency and at a single position are considerably more computationally tractable than their multi-frequency, bulk medium counterparts. Here we concentrate on a simple and universal phenomenon which is well-approximated by radiation of a single frequency interacting with a point-like atom-resonant energy transfer (RET). Extension to bulk media would not involve too much extra computational overhead since this would still only require two Green's tensors, but calculations of quantities dependent on a continuous spectrum (e.g. ground-state Casimir-Polder forces) would require re-calculation of G at very many frequencies.
We will work in the dipole approximation and aim to optimise the RET rate Γ between dipole moments d A and d D , meaning we take [34]: We then have simply; which is the well-known expression of the resonance energy transfer rate Γ. Using this in (8), after some algebra one finds which is the equation we will work with from here on. To simplify our presentation of the main features of the method we restrict ourselves to systems with translational invariance along one axis, in other words those which can be considered as effectively two-dimensional. In order to validate the two-dimensional RET results that we will calculate (as well as the general FDTD approach), it is necessary to have an analytic expression for RET in two dimensions. Formally, 2D-RET is equivalent to taking a pair of 'line dipoles' each consisting of two infinitely extended parallel oppositely-charged wires in three dimensions, as discussed in detail in [38]. The Green's tensor from [38] can be directly substituted into (10), resulting in a lengthy expression, which can be simplified by noting that in situations of practical interest the dipoles are often randomly oriented necessitating an isotropic average (see appendix B), which gives; where H (1) n and H (2) n are Hankel functions of the first and second kind respectively, ζ = ω D ρ/c, ρ = |r A − r D |, d A = |d A | and d D = |d D |. The Green's tensor from [38] can also be substituted into equation (11) to give an initial indication of where material should be placed in order to optimise the RET process, this is shown in figure 2 Equation (12) can be used to validate our general numerical approach, in which we used the free FDTD software Meep [39] to calculate the Green's tensors entering (10) and (11). In order to do this we note that the ij component of the Green's tensor G r, r , ω describes the ith component of an electric field at observed at r due to the j component of a current source at r 0 ω d 3 r G r, r , ω · j(r , ω).
(13) Figure 2. Merit function gradient (11) using the Green's tensor from [38] and dipoles aligned along the x axis with donor transition wavelength of π μm. Only the functional form and sign (not the magnitude) of the merit function gradient are of relevance to an optimisation, so the values here are shown normalised to the magnitude of δF at the midpoint between the two dipoles. For a point source at r s j(r , ω) = j(ω)δ(r − r s ) (14) so equation (13) becomes E (r, ω) = iμ 0 ωG(r, r s , ω) · j(ω) (15) where j(ω) ≡ j(r s , ω) is the source current in the frequency domain. The FDTD implementation found in Meep (or in other FDTD tools) can calculate the electric field of a time-domain source j(t), which gives us everything we need to deduce G. To do this, we introduce a current source in a similar way to reference [40], namely a short Gaussian pulse (the built-in Meep function Gaussian source) polarised in the j direction. The corresponding time-domain current is given by; where A is an arbitrary amplitude [appearing on both sides of equation (15)], t 0 is the time at which the maximum is reached and w is the temporal width of the Gaussian. The Fourier transform of this is;  The time-domain simulation is allowed to run for long enough that all fields have decayed away, which is taken here as being 100 times the temporal width of the Gaussian (no change is observed in the results when the simulation time is varied, as long as it is more than roughly 10 times with width of the Gaussian). The result is a set of time-domain electric fields, which are then Fourier-transformed. Dividing these transformed fields component-wise by iμ 0 ωj(ω), one row of the Green's tensor is obtained, corresponding to a particular source polarisation direction j. This can be done in parallel for all three source polarisations, giving all nine components of G. In most situations the Green's tensor is symmetric so that only six of these components are independent, allowing further increases in efficiency. The RET rate obtained by calculating the FDTD Green's tensor in this way and substituting it into equation (10) is compared with the analytic expression (12) in figure 3. Close agreement is found until separations are comparable to pixel size, as is to be expected in any numerical scheme. We can now calculate the effect of arbitrary 2D geometries on RET. To do this we examine an analogue of the well-known Purcell factor for spontaneous decay [21], by investigating a Purcell-type factor F p = Γ/Γ 0 , where Γ 0 is the RET rate in vacuum. Broadly analogous quantities for laser cavities have previously been optimised via algorithmic design (e.g. Q factors in [41]), but without coupling to atoms and without the power and flexibility of the Green's tensor approach. We begin by choosing some geometries which are expected to enhance RET, these are shown in figure 4, giving a maximum F p in the hundreds. Iterative optimisation techniques can dramatically improve on this performance. In order to demonstrate this, we use the additive approach shown in figure 1, the results of which are shown in figure 5. An extremely large enhancement is found, reaching a factor of approximately 10 5 after 250 iterations-orders of magnitude higher than any enhancement found in the traditional designs shown in figure 4. It is worth noting that this extraordinarily high enhancement is achieved with a much smaller amount of dielectric material than in the traditional designs, even though the dielectric constant is identical. While a full eigenmode analysis of the structures shown in figure 5 is not the focus of the current work, some sense can be made of the apparently chaotic-looking designs by noting the waveguide-like structures in the regions into which dipole radiation is strongest (perpendicularly to the dipole moment). These consist of thin lines of dielectric material spaced approximately half a wavelength apart, acting in some sense like highly-reflective Bragg mirrors for the excitations travelling from donor to acceptor. Finally we note that while all our results are scale-invariant (depending only on the ratio of the transition wavelength to the interatomic separation), the minimum feature size in our simulations is 0.1 μm-this is broadly consistent with the level to which complex structures can be manufactured (see, e.g. [43]).

Conclusions
In this article we have presented a convenient and system-agnostic version of adjoint optimisation of electromagnetism based entirely on the electromagnetic dyadic Green's tensor. This allows the techniques of inverse design to be applied to any of the vast number of interactions and processes which can be expressed in terms of this tensor [19,20,22,24,26,28,29,31,33]. As a simple example we chose resonant energy transfer in two dimensions, showing orders of magnitude improvement in engineering potential compared to hand-made designs. The methods of inverse design presented here could also be used to investigate exceptional points in parameter space found by minimising some distance functional. One major extension of our work would be to non-reciprocal media, which would be taken into account by simply by leaving the Green's tensors in equation (6) untransposed, although the resulting simulations would no longer necessarily be able to take advantage of the speedup provided by the adjoint method. Further extensions of our work will include exploitation of the general level-set optimisation equation presented here, three-dimensional simulations and consideration of other observables including for example quantum yield of fluorescence processes and environment-induced coherence. and ∇ × ∇ × G r, r , ω − ω 2 c 2 ε (r, ω) G r, r , ω = Iδ r − r . If the background medium is vacuum, then the permittivity change χ (r, ω) is which means it is equal to the susceptibility, and therefore can be converted to a polarisability α using the Clausius-Mosotti law in the dilute limit [45]: which is equation (5) in the main text.