General point dipole theory for periodic metasurfaces: magnetoelectric scattering lattices coupled to planar photonic structures

We study semi-analytically the light emission and absorption properties of arbitrary stratified photonic structures with embedded two-dimensional magnetoelectric point scattering lattices, as used in recent plasmon-enhanced LEDs and solar cells. By employing dyadic Green's function for the layered structure in combination with Ewald lattice summation to deal with the particle lattice, we develop an efficient method to study the coupling between planar 2D scattering lattices of plasmonic, or metamaterial point particles, coupled to layered structures. Using the `array scanning method' we deal with localized sources. Firstly, we apply our method to light emission enhancement of dipole emitters in slab waveguides, mediated by plasmonic lattices. We benchmark the array scanning method against a reciprocity-based approach to find that the calculated radiative rate enhancement in k-space below the light cone shows excellent agreement. Secondly, we apply our method to study absorption-enhancement in thin-film solar cells mediated by periodic Ag nanoparticle arrays. Lastly, we study the emission distribution in k-space of a coupled waveguide-lattice system. In particular, we explore the dark mode excitation on the plasmonic lattice using the so-called Array Scanning Method. Our method could be useful for simulating a broad range of complex nanophotonic structures, i.e., metasurfaces, plasmon-enhanced light emitting systems and photovoltaics.


I. 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][3][4][5][6][7][8][9] , and photovoltaics 10, 11 . Current developments in fundamental electromagnetics research particularly focus on plasmonic lattices [9][10][11][12][13][14][15] , metamaterials 16,17 , and metasurfaces [18][19][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,21 . 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 34,35 , 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 22 . 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 23 , for asymmetric transmission 24,25 , 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][7][8][9]26,27 , 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 x Ga 1−x N 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. 1a). 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. 1b). 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,28,29 . Thereby absorption can be enhanced, and cell thicknesses can be reduced (Fig. 1c). In all these scenarios, the periodic plasmon structure is placed in, or on, a nontrivial stratified dielectric system with a complex mode structure.
While the relevance of periodic plasmonic lattices and metasurfaces embedded in stratified systems is large, no simple design or analysis tool is available. Results usually are compared to brute force numerical calculations such as FDTD or FEM modeling [30][31][32][33] . This method provides accurate full-wave simulation results, but involves 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. Considerable efforts have been devoted to developing dipole lattice models 22,34,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 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.

II. 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 function. We focus on stratified systems, meaning that the Green 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 rahter 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 II A introduces the Green function, section II B sets up the infinite-lattice dipole problem in momentum space, and section II C introduces the array scanning method. In section II D we discuss an asymptotic expression for evaluating far field radiation, and for calculating the local density of states via reciprocity.
A. 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, whereḠ (e/h)e denotes the electric/magnetic field Green's tensor from an electric point source, andḠ (e/h)h denotes electric/magnetic field Green's tensor from a magnetic point source. The fields provided by an object with electric current j e (r) and magnetic current j h (r) are given by In this work we focus on scatterers that can be modelled as point dipoles, requiring that they are moderate in size (order λ/2π) and sufficiently spaced (spacing at least the particle radius).
Although the examples that we present all deal with plasmonic scatterers that only have an electric polarizability, for completeness we also incorporate a magnetic dipole response in the model. Thereby, the model also applies to many metasurfaces. For a point scatterer, the induced electric and magnetic dipole moments P and M can be expressed in terms of induced currents through j e (r) = −iω 0 P δ(r − r 0 ), j h (r) = −iω 0 M δ(r − r 0 In this work we focus on lattices embedded in planar photonic structures, consisting of an arbitrary multilayer stack. In this case, to construct the multilayer Green function one uses plane waves, decomposed into TM and TE polarization for each given propagation direction, as specified by parallel momentum k || (momentum component in the plane parallel to the interfaces) as a complete basis to describe any field in space. By convention, this set of plane waves is expressed through a set of basis functions M (k || ) defined as describing the electric field of TE waves, and a dual set (N (k || )) that describes the electric field of TM waves. One set of basis functions can be derived from the other sinceḠ ee =Ḡ hh , the M (k || ) functions not only describe the E-field of TE waves, but also describe the TM-wave H-fields. From the relation between E-field and H-field in source-free regions, G eh (r, r ) = −∇×Ḡ hh (r, r ), one finds ∇×M (k || ) = k 0 N (k || ). Conversely, one can also derive ∇×N (k || ) = k 0 M (k || ), where curl operates on the coordinates of detection point r. This relation facilitates a straightforward computation of all subtensors ofḠ in terms of M (k || ) and N (k || ). The appendices report the functions M (k || ) and N (k || ) for free space (Appendix A, essentially plane waves), and for arbitrary stratified systems. In any multilayer structure, the functions M (k || ) and N (k || ) can be found by applying the transfer matrix formalism, that introduces the wave-vector dependent fresnel coefficients into the problem. Appendix C generalizes the aforementioned reported recipe [41][42][43] for Green's tensorḠ ee to the complete description of the E and H fields induced by a single magneto-electric point scatter.

B. 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 1 + na 2 ) 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 whereᾱ is the dynamic polarizability tensor,ḡ(R l , R l ) = 4π , and where p l (m l ) represents the induced electric (magnetic) dipole moment in the point scatter seated at R l . By using the stratified medium dyadic Green's functionḡ(r, r ), given in Appendix Eq. (B2), we ensure that all multiple scattering interactions between the induced dipole moments p l that are mediated by both propagating and guided modes are included. The driving field e k || (R l ) (h k || (R l )) is understood to be the field solution for the full stratified system in absence of the nanoscatterers in response to an incident plane wave E 0 in of parallel momentum k || , as specified in 44,45 . According to de Vries et al. 46 , the dynamical polarizability tensorᾱ of a point scatter can be written asᾱ = ᾱ −1 static − 4πk 2Ḡ b (r 0 , r 0 ) −1 , whereᾱ static is a quasistatic polarizability (no radiation damping), and whereḠ b is the background Green's tensor at the position of the point scatter itself. A regularization procedure needs to be carried out to remove the singularity inḠ b . Since this regularization is not specific to the stratified system, we refer to Ref. 46 for its derivation. For a point scatterer in free space, the routine boils down to taking solely the imaginary part ofḠ b , which accounts for the radiation damping to fulfill the optical theorem. In a complex system, one corrects the polarizability tensor of a point scatter by adding to the inverse static polarizability 4πk 2 iImḠ f ree (referring to the Green function of homogeneous space filled with the medium directly surrounding the scattering), plus the scattered part 4πk 2Ḡ s (r 0 , r 0 ) of the Green's tensor that contains all the multiple reflections in the multilayer stack, and is regular. Mathematically, one adds to the inverse of static polarizability, i.e., . Returning to Eq. (4), the translation invariance of the coupled lattice-waveguide system implies that a Bloch wave form T can be used to simplify the infinite set of coupled equations to a single equation whereḡ = is the lattice summation of the stratified system Green's function without self-term, given byḡ = = l =0ḡ R 0 is the coordinates of the particle labeled by l = [0, 0] on the lattice. Even for dipole interactions in free space, it is a challenge to calculate the lattice sumḡ = , which can be tackled using Ewald summation techniques. We have adapted the recipe of Linton 47 for summations involving just scalar Green functions in free space to deal also with dyadic Green functions of arbitrary stratified systems. To this end, we reformulate the total Green's tensor of stratified structures asḡ(r, r ) =ḡ 0 (r, r ) +ḡ s (r, r ), wherē g 0 andḡ s are the free space term and the scattered term due to the interfaces, respectively. The explicit form ofḡ 0 andḡ s can be seen in Appendix(A-C). This separation ofḡ(r, r ) leads tō The 'free' term (entire term in {.}) in Eq. (6) can be efficiently calculated using Ewald summation. We use the recipe reported in Ref 36 and 38. At first sight, the summation overḡ s in the third term in Eq. (6) must be much more daunting, since the scattered part of the Green function is formulated as a Sommerfeld integral over parallel wavector containing the M and N modes of the multilayer. Remarkably, the sum of poorly convergent complex integrals simplifies to a rapidly converging sum of just the integrand at select k || that rapidly converges due to the evanescent nature of the higher diffraction orders. Taking one block in g s for an interface as an example, we explicitly formulate lḠ ee s (R 0 , R l )e ik || ·R l for the example of a lattice placed on top of a slab waveguide as follows, Here we refer to the Appendix for a listing of the quantities involved -essentially γ 1,23 s,p accounts for the polarization-dependent fresnel coefficients of the multilayer system, generalized to extend to large k || . The crucial realisation is that in this complicated integrand there is no dependence on R l except in the very last term, i.e., the phase factor e i(k || −q || )·R l . Hence one can pull the summation over lattice sites into the integrand and apply the completeness relation (where G l = mG x + nG y , G x and G y are reciprocal lattice vectors and A denotes the real-space unit cell area) to convert the real space summation of integrals in Eq. (7) into a reciprocal space sum with no complex integrals as follows, In Eq. (8), the sum can be easily truncated (despite the existence of guided modes with a very long in-plane real-space range), due to the fact that each term carries a e i2kzd (d denotes the distance of particle center to its nearest neighboring interface), which exponentially decreases as |k || + G l | is large enough. The same rationale holds for all blocks inḡ s . Once the induced moment p 0 and m 0 is obtained, it is straightforward to calculate the total field at any location as C. 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 k || = J δ(r − (r + R l ))e ik || ·R 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., 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Ḡ whereḡ b (r, r ) is the background dyadic green function (i.e., that of the stratified system in absence of the lattice), the integrand of the second termĪ(k || , r, r 0 ) is given byĪ(k || , r, r 0 ) = g = (r, r 0 )Q(r 0 , r 0 )g = (r 0 , r ), and g = (r 0 , r ) = lḡ (r 0 , r )e ik || ·R l is again a lattice-summed Green function, though now including the self-term, i.e., l = 0 term. The interpretation ofĪ(k || , r, r 0 ) is that the first term gives direct radiation from the source without interaction with the lattice. The second term contains the emission scattered via the lattice, and is the term containing the array scanning. The integrand of second term has three factors, that should be read from right to left: gives the field strength at the particle lattice site r 0 excited by the array of discrete current source J e ik || ·R l located at r + R l , ] −1 quantifies how much the lattice is polarized by the phased array source J e ik || ·R l , 3. g = (r, r 0 ) = lḡ (r, r 0 )e ik || ·R l , describes the field strength at r of field radiated by the polarized lattice.
The calculation ofḡ = (r, r 0 ) andḡ = (r 0 , r ) can be performed by applying the completeness relation as discussed for g = (r 0 , r 0 ). A major difference, however, is that now the source point r 0 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 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.,Ḡ EE,21 slab (r 0 , r ), to illustrate this useful cancellation result, which is given bȳ where T s,ij = 2 kzi+kzj , and T p,ij = 2 √ i j kzi j +kzj i , ψ = (k 2z − k 1z )d. Importantly, such cancellations carry over to complex layered structures, meaning that numerical calculation greatly benefits from judicious ordering of terms.
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.

D. 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, where s = [s x , s y , s z ] = [ x r , y r , z r ] = [ kx k0 , ky k0 , kz k0 ] parametrizes the far-field 'viewing direction', and where E(k x , k y ; z = z r ) is the angular spectrum of the near field at an arbitrarily chosen detection plane z = z r in the half-space that extends to the far field detector. For plane wave incidence with in plane wave vector k || , the only Fourier components of E(k x , k y ; z = z r ) contributing to the far field are at [k x , k y ] = [k || + G l ] ,where l labels the discrete diffraction orders.
In this work we are particularly concerned with mapping the fractional, i.e., parallel wave-vector 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 ,ḡ for radiation at an output angle given by [k x /k 0 , k y /k 0 ] = cosθ[cosφ, sinφ]. Here the termḡ b ([k x , k y , z = z r ], r ) represents the far field emitted directly by the source (no interaction with the lattice). The second term accounts for radiation that appears via interaction with the lattice. It contains three sub-terms multiplicating with each other: the first sub-termḡ = (r 0 , r ) accounts for how strongly the source drives the lattice, the second sub-termQ plays the role of polarizability, and the final sub-termḡ = l=0 quantifies how each particle in the lattice scatters out. Explicit expressions of g = (r 0 , r ) andQ(r 0 , r 0 ) are given in Eq. (10), while the term g = l=0 ([k x , k y , z = z r ], r 0 ) only takes the zero diffraction order, i.e., l = 0. We can use Eq. (13) to calculate the radiative local density of states by taking the sampling points of k || within the light cone, without the need of unfolding k || into corresponding diffraction orders for large pitches (larger than half wavelength), compared to sampling points taken from Brillouin zone.
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 1 kz 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| 2 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 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.

A. 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.
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ᾱ static = (3.78 + 0.135i) · 10 −3 λ 3 0 (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 Fig. 2 (a,b,c,d), and calculated according to reciprocity, i.e., as a near-field enhancement upon far-field TM-wave illumination, as shown in Fig. 2 (e,f,g,h). Results are shown for four different lattice constants, chosen as λ 0 /4 (subdiffractive), λ 0 /2 (diffractive), 3λ 0 /4 (close to 2nd order Bragg diffraction for the waveguide mode) and λ 0 .
In Fig. 2 (a,e), the lattice is effectively a dense, non-diffractive metasurface, as the pitch is chosen so small (λ 0 /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, 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λ0) 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λ0, 1/2λ0, 3/4λ0, λ0 respectively. λ0 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) · 10 −3 λ 3 0 (typical for, e.g., a 80 nm aluminum particle in the green). The calculation is done for a single frequency determined by 2π/λ0, and taking a source height Z0 = 1/10λ0. 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 Fig. 2 (b,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 way 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 Fig. 2 (c,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 III C, 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 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. of the calculation in Fig. 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.

B. 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, 53 . Figure 3 (a,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 amorhous 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.
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 Fig. 3 (a,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.

C. Excitation of dark modes by point sources
For many interesting quantities, plane wave calculations, together with the reciprocity technique suffice. These include angledependent 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 = 116nm) 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 53 ). 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 γ 0 . 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.
To further study the dark mode excitation, we plot the integrands of the scattering part of the Green function, i.e.,Ī(k || , r, r 0 ) in Eq. (10), in one quadrant of the full Brillouin zone in Fig. 4 (d-k) 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 (b-h), blue (black) dotted lines indicate the light line in air (GaAs). Due to C 4 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  TABLE I. 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). /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 (c), the contribution of the waveguide mode to the integrand of scattering part of Green function is weak, as shown in the contour plot of log 10 |Im[G zz ]|. For larger in-plane momentum, the integrand of the scattering part of Green 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 Fig. 4 (d-h) 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 Fig. 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 0 ) 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 0 ), 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 Fig. 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. (I). As evident in Fig. 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 0 , 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 0 , starts to contribute due to the stronger plasmonic hybridization. As the wavelength increases further, the emission with |k || | ≥ n GaAs k 0 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 0 , k 0 < |k || | ≥ n GaAs k 0 , and |k || | ≥ n GaAs k 0 , that ultimately determines how emission occurs in the complex photonic environments 54 . Our method can aid evaluating different designs on basis of the fact that the integrand in the ASM method can be separated in different contributions. To conclude this example, 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.

IV. CONCLUSION
In closing, we presented a general point dipole theory for periodic metasurfaces that are embedded in arbitrary stratified layer systems. With the assistance of Ewald's summation technique and dyadic Green's function of stratified layers, we study magnetoelectric scattering lattices coupled to layered photonic structures. One can with this method rapidly understand the physics of complicated structures by approximating 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.
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