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.

Spontaneous emission is the irreversible decay of a quantum emitter upon emission of a photon that arises due to interaction with the fluctuating electromagnetic vacuum field. If the local density of available photon states (LDOS), and thereby of vacuum fluctuations, is modified the spontaneous emission rate may be inhibited or enhanced. This effect was discussed first for microcavities by Purcell in 19461. The effect was clearly demonstrated by Drexhage in 1966 in a fluorescence experiment using a rare earth ion placed in front of a mirror2. Over the past decades technological advances have made precise fabrication of nanostructered materials possible, allowing for tailoring the LDOS3456. Engineering the LDOS is attractive since it controls light-matter interaction such as thermal emission, absorption, and spontaneous emission.

Metamaterials and metasurfaces are nanostructured three- and two dimensional materials that aim to mimic homogeneous materials and interfaces, but with unconventional material properties. Especially metasurfaces and plasmonic lattices have recently attracted interest in the framework of spontaneous emission control due to their guiding properties, broad optical resonances and high field enhancements78910, the possibility of diverging LDOS in hyperbolic metamaterials11, and the aspect of controlling magnetic and chiral transitions1213. Already soon after the first metamaterials were made, Ruppin and Martin14 and Kästel and Fleischhauer15 analyzed the classical Drexhage experiment, but envisioning magnetic, and negative index continuous metamaterials to modulate the LDOS. This thought experiment is interesting for a few reasons: First, metamaterials are often proposed with near-field applications in mind, such as super-resolution imaging, or cloaking. Yet, whether a metamaterial medium built out of, for instance, magnetically polarizable scatterers, acts as a magnetic medium, is usually tested from the far field with just a single input wave vectors. Instead, a more comprehensive test would be to measure the local density of states, since Drexhage’s effect incorporates the different phase upon reflection1415, and sums over all wave vectors. Second, since an emitter is a point-like probe, approaching it to a metamaterial while measuring the lifetime is a direct method to probe at which source-material separations the effective medium approximation holds despite the inherently discrete geometry of metamaterials. These questions have, to the best of our knowledge, not been addressed previously.

To address these issues we present a semi-analytical point-dipole model, that allows for swift calculations of the LDOS of lattices of arbitrary electric, magnetic and bianisotropic dipolar scatterers. We utilize the method on two types of lattices consisting of isotropically scattering particles with an electric and a magnetic response, respectively. By comparing the Drexhage effect of the two lattices we explore the validity of treating a surface of 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 Abajo16. An extension to the full magneto-electric case was presented in71718. In the following we shall use results derived in7 to which we refer the reader for further details. We consider a 2D periodic lattice of point scatterers in the dipole approximation positioned at  = m1 + n2, where m and n are integers, and 1 and 2 are the real space basis vectors. Previous work, based on finite difference time domain simulations19 and quasistatic multipole theory20, has shown that the dipole approximation is warranted for d > 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 and according to72122

For ease of notation we use a rationalized unit system as described in ref.22 where e.g. ||/|| = 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 which can for instance be derived from an LC model. Here denotes matrix inversion, k denotes the wave number, is the 6-dimensional identity tensor22. The magnetoelectric static polarizability is decomposed as

where each matrix element is a 3 × 3 dimensionless matrix. The diagonals reflect a purely electric (magnetic) response, whereas the off-diagonal tensors , 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 00 is set by the sum of the incident field and the field of all other dipoles in the lattice7

where || is the parallel momentum of the incident plane wave,

and 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 ||23. However, the resulting ||-integral unfortunately converges poorly, especially for small distances to the dipole source24. Instead we shall use a technique referred to as the array scanning method24. We consider a single point source dipole with a point current at position 0, () = −iωδ( − 0). Here we use the notion of point current in a more generalized magnetoelectric context, where is a 6-element vector describing both the electric and magnetic dipole, i.e.  = −iω(, )T. We may synthesize this single point source by summing infinite phased arrays of point sources:

The original single source current is recovered from the phased array as

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

where we defined 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).

Similar to equation (8), we may evaluate the scattered field at a position by multiplying the induced dipole with giving

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:

where

Using equation (11), the decay rate, γ(0), of an emitter relative to the decay rate, γvac(0), in vacuum is calculated as23:

where † is the conjugate transpose and is the normalized electric (magnetic) dipole moment. In this work we will solely consider electric dipole transitions as source . We note that the computation of the summation in equation (5) is carried out using Ewald summation25 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 resonances18. Also, more advanced methods for retrieving the polarizability, e.g. surface integral equations26, 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 , and (2) Scatterers with an isotropic magnetic response by setting , and . 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 μm2), as studied in depth by Husnik et al.2728.

Table 1

Used parameters for the calculations.

Parameter	Value	Description
d1 = d2 ≡ d	300 nm	Lattice constant.
λ0 = 2πc/ω0	1.5 μm	Res. wavelength of particles.
Γ	83 THz	Ohmic damping in particles17.
V	(90 nm)3	Volume of scatterer.

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 || = (0, 0)d/2, || = (1, 0)d/2, || = (1, 1)d/2, and || = (1, 1)d/2, oriented parallel to the lattice plane along , 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 experiments232930. 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 Martin14 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.

Figure 1

Calculated relative lifetime of an electric dipole as a function of distance for the case of a lattice with isotropic magnetic scatterers (dashed) and electric scatterers (solid), at the resonance frequency ω0 for the four in-plane positions (x, y) = (d/2, 0) (blue), (x, y) = (0, d/2) (red), (x, y) = (d/2, d/2) (yellow), and (x, y) = (d/2, 0) (purple) as depicted in the inset.

(a) Dipoles oriented along parallel to the lattice. (b) Dipoles oriented perpendicular to the surface. Thin black lines are calculated lifetimes assuming a homogeneous planar interface.

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 assumed111415. 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 interface2329, which for electric sources perpendicular, respectively parallel to an interface read

and

Here and r 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 lattice7. 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/r3 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 . This is confirmed from the calculated dispersion of the lattice mode with induced dipoles parallel to , presented in Fig. 2(a) (for details on the calculation of the modal dispersion we refer to ref.7). At points X ≡ || = (π/d, 0) and M ≡ || = (π/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 || = 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 || = Γ ≡ (0, 0), remains.

Figure 2

Spectral dependence of the relative lifetime of a dipole emitter placed in the plane (z = 0) of an electric isotropically scattering lattice, as a function of emission frequency and in-plane position along straight paths as illustrated by the colored arrows and the inset.

(a) Dipole along parallel to the lattice. (b) Dipole along perpendicular to the lattice. Markers illustrate different positions of the emitter. The inferred coupling to guided modes, associated with these positions, are marked on the calculated dispersion of the two modes with (c) induced dipoles parallelle to the plane of the lattice along and (d) induced dipoles perpendicular to the lattice plane.

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 Martin14 and Kästel and Fleischhauer15. 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 dipoles18.

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 character3233, 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 cells18, multilayered unit cells, and bi-anisotropic or hyperbolic metasurfaces.

Methods

The integrant in equation (12) typically contains sharp features over the ||-plane, so an adaptive numerical integration is crucial. Furthermore, since

there is a branchpoint at |||| = k. Since the routine method of performing the integration over || into the complex plane34 would require complex ||-lattice sums, we avoid it. Instead we split the integration into two different domains. For k|| within the light cone we use polar coordinates

Outside the light cone cartesian coordinates (k, k) are used for

were used. Rather than computing all 36 tensor elements we directly calculated

where

Additional Information

How to cite this article: Lunnemann, P. and Koenderink, A. F. The local density of optical states of a metasurface. Sci. Rep. 6, 20655; doi: 10.1038/srep20655 (2016).