Nonlinear optics of graphene and other 2D materials in layered structures

We present a theoretical framework for nonlinear optics of graphene and other 2D materials in layered structures. We derive a key equation to find the effective electric field and the sheet current density in the 2D material for given incident light beams. Our approach takes into account the effect of the surrounding environment and characterizes its contribution as a structure factor. We apply our approach to two experimental setups, and discuss the structure factors for several nonlinear optical processes including second harmonic generation, third harmonic generation, and parametric frequency conversion. Our systematic study gives a strict extraction method for the nonlinear coefficients, and provides new insights in how layered structures influence the nonlinear signal observed from 2D materials.


I. INTRODUCTION
Graphene possesses extraordinary optical properties, 1 including broadband absorption, 2 the existence of tightly confined plasmon modes, [3][4][5] and extremely strong optical nonlinearities. [6][7][8] Most of these properties strongly depend on the chemical potential, which can be controlled chemically, 9,10 optically, 11 or with the use of an external gate voltage. 12 Since the first experimental demonstration of strong parametric frequency conversion (PFC) in graphene, 13 much attention has been paid to its optical nonlinearities. With the integration of graphene onto photonic chips becoming a mature technology, 1,14,15 graphene is recognized as a potential resource for many photonic devices that exploit optical nonlinearities, including saturable absorbers, 16 broadband optical modulators, [17][18][19] optical switches, 20 and wavelength converters. 15 For these applications to be developed, a full understanding of the nonlinear radiation of graphene is necessary. In general, this is determined by both the intrinsic optical nonlinearity of graphene, and the effects of the surrounding environment. Neither is sufficiently understood at this point.
The intrinsic optical nonlinearity of graphene has been studied extensively. Various experiments with graphene in free-space configurations, on waveguides, and on photonic crystals have been performed to investigate different nonlinear phenomena, including PFC, 13,21 third harmonic generation (THG), [21][22][23][24][25] second harmonic generation (SHG), [26][27][28][29][30][31] Kerr effects (or self-phase modulation, SPM) and two photon absorption, [32][33][34][35][36][37][38][39] and coherent current injection. 40,41 The extracted effective bulk third order susceptibilities χ (3) eff are in the range of 10 −13 − 10 −19 m 2 /V 2 , while the values of the effective bulk second order susceptibilities χ (2) eff are seldom extracted. Recent experiments have shown the unusual negative sign of graphene's nonlinear refractive index, 37,38,42 the ability to tune the nonlinearity by varying the chemical potential, 21,25 and the possibility of having quasiexponentially growing SPM in graphene on waveguides. 39 Most of the existing microscopic theories [43][44][45][46][47][48][49][50][51][52][53][54] focus on the dependence of the sheet conductivities σ (n) (∝ −iωχ (n) eff ) on incident light frequencies, temperature, and doping levels. When scattering is described within a phenomenological relaxation time approximation, the calculated χ (3) eff are about two orders of magnitude smaller than most values extracted from earlier experiments for lightly doped graphene. 43,44,46 Several effects have been considered that might bring theory in better agreement with these experiments, including saturation of the optical nonlinearity, 45,52 the influence of cascaded second order processes, 53,55 and novel plasmonic effects, 53,56 but none of them can sufficiently enhance the calculated conductivities. Recent experiment on light propagating in graphene-covered waveguides 39 has shown that freecarrier refraction plays an important role, and that a treatment beyond perturbation theory is required there. This may be necessary in other excitation configurations as well. However, agreement with perturbative calculations is satisfactory in recent free-space experiments, where the doping level of graphene is tuned over a large range. 21,25 Besides the insufficient understanding of the physical mechanism for the intrinsic optical nonlinearity, the accurate extraction of the effective susceptibility itself from experimental data is also an outstanding problem, especially considering the one-atom thickness. Experimentally, two methods are widely used, both treating graphene as a thin film 13 with an effective thickness d gr ≈ 3.3Å. One is simply to compare the radiation signal to a film with known susceptibility, and the other is to utilize the well-known results for a thin bulk sample in the limit of an undepleted fundamental. 22,57 Neither strategy is based on a strict derivation of the response of a monolayer sample, and the effective susceptibility extracted from either method depends on the artificial thickness assigned to the graphene layer. For usual bulk materials, the generated nonlinear radiation can be worked out using coupled-mode theory, 57 in which a coherence length associated with the phase-matching condition arises in the derivation of the equations for the propagating fields. It has been pointed out that for weakly coupled 2D layers the atomic spacing between the layers should guarantee coherent radiation from the layers. 8,58 Despite many investigations using standard software, 22,59,60 and nonlinear boundary conditions, [61][62][63][64] there is still no systematic treatment of the response of graphene, taking into account its 2D nature, which could later be generalized to treat multilayer samples.
In this work we focus on how the surrounding environment affects the nonlinear radiation of general 2D materials inside layered structures. We model the current density in graphene as a "current sheet" described by a Dirac δ-function, and derive a key equation for graphene nonlinear optics, from which the effective electric fields and current density inside 2D materials can be determined as a response to incident laser beams. These results are easily combined with the transfer matrix formalism, and are used to analyze the nonlinear fields generated by a graphene sheet on multilayered structures. The output fields E out are connected to incident fields E in in a form E out ∝ σ (n) β (n) E n in , with the structure factor β (n) describing the environmental effects. In certain structures, the contribution from the structure factor can be extremely large, providing a controllable way to enhance the generated nonlinear signals. 55,59,60,65,66 The equations we establish can be readily applied to the entire family of 2D materials, including bi-layer graphene, functionalized graphene, monolayer transition-metal dichalcogenides, black phosphorene, silicene, stanene, and so on. We also discuss whether or not a cascaded second order response can modify extracted effective third order nonlinear response coefficients.
We organize this paper as follows: In Sec. II we work out the equations for a suspended 2D layer, and obtain a consistent set of equations for the second and third order nonlinear response; in Sec. III we extend these results to a system with 2D materials embedded inside a multilayered structure; in Sec. IV we apply these results to a graphenecovered multilayered structure. We conclude in Sec. V.

II. SUSPENDED 2D LAYER
We first consider the optics of a suspended monolayer (or nearly-monolayer) structure, which is centered Here E α j;λ (k) are the amplitudes of the α-polarized components of the electric field propagating (or evanescent) in the jth region in the λ direction;ŝ andp± are polarization vectors, with a compact notation k = (κ, ω). The sketch of the polarization vectors and the propagation directions is schematic, since in generalp±can have complex Cartesian components and the fields can be evanescent (See Appendix A). The 2D layer experiences an electric field E eff (k) and possesses a sheet current J(k).
at z = 0, as shown in Fig. 1. In the neighborhood of the monolayer the microscopic electric field and current density have a very complicated spatial dependence, with variations on an atomic scale, and their full description requires a very careful treatment. 67 However, for the calculation of radiation at optical frequencies with wavelengths much larger than these dimensions, a "δ-model" for the current density can be introduced, where k = (κ, ω) indicates the in-plane wave vector κ and frequency ω; for details of the notation see Appendix A.
Here J(k) is a sheet current density. We will see that it mostly appears in the form which has the dimension of an electric field. Within this δ-model all of space is split into three regions, z > 0, z = 0, and z < 0, and both the regions z > 0 and z < 0 are vacuum. The electric field in the latter two regions are described by amplitudes E α j;λ (k), where j = 1, 2 is an index identifying the region (z > 0 and z < 0 respectively), λ = + (−) stands for the upward (downward) propagating or evanescent direction, and α = s (p) identifies the s (p)-polarized component. Using a Green function strategy 68 (see Appendix B), we can connect the fields in the 1st and 2nd regions by with The vectorsê α 0;λ (k) are the polarization vectors for the indicated polarizations α and propagation (or evanescent) directions λ, which are defined in Appendix A; ω ≡ ω/c and w 0 (k) = √ ω 2 − κ 2 . The fields in Eqs. (3) and (4) include any incident fields that are present, from sources or media above (amplitude E α 1;− (k)) or below (amplitude E α 2;+ (k)) the monolayer.
To determine J (k) and thus find the total fields above and below the monolayer one needs to specify the response of the medium to the full field, including that from the monolayer itself. Dealing with theẑ component of the response is here particularly difficult, because thê z component of the field would vary strongly over that microscopic thickness. Largely this strongly varying field would be included in a standard calculation of the response as done in condensed matter physics, and thus would be included in such a calculation of the linear or nonlinear response coefficients; we need to identify an effective driving field E eff (k) that would lead to a calculation of these response coefficients. We follow the standard approach (see, e.g. 63 ) and take E eff (k) to be the average of the fields above and below the 2D layer, The task of determining the dependence of J (k) on E eff (k) is then assigned to condensed matter physics; it has been widely studied both perturbatively 43,44,47,48,53 and numerically. 45,54 Other strategies and their associated definitions of E eff (k) could be considered, but would just shift how much of the problem was relegated to the task of determining J [k; E eff (k)], and how much was relegated to the task of constructing E eff (k) from the incident field and the field from the microscopic current density itself. The advantage of the use of Eq. (6) is that it is a simple expression but still completely contains the radiation reaction field due to the current sheet, so that the effect of E eff (k) on the sheet includes the consequences of its radiation carrying energy away; energy conservation is thus respected.

A. Effective fields and sheet current density inside the 2D layer
Once the functional J [k; E eff ] is known, Eqs. Here is the incident field at the 2D layer (the field without the presence of the 2D layer), and Q 0 (k) = − ω w0ŝ (k)ŝ(k) − w0 ωκκ − κ 2 w0 ωẑẑ is a radiation coefficient tensor matrix associated with the 2D layer. The second term at the right hand side of Eq. (7) gives the onsite fields induced by the sheet current. Equation (7) is our key result; what is also needed to solve for the radiated fields is the identification of the functional J [k; E eff ] from a response theory in condensed matter physics. Using J [k; E eff ] together with Eq. (7) then leads to a self-consistent solution for the fields and the sheet current.
For weak electric fields, the sheet current density can be written into a power expansion of the effective fields. Up to the third order terms, it is Here the tensors η (1) are associated respectively with the linear, second order, and third order conductivities.
The linear conductivity σ (1) (k) has been widely discussed. 69,70 In our calculations it is a good approximation to set σ (1) (1) (ω) because in the microscopic calculation of the conductivity the light wave vector is much smaller than the involved electron wave vectors, and is usually ignored. Furthermore, for most 2D materials, the linear conductivity tensor only has nonzero components σ (1);xx = σ (1);yy = σ (1) (ω) and σ (1);zz = σ (1) is usually ignored, its value is not zero and it is important for certain applications; 5 as well, including it would be useful in characterizing samples. As with σ (1) (k), it is also a good approximation to neglect the dependence of σ (3) on the wave vector and set all κ = 0. While for σ (2) this would hold if the material lacked inversion symmetry, in the presence of such symmetry (such as in graphene), at least the dependence on κ to first order must be kept, which corresponds to keeping electric-quadrupolelike and magnetic-dipole-like contributions to the second order nonlinear response. 49,50,53 Usually the nonlinear term in Eq. (8) is much smaller than the linear term, so it is convenient to isolate the nonlinear response. From Eqs. (7) and (8) the effective field at the 2D layer is formally given by and the sheet current density is withÎ =ŝ(k)ŝ(k) +kk +ẑẑ being the unit dyadic.
In this paper, we focus on the perturbative solutions of Eqs. (10a) and (9). With respect to the effective fields, and The third order term J nl (k) has two contributions: One is from the direct third order processes, which is widely discussed in literature; the other is from the cascaded second order processes. The results of a carefully designed experiment 55 suggest that cascaded second order processes could lead to a significant change in the reflection spectrum due to the generation of surface plasmons.

B. Outgoing fields
Substituting Eq. (10b) into Eqs. (3) and (4), we get the outgoing fields as Here r α (k) and t α (k) are Fresnel coefficients of a suspended 2D layer, and are given by The radiation term from the nonlinear current is Note that in the absence of a nonlinear current the results in Eqs. (12a) and (12b) in terms of the Fresnel coefficients are just what one would expect. For use in the next section it is convenient to write these results in a transfer matrix formalism. We can rewrite Eqs. (12a) and (12b) to give where M α 2D (k) is a transfer matrix for the 2D material layer, and it is given from the Fresnel coefficients as For a s-or p-polarized light beam incident from above, the absorption induced by this suspended layer is given . At normal incidence, where the difference between s− and p−polarized light vanishes, it is 2η In this section we consider the effects of the surrounding environment on the nonlinear optics of a 2D material, with a structure shown in Fig. 2 (a). The incident fields are E α cl;− (k) from the cladding and E α sb;+ (k) from the substrate, while the outgoing fields are E α cl;+ (k) and E α sb;− (k). We only consider the nonlinear optics from the 2D layer; for other layers, nonlinear radiation can be obtained from the coupled mode theory 57 or, in an undepleted pump approximation, by the bulk medium version of these equations 68 . To utilize the transfer matrix formalism, we slightly deform the structure into Fig. 2 (b), anticipating that we will treat the thickness of the vacuum region containing the 2D layer as vanishingly small. Now the expressions in Sec. II can be straightforwardly extended. We denote the transfer matrix for the multilayered structure from the cladding layer to the upper surface of the 2D layer region as M α cu (k), and the transfer matrix for the multilayered structure from the lower surface of the 2D layer region to the substrate as M α ls (k). Then we have Eliminating E α j;± (k) from Eqs. (15), (17a), and (17b), the fields at the cladding connect with the fields at the substrate by a transfer matrix formalism as with the total transfer matrix and a transfer matrix for the nonlinear radiation Converting to a scattering matrix form, we obtain the outgoing fields as where the Fresnel coefficients appearing here are associated with the total transfer matrix M α cs (k) by In the absence of the nonlinear terms the results from Eq. (21) for E α cl;+ (k) and E α sb;− (k) in terms of the Fresnel coefficients are just what one would expect. For the nonlinear radiation, the key quantity is still G α 2D;λ (k), or equivalently the nonlinear current J nl (k), which strongly depend on the total field E eff (k) at the 2D layer. Again, similar to the strategy for a suspended 2D layer, it can be found from Eqs. (3), (4), (6), (17a), and (17b) as where the incident fields E 2D (k) at the 2D layer, and the radiation coefficients Q(k) are given in Appendix C. Equation (23) has the same structure as Eq. (7), but each term is modified to include the effects from the environment, i.e., other layers. All the results from Eq. (10a) to Eq. (11d) can be simply transferred here after replacing E 0 2D (k) and Q 0 (k) by E 2D (k) and Q(k), respectively. For example, comparing to Eq. (10a), the effective field depends on the nonlinear current according to We illustrate our approach by considering a graphenecovered multilayered structure shown in Fig. 3. We assume there is incident light only from the cladding layer, so E α sb;+ (k) = 0, and we are interested in the output field E α cl;+ (k). We ignore any nonlinear sources in the multilayer, cladding, and substrate. In the widely accepted treatment for graphene the optical response in the z direction is ignored, so we put η (1) ⊥ (ω) → 0, and assume as well that the nonlinear current has only inplane components and responds only to the in-plane electric field components. We then require only the inplane components of the effective fields, which we write as E eff; (k) = E eff;s (k)ŝ(k) + E eff;κ (k)κ, and similarly will have only an in-plane nonlinear current J nl; (k). With the z-component of the current sheet thus ignored and the z-component of the effective field not required, Eq. (24) is greatly simplified to give for j = s or κ. Using Eq. (C6a) and (C6b) for E 2D;j (k), the in-plane effective electric fields are found to be and from the transfer matrices the total reflectivity R α cs (k) is found to be Here R α ls (k) is the reflection coefficient without the presence of graphene, τ s = 1, and τ p = −1; in the second expression the result is displayed explicitly in terms of the linear conductivity of the 2D layer. Finally, from Eq. (21) the output fields can then be written into terms of the nonlinear current as According to earlier results, 44,50 the nonlinear currents up to the third order terms are with k a = k − k 1 and k b = k − k 1 − k 2 . For the second order response, S xxyy Q and S xyxy Q are response coefficients associated with the quadrupole-like contribution, and S xxyy M is a response coefficient associated with the magnetic dipole-like contribution; they are related to the parameters S dabc Q/M defined earlier 50 by S = S/(2cǫ 0 ). For later use, the components along the x-direction of the conductivity tensors are given by with S xxxx (ω 1 , ω 2 ) = 2S xxyy Q (ω 1 , ω 2 ) + S xyxy Q (ω 1 , ω 2 ). We note that the argument k = (κ, ω) of E α cl;− (k) varies continuously, and thus the equations can describe a pulse of light with arbitrary polarization incident on only a finite region of the structure shown in Fig. 3. Here we consider simple cases with incident plane waves at discrete k j , which is performed by the transformation All new k j generated by the nonlinearity are then also discrete. In the following two subsections, we explicitly give the perturbation formulas for SHG, THG, and PFC.

A. SHG and THG
In this section we derive the formulas for the SHG and THG signals, induced by a single p−polarized incident laser beam with amplitude E p cl;− (k), as shown in Fig. 3. We list all perturbative quantities with respect to the incident fields. The perturbative effective fields at the graphene layer are given as The perturbative current responses for SHG and THG are Substituting them into Eq. (29), the radiation fields for SHG and THG only include the p-polarized components as with η SHG (k) = η 2 (k, k) and Here all quantities indicated by β are dimensionless structure factors, which include the effects from the lower layers of the structure. They are given by As an example, we consider these structure related quantities for a THG experimental setup as used by Kumar et al.. 22 The layer structure is shown in Fig. 4(a); it is composed of air cladding/graphene/300 nm thick SiO 2 /silicon substrate. The incident light is at ω = 0.72 eV, with an incident angle θ or in-plane wavevector κ = ω sin θx. The dielectric constants we used are listed in Table I, and the linear conductivity of graphene is taken as σ 0 .
unity, the layers underneath reduce the harmonic radiation significantly. The dependence of these structure factors on the substrate reflection coefficient is clear. Equation (36) includes the contribution to the THG coefficients from the cascaded processes, with the values of |P(2k)| smaller than 0.5. By checking graphene's con-ductivities listed in Table I, we find the contribution from the cascaded second order processes to the THG process are about 9 orders of magnitude smaller than the direct third order response. This is not surprising, because the second order response in graphene is induced by theoretically forbidden optical processes and is a weak effect, which gives very small η 2 (k, k) and η 2 (k, 2k). And P(2k) does not undergo any resonance for the structure and parameters used in the experiment.
Following a similar strategy, in the next section we look at parametric frequency conversion. The PFC process occurs when the structure is excited by two laser beams: one is a pump beam at k p , the other is a signal beam at k s . We are interested in the output light at a difference frequency k d = k p − k s , which is generated by the second order nonlinearity, and in the idler light at k i = 2k p − k s , which is generated by the third order nonlinearity. In the experiment of Constant et al., 55 due to a carefully designed structure and appropriate incident angles the light at k d can be evanescent and on resonance with the plasmon modes of the whole structure. As shown in Fig. 5, the structure is composed of graphene covered SiO 2 . For simplicity, we consider that the p-polarized incident beams are in the same incident plane, with in-plane wave vectors κ p = κ pκp and κ s = κ s (−κ p ), and with frequencies ω p > ω s .
In a calculation similar to that of the previous section, we find the output fields at k d and k i are The prefactors 2 and 3 come from the permutation of the incident field components, and the minus signs appear due to the opposite incident directions of these two beams. The terms related to the conductivities are η DFG (k p , k s ) = η 2 (k p , −k s ) and η PFC (k p , k s ) = η 3 (ω p , ω p , −ω s ) + η 3c (k p , k s ). The term η 3c includes the contribution from the cascaded second order processes to the third order response, and is given as Here only the cascaded process involving the light at k d is taken into account. The dimensionless structure factors are Besides involving a product of two second order conductivities, the cascaded process coefficient η 3c (k p , k s ) in Eq. (39) is also proportional to the term P(k d ). This illustrates how the structure underneath the graphene can be used to tune the cascaded process: In the experiment by Constant et al.,55 this term is maximized by matching k d with the surface plasmon polariton (SPP) resonance of the entire structure. In an ideal case, when all losses are ignored, R p cs (k d ) diverges at the resonance; with losses, its value becomes finite, but still very large. As an illustration, this can be clearly seen in Fig. 6(a). For the parameters of graphene we choose the chemical potential as µ = 1.11 eV, the relaxation parameters Γ i = Γ e = 1 meV, and assume zero temperature. The very high chemical potential, which can be realized by a gate voltage, 21 is chosen to satisfy resonant conditions in calculating the second order and third order conductivities. The calculation is made for signal light at wavelength 615 nm ( ω s = 10.22 µm −1 ) with an incident angle θ s = 125 • , and pump light at wavelength λ p varying from 547 nm to 615 nm ( ω p varying from 11.49 µm −1 to 10.22 µm −1 ) with an incident angle θ p = 15 • . As the pump light wavelength varies, the in-plane wave vector and frequency of the light at k d = k p − k s vary along the black line in Fig. 6(a). The linear conductivities for graphene are shown in Fig. 6(c). The maximal value of |R p cs (k d )| along the black line in Fig. 6(a) is about 41, and occurs for ω p = 11.22 µm −1 where the SPP is excited.
We first check the structure factor for PFC, which is shown in Fig. 6(b). As λ p varies, |β PFC | is at the order of 0.05 and changes only slightly. This small value results mostly because of the oblique incidence of the light, as captured by the term w 0 (k)/ ω in the expression of Eq. (33) for P(k). For the difference frequency generation, the structure factor |β DFG | can be as large as 10 at the SPP resonance at ω p = 11.22µm −1 . Obviously, the light at k d is evanescent in air; however, it is could be detected by the use of near field techniques.
Interestingly, the cascaded second order processes can exceed the direct third order process at certain frequencies, mainly because in this case we have very large second order conductivities as well as a SPP resonance, as shown in Fig. 6  (e) ηPFC(kp, ks), η3c(kp, ks), η3(ωp, ωp, −ωs), and P(k d ). All curves indicated by horizontal arrows use the right y-axis. In the calculation, the parameters for graphene are taken as µ = 1.11 meV, Γi = Γe = 1 meV, and T = 0 K. ω p , the first is at ω A = [ ω s + 2µ/( c)]/2 = 10.73µm −1 , corresponding to the condition 2 ω A − ω s = 2µ; the second is at ω B = ω s ; and the third is at ω C = 2µ/( c) = 11.25µm −1 , corresponding to the condition ω C = 2µ. Our calculations show that η 2 (k p , −k s ) shows peaks at ω B and ω C ; η 2 (k p , k d ) shows peaks at ω A , ω B , and ω C ; and η 3 (ω p , ω p , −ω s ) shows peaks at ω A and ω c . For ω p around ω B and ω C , the cascaded second order processes dominate because of the very large second order conductivities. But there is another peak for ω p ∼ 11.22 µm, located at the same peak position of P(k d ), which is induced by the SPP resonance.
All these calculations are for ideal parameters. For higher temperature, larger relaxation parameters, and smaller chemical potential, the second order conductivities can be reduced greatly, as well the value of η 3c (k p , k s ). Because the graphene layer is put in an asymmetric environment, an additional contribution to the second order nonlinearity can arise from the interface effect. However, it contributes only to the out-of-plane tensor components, which are ignored in the usual model for graphene.

V. CONCLUSION
We have derived the basic equations for analyzing the linear and nonlinear optics of graphene and other 2D materials in layered structures. For linear optics, we obtained the Fresnel coefficients of a suspended 2D layer. Complementing existing work in the literature, these coefficients take into account the optical response along the out-of-plane direction of the 2D layers, which is usually ignored, so they can be used to extract the out-of-plane linear conductivity components from experimental data. In the nonlinear regime, we derived an equation for the total electric field at the 2D material including the feedback radiation of its own sheet current density. After combining this with current response theory, which is widely adopted for calculating the conductivities, a set of consistent equations for the total field and the current can be obtained; after solving them either analytically or numerically, the nonlinear radiation can be found. We also derived the perturbative expressions at weak incident light beams. We discussed how the structures themselves affect the output field intensity. For some exper-imental scenarios the structures can reduce the radiation by several orders of magnitude, compared to that for a suspended sample. Our results can serve as a starting point for further nonlinear optical investigations of graphene or 2D materials in layered structures. In a homogeneous medium with dielectric constant ε(ω), the electric fields can be written as with R = xx + yŷ. Throughout this work, we use a compact notation k = (κ, ω), and so E(z; κ, ω) = E(z; k). The fields can be decomposed into upward and downward propagating (or evanescent) fields, with s-and p-polarized components The polarization vectors,ê α λ (k) with α = s, p for the field polarization and λ = +, − for the propagation direction, are defined bŷ e s ± (k) ≡ŝ(k) =κ ×ẑ , where w(k) = ε(ω) ω 2 − κ 2 with ω = ω/c, is taken such that Im[w] ≥ 0, with Re[w] > 0 if Im[w] = 0. The polarization vectors are in general complex quantities; for propagating fields in lossless materials, their directions are shown in Fig. 1. Although these vectors satisfŷ s(k) ·ŝ(k) = 1 andp ± (k) ·p ± (k) = 1, the "length" p ± (k) · [p ± (k)] * can be far larger than unity for an evanescent field.
To identify different regions or layers, we will use an extra subscript j in E α j;λ (z; k), w j (k),ê α j;λ (k), and so on.
Appendix B: Green function method for a sheet current density For a sheet current density J(z; k) embedded in a homogeneous medium with dielectric constant ε(ω), the radiation field can be calculated using a Green function strategy 68 as with In the region z > 0 of the setup in Fig. 1, the induced field E ind (z > 0; k) only propagate upwards, and it can be written as E ind (z > 0; k) = α=s,p F α + (k)e iw(k)zêα + (k) (B3) with F α + given in Eq. (5). Then the upward fields E 1;+ (k) include two contributions: one is from the incident fields E 2;+ (k) that is propagated from the region z < 0, as if the sheet current density were absent; the other is the induced field F α + (k). Then we get Eq. (3). Similarly, we can get Eq. (4).
Appendix C: Incident field E eff (k) at the 2D layer and the radiation coefficients Q(k) For a 2D layer inside a multilayered structure, it is convenient to write the incident field E 2D at the 2D layer in a coordinate formed by {ŝ,κ,ẑ} as E 2D = E 2D;sŝ + E 2D;κκ + E 2D;zẑ . (C1) Without causing any confusion, the k-dependence of each quantity is implicitly shown. From Eqs. (3), (4), (6), (17a), and (17b), the effective field E eff at the 2D layer can be written into Eq. (23) with