1. Introduction

In nanophotonics, periodically structured media and periodic texturing of devices are of large interest both for fundamentally new approaches to controlling light, and for applications in spectroscopy [1], solid-state lighting [2–9], and photovoltaics [10,11]. Current developments in fundamental electromagnetics research particularly focus on plasmonic lattices [9–15], metamaterials [16,17], and metasurfaces [18–20]. Plasmonic lattices are two-dimensional subwavelength or diffractive periodic structures of metal nanoparticles, such as spheres, cylinders, bars or pyramids, in which the large scattering strength of plasmonic localized resonances is combined with collective scattering effects to give strong local fields [13]. Currently these are of large interest due to the possibility of strong coupling between localized modes and diffraction conditions [14], and the possibility of strong coupling of the resulting hybrid photonic modes with embedded emitters. In metamaterials the individual inclusions are generalized from the plasmonic case to also support magnetic, or bi-anisotropic resonances, as is the case in split rings [21, 22], and helices [17]. Metasurfaces, finally, are 2D optically thin polarizable particle arrays, in which individual magneto-electric inclusions may be rationally designed to be non-identical, and if periodic, often take the form of frequency selective surfaces, where a single unit cell may contain many magneto-electric inclusions [23]. For gradient structuring of the meta-atoms in a unit cell, such type of metasurfaces can be useful in the wavefront engineering for creating ultrathin optical components [18,19] such as lenses, wave plates and holograms. Metasurface have further been studied extensively as perfect absorbers [24], for asymmetric transmission [25,26], and interfacing bound modes and free propagation photons [2–8].

Our paper is motivated by a host of applications of periodic plasmonic lattices and metasurfaces that have been demonstrated in recent experiments. In particular, periodic texturing of devices with metal structures has been proven to be highly successful in solid-state lighting applications [6–9, 27, 28], and for improved photovoltaics [10, 11]. Figure 1 summarizes several of the main scenarios. In solid state lighting, emission is typically generated in a high index GaN or In_xGa_1−xN semiconductor light emitting device (LED) that is grown on a sapphire substrate. Light that can not be extracted due to total internal reflection and guided modes of the high index layer, can be collected if the LED is textured with a dielectric or plasmonic pattern (Fig. 1(a)). A second scenario appears in current white light LEDs, where emission is generated by using an efficient blue LED as pump for a phosphorous wavelength conversion layer [9]. Here the challenge is to minimize the use of phosphorous material by using plasmonic texturing as a means to enhance blue light absorption, and to optimize (directional) phosphor emission extraction (Fig. 1(b)). A third application scenario is in silicon photovoltaics. Crystalline and amorphous silicon solar cells suffer from long absorption lengths, and high interface reflection. Periodic structuring by plasmon particles can be used to reduce reflections, and to cause diffractive coupling to guided modes in the solar cell [10,11,29,30]. Thereby absorption can be enhanced, and cell thicknesses can be reduced (Fig. 1(c)). In all these scenarios, the periodic plasmon structure is placed in, or on, a nontrivial stratified dielectric system with a complex mode structure.

 

Fig. 1 Sketches of canonical scenarios that our theory can address. (a) Plasmonic lattice patterned on a GaN material based LED structure to improve outcoupling; (b) Remote phosphor for white light emission from LED systems, pumped by a blue LED. Plasmonic texturing can improve absorption, and redirect emission. (c) Thin film solar cell with texturing to improve light incoupling.

Download Full Size | PDF

While the relevance of periodic plasmonic lattices and metasurfaces embedded in stratified systems is large, no simple design or analysis tool is available. Experimental results usually are compared to brute force numerical calculations such as FDTD or FEM modeling [31–34]. These methods can provide accurate full-wave simulation results, but involve extensive computational efforts. Since for a given structure, cataloguing the angle- and frequency dependent response usually takes of order 100 seconds per angle and frequency, the challenge of testing many designs is herculean. Moreover, the mechanism of the modes hybridization, or the interplay between the particle resonances, lattice resonances, and guided stack modes is usually very difficult to retrieve from full wave numerical calculations. It is desirable to have a physically transparent and conceptually simple model that can be applied to plasmonic and metamaterial lattices embedded in arbitrary stratified systems. In this work we present a dipole lattice model that can deal with arbitrary stratified systems, in which the planar scatters are embedded and approximated as point scattering dipoles. Considerable efforts have been devoted to studying the origin of magneto-electric coupling of single scatter, i.e., induced by breaking the symmetry of the scattering particle itself or the hosting environment [21,22,35], as well as in developing dipole lattice models [21,23,36] for plasmonic and magneto-electric scatterer systems, where retardation as well as the cross coupling between electric dipole and magnetic dipole moments are taken into account through Ewald lattice summation. We recently extended this to lattices in front of a simple interface [37,38]. In this paper, we treat dipolar lattices arbitrarily positioned in arbitrarily stratified structures, as sketched in the possible scenarios in Fig. 1. Our model is not limited to any specific range of periodicities or scattering strengths of the lattice, and can thus handle metamaterials, diffractive gratings, and any regime in between. Moreover, using the “array scanning method” we show how to deal with single point emitters driving the periodic systems.

In section II, we give a comprehensive description of the model, which includes the construction of the dyadic Green’s function of magnetoelectric point scatter lattices coupled to planar stratified stacks, the array scanning method, the local density of optical states, and far field emissions patterns. In section III, we present results for solid-state lighting and photovoltaic scenarios. Section IV presents the conclusions of the paper.

2. Theory

Our overall strategy to build a theoretical model is the following. First, we set up a general theory for interaction between dipolar scatterers embedded in a complex dielectric environment. The scatterers are captured in a polarizability, while the interactions are quantified by a Green’s function. We focus on stratified systems, meaning that the Green’s function can be built in an angular spectrum representation, separating TE and TM waves. As next step, we determine how to solve for interactions in infinite lattices using Ewald summation, on the proviso that the excitation field is not localized, but rather has some definite parallel momentum, as would be provided in a plane-wave scattering experiment. Finally we discuss how to obtain the response to localized driving from superposition of the momentum-space solution by the Array Scanning Method (ASM). Section 2.1 introduces the Green’s function, section 2.2 sets up the infinite-lattice dipole problem in momentum space, and section 2.3 introduces the array scanning method. In section 2.4 we discuss an asymptotic expression for evaluating far field radiation, and for calculating the local density of states via reciprocity.

2.1. Dyadic Green’s function for a single magnetoelectric point dipole

For self-consistency, we start from the definition of the Green’s tensor [39,40], which reads,

2.2. 2D lattice and plane wave excitation

We now turn to the case of a 2D periodic magneto-electric point scatter lattice coupled to a slab waveguide under a plane wave illumination. We consider the response of a 2D lattice (lattice vectors R_l = ma₁ + na₂) of point dipoles of polarizability ᾱ to a driving field of parallel wave vector k_‖, where we index lattice sites with the label l = [m, n]. The response is set by solving a self-consistent coupled equation given by

Returning to Eq. (4), the translation invariance of the coupled lattice-waveguide system implies that a Bloch wave form [p_l, m_l]^T = e^ik_‖·R_l [p₀, m₀]^T can be used to simplify the infinite set of coupled equations to a single equation

2.3. Point dipole excitation

We proceed to discuss the coupled lattice-waveguide under a point source excitation by employing the array scanning method (ASM) reviewed, for instance, by Capolino et al. in [48]. The use of Bloch’s theorem, and the Ewald lattice summation technique, as we presented sofar, restricts one to finding the response to extended driving field of definite parallel momentum, as it requires periodicity of the entire system (geometry plus driving) under study. In principle, one can decompose the driving by a single point source in plane waves, and thereby obtain the response of a lattice to local driving by superposition. More practical for implementation is the Array Scanning Method which states that a single point source can be constructed from integration over phased arrays of sources [48]. Provided one can solve for any well-defined phased array $J_{{\mathit{k}}_{\Vert }}=\mathit{J}\sum \delta (\mathit{r}-({\mathit{r}}^{\prime }+{\mathit{R}}_{\mathit{l}}))e^{i{\mathit{k}}_{\Vert }\cdot {\mathit{R}}_{\mathit{l}}}$ of definite parallel momentum, and the same periodicity as the lattice, one can retrieve the single point dipole source by integrating J_k‖ over the Brillouin zone, i.e., $\mathit{J}({\mathit{r}}^{\prime })=\frac{A}{{\left(2\pi \right)}^2}\underset{BZ}{\int }{\mathit{J}}_{{\mathit{k}}_{\Vert }}d{\mathit{k}}_{\Vert }$. For each individual k_‖, the problem is directly compatible with the formulation that was set up for plane wave excitation. When summing the results over all k_‖, one has to ensure each k_‖ receives its own weighting, as given by the background Green’s tensor. The central result of this approach is that the field due to a point source is

The calculation of ḡ⁼(r, r₀) and ḡ⁼(r₀, r′) can be performed by applying the completeness relation as discussed for ḡ^≠(r₀, r₀). A major difference, however, is that now the source point r₀ might be in a different layer than the lattice, as opposed to the lattice sum appearing for particle-particle interactions. An apparent problem is that the Green’s tensor introduces 1/k_nz terms, which introduce a singularity at the light line for medium n. Detailed analysis shows that the 1/k_1z in the Green’s tensor cancels out with the k_1z proportionality of the fresnel transmission coefficient, which removes the numerical singularity. As example, consider a three-layer system, with an emitter in layer 1 (r′ = [0, 0, z′]), while the particle lattice is in layer 2, so that we need a corresponding dyadic Green’s tensor that accounts for the field propagation across the interface through a generalized fresnel transmission coefficient. We take only one block in ḡ_b (r, r′), e.g., ${\overline{\mathit{G}}}_{\mathit{slab}}^{EE,21}({\mathit{r}}_0,{\mathit{r}}^{\prime })$, to illustrate this useful cancellation result, which is given by

As regards to performing the integrations over the Brillouin zone of the type above, we note that even if the 1/k_nz singularities are eliminated, singularities still can occur due to poles in the transmission coefficients at k_‖ that match guided modes of the stratified system. We implemented a finite element scheme to efficiently carry out the numerical integration. The basic idea is to discretize the integration domain into triangles on which a simple quadrature rule is evaluated, and to adaptively refine the triangle density to obtain convergence also at guided mode contributions.

2.4. Far field asymptotic form and radiative LDOS

So far we discussed how to solve for the response of scattering particle lattices to extended and local excitation, without going into the question how to calculate far field observables. Using far field asymptotic forms of Eq. (9) enables one to calculate transmission, reflection, diffraction efficiency, and absorption with ease for any incident parallel momentum. The asymptotic form of Eq. (9) can be obtained [40] by,

In this work we are particularly concerned with mapping the fractional, i.e., parallel wavevector resolved, radiative local density of states, i.e., the radiated flux per solid angle in the far field that is due to a localized source. In the ASM scheme, the asymptotic form of the far field Green’s function can be approximated from Eq. (10) by taking the angular spectrum at a reference plane z = z_r,

An independent method to obtain the radiative local density of states is to use the reciprocity theorem. Reciprocity implies that the local field enhancement inside the waveguide upon far field plane wave illumination is equivalent to calculating the plane wave strength in the far field due to a localized source [49]. Strictly, one needs to identify the local electric field enhancement in the waveguide, upon far field plane wave illumination weighted by a factor of $\frac{1}{k_z}$ with the fractional radiative local density of states enhancement. Hence, one can calculate the far field emission pattern of a localized source coupled to the lattice-waveguide simply by calculating the local field at the emitter position due to a plane wave scattered by the coupled lattice-waveguide, incident from the direction of interest. The reciprocity technique can be easily extended to a bulk light source, i.e., emitters randomly distributed inside the waveguide in a certain region, by integrating |E|² over the active area [50] to efficiently obtain the far field emission pattern. However, the disadvantage of the reciprocity method is that it will not allow to determine the guided-mode contribution to LDOS, nor of absorption induced by the lattice. Any contribution to LDOS (imaginary part of the full system Green’s dyadic) that does not appear as far field power, is not accessible to the reciprocity method. Thereby both the ASM and reciprocity are useful, and provide complementary information.

3. Results and discussion

In this section we discuss results of our formalism with the aim to illustrate the applicability of the model to different scenarios. These are (1) plasmon-enhanced directional emission in solid-state lighting, (2) light trapping in textured thin-film solar cells, and (3) unraveling the Purcell factor provided by a complex plasmonic system into quenching, guided mode contributions, and free space contributions.

3.1. Comparison of the far field pattern between reciprocity method and ASM scanning method

We consider the fractional radiative local density of optical states (FLDOS) in a scenario similar to that recently explored by Lozano et al. [9] and Rodriguez et al. [50] for improvement of white light generation in thin phosphors that are excited by blue LEDs. In this scenario, a random ensemble of emitters in a thin layer of moderate index (about 1.5) is coupled to a lattice of plasmon particles to benefit from “diffractive “surface lattice resonances, and diffractive outcoupling of light emitted into waveguide modes. Experiments [9, 50] in particular showed highly directional emission (several degrees wide emission cone) into angles matching grating outcoupling of waveguided modes. Here we present calculations of the fractional radiative LDOS for four different lattice pitches, to at the same time demonstrate the applicability of the point dipole model, and compare the ASM method to the reciprocity method. In Fig. 2 we consider the angle-dependent far-field emission enhancement for a single emitter coupled to this waveguide-lattice structure, divided by the emission of a bare emitter in the same slab structure without particles.

 

Fig. 2 Angle-dependent far-field emission enhancement of an electric dipole located insidea slab waveguide coupled to a plasmonic lattice, as sketched in Fig. 1 (a). The host medium of the lattice is air, below which a waveguiding layer (index 1.6, thickness 0.2λ₀) is placed on a perfect reflector. a–d (e–h) shows emission pattern calculated from ASM method (reciprocity method) at four lattice constants, i.e., 1/4λ₀, 1/2λ₀, 3/4λ₀, λ₀ respectively. λ₀ is the vacuum wavelength. White lines are the folded slab waveguide inside the light cone. In (b,c,d), the white number pairs denote the diffraction order. The nanoparticle is taken as metallic sphere defined through an electrostatic polarizability of $(3.78+0.135i)\cdot {10}^{-3}{\lambda }_0^3$ (typical for, e.g., a 80 nm aluminum particle in the green). The calculation is done for a single frequency determined by 2π/λ₀, and taking a source height Z₀ = 1/10λ₀.

Download Full Size | PDF

As stratified system we take a n=1.6 slab (the phosphor) positioned on a perfect electric conductor, and a thickness that is 1/5th of the free space emission wavelength (550 nm, green emission), i.e., 110 nm. As particles we take metal nanoparticles placed on the slab, where the particles are defined through a static polarizability ${\overline{\alpha }}_{\mathit{static}}=(3.78+0.135i)\cdot {10}^{-3}{\lambda }_0^3$ (typical for, e.g., 80 nm Aluminum particles in the green). For this example we consider a single source position and orientation, although for actual phosphors one would need to perform incoherent averaging over all source positions and orientations. For this example, we consider the source as an electric dipole moment perpendicular to the interfaces, so that only TM modes are excited. Fig. 2 shows both the angle-dependent far-field emission enhancement calculated according to the ASM method, with assistance of Eq. (13), as shown in Figs. 2 (a)–(d), and calculated according to reciprocity, i.e., as a near-field enhancement upon far-field TM-wave illumination, as shown in Figs. 2 (e)–(h). Results are shown for four different lattice constants, chosen as λ₀/4 (subdiffractive), λ₀/2 (diffractive), 3λ₀/4 (close to 2nd order Bragg diffraction for the waveguide mode) and λ₀.

In Fig. 2 (a,e), the lattice is effectively a dense, non-diffractive metasurface, as the pitch is chosen so small (λ₀/4) that no diffraction contributes to the far field emission pattern. For this lattice pitch, the directional radiation enhancement appears as a cone of about 30 degree opening angle. Here the reader should be warned that we plot enhancement, and when dividing to obtain photoluminescence enhancement, one divides by the bare dipole emission, which is dark near the vertical for the chosen vertical dipole orientation. As the lattice constant increases, the waveguide mode, with dispersion ω = |k_‖|n_mode diffracts, and by diffraction is partially folded inside the light cone, thus coupling out into the far field. For the larger pitch, emission from the dipole that is emitted into the waveguide mode can hence be collected, and is manifested as bright arcs with a radius of curvature given by the waveguide mode index in Figs. 2 (b) and (f), and the panel pairs (c,g) and (d,h). For larger pitch, the arcs approach, and show a anticrossing, as the condition of second order Bragg diffraction for the TM waveguide mode is passed. The bright arcs in enhancement follow essentially the repeated-zone scheme folding of the waveguide slab mode dispersion by the plasmonic lattice (waveguide dispersion replicated every reciprocal lattice vector). Compared with the pure waveguide folding as indicated by the white dash lines (unperturbed waveguide dispersion), plotted in repeated zone scheme dispersion, the bright arcs are broadened and slightly shifted. The broadening reflects a finite propagation length of the folded wave guide mode, due to the intrinsic Ohmic losses from the point scatters and the inclusion of radiation damping. As regards the shift of the mode to smaller mode index, we note that the particles push the mode profile out of the high-index material thereby lowering the mode index, shifting features away from the drawn circles. Qualitatively, i.e., in terms of showing sharp features that reflect the hybridized plasmonic-waveguide modes in a repeated zone scheme, the calculations are in good accord with observations made by Lozano et al. and Schokker et al. [9,12]. In terms of magnitude of the enhancement, as shown in Figs. 2 (c) and (g), the predicted FLDOS enhancements are large, of order 100 (panel d) to 1000 (panel (c)). This large value is due to the fact that we use a perpendicular dipole for excitation, which in the bare layer case hardly radiates at small emission angles. Since our results are normalized to the bare dipole donut-shaped emission pattern, this emphasizes emission at small k_‖ as a pronounced enhancement.

A significant advantage of the reciprocity-based approach is that it requires only plane wave calculations, and that results can be obtained for all source positions and orientations in the unit cell in one calculation. This should be contrasted to the ASM method, that needs a full calculation per dipole position and orientation. However, the reciprocity based method can not separate out contributions to the local density of states that do not result in far field radiation, such as quenching and emission that remains trapped in waveguide modes. As we describe in section 3.3, the ASM does provide means to unravel these phenomena. To highlight the calculation speed advantage of the ASM (the slowest of our two methods), over full wave simulations, we give a brief comparison of the computational time of our method and the full-wave simulation package Comsol. For ASM, each frame of the calculation in Figs. 2 (a)–(d) with 200 ×100 computation points takes roughly 1 minute. For a similarly complex structure, Comsol will take roughly 1 minute for a single computation point, on similar computer hardware. This significant difference in computational time signifies the efficiency of our method, as it allows to rapidly screen different geometries. Beyond being faster by a factor over 10⁴, our method can also be efficiently parallelized to screen different particle types populating the same lattice and same stratified system. The reason is that almost all the computational effort is invested in calculating the stratified system lattice Green’s function. This result does not need to be recalculated if the lattice is populated by different polarizability tensors.

3.2. Photovoltaic system

As second example of the utility of our method we apply our dipole theory to a realistic photovoltaic (PV) structure, in which the plasmonic lattice is coupled to thin-film silicon solar cells containing several layers, as studied by van Lare [10] et. al. In particular, the application area is to amorphous silicon (a-Si:H) photovoltaics, in which 350 nm of silicon is sandwiched between an ITO contact, on one side, and a 80 nm ZnO buffer layer to a 200 nm thick silver back contact on the other side. One challenge of such cells is that the silicon is much thinner than the absorption length, especially in the red. Consequently, these cells have up to 80% reflectivity for wavelengths above 650 nm. The aim of plasmonic structuring is to improve absorption in the amorphous silicon layer, despite the fact that it is far thinner than the bulk absorption length, by coupling into guided modes. This example illustrates the potential application of our method to rapidly screen complex photonic structures in thin-film photovoltaics. In the construction of our dyadic Green’s function, we implement the full layer-structure of the thin-film silicon solar cells. However, to a good approximation, we take the 200 nm thick Ag back contact as infinite. We use tabulated optical constants for all materials, taken from [10,51].

Figures 3(a) and (b) shows the calculated reflection, absorption in a-Si:H i-layer for the bare layer geometry, and geometries that include embedded plasmonic particle lattices of pitch 500 nm, and two different particle sizes (radius of 70 nm, and 120 nm). In Fig. 3(a), it is evident that without the plasmonic lattice, the cell shows minimum reflectivity at 550 nm, and a large 80% reflectivity for wavelengths above 650 nm. This is in good accord with the full wave simulations of van Lare, and highlights a main problem of thin film amorphous solar cell, namely that light is reflected at the back contact after a single pass, without being absorbed. The plasmonic lattices strongly reduce the cell reflection in a set of narrow spectral lines, where the introduced reflection minima are most pronounced for the larger particles, again in good accord with the work of van Lare [10]. Fig. 3 (b) shows the absorption inside the amorphous silicon material only. The plasmonic structure introduces pronounced absorption peaks at wavelengths near the reflection minima, that are due to diffractive coupling into waveguide modes.

 

Fig. 3 (a) Calculated reflection from the layer geometry shown in the inset, which reflects the cell geometry of a realistic photovoltaic cell based on amorphous silicon [10]. (b) Calculated absorption in the a-Si:H i-layer. Red line for the bare layer geometry without plasmonic lattice, blue line for the plasmonic lattice with particle radius of 70 nm, black line for the plasmonic lattice with particle radius of 120 nm. (c) Absorption as function of incident angle and photon energy. In (a–b), light is illuminated at normal incidence. In (c), the incident light has TE polarization.

Download Full Size | PDF

To further quantify the role of the diffractive coupling between the guided modes and the free propagating modes in enhancing absorption, we study an ω-k diagram of the absorption [11]. Figure 3(c) shows the calculated absorption in the a-Si:H i-layer for the case where a plasmonic lattice with particles of radius 120 nm is embedded in the system, as function of incident angle (with TE polarization) and photon energy. As in Figs. 3 (a) and (b), the absorption diagram as function of frequency and wave vector shows sharp and dispersive features on top of a background that is independent of angle, that shows a sharp step in absorption at 0.6 μm wavelength. This sharp step is due to a sharp change in absorption coefficients of a-Si:H layer at 0.6 μm. The very sharp and dispersive features are due to diffractive resonances, i.e, diffractive coupling to the guided modes in the a-Si:H layer, as well as lattice resonances given by the sharp bright lines. Very similar wavelength and angle-dependent features in photocurrent have been measured for actual cells by Ferry et al. [11], and such measurements can be an important tool to unravel the physics behind photovoltaic performance improvements. Beyond the relevance of this particular amorphous silicon cell example, the calculations evidence that our method can function as a rapid means to screen different light trapping strategies without the need for full wave simulations.

3.3. Excitation of dark modes by point sources

For many interesting quantities, plane wave calculations, together with the reciprocity technique suffice. These include angle-dependent extinction, absorption, reflection, refraction and diffraction and the differential radiated power for a localized source to any far field angle. However, for some highly interesting quantities, the array scanning method is indispensable. In general, a photonic structure may offer guided modes as efficient channels for emission, or in case of a plasmonic structure, may cause absorption. Generally, the total decay rate enhancement is due to the sum of induced quenching, emission into guided modes, and an increased rate of emission into far field modes. Only the last term is accessible to the reciprocity technique. The array scanning method, on the other hand, can be used to unravel quenching, and guided mode contributions. The full Brillouin zone integral in the array scanning method provides the sum of all three contributions, i.e., the full local density of states. In addition, as it gives the field at any position due to a localized source, one can calculate absorption. Finally, analyzing the integrand of Eq. (10) provides insight in the specific wave vectors that contribute to the local density of states, for instance pointing at guided mode contributions to the emission.

As third example to illustrate our method, we discuss a fluorescence problem tackled by the ASM method in which part of the emission enhancement of a dipole source is funneled into guided modes. As geometry we take a thin (d = 116 nm) layer of GaAs in air, and assume that a square lattice (lattice constant p) of aluminum particles (radius 0.25p) is embedded in the slab at about 2/3rds of the slab thickness (all dielectric constants taken from Palik [51]). We assume a source located just above the slab top surface (25 nm separation). This geometry is not altogether realistic from an experimental stand point but has many illustrative challenges from a computational point of view. These are (A) as opposed to the previous examples the particles are inside a layer while the source is on top, (B) strong coupling of the particles to the slab waveguide mode, (C) no mirror symmetry in the geometry, so that particles couple to TE and TM modes, and (D) evanescent coupling of the source to the slab wave guide mode. We consider this computationally tough problem in a range of near- to mid-IR wavelengths, and for two pitches.

As shown in Fig. 4 (a), we study the total decay rate enhancement (in current literature often referred to as Purcell factor (F)) of a dipolar emitter, coupled to the scatterer/waveguide system. The total decay rate enhancement is given by the total emission rate γ_tot divided by the vacuum rate γ₀. In Fig. 4(a), results are shown for the emitter oriented in plane (along x-direction, black lines), and normal to the plane (along z-direction, blue lines). In this diagram solid lines refer to p=145 nm, while dashed lines refer to p =173 nm. For x-dipoles, the Purcell factor increases moderately, towards long wavelengths, which is due to increased material losses in aluminum particles at long wavelength. In contrast, the increase of the Purcell factor for a z-dipole coupled to the same configuration of lattice-waveguide system is pronounced, attaining values up to 40 or 80 for different pitches, as evident by the sharp peaks denoted by F (1.37 μm) and G (1.57 μm) in Fig. 4 (a). The large enhancement of γ_tot results from the excitation of dark modes, i.e., guided modes, of the plasmonic lattice for long wavelengths, rather than merely from increased material losses. This is evident from two facts: (1) the peaks F, G are not close to the wavelength of maximum material losses, (2) the enhanced emission rate, up to 80, for z-dipoles far exceeds that for x-dipoles, which is less than 8.

 

Fig. 4 Emission rate of a dipole coupled to an aluminum sphere array embedded in a GaAs slab waveguide for an x-oriented dipole (black curves, left y-axis), and for a z-oriented dipole (blue curves, righthand y-axis). Solid (dashed) curves are for a lattice pitch of 145 (173) nm. (b) Contour plot of the integrand of a bare slab Green’s function ${\text{log}}_{10}\left|\mathit{Im}\left[G_{zz}^s\right]\right|$. (c–g) Contour plot of Im[G_zz], i.e., Ī(k_‖, r, r₀) in Eq. (10) for the special points, as marked by C,D, E, F, G in (a). In (b–g), the blue (black) dotted lines indicate the folded light line in air (GaAs), the white and magenta (red) dotted lines for TE (TM) bound modes.

Download Full Size | PDF

To further study the dark mode excitation, we plot the integrands of the scattering part of the Green’s function, i.e., Ī(k_‖, r, r₀) in Eq. (10), in one quadrant of the full Brillouin zone in Figs. 4(c)–(g) for special points C,D, E, F, G marked in Fig. 4(a), in all cases focusing on the z-oriented dipole case. For reference, one can also evaluate the integrand for a ‘bare’ slab, containing a lattice of particles of zero polarizability, shown in Fig. 4 (b) at wavelength of 1 μm. In all panels, the plot range encompasses exactly the positive quadrant of the first Brillouin zone. In Figs. 4 (b)–(g), blue (black) dotted lines indicate the light line in air (GaAs). Due to C₄ symmetry of the square lattice, the integrands in the full Brillouin zone can be obtained by mirroring along k_x = 0, and k_y = 0. In order to directly compare the k-space spectrum content, we use the same colormap for all contour plots in (c–g). The bare slab system supports three waveguides modes at short wavelength as shown in (b), i.e., two TE modes indicated by white and magenta dashed lines, plus one TM mode indicated by a red dashed line. At longer wavelength the fundamental TM /TE mode indicated by the magenta/red dotted line remains, while the higher order TE mode (white dotted line) vanishes.

In the bare slab waveguide in Fig. 4 (b), the contribution of the waveguide mode to the integrand of scattering part of Green’s function is weak, as shown in the contour plot of log₁₀ |Im[G_zz]|. For larger in-plane momentum, the integrand of the scattering part of Green’s function approaches 0. However, with the lattice of scatterers in place, the integrands of the second term in Eq. (10) are strongly mediated by the periodic scatterers. At short wavelengths, i.e., (c,d), there are three major features: (1) the bare slab waveguide supports well-defined waveguide modes, indicated by the dark arcs; (2) diffraction occurs, where the guided modes of the slab waveguide folds back into the Brillouin zone; (3) dark plasmonic modes are excited on the lattice in Figs. 4(d)–(g) in a limited region of k-space close to the folded guided mode and is evident as very large integrand (bright features). At longer wavelengths (e–g), the dipole lattices become denser with respect to vacuum wavelength, so that back-folding no longer leads to intersections appearing with back-folded waveguide modes. Plasmon hybridization becomes stronger, and hence lattice resonances may play a significant role in the enhancement of F. Indeed, the integrand in Figs. 4(e)–(g) shows that there are considerable plasmonic contributions (broad bright features) in the Brillouin zone within the (folded) light cone (roughly 3.35k₀) given by GaAs, besides the narrow contribution from guided modes, which are only weakly confined at these long wavelengths. In Fig. 4 (d), a major part of the emission is inside the folded light cone given by the slab medium (GaAs), and includes a considerable amount inside the vacuum light cone (k₀), which indicates that a significant fraction of the total emission can be delivered into the far field. In contrast, one notices that in the integrand in Figs. 4(e)–(g), the major contribution to γ_tot lies outside the vacuum light cone. Hence, most of the emission is funneled into the dark modes of the lattice, and is eventually dissipated due to material losses.

The emission distribution in k-space, i.e., the fraction of the emission fuelled into far-field, as well as into the dark modes, for points C–G are tabulated in Table. (1). As evident in Figs. 4(c)–(e), the the major contribution to the integrand appears at parallel wave vectors beyond the light-line of the host material of the lattice. The fraction of emission into the far field, i.e., |k_‖| < k₀, reaches its maximum value (0.35%) at the point marked D. Further to the red, i.e, at point E, we see that emission with large parallel k-component, i.e., |k_‖| ≥ n_GaAs k₀, starts to contribute due to the stronger plasmonic hybridization. As the wavelength increases further, the emission with |k_‖| ≥ n_GaAs k₀ quickly dominates, which yields nearly 100% emission funnelled into the dark mode of the plasmon lattice. Due to the large |k_‖|, those modes are tightly confined by the plasmonic lattice, and quickly dissipate due to material losses. It is important to realize that designing the k-space emission distribution is highly non-trivial, and the rules behind such type of design are beyond the scope of this paper. Essentially, it is the competition among channels, i.e., the emission with |k_‖| < k₀, k₀ < |k_‖| ≤ n_GaAs k₀, and |k_‖| ≥ n_GaAs k₀, that ultimately determines how emission occurs in the complex photonic environments [52].

Table 1. Fraction of emission of z-dipole coupled to waveguide-lattice system in k-space. C, D, E, F, G denote the special points marked in Fig. 4 (a).

View Table

To conclude, our method can aid evaluating different photonic designs for enhancing light-matter interaction strength on basis of the fact that the integrand in the ASM method can be separated in different contributions. This should be contrasted to full-wave simulation tools that, first of all, struggle with single sources in infinitely periodic systems, and secondly, do not provide an explicit separation in poles. The ASM-method provides key insights in the emission distribution in k-space, which could be utilized in many problems that require engineering the emission distribution in k-space as well as the separation of emission into far field, guided modes, and absorption.

4. Conclusion

In closing, we presented a general point dipole theory for periodic metasurfaces that are embedded in arbitrary stratified layer systems, and are driven either by plane waves, or by single point emitters. With the assistance of Ewald’s summation technique and dyadic Green’s function of stratified layers, we can access complex magnetoelectric scattering lattices coupled to layered photonic structures. One can with this method rapidly understand the physics of complicated structures as long as one can approximate meta-atoms as point dipoles with magnetoelectric responses, while accounting for all particle-particle interactions self-consistently via freely propagating photons as well as bound modes provided by the stratified layers. We further extended the plane wave excitation case to excitation by a single, localized, dipolar source, through the array scanning method. Thereby, all the relevant optical properties, such as LDOS and far field emission patterns, of the coupled waveguide lattice, can be extracted straightforwardly. Due to the generic nature of the multiple scattering framework that is adapted in this paper, our approach to study periodic lattices that contain just one one scatter per unit cell can be easily extended to complex periodic lattices which contain more than one scatters in each unit cell, [53]. In this extension, one decomposes the complex lattice as the sum of simple lattices, the number of which equals to the number of meta-atoms per primitive unit cell; next one extracts the the effective polarizability tensor for each simple lattice separately; lastly, one builds up a coupled dipole equation for the induced dipole in each scatterer in the unit cell, with the polarizability tensor for each meta-atom being the effective polarizability tensor dressed by the lattice Green’s function. With the aforementioned procedure, our general point dipole theory can be readily extended to very complex lattices, meaning that many different gradient meta-surfaces can be studied.

As applications of the proposed method, we demonstated that it applies readily to solid-state lighting scenarios where point emitters are coupled to periodic plasmonic lattices, in turn embedded in high index layered systems. Also, we applied our method to a plasmon enhanced thin-film solar cell structure containing a plasmonic nanoantenna array. We find that our method can efficiently handle such structures, capturing salient features due to diffractive coupling in terms of extinction, absorption and reflection. Due to high efficiency and flexibility of the dipole theory, it can be used for both diffractive, dilute lattices, and for dense lattices with complex clusters of inclusions that form a metamaterial or metasurface. While sacrificing accuracy in approximating geometry in the point dipole approximation, the method has a crucial speed advantage of about 10⁴, and can deal with stratified systems exactly, instead of relying on a finite volume truncation. We thus envision that our method is useful for studying metasurfaces, as well as related applications such as light-emitting diodes and solar cells.

Appendix 1: Validity check on the transmission and reflection of the point dipole model in layered structures

In this appendix we study the validity of the point dipole approximation by comparing reflectivity and transmissivity as calculated by our theory with simulated results from Multiphysics COMSOL 5.2 for an example stratified system. In this model we take a three-layer system, (superstrate, thin layer, substrate). The thin layer contains particles with dielectric constant specified by a Drude model $∊(\omega )=1-{\omega }_p^2/({\omega }^2+i\omega \gamma )$. In this particular test [54] we take the thin layer (thickness d) as air, and vary the superstrate and substrate index (∊₁ and ∊₃), while the silver scatterers are placed in a square grid at height z₀ in the slab. We consider normal incidence light with polarization along x, and analyze reflectivity (R) and transmission (T).

For a fair comparison, it is important that the point dipole model is applied with a correct polarizability tensor ᾱ₀ of the single scatterer. We extract ᾱ₀ of the single scatter from a single scatterer full-wave extinction calculation. For a spherical scatterer, the polarizability tensor is diagonal ᾱ₀ = α₀Ī and in the dipole approximation α₀ directly links to to extinction via Im(α₀) = σ_ext/4πk₀. We have verified that a Lorentz line-shape to the calculated extinction (${\sigma }_{\mathit{ext}}=L\frac{{\omega }_0^2}{{\omega }_0^2-{\omega }^2-j\gamma \omega }$) is a good fit to the full wave simulation.

In our test, we take the distance between superstrate and substrate d as 75nm, with the particles midway. Each scatter has a 50 nm diameter, small compared to the wavelength. We consider the lattice in free space, ∊₁ = ∊₂ = ∊₃ = 1, next to a single interface (∊₃ = 9), and in between two interfaces (∊₁ = 9, ∊₂ = 1, ∊₃ = 9) that cause multiple reflections. The explicit form of Green’s function can be seen in Appendix(B–C). Figure 5 (a) shows the α₀ and α_eff for each structure mentioned above, where it should be noted that the polarizabilities shifts as a function of the choice of environment. Figs. 5 (b)–(d) shows R and T for the three aforementioned structures using our point dipole approached compared with fullwave simulation calculated using COMSOL, which shows excellent agreements. Thus, we conclude that the point dipole theory has a good accuracy in the analysis of lattice scattering in stratified structures, provided one starts with an accurate polarizability tensor.

 

Fig. 5 (a) Polarizability α₀ and α_eff for free space(1), single interface(2) and slab(3) respectively. (b)–(d) Reflectivity and transmissivity of point dipole theory and COMSOL. Red and blue line indicate the reflectivity and transmissivity calculated using point theory, while red and blue circles indicate the reflectivity and transmissivity of Multiphysics COMSOL 5.2 simulation. (b), (c) and (d) indicate the structure of free space, single interface and slab respectively.

Download Full Size | PDF

Appendix 2: Construction of free space dyadic Green’s function based on plane wave expansion

For a given wave vector k = [k_‖, k_z], the TE and TM polarization modes associated with a plane wave are given by E_TE (r) = M(k_‖, ±₁, ±₂)e^{ik_‖·r_‖}, and E_TM (r) = N(k_‖, ±₁, ±₂)e^{ik_‖·r_‖} respectively, where M and N are defined as follows,

(14a)$$\mathit{M}\left({\mathit{k}}_{\Vert ,{\pm }_1,{\pm }_2}\right)=e^{i_{{\pm }_2}{k}_{z}z}\left(\widehat{x}\frac{-k_y}{k_{\Vert }}+\widehat{y}\frac{k_x}{k_{\Vert }}\right),$$

(14b)$$\mathit{N}\left({\mathit{k}}_{\Vert ,{\pm }_1,{\pm }_2}\right)=e^{i{\pm }_2{k}_{z}z}\left[\widehat{z}\frac{i{k}_{\Vert }}{k_0}-\widehat{x}\frac{i({\pm }_1)k_z{k}_x}{k_{\Vert }{k}_0}-\widehat{y}\frac{i({\pm }_1)k_z{k}_y}{k_{\Vert }{k}_0}\right],$$

where r = [r_‖, z], k_‖ = |k_‖|, $k_0=\sqrt{k_{\Vert }^2+k_z^2}$. The first sign factor describes the propagation direction along z-axis, while the second sign factor denotes the sign of k_z in the construction of the basis function. Throughout the paper, we will use M(r, ±, ±) and N(r, ±, ±) for short. In the above formulation, the prefactors of TE and TM waves are arranged to fulfill the relation of

(15)$$\begin{array}{l}\nabla \times \left[e^{i{\mathit{k}}_{\Vert }\cdot {\mathit{r}}_{\Vert }}M({\mathit{k}}_{\Vert },\pm ,\pm )\right]=e^{i{\mathit{k}}_{\Vert }\cdot {\mathit{r}}_{\Vert }}{k}_{0}N({\mathit{k}}_{\Vert },\pm ,\pm ).\\ \nabla \times \left[e^{i{\mathit{k}}_{\Vert }\cdot {\mathit{r}}_{\Vert }}N({\mathit{k}}_{\Vert },\pm ,\pm )\right]=e^{i{\mathit{k}}_{\Vert }\cdot {\mathit{r}}_{\Vert }}{k}_{0}M({\mathit{k}}_{\Vert },\pm ,\pm ).\end{array}$$

In free space, the dyadic Green’s function of the electric field can be expanded using the plane waves as follows (upper sign applies for z >= z′),

(16)$$\begin{array}{l}{\overline{\mathit{G}}}_0^{ee}(\mathit{r},{\mathit{r}}^{\prime })=\frac{i}{8{\pi }^2}\int d{k}_{x}d{k}_y\frac{e^{i\theta }}{k_z}\left[{\mathit{M}}_{+}\otimes {\mathit{M}}_{-}^{\prime }-{\mathit{N}}_{+}\otimes {\mathit{N}}_{-}^{\prime }\right]\\ {\overline{\mathit{G}}}_0^{he}(\mathit{r},{\mathit{r}}^{\prime })=\frac{i{k}_0}{8{\pi }^2}\int d{k}_{x}d{k}_y\frac{e^{i\theta }}{k_z}\left[{\mathit{N}}_{+}\otimes {\mathit{M}}_{-}^{\prime }-{\mathit{M}}_{+}\otimes {\mathit{N}}_{-}^{\prime }\right],\\ {\overline{\mathit{G}}}_0^{hh}(\mathit{r},{\mathit{r}}^{\prime })=\frac{i}{8{\pi }^2}\int d{k}_{x}d{k}_y\frac{e^{i\theta }}{k_z}\left[-{\mathit{N}}_{+}\otimes {\mathit{N}}_{-}^{\prime }+{\mathit{M}}_{+}\otimes {\mathit{M}}_{-}^{\prime }\right],\\ {\overline{\mathit{G}}}_0^{he}(\mathit{r},{\mathit{r}}^{\prime })=\frac{i{k}_0}{8{\pi }^2}\int d{k}_{x}d{k}_y\frac{e^{i\theta }}{k_z}\left[-{\mathit{M}}_{+}\otimes {\mathit{N}}_{-}^{\prime }+{\mathit{N}}_{+}\otimes {\mathit{M}}_{-}^{\prime }\right].\end{array}$$

where $\theta ={\mathit{k}}_{\Vert }\cdot ({\mathit{r}}_{\Vert }-{\mathit{r}}_{\Vert }^{\prime })$, M_(±)/N_(±) = M/N(k_‖, ±, ±), ${\mathit{M}}_{(\mp )}^{\prime }/{\mathit{N}}_{(\mp )}^{\prime }=\mathit{M}/\mathit{N}({\mathit{k}}_{\Vert },\pm ,\mp )$.

Appendix 3: Dyadic Green’s function for stratified structures

Dyadic Green’s function for single interface

The dyadic Green’s function for an interface or stratified structure can be expanded using the plane waves, i.e., ${\overline{\mathit{G}}}_{\mathit{if}}(\mathit{r},{\mathit{r}}^{\prime })={\overline{\mathit{G}}}_0(\mathit{r},{\mathit{r}}^{\prime })+{\overline{\mathit{G}}}_s^{\mathit{if}}(\mathit{r},{\mathit{r}}^{\prime })$, where Ḡ₀(r, r′) is the free space term given by Eq. (16).

(17)$$\begin{array}{l}{\overline{\mathit{G}}}_{s,ee}^{\mathit{if}}=\frac{i}{8{\pi }^2}\int d{\mathit{k}}_{\Vert }\frac{e^{i\theta }}{k_z}\left\{{\gamma }_s{\mathit{M}}_{+}\otimes {\mathit{M}}_{+}^{\prime }-{\gamma }_p{\mathit{N}}_{+}\otimes {\mathit{N}}_{+}^{\prime }\right\}\\ {\overline{\mathit{G}}}_{s,he}^{\mathit{if}}=\frac{i{k}_0}{8{\pi }^2}\int d{\mathit{k}}_{\Vert }\frac{e^{i\theta }}{k_z}\left\{{\gamma }_s{\mathit{N}}_{+}\otimes {\mathit{M}}_{+}^{\prime }-{\gamma }_p{\mathit{M}}_{+}\otimes {\mathit{N}}_{+}^{\prime }\right\},\\ {\overline{\mathit{G}}}_{s,hh}^{\mathit{if}}=\frac{i}{8{\pi }^2}\int d{\mathit{k}}_{\Vert }\frac{e^{i\theta }}{k_z}\left\{{\gamma }_p{\mathit{M}}_{+}\otimes {\mathit{M}}_{+}^{\prime }-{\gamma }_s{\mathit{N}}_{+}\otimes {\mathit{N}}_{+}^{\prime }\right\},\\ {\overline{\mathit{G}}}_{s,eh}^{\mathit{if}}=\frac{i{k}_0}{8{\pi }^2}\int d{\mathit{k}}_{\Vert }\frac{e^{i\theta }}{k_z}\left\{{\gamma }_p{\mathit{N}}_{+}\otimes {\mathit{M}}_{+}^{\prime }-{\gamma }_s{\mathit{M}}_{+}\otimes {\mathit{N}}_{+}^{\prime }\right\}.\end{array}$$

Dyadic Green’s function for a slab structure

The dyadic Green’s function for a slab structure can be expanded using the plane waves as follows [41–43], see the geometry in Fig. 1 (a), in which both the emitter and the detection point are located in the layer 2,

(18)$$\begin{array}{ll}{\overline{\mathit{G}}}_{2,2}^{ee}(\mathit{r},{\mathit{r}}^{\prime })=\hfill & \hfill \int d{\mathit{k}}_{\Vert }\frac{i{e}^{i\theta }}{8{\pi }^2{k}_{z2}}\{\frac{1}{{\mathrm{\Delta }}_s}\left[{\mathit{M}}_{+}+{\gamma }_{s,21}(k_{\Vert }){\mathit{M}}_{(-)}\eta \right]\otimes \left[{\mathit{M}}_{(-)}^{\prime }+{\gamma }_{s,23}(k_{\Vert }){\mathit{M}}_{(+)}^{\prime }\right]-\\ \hfill & \hfill \frac{1}{{\mathrm{\Delta }}_p}\left[{\mathit{N}}_{(+)}+{\gamma }_{p,21}(k_{\Vert }){\mathit{N}}_{(-)}\eta \right]\otimes \left[{\mathit{N}}_{(-)}^{\prime }+{\gamma }_{p,23}(k_{\Vert }){\mathit{N}}_{(+)}^{\prime }\right]\},\\ {\overline{\mathit{G}}}_{2,2}^{he}(\mathit{r},{\mathit{r}}^{\prime })=\hfill & \hfill \int d{\mathit{k}}_{\Vert }\frac{i{e}^{i\theta }{k}_0}{8{\pi }^2{k}_{z2}}\{\frac{1}{{\mathrm{\Delta }}_s}\left[{\mathit{N}}_{+}+{\gamma }_{s,21}(k_{\Vert }){\mathit{N}}_{(-)}\eta \right]\otimes \left[{\mathit{M}}_{(-)}^{\prime }+{\gamma }_{s,23}(k_{\Vert }){\mathit{M}}_{(+)}^{\prime }\right]-\\ \hfill & \hfill \frac{1}{{\mathrm{\Delta }}_p}\left[{\mathit{M}}_{(+)}+{\gamma }_{p,21}(k_{\Vert }){\mathit{M}}_{(-)}\eta \right]\otimes \left[{\mathit{N}}_{(-)}^{\prime }+{\gamma }_{p,23}(k_{\Vert }){\mathit{N}}_{(+)}^{\prime }\right]\},\\ {\overline{\mathit{G}}}_{2,2}^{hh}(\mathit{r},{\mathit{r}}^{\prime })=\hfill & \hfill \int d{\mathit{k}}_{\Vert }\frac{i{e}^{i\theta }}{8{\pi }^2{k}_{z2}}\{\frac{1}{{\mathrm{\Delta }}_s}\left[{\mathit{N}}_{+}+{\gamma }_{s,21}(k_{\Vert }){\mathit{N}}_{(-)}\eta \right]\otimes \left[{\mathit{N}}_{(-)}^{\prime }+{\gamma }_{s,23}(k_{\Vert }){\mathit{N}}_{(+)}^{\prime }\right]-\\ \hfill & \hfill \frac{1}{{\mathrm{\Delta }}_p}\left[{\mathit{M}}_{(+)}+{\gamma }_{p,21}(k_{\Vert }){\mathit{M}}_{(-)}\eta \right]\otimes \left[{\mathit{M}}_{(+)}^{\prime }+{\gamma }_{p,23}(k_{\Vert }){\mathit{M}}_{(-)}^{\prime }\right]\},\\ {\overline{\mathit{G}}}_{2,2}^{eh}(\mathit{r},{\mathit{r}}^{\prime })=\hfill & \hfill \int d{\mathit{k}}_{\Vert }\frac{i{e}^{i\theta }{k}_0}{8{\pi }^2{k}_{z2}}\{\frac{1}{{\mathrm{\Delta }}_s}\left[{\mathit{M}}_{(+)}+{\gamma }_{s,21}(k_{\Vert }){\mathit{M}}_{(-)}\eta \right]\otimes \left[{\mathit{N}}_{(-)}^{\prime }+{\gamma }_{s,23}(k_{\Vert }){\mathit{N}}_{(+)}^{\prime }\right]-\\ \hfill & \hfill \frac{1}{{\mathrm{\Delta }}_p}\left[{\mathit{N}}_{(+)}+{\gamma }_{p,21}(k_{\Vert }){\mathit{N}}_{(-)}\eta \right]\otimes \left[{\mathit{M}}_{(-)}^{\prime }+{\gamma }_{p,23}(k_{\Vert }){\mathit{M}}_{(+)}^{\prime }\right]\}.\end{array}$$

where Δ_s/p = 1 − γ_s/p,21(k_‖)γ_s/p,23(k_‖)e^i2k_zd, and ${\gamma }_{s,ij}=\frac{k_{zi}-k_{zj}}{k_{zi}+k_{zj}}$, ${\gamma }_{p,ij}=\frac{{∊}_j{k}_{zi}-{∊}_i{k}_{zj}}{k_{zi}{∊}_j+k_{zj}{∊}_i}$, η = e^i2k_zd.

Funding

National Natural Science Foundation of China (Grant No. 61405067); Fundamental Research Funds for the Central Universities, (HUST: 2017KFYXJJ027); The Netherlands Organization for Scientific Research.

Acknowledgment

Y.T. Chen would like to thank Prof. S. A. Tretyakov for fruitful discussions. Y.T. Chen acknowledges financial support from the National Natural Science Foundation of China (Grant No. 61405067), and the Fundamental Research Funds for the Central Universities, HUST: 2017KFYXJJ027. This work is part of the research program of “The Netherlands Organization for Scientific Research (NWO)”. A.F.K. gratefully acknowledges an NWO-Vidi and Vici grant for financial support.

References and links

1. K. A. Willets and R. P. Van Duyne, “Localized surface plasmon resonance spectroscopy and sensing,” Annu. Rev. Phys. Chem. 58, 267–297 (2007). [CrossRef]  

2. I. Schnitzer, E. Yablonovitch, C. Caneau, T. J. Gmitter, and A. Scherer, “30% external quantum efficiency from surface textured, thin-film light-emitting diodes,” Appl. Phys. Lett. 63(16), 2174–2176 (1993). [CrossRef]  

3. J. Vuckovic, M. Lončar, and A. Scherer, “Surface plasmon enhanced light-emitting diode,” IEEE J. Quantum Electron. 36(10), 1131–1144 (2000). [CrossRef]  

4. K. Okamoto, I. Niki, A. Shvartser, Y. Narukawa, T. Mukai, and A. Scherer, “Surface-plasmon-enhanced light emitters based on InGaN quantum wells,” Nat. Mater. 3(9), 601–605 (2004). [CrossRef]   [PubMed]  

5. D. M. Yeh, C. F. Huang, C. Y. Chen, Y. C. Lu, and C. C. Yang, “Localized surface plasmon-induced emission enhancement of a green light-emitting diode,” Nanotechnology 19(34), 345201 (2008). [CrossRef]   [PubMed]  

6. J. Henson, J. DiMaria, E. Dimakis, T. D. Moustakas, and R. Paiella, “Plasmon-enhanced light emission based on lattice resonances of silver nanocylinder arrays,” Opt. Lett. 37(1), 79–81 (2012). [CrossRef]   [PubMed]  

7. W. L. Barnes, “Electromagnetic crystals for surface plasmon polaritons and the extraction of light from emissive devices,” J. Lightwave Technol. 17(11), 2170 (1999). [CrossRef]  

8. H. Rigneault, F. Lemarchand, A. Sentenac, and H. Giovanini, “Extraction of light from sources located inside waveguide grating structures,” Opt. Lett. 24(3), 148–150 (1999). [CrossRef]  

9. G. Lozano, D. J. Louwers, S. R. K. Rodríguez, S. Murai, O. T. Jansen, M. A. Verschuuren, and J. Gómez Rivas, “Plasmonics for solid-state lighting: enhanced excitation and directional emission of highly efficient light sources,” Light Sci. Appl. 2(5), e66 (2013). [CrossRef]  

10. M. van Lare, F. Lenzmann, M. A. Verschuuren, and A. Polman, “Mode coupling by plasmonic surface scatterers in thin-film silicon solar cells,” Appl. Phys. Lett. 101(22), 221110 (2012). [CrossRef]  

11. V. E. Ferry, M. A. Verschuuren, H. B. T. Li, E. Verhagen, R. J. Walters, R. E. I. Schropp, H. A. Atwater, and A. Polman, “Light trapping in ultrathin plasmonic solar cells,” Opt. Express 18(102), A237–A245 (2010). [CrossRef]   [PubMed]  

12. A. H. Schokker and A. F. Koenderink, “Lasing at the band edges of plasmonic lattices,” Phys. Rev. B 90(15), 155452 (2014). [CrossRef]  

13. M. Auguié and W. L. Barnes, “Collective resonances in gold nanoparticle arrays,” Phys. Rev. Lett. 101(14), 143902 (2008). [CrossRef]   [PubMed]  

14. P. Törmä and W. L. Barnes, “Strong coupling between surface plasmon polaritons and emitters: a review,” Rep. Prog. Phys. 78(1), 013901 (2015). [CrossRef]  

15. V. G. Kravets, F. Schedin, and A. N. Grigorenko, “Extremely narrow plasmon resonances based on diffraction coupling of localized plasmons in arrays of metallic nanoparticles,” Phys. Rev. Lett. 101(8), 087403 (2008). [CrossRef]   [PubMed]  

16. I. Sersic, M. Frimmer, E. Verhagen, and A. F. Koenderink, “Electric and magnetic dipole coupling in near-infrared split-ring metamaterial arrays,” Phys. Rev. Lett. 103(21), 213902 (2009). [CrossRef]  

17. J. K. Gansel, M. Thiel, M. S. Rill, M. Decker, K. Bade, K V. Saile, G. V. Freymann, S. Linden, and M. Wegener, “Gold helix photonic metamaterial as broadband circular polarizer,” Science 325(5947), 1513–1515 (2009). [CrossRef]   [PubMed]  

18. N. Yu, P. Genevet, M. A. Kats, F. Aieta, J. P. Tetienne, F. Capasso, and Z. Gaburro, “Light propagation with phase discontinuities: generalized laws of reflection and refraction,” Science 334(6054), 333–337 (2011). [CrossRef]   [PubMed]  

19. N. Yu and F. Capasso, “Flat optics with designer metasurfaces,” Nat. Mater. 13(2), 139–150 (2014). [CrossRef]   [PubMed]  

20. A. E. Minovich, A. E. Miroshnichenko, A. Y. Bykov, T. V. Murzina, D. N. Neshev, and Y. S. Kivshar, “Functional and nonlinear optical metasurfaces,” Laser Photon. Rev. 9(2), 195–213 (2015). [CrossRef]  

21. I. Sersic, C. Tuambilangana, T. Kampfrath, and A. F. Koenderink, “Magnetoelectric point scattering theory for metamaterial scatterers,” Phys. Rev. B 83(24), 245102 (2011). [CrossRef]  

22. J. Xu, B. B. Wu, and Y. T. Chen, “Elimination of polarization degeneracy in circularly symmetric bianisotropic waveguides: a decoupled case,” Opt. Express 23(9), 11566–11575 (2015). [CrossRef]   [PubMed]  

23. Y. Ra’di, C. R. Simovski, and S. A. Tretyakov, “Thin perfect absorbers for electromagnetic waves: theory, design, and realizations,” Phys. Rev. Appl. 3(3), 037001 (2015). [CrossRef]  

24. N. I. Landy, S. Sajuyigbe, J. J. Mock, D. R. Smith, and W. J. Padilla, “Perfect metamaterial absorber,” Phys. Rev. Lett. 100(20), 207402 (2008). [CrossRef]   [PubMed]  

25. V. A. Fedotov, P. L. Mladyonov, S. L. Prosvirnin, A. V. Rogacheva, Y. Chen, and N. I. Zheludev, “Asymmetric propagation of electromagnetic waves through a planar chiral structure,” Phys. Rev. Lett. 97(16), 167401 (2006). [CrossRef]   [PubMed]  

26. C. Menzel, C. Helgert, C. Rockstuhl, E. B. Kley, A. Tünnermann, T. Pertsch, and F. Lederer, “Asymmetric transmission of linearly polarized light at optical metamaterials,” Phys. Rev. Lett. 104(25), 253902 (2010). [CrossRef]   [PubMed]  

27. A. David, H. Benisty, and C. Weisbuch, “Photonic crystal light-emitting sources,” Rep. Prog. Phys. 75(12), 126501 (2012). [CrossRef]   [PubMed]  

28. D. Iida, A. Fadil, Y. T. Chen, Y. Y. Ou, O. Kopylov, M. Iwaya, T. Takeuchi, S. Kamiyama, I. Akasaki, and H. Y. Ou, “Internal quantum efficiency enhancement of GaInN/GaN quantum-well structures using Ag nanoparticles,” AIP Adv. 5(9), 097169 (2015). [CrossRef]  

29. K. R. Catchpole and A. Polman, “Design principles for particle plasmon enhanced solar cells,” Appl. Phys. Lett. 93(19), 191113 (2008). [CrossRef]  

30. K. R. Catchpole and A. Polman, “Plasmonic solar cells,” Opt. Express 16(26), 21793–21800 (2008). [CrossRef]   [PubMed]  

31. V. A. Fedotov, J. Wallauer, M. Walther, M. Perino, N. Papasimakis, and N. I. Zheludev, “Wavevector selective metasurfaces and tunnel vision filters,” Light Sci. Appl. 4(7), e306 (2015). [CrossRef]  

32. A. Pors and S. I. Bozhevolnyi, “Plasmonic metasurfaces for efficient phase control in reflection,” Opt. Express 21(22), 27438–27451 (2013). [CrossRef]   [PubMed]  

33. Y. T. Chen, T. R. Nielsen, N. Gregersen, P. Lodahl, and J. Mørk, “Finite-element modeling of spontaneous emission of a quantum emitter at nanoscale proximity to plasmonic waveguides,” Phys. Rev. B 81(12), 125431 (2010). [CrossRef]  

34. Y. T. Chen, N. Gregersen, T. R. Nielsen, J. Mørk, and P. Lodahl, “Spontaneous decay of a single quantum dot coupled to a metallic slot waveguide in the presence of leaky plasmonic modes,” Opt. Express 18(12), 12489–12498 (2010). [CrossRef]   [PubMed]  

35. A. E. Miroshnichenko, A. B. Evlyukhin, Y. S. Kivshar, and B. N. Chichkov, “Substrated-induced resonant magneto-electric effects for dielectric nanoparticles,” ACS Photonics , 2(10), 1423–1428 (2015). [CrossRef]  

36. P. Lunnemann, I. Sersic, and A. F. Koenderink, “Optical properties of two-dimensional magnetoelectric point scattering lattices,” Phys. Rev. B 88(24), 245109 (2013). [CrossRef]  

37. A. Kwadrin and A. F. Koenderink, “Probing the electrodynamic local density of states with magnetoelectric point scatterers,” Phys. Rev. B 87(12), 125123 (2013). [CrossRef]  

38. A. Kwadrin, C. I. Osorio, and A. F. Koenderink, “Backaction in metasurface etalons,” Phys. Rev. B 93(10), 104301 (2016). [CrossRef]  

39. A. Eroglu and J. K. Lee, “Simplified formulation of dyadic Green’s functions and their duality relations for general anisotropic media,” Progress in Electromagnetic Research, PIER 77, 391–408 (2007). [CrossRef]  

40. L. Novotny and B. Hecht, Principles of Nano-optics (Cambridge University Press, 2006). [CrossRef]  

41. V. V. Klimov and M. Ducloy, “Quadrupole transitions near an interface: general theory and application to an atom inside a planar cavity,” Phys. Rev. A 72(4), 043809 (2005). [CrossRef]  

42. R. L. Hartman, S. M. Cohen, and P. T. Leung, “A note on the Green dyadic calculation of the decay rates for admolecules at multiple planar interfaces,” J. Chem. Phys. 110(4), 2189–2194 (1999). [CrossRef]  

43. D. H. S. Cheng, “On the formulation of the dyadic Green’s function in a layered medium,” Electromagnetics 6(2), 171–182 (1986). [CrossRef]  

44. H. P. Urbach and G. L. J. A. Rikken, “Spontaneous emission from a dielectric slab,” Phys. Rev. A 57(5), 3913 (1998). [CrossRef]  

45. J. Jung, T. Søndergaard, T. G. Pedersen, K. Pedersen, A. N. Larsen, and B. B. Nielsen, “Dyadic Green’s functions of thin films: applications within plasmonic solar cells,” Phys. Rev. B 83(8), 085419 (2011). [CrossRef]  

46. P. de Vries, D. V. Coevorden, and A. Lagendijk, “Point scatterers for classical waves,” Rev. Mod. Phys. 70(2), 447 (1998). [CrossRef]  

47. C. M. Linton, “Lattice sums for the Helmholtz equation,” SIAM Rev. 52(4), 630–674 (2010). [CrossRef]  

48. F. Capolino, D. R. Jackson, R. W. Donald, and B. F. Leopold, “Comparison of methods for calculating the field excited by a dipole near a 2-D periodic material,” IEEE Trans. Antennas Propagat. 55(6), 1644–1655 (2007). [CrossRef]  

49. O. T. A. Janssen, A. J. H. Wachters, and H. P. Urbach, “Efficient optimization method for the light extraction from periodically modulated LEDs using reciprocity,” Opt. Express 18(24), 24522–24535 (2010). [CrossRef]   [PubMed]  

50. S. R. K. Rodriguez, Y. T. Chen, T. P. Steinbusch, M. A. Verschuuren, A. F. Koenderink, and J. Gómez Rivas, “From weak to strong coupling of localized surface plasmons to guided modes in a luminescent slab,” Phys. Rev. B 90(23), 235406 (2014). [CrossRef]  

51. E. D. Palik, Handbook of Optical Constants of Solids (Academic Press, 1998).

52. C. Qu, S. J. Ma, J. M. Hao, M. Qiu, X. Li, S. Y. Xiao, Z. Q. Miao, N. Dai, Q. He, S. L. Sun, and L. Zhou, “Tailor the functionalities of metasurfaces based on a complete phase diagram,” Phys. Rev. Lett. 115(23), 235503 (2015). [CrossRef]   [PubMed]  

53. A. Kwadrin and A. F. Koenderink, “Diffractive stacks of metamaterial lattices with a complex unit cell: Self-consistent long-range bianisotropic interactions in experiment and theory,” Phys. Rev. B 89(4) 045120 (2014). [CrossRef]  

54. The matlab scripts of calculating T and R for the slab structure with embedded plasmonic sphere lattice in Appendix A can be requested via www.koenderink.info or yuntian@hust.edu.cn