The local density of optical states of a metasurface

While metamaterials are often desirable for near-field functions, such as perfect lensing, or cloaking, they are often quantified by their response to plane waves from the far field. Here, we present a theoretical analysis of the local density of states near lattices of discrete magnetic scatterers, i.e., the response to near field excitation by a point source. Based on a pointdipole theory using Ewald summation and an array scanning method, we can swiftly and semi-analytically evaluate the local density of states (LDOS) for magnetoelectric point sources in front of an infinite two-dimensional (2D) lattice composed of arbitrary magnetoelectric dipole scatterers. The method takes into account radiation damping as well as all retarded electrodynamic interactions in a self-consistent manner. We show that a lattice of magnetic scatterers evidences characteristic Drexhage oscillations. However, the oscillations are phase shifted relative to the electrically scattering lattice consistent with the difference expected for reflection off homogeneous magnetic respectively electric mirrors. Furthermore, we identify in which source-surface separation regimes the metasurface may be treated as a homogeneous interface, and in which homogenization fails. A strong frequency and in-plane position dependence of the LDOS close to the lattice reveals coupling to guided modes supported by the lattice.

subdiffractive pitch, and composed of strong electric and magnetic scatterers as a homogenized electric or magnetic mirror. Furthermore, regimes in which the materials may suitably be treated as a homogeneous material are identified. We spectrally resolve the LDOS in regions in the lattice plane, revealing an increased LDOS by coupling to guided lattice modes.

Results
Theoretical framework. The optical response to plane wave excitation of 2D periodic lattices of electric polarizabilities has previously been reviewed by de Abajo 16 . An extension to the full magneto-electric case was presented in 7,17,18 . In the following we shall use results derived in 7 to which we refer the reader for further details. We consider a 2D periodic lattice of point scatterers in the dipole approximation positioned at R mn = md 1 + nd 2 , where m and n are integers, and d 1 and d 2 are the real space basis vectors. Previous work, based on finite difference time domain simulations 19 and quasistatic multipole theory 20 , has shown that the dipole approximation is warranted for d i > 3b, where b is the radius of the spheres. Each particle is described by a polarizability tensor, α ↔ , that relates the induced electric and magnetic dipole moment, p and m, to a driving electric and magnetic field E and H according to 7,21,22 ( ) For ease of notation we use a rationalized unit system as described in ref. 22 where e.g. |E|/|H| = 1 for a plane wave. We note that α ↔ is subject to symmetry constraints and must be made electrodynamically consistent, bound by the optical theorem. This is achieved by addition of radiation damping, , to the electrostatic polarizability α ↔ 0 which can for instance be derived from an LC model. Here . −1 denotes matrix inversion, k denotes the wave number,  is the 6-dimensional identity tensor 22 . The magnetoelectric static polarizability is decomposed as , describe bianisotropy, such as the electric response to magnetic fields (and vice versa).  ω ( ) is a Lorenzian prefactor, typical for a plasmon resonance, with resonance frequency ω 0 , Ohmic damping Γ and amplitude governed by the volume of the scatterer V. The induced dipole moment on a scatterer at the origin R 00 is set by the sum of the incident field and the field of all other dipoles in the lattice 7 where k || is the parallel momentum of the incident plane wave, is the 6 × 6 dyadic Green function of the medium surrounding the lattice. For our case, we shall assume the surrounding medium to be vacuum.
Calculating the LDOS in front of the lattice requires evaluating the scattered field arising from a single point source, instead of from a plane wave of definite parallel momentum. One approach would be to expand the field of the dipole in its parallel plane waves k || 23 . However, the resulting k || -integral unfortunately converges poorly, especially for small distances to the dipole source 24 . Instead we shall use a technique referred to as the array scanning method 24 . We consider a single point source dipole with a point current j at position r 0 , j(r) = − iωδ(r − r 0 )j. Here we use the notion of point current in a more generalized magnetoelectric context, where j is a 6-element vector describing both the electric and magnetic dipole, i.e. j = − iω(μ e , μ m ) T . We may synthesize this single point source by summing infinite phased arrays of point sources: where BZ denotes the Brillouin zone and  is the real-space unit cell area. We denote all quantities related to the phased array with a tilde. The incident field at the origin, generated by the phased array, is found by propagating the fields from each dipole in the phased array. We get that acts as a field propagator of an array of dipole source.
The induced dipole moment of the scatterer at the origin, driven by the phased array, is found using equation (4) and equation (8).  (8), we may evaluate the scattered field at a position r by multiplying the induced dipole with The scattered field from the original single dipole source is found by integrating the scattered field, generated by the phased array, over the entire Brillouin zone: Using equation (11), the decay rate, γ(r 0 ), of an emitter relative to the decay rate, γ vac (r 0 ), in vacuum is calculated as 23 : where † is the conjugate transpose and μ e µ ( ) m is the normalized electric (magnetic) dipole moment. In this work we will solely consider electric dipole transitions as source µ ( = ) 0 m . We note that the computation of the summation in equation (5) is carried out using Ewald summation 25 described in Supplementary material, and details of the integral in equation (12) is computed in practice are described in Methods. Moreover, while we only consider a single magneto-electric dipole mode of the scatterers, the model may easily be extended to treat stacked lattices as well as complex unit cells consisting of different scatterers, to mimic multipolar resonances 18 . Also, more advanced methods for retrieving the polarizability, e.g. surface integral equations 26 , may be used as input.

Numerical examples.
As examples we shall consider non-diffractive square lattices of strong scatterers. We calculate the LDOS near lattices of two types of scatterers: (1) Scatterers with an isotropic electric response (i.e. plasmonic spheres) by setting , and (2) Scatterers with an isotropic magnetic response by setting All parameter values used are presented in Table 1. The parameters are chosen so that the electric scatterers match the polarizability, extinction cross section and albedo found experimentally for plasmonic scatterers at telecom frequencies (extinction cross section 0.38 μ m 2 ), as studied in depth by Husnik et al. 27,28 .
The calculated LDOS modulation (plotted as predicted fluorescence lifetime normalized to lifetime in vacuum) as a function of distance is presented in Fig. 1, for an electric dipole source positioned at the four symmetry  points r || = (0, 0)d/2, r || = (1, 0)d/2, r || = (1, 1)d/2, and r || = (1, 1)d/2, oriented parallel to the lattice plane along x, Fig. 1(a), and perpendicular to the lattice plane along ẑ, Fig. 1(b). The relative lifetimes oscillate as a function of distance with a periodicity of about λ/2, as encountered in typical Drexhage-type experiments 2,3,29,30 . Comparing the electric versus magnetic lattices, we note that the oscillations in lifetime are π out of phase. A similar effect was predicted by Ruppin and Martin 14 for hypothetical 'magnetic mirrors' , i.e., for reflection at a medium that presents μ = − ∞, ε = 1, as opposed to ε = − ∞, μ = 1 for a normal electric mirror. In their work, the difference is associated with a π difference in Fresnel reflection coefficients that appears when interchanging magnetic permeability and electric permittivity. The calculated Drexhage oscillations, and their reversal in phase with exchanging the nature of the scatterers hence confirms that electric (magnetic) particle lattices act as effective electric (magnetic) reflective interfaces. Considering the case of an electrically scattering lattice, solid lines in Fig. 1, we note, that for distances beyond ~2d, or equivalently about λ/3, the lifetimes at the four different positions are indiscernible. Above this distance, the lattice is well approximated as an effective homogeneous material, as often assumed 11,14,15 . To qualify this statement further, we calculated the angle-dependent far field reflection coefficients (using equation (4) and equation (15) in ref. 17). These reflection constants can be used as input to textbook expressions for the LDOS near a homogeneous interface 23,29 , which for electric sources perpendicular, respectively parallel to an interface read Im Here ( − ) / k k k z 2 2 1 2 and r s,p represent the s-and p-reflection coefficient, and integrating up to k || = k accounts for all far-field reflection effects. We find excellent agreement for distances beyond a few lattice constants. This delineates the validity of using far field measurements to obtain effective material parameters. Furthermore, the notion of "effective material parameter" should be read as meaning that the medium is fully quantified by its far-field reflection for all angles, irrespective of the question if these reflection constants are consistent with any ε and μ.
For closer distances than ~2d, the discrete nature of the lattice is revealed in the position dependence of the decay. For all four positions, the lifetime rapidly decreases for short decreasing distances. Naturally, very close to a scattering sphere we expect a decrease associated with the near field of a single sphere. This should occur for ranges of order 50 nm (a/5) 31 .
For intermediate distances we identify a third effect, namely coupling to guided modes in the lattice 7 . To investigate contributions from guided modes we calculated the relative lifetime as a function of its emission frequency and in-plane position for a parallel and perpendicular dipole positioned in the plane of an electric isotropically scattering lattice, presented in Fig. 2(c,d). Firstly, we note, that for a perpendicular (parallel) dipole positioned in the plane of the lattice, all electric field components in the plane are perpendicular (parallel) to the lattice plane. Hence we expect coupling to modes with induced dipoles being purely perpendicular (parallel) to the plane. Considering the case of a parallel dipole (Fig. 2(d)) we firstly notice that close to the scattering element at (x, y) = (0, 0), the lifetime drastically decreases owing to the 1/r 3 scaling of the near field of the scatterer. Elsewhere, distinct bands of reduced lifetimes are resolved for frequencies different from the resonance frequency of the individual scatterer. This indicates that the source dipole couples not simply to the individual scattering elements, rather it couples to a guided lattice mode that is frequency dispersive. E.g. near (x, y) = (0, d/2), marked with a red circle, a significant reduction of the lifetime occurs for blue shifted frequencies relative to the single particle resonance frequency (ω 0 ). Symmetry of the lattice and the field lines of a dipole imply that the band arises from coupling to a longitudinal in-plane mode (LI) where the induced dipoles are arranged in a head to head configuration along x. This is confirmed from the calculated dispersion of the lattice mode with induced dipoles parallel to x, presented in Fig. 2(a) (for details on the calculation of the modal dispersion we refer to ref. 7). At points X ≡ k || = (π/d, 0) and M ≡ k || = (π/d, π/d), the mode is blueshifted with a flat slope thus giving rise to a large LDOS. Similarly, near (x, y) = (d/2, 0), marked with a triangle, a reduction is seen to occur for red shifted frequencies corresponding to a transverse in-plane mode with k || = Y ≡ (0, π/2). In the case of a dipole perpendicular to the lattice, (Fig. 2(d)), the calculated lifetime is symmetric about (0, 0) owing to the four-fold rotational symmetry of the lattice. Two bands appear near (0, d/2) and (d/2, 0) with one being slightly red shifted, the other blue shifted relative to the resonance frequency ω 0 . Comparing with the calculated mode with induced dipole momements perpendicular to the lattice plane, shown in Fig. 2(b), we conclude that the red shifted band is associated with coupling to a transverse guided mode with the induced dipoles perpendicular to the lattice, while the blue shifted resonance arise from coupling to a non-guided mode with wavevectors near above the light line. Coupling to this leaky mode is only achieved close to the lattice, since only in the near field of a radiating dipole does it contain wave vectors parallel to its dipole moment. Due to symmetry at (1, 1)d/2, marked by + , only coupling to the blue shifted non-guided mode with vanishing in-plane wavevectors near the point k || = Γ ≡ (0, 0), remains.

Discussion and Conclusion
In conclusion, we have presented a simple point dipole method using the array scanning method for calculating the LDOS of an arbitrary magnetoelectric infinite 2D lattice. The primary motivation to tackle this problem was to assess in how far analyzing a metamaterial as effectively homogeneous is reasonable in an actual scenario where it interacts with a localized object in its near field. As example, we calculated the lifetime of a dipole in front of electric and magnetic isotropically scattering spheres. We found that a lattice of magnetic scatterers shows characteristic oscillations of the LDOS as a function of distance, shifted in phase compared to those at an electric scattering lattice. This confirms that a metamaterial can appear as a magnetic mirror also in "Drexhage" experiments that are not limited to probing by a single far field incidence angle, as was first proposed by Ruppin and Martin 14 and Kästel and Fleischhauer 15 . Our results reveal that for distances beyond 2d ~ λ/3, the surfaces can be well approximated as an effective homogenous interface, with electric and magnetic properties taken from far field reflection constants. For somewhat shorter distances the lifetime shows a dependence on both in-plane position and frequency that is due to the discrete nature of the lattice, and coupling to lattice guided modes, which is not captured by far field reflection constants. At even shorter distances comparable to feature sizes of the scatterer, where microscopic detail matters, equation (7) of our work remains valid, however, the dipole approximation breaks down. Microscopically, one could use a full-wave solver (FDTD, COMSOL) for every wave vector in the integral in equation (7). In practice, however, this leads to an impractical computational burden. As an intermediate, and more tractable, approach we propose to improve the microscopic detail captured by our model by using multiple dipoles to describe a single scatterer, instead of using single dipoles 18 .
These results are of fundamental interest to the question how one probes the range of validity of effective medium parameters in near field geometries. Furthermore, our method is excellently suited for emitters with an excited state subject to competing radiative decay pathways with electric, magnetic, and mixed character 32,33 , where the calculated LDOS for the magnetic and electric transitions may be used as coefficients in the rate equations for the density of states of the emitter. Finally, our method can be easily extended to diffractive plasmonic systems, arbitrarily complex unit cells 18 , multilayered unit cells, and bi-anisotropic or hyperbolic metasurfaces.