Quantized fluctuational electrodynamics for three-dimensional plasmonic structures

TL;DR

This paper extends the quantized fluctuational electrodynamics (QFED) formalism from one dimension to three dimensions, enabling the analysis of complex plasmonic structures and their local photon properties.

Contribution

The authors develop a 3D QFED formalism that overcomes previous limitations, allowing detailed study of plasmonic and other three-dimensional nanostructures.

Findings

Demonstrated the formalism on a quantum-well with metallic coating

Calculated local photon numbers and field temperatures

Showed applicability to plasmonic structures

Abstract

We recently introduced a quantized fluctuational electrodynamics (QFED) formalism that provides a physically insightful definition of an effective position-dependent photon-number operator and the associated ladder operators. However, this far the formalism has been applicable only for the normal incidence of the electromagnetic field in planar structures. In this work, we overcome the main limitation of the one-dimensional QFED formalism by extending the model to three dimensions, allowing us to use the QFED method to study, e.g., plasmonic structures. To demonstrate the benefits of the developed formalism, we apply it to study the local steady-state photon numbers and field temperatures in a light-emitting near-surface InGaN quantum-well structure with a metallic coating supporting surface plasmons.

Equations

$$\nabla \cdot \mathbf{D}$$

$$\nabla \cdot \mathbf{B}$$

$$\nabla \times \mathbf{E}$$

$$\nabla \times \mathbf{H}$$

$$\displaystyle\nabla\times\Big{(}\frac{\nabla\times\mathbf{E}}{\mu_{0}\mu}\Big{)}-\omega^{2}\varepsilon_{0}\varepsilon\mathbf{E}$$

$${\hat{\mathbf{E}}}^{+}(\mathbf{r},\omega )$$

$$+k_0\int {\overleftrightarrow{\mathbf{G}}}_{\mathrm{em}}(\mathbf{r},\omega ,{\mathbf{r}}^{\prime })\cdot {\hat{\mathbf{J}}}_{\mathrm{m}}^{+}({\mathbf{r}}^{\prime },\omega )d^3{r}^{\prime },$$

$${\hat{\mathbf{H}}}^{+}(\mathbf{r},\omega )$$

$$+i\omega {\epsilon }_0\int {\overleftrightarrow{\mathbf{G}}}_{\mathrm{mm}}(\mathbf{r},\omega ,{\mathbf{r}}^{\prime })\cdot {\hat{\mathbf{J}}}_{\mathrm{m}}^{+}({\mathbf{r}}^{\prime },\omega )d^3{r}^{\prime }.$$

$${\hat{a}}_j(\mathbf{r},\omega )$$

$$\displaystyle\hskip 14.22636pt\times\int\Big{[}\sqrt{\rho_{\mathrm{NL},j,\mathrm{e}}(\mathbf{r},\omega,\mathbf{r}^{\prime})}\hat{f}_{\mathrm{e}}(\mathbf{r}^{\prime},\omega)$$

$$\displaystyle\hskip 14.22636pt+\sqrt{\rho_{\mathrm{NL},j,\mathrm{m}}(\mathbf{r},\omega,\mathbf{r}^{\prime})}\hat{f}_{\mathrm{m}}(\mathbf{r}^{\prime},\omega)\Big{]}d^{3}r^{\prime},$$

$$⟨{\hat{n}}_j(\mathbf{r},\omega )⟩=\frac{\int {\rho }_{\mathrm{NL},j}(\mathbf{r},\omega ,{\mathbf{r}}^{\prime })⟨\hat{\eta }({\mathbf{r}}^{\prime },\omega )⟩d^3{r}^{\prime }}{\int {\rho }_{\mathrm{NL},j}(\mathbf{r},\omega ,{\mathbf{r}}^{\prime })d^3{r}^{\prime }},$$

$${\rho }_{\mathrm{NL},\mathrm{e}}(\mathbf{r},\omega ,{\mathbf{r}}^{\prime })$$

$$\displaystyle=\frac{2\omega^{3}}{\pi c^{4}}\Big{(}\varepsilon_{\mathrm{i}}(\mathbf{r}^{\prime},\omega)\mathrm{Tr}[\overset{\text{\tiny$\leftrightarrow$}}{\mathbf{G}}_{\mathrm{ee}}(\mathbf{r},\omega,\mathbf{r}^{\prime})\cdot\overset{\text{\tiny$\leftrightarrow$}}{\mathbf{G}}_{\mathrm{ee}}^{\dagger}(\mathbf{r},\omega,\mathbf{r}^{\prime})]$$

$$\displaystyle\hskip 14.22636pt+\mu_{\mathrm{i}}(\mathbf{r}^{\prime},\omega)\mathrm{Tr}[\overset{\text{\tiny$\leftrightarrow$}}{\mathbf{G}}_{\mathrm{em}}(\mathbf{r},\omega,\mathbf{r}^{\prime})\cdot\overset{\text{\tiny$\leftrightarrow$}}{\mathbf{G}}_{\mathrm{em}}^{\dagger}(\mathbf{r},\omega,\mathbf{r}^{\prime})]\Big{)},$$

$${\rho }_{\mathrm{NL},\mathrm{m}}(\mathbf{r},\omega ,{\mathbf{r}}^{\prime })$$

$$\displaystyle=\frac{2\omega^{3}}{\pi c^{4}}\Big{(}\varepsilon_{\mathrm{i}}(\mathbf{r}^{\prime},\omega)\mathrm{Tr}[\overset{\text{\tiny$\leftrightarrow$}}{\mathbf{G}}_{\mathrm{me}}(\mathbf{r},\omega,\mathbf{r}^{\prime})\cdot\overset{\text{\tiny$\leftrightarrow$}}{\mathbf{G}}_{\mathrm{me}}^{\dagger}(\mathbf{r},\omega,\mathbf{r}^{\prime})]$$

$$\displaystyle\hskip 14.22636pt+\mu_{\mathrm{i}}(\mathbf{r}^{\prime},\omega)\mathrm{Tr}[\overset{\text{\tiny$\leftrightarrow$}}{\mathbf{G}}_{\mathrm{mm}}(\mathbf{r},\omega,\mathbf{r}^{\prime})\cdot\overset{\text{\tiny$\leftrightarrow$}}{\mathbf{G}}_{\mathrm{mm}}^{\dagger}(\mathbf{r},\omega,\mathbf{r}^{\prime})]\Big{)},$$

$${\rho }_{\mathrm{NL},\mathrm{tot}}(\mathbf{r},\omega ,{\mathbf{r}}^{\prime })$$

$$=\frac{|\epsilon (\mathbf{r},\omega )|}{2}{\rho }_{\mathrm{NL},\mathrm{e}}(\mathbf{r},\omega ,{\mathbf{r}}^{\prime })+\frac{|\mu (\mathbf{r},\omega )|}{2}{\rho }_{\mathrm{NL},\mathrm{m}}(\mathbf{r},\omega ,{\mathbf{r}}^{\prime }).$$

$${\rho }_j(\mathbf{r},\omega )=\int {\rho }_{\mathrm{NL},j}(\mathbf{r},\omega ,{\mathbf{r}}^{\prime })d^3{r}^{\prime }.$$

$${\rho }_j(\mathbf{r},\omega )$$

$$⟨\hat{E}(\mathbf{r},t{)}^2{⟩}_{\omega }$$

$$⟨\hat{H}(\mathbf{r},t{)}^2{⟩}_{\omega }$$

$$⟨\hat{u}(\mathbf{r},t){⟩}_{\omega }$$

$${\mathbf{\rho }}_{\mathrm{IF}}(\mathbf{r},\omega ,{\mathbf{r}}^{\prime })$$

$$=\frac{2{\omega }^3{n}_{\mathrm{r}}(\mathbf{r},\omega )}{\pi c^4}$$

$$\displaystyle\hskip 14.22636pt\times\Big{(}\mu_{\mathrm{i}}(\mathbf{r}^{\prime},\omega)\mathrm{Im}\Big{[}\mathrm{Tr}[\overset{\text{\tiny$\leftrightarrow$}}{\mathbf{G}}_{\mathrm{mm}}(\mathbf{r},\omega,\mathbf{r}^{\prime})\times\overset{\text{\tiny$\leftrightarrow$}}{\mathbf{G}}_{\mathrm{em}}^{\dagger}(\mathbf{r},\omega,\mathbf{r}^{\prime})]\Big{]}$$

$$\displaystyle\hskip 14.22636pt-\varepsilon_{\mathrm{i}}(\mathbf{r}^{\prime},\omega)\mathrm{Im}\Big{[}\mathrm{Tr}[\overset{\text{\tiny$\leftrightarrow$}}{\mathbf{G}}_{\mathrm{ee}}(\mathbf{r},\omega,\mathbf{r}^{\prime})\times\overset{\text{\tiny$\leftrightarrow$}}{\mathbf{G}}_{\mathrm{me}}^{\dagger}(\mathbf{r},\omega,\mathbf{r}^{\prime})]\Big{]}\Big{)},$$

$$⟨\hat{\mathbf{S}}(\mathbf{r},t){⟩}_{\omega }=\hslash \omega v(\mathbf{r},\omega )\int {\mathbf{\rho }}_{\mathrm{IF}}(\mathbf{r},\omega ,{\mathbf{r}}^{\prime })⟨\hat{\eta }({\mathbf{r}}^{\prime },\omega )⟩d^3{r}^{\prime },$$

$${\hat{Q}}_{\mathrm{em}}(\mathbf{r},t)$$

$${\hat{Q}}_{\mathrm{abs}}(\mathbf{r},t)$$

$$\hat{Q}(\mathbf{r},t)=$$

$$⟨\hat{Q}(\mathbf{r},t){⟩}_{\omega }$$

$$=\hslash {\omega }^2{\epsilon }_{\mathrm{i}}(\mathbf{r},\omega ){\rho }_{\mathrm{e}}(\mathbf{r},\omega )[⟨\hat{\eta }(\mathbf{r},\omega )⟩-⟨{\hat{n}}_{\mathrm{e}}(\mathbf{r},\omega )⟩]$$

$$+\hslash {\omega }^2{\mu }_{\mathrm{i}}(\mathbf{r},\omega ){\rho }_{\mathrm{m}}(\mathbf{r},\omega )[⟨\hat{\eta }(\mathbf{r},\omega )⟩-⟨{\hat{n}}_{\mathrm{m}}(\mathbf{r},\omega )⟩].$$

$$\displaystyle\nabla_{\mathbf{r}}\times\Big{(}\frac{\nabla_{\mathbf{r}}\times\overset{\text{\tiny$\leftrightarrow$}}{\mathbf{G}}_{\mathrm{ee}}(\mathbf{r},\omega,\mathbf{r}^{\prime})}{\mu(\mathbf{r},\omega)}\Big{)}-k_{0}^{2}\varepsilon(\mathbf{r},\omega)\overset{\text{\tiny$\leftrightarrow$}}{\mathbf{G}}_{\mathrm{ee}}(\mathbf{r},\omega,\mathbf{r}^{\prime})$$

$$=\overleftrightarrow{\mathbf{I}}\delta (\mathbf{r}-{\mathbf{r}}^{\prime }),$$

$${\hat{\mathbf{E}}}^{+}(\mathbf{r},\omega )$$

$$\displaystyle\hskip 14.22636pt-\nabla_{\mathbf{r}^{\prime}}\times\Big{(}\frac{\hat{\mathbf{J}}_{\mathrm{m}}^{+}(\mathbf{r}^{\prime},\omega)}{\mu_{0}\mu(\mathbf{r}^{\prime},\omega)}\Big{)}\Big{]}d^{3}r^{\prime}$$

$$=i\omega {\mu }_0\int {\overleftrightarrow{\mathbf{G}}}_{\mathrm{ee}}(\mathbf{r},\omega ,{\mathbf{r}}^{\prime })\cdot {\hat{\mathbf{J}}}_{\mathrm{e}}^{+}({\mathbf{r}}^{\prime },\omega )d^3{r}^{\prime }$$

$$+k_0\int {\overleftrightarrow{\mathbf{G}}}_{\mathrm{em}}(\mathbf{r},\omega ,{\mathbf{r}}^{\prime })\cdot {\hat{\mathbf{J}}}_{\mathrm{m}}^{+}({\mathbf{r}}^{\prime },\omega )d^3{r}^{\prime },$$

Full text

I Introduction

The quantum optical processes in lossy and lossless material systems have been widely studied during the last few decades. This has led to advances, e.g., in nanoplasmonics Sorger2011 ; Oulton2009 ; Huang2013 ; Sadi2013a , near-field microscopy Taubner2006 ; Hillenbrand2002 , thin-film light-emitting diodes Nakamura2013 ; Heikkila2013 , photonic crystals Russell2003 ; Akahane2003 , and metamaterials Tanaka2010 ; Mattiucci2013 . For describing spatial field evolution in resonant structures, one of the most widely used quantization approaches has been the input-output formalism of the photon creation and annihilation operators. The formalism was originally developed for describing lossless and dispersionless dielectrics Knoll1987 and was later extended for lossy and dispersive media Knoll1991 ; Allen1992 ; Huttner1992 ; Barnett1995 ; Matloob1995 ; Matloob1996 . The early studies clearly highlighted that the noise and field operators in nonuniform systems are position dependent and that the vector potential and electric-field operators obey the well-known canonical commutation relation as expected Barnett1995 ; Matloob1995 . However, the canonical commutation relations did not extend to the photon creation and annihilation operators, which were found to exhibit anomalies in resonant structures Ueda1994 ; Raymer2013 ; Barnett1996 ; Aiello2000 ; Stefano2000 . It was first concluded that these anomalies have no physical significance. Although the formalism was later successfully used to study, e.g., amplifying media and spontaneous decay in left-handed media Dung2003 ; Khanbekyan2005 ; Raabe2007 ; Raabe2008 , it was also recently shown that the anomalous commutation relations should, e.g., lead to the existence of a threshold for second-harmonic generation inside microcavities Gauvin2014 ; Collette2013 . The anomalous commutation relations have also been found to prevent a systematic description of local thermal balance between the field and interacting media Partanen2014a ; Partanen2015a .

We recently solved the cavity commutation relation anomaly and photon-number problem by introducing a quantized fluctuational electrodynamics (QFED) model to describe photon number and showed that the expectation values of the properly normalized annihilation and creation operators result in a meaningful photon-number model and thermal balance conditions Partanen2014a ; Partanen2014c ; Partanen2015a ; Partanen2016b . This far, our models have been strictly one-dimensional and limited to normal incidence in planar structures, which has provided an adequate framework for describing the fundamental properties of cavity fields. However, considering the associated transparent description of the photon number and field temperature, it becomes reasonable to ask how the description can be expanded to more complex systems involving, e.g., plasmons that have been of great topical interest Barnes2003 ; Okamoto2004 ; Yeh2008 ; Bonnand2006a ; Bonnand2006b ; Tanaka2010 ; Pitarke2007 and whose description could benefit from the new methodology clearly separating the local density of states (LDOS) and the photon number. Here we therefore present a generalized QFED model to account for fully three-dimensional propagation as well as the associated spectral expansion for planar structures. We also demonstrate the benefits of the formalism by applying it to study the local steady-state field properties and plasmonic interactions in a light-emitting near-surface InGaN quantum-well (QW) structure with a metallic coating supporting surface plasmons (SPs).

This paper is organized as follows: The theory of the QFED method is presented in Sec. II. As a background for QFED, we first review the general three-dimensional noise operator formalism and the use of Green’s functions to obtain the solutions of the electromagnetic (EM) fields. This is followed by a presentation of the new contribution to the theory: the properly normalized position- and frequency-dependent photon ladder operators, the related photon-number presentation, and the generalized forms of the densities of states in the QFED method. After introducing the ladder and number operators, we also briefly review how the operators can be used to present the associated models for the field fluctuations, Poynting vector, and absorption and emission operators. Note that the expectation values of these macroscopic field quantities are equivalent to the values obtained by using the conventional fluctuational electrodynamics. In Sec. III, we demonstrate the applicability and study the physical implications of the introduced QFED method in an example InGaN QW geometry.

II Quantized fluctuational electrodynamics method

In this section, we outline the derivation of the three-dimensional QFED theory. Detailed derivations are given in the appendixes. We start by introducing the fundamental equations of the conventional fluctuational electrodynamics theory and its quantization in Sec. II.1 and the solution of fields using the dyadic Green’s functions in Sec. II.2. Then, in Sec. II.3, we present the properly normalized photon ladder operators and the related photon-number and density-of-states concepts that expand the classical and previously used quantized versions of fluctuational electrodynamics to enable an unambiguous photon-level description of the three-dimensional system. In Secs. II.4 and II.5, we focus on calculating the Poynting vector operator and the thermal balance predicted by the quantized theory using the newly established operators.

II.1 Noise operator formalism

Maxwell’s equations describe electric and magnetic fields generated by currents and charges in matter. They relate the electric field strength $\mathbf{E}$, the magnetic field strength $\mathbf{H}$, the electric flux density $\mathbf{D}$, and the magnetic flux density $\mathbf{B}$ to the free electric charge density $\rho_{\mathrm{f}}$ and the free electric current density $\mathbf{J}_{\mathrm{f}}$. In the frequency domain, Maxwell’s equations are written for positive frequencies as Partanen2016b

[TABLE]

Here we have additionally used the constitutive relations $\mathbf{D}=\varepsilon_{0}\varepsilon\mathbf{E}+\delta\mathbf{P}$ and $\mathbf{B}=\mu_{0}(\mu\mathbf{H}+\delta\mathbf{M})$, where $\varepsilon_{0}$ and $\mu_{0}$ are the permittivity and permeability of vacuum, $\varepsilon=\varepsilon_{\mathrm{r}}+i\varepsilon_{\mathrm{i}}$ and $\mu=\mu_{\mathrm{r}}+i\mu_{\mathrm{i}}$ are the relative permittivity and permeability of the medium with real and imaginary parts denoted by subscripts $\mathrm{r}$ and $\mathrm{i}$, and the polarization and magnetization fields $\delta\mathbf{P}$ and $\delta\mathbf{M}$ denote the polarization and magnetization that are not linearly proportional to the respective field strengths Sipe1987 . In the context of the fluctuational electrodynamics and the present work, $\delta\mathbf{P}$ and $\delta\mathbf{M}$ describe the small thermal fluctuations of the linear polarization and magnetization fields Partanen2016b . For the remainder of this work, the current density of free charges $\mathbf{J}_{\mathrm{f}}$ is also included in the electric permittivity for notational simplicity.

From Maxwell’s equations in Eqs. (1)–(4) it follows that the electric field obeys the well-known equation Partanen2016b

[TABLE]

where the terms $\mathbf{J}_{\mathrm{e}}=-i\omega\delta\mathbf{P}$ and $\mathbf{J}_{\mathrm{m}}=-i\omega\mu_{0}\delta\mathbf{M}$ represent the polarization and magnetization currents that act as field sources in the noise operator theory Partanen2016b and in the classical fluctuational electrodynamics Narayanaswamy2014 ; Polimeridis2015 . The electric term $\mathbf{J}_{\mathrm{e}}$ includes contributions from both the electric currents due to free charges and polarization terms associated with dipole currents and thermal dipole fluctuations. For the magnetic current term $\mathbf{J}_{\mathrm{m}}$, the only contribution arises from the magnetic dipoles. Note that, after solving the electric field from Eq. (5), the calculation of the magnetic field is straightforward using Faraday’s law in Eq. (3).

In the previously known noise operator framework, we use the canonical quantization of the above classical equations where the classical field quantities are replaced by corresponding quantum operators Dung2003 ; Matloob1995 ; Partanen2014a ; Partanen2016b . The electric and magnetic noise current operators $\hat{\mathbf{J}}_{\mathrm{e}}^{+}(\mathbf{r},\omega)$ and $\hat{\mathbf{J}}_{\mathrm{m}}^{+}(\mathbf{r},\omega)$ are written in terms of bosonic source field operators $\hat{f}_{\mathrm{e}}(\mathbf{r},\omega)$ and $\hat{f}_{\mathrm{m}}(\mathbf{r},\omega)$ as $\hat{\mathbf{J}}_{\mathrm{e}}^{+}(\mathbf{r},\omega)=\sum_{\alpha}j_{\mathrm{0,e}}(\mathbf{r},\omega)\hat{\mathbf{e}}_{\alpha}\hat{f}_{\mathrm{e}}(\mathbf{r},\omega)$ and $\hat{\mathbf{J}}_{\mathrm{m}}^{+}(\mathbf{r},\omega)=\sum_{\alpha}j_{\mathrm{0,m}}(\mathbf{r},\omega)\hat{\mathbf{e}}_{\alpha}\hat{f}_{\mathrm{m}}(\mathbf{r},\omega)$, where $\hat{\mathbf{e}}_{\alpha}$ are unit vectors for the three coordinate directions $\alpha\in\{x,y,z\}$ Partanen2016b ; Partanen2014a . The operators $\hat{f}_{\mathrm{e}}(\mathbf{r},\omega)$ and $\hat{f}_{\mathrm{m}}(\mathbf{r},\omega)$ obey the canonical commutation relation $[\hat{f}_{j}(\mathbf{r},\omega),\hat{f}_{k}^{\dagger}(\mathbf{r}^{\prime},\omega^{\prime})]=\delta_{jk}\delta(\mathbf{r}-\mathbf{r}^{\prime})\delta(\omega-\omega^{\prime})$, with $j,k\in\{\mathrm{e,m}\}$. The normalization factors $j_{\mathrm{0,e}}(\mathbf{r},\omega)$ and $j_{\mathrm{0,m}}(\mathbf{r},\omega)$ have been determined to be $j_{\mathrm{0,e}}(\mathbf{r},\omega)=\sqrt{4\pi\hbar\omega^{2}\varepsilon_{0}\varepsilon_{\mathrm{i}}(\mathbf{r},\omega)}$ Partanen2014a ; Khanbekyan2005 and $j_{\mathrm{0,m}}(\mathbf{r},\omega)=\sqrt{4\pi\hbar\omega^{2}\mu_{0}\mu_{\mathrm{i}}(\mathbf{r},\omega)}$ Partanen2016b ; Dung2003 .

II.2 Green’s functions

In order to write the solution of Eq. (5) in a general form, we apply the conventional dyadic Green’s function formalism Eckhardt1984 ; Paulus2000 , where the field solutions are written as

[TABLE]

[TABLE]

Here \overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{\mathrm{ee}}(\mathbf{r},\omega,\mathbf{r}^{\prime}) is the electric Green’s function, \overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{\mathrm{mm}}(\mathbf{r},\omega,\mathbf{r}^{\prime}) is the magnetic Green’s function, and \overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{\mathrm{em}}(\mathbf{r},\omega,\mathbf{r}^{\prime}) and \overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{\mathrm{me}}(\mathbf{r},\omega,\mathbf{r}^{\prime}) are the exchange Green’s functions. For completeness, the relations between these Green’s functions are shown in Appendix A, and the Green’s functions are explicitly presented for stratified media in Appendix B.

II.3 Photon numbers and densities of states

In analogy with the one-dimensional QFED formalism Partanen2016b , we can define the position- and frequency-dependent effective photon ladder operators $\hat{a}_{j}(\mathbf{r},\omega)$, which obey the canonical commutation relation $[\hat{a}_{j}(\mathbf{r},\omega),\hat{a}_{j}^{\dagger}(\mathbf{r},\omega)]=\delta(\omega-\omega^{\prime})$, for the electric, magnetic, and total EM fields, $j\in\{\mathrm{e,m,tot}\}$. These operators and the corresponding effective photon-number expectation values $\langle\hat{n}_{\mathrm{j}}(\mathbf{r},\omega^{\prime})\rangle$ are given by

[TABLE]

[TABLE]

where $\langle\hat{\eta}(\mathbf{r}^{\prime},\omega)\rangle$ is the source-field photon-number expectation value related to the bosonic noise operators as $\langle\hat{\eta}(\mathbf{r},\omega)\rangle=\int\langle\hat{f}_{\mathrm{e}}^{\dagger}(\mathbf{r},\omega)\hat{f}_{\mathrm{e}}(\mathbf{r}^{\prime},\omega^{\prime})\rangle d^{3}r^{\prime}d\omega^{\prime}=\int\langle\hat{f}_{\mathrm{m}}^{\dagger}(\mathbf{r},\omega)\hat{f}_{\mathrm{m}}(\mathbf{r}^{\prime},\omega^{\prime})\rangle d^{3}r^{\prime}d\omega^{\prime}$ and $\rho_{\mathrm{NL},j}(\mathbf{r},\omega,\mathbf{r}^{\prime})$ are the nonlocal densities of states (NLDOSs) for the electric, magnetic, and total EM fields, written as

[TABLE]

The NLDOS components $\rho_{\mathrm{NL},j,\mathrm{e}}(\mathbf{r},\omega,\mathbf{r}^{\prime})$ and $\rho_{\mathrm{NL},j,\mathrm{m}}(\mathbf{r},\omega,\mathbf{r}^{\prime})$, with $j\in\{\mathrm{e,m}\}$, in Eq. (8) denote, respectively, the first and the second terms of Eqs. (10) and (11). The total NLDOS terms $\rho_{\mathrm{NL,tot,e}}(\mathbf{r},\omega,\mathbf{r}^{\prime})$ and $\rho_{\mathrm{NL,j,m}}(\mathbf{r},\omega,\mathbf{r}^{\prime})$ are calculated by using Eq. (12) with the corresponding terms in the electric and magnetic NLDOSs.

Note that the expressions for the photon ladder operators and the photon numbers in Eqs. (8) and (9) are the same as the expression in the one-dimensional formalism Partanen2016b , but the NLDOSs are different. The derivation of these NLDOSs is presented in Appendix C, and for general stratified media, the densities of states are presented in Appendix D. The LDOSs $\rho_{j}(\mathbf{r},\omega)$ are given in terms of the NLDOSs by

[TABLE]

It is well-known that, in vacuum, the imaginary parts of the traces of the dyadic Green’s functions give the electric and magnetic LDOSs $\rho_{\mathrm{e}}(\mathbf{r},\omega)$ and $\rho_{\mathrm{m}}(\mathbf{r},\omega)$ as Joulain2003 ; Joulain2005

[TABLE]

where $j\in\{\mathrm{e,m}\}$. A similar relation also applies for the normal components of the Fourier-transformed quantities in layered media Partanen2014c ; Partanen2016b , and typically, also, the spatially resolved form in Eq. (14) is expected to be valid inside lossy media. However, in lossy media, these LDOSs are generally known to become infinite due to the contribution of evanescent waves Joulain2003 ; Joulain2005 .

In terms of the photon-number expectation values in Eq. (9) and the LDOSs in Eq. (13), the spectral electric and magnetic field fluctuations and the energy density are given by Partanen2014c

[TABLE]

Here the subscript $\omega$ denotes the contribution of $\omega$ to the total quantities which are obtained as integrals over positive frequencies.

II.4 Quantized Poynting vector operator

To conform with our earlier works and to enable describing energy flow in detail, we also find the three-dimensional generalized expression for the Poynting vector. For an optical mode, the quantum optical Poynting vector is defined as a normal-ordered operator in terms of the positive- and negative-frequency parts of the electric and magnetic field operators as $\hat{\mathbf{S}}(\mathbf{r},t)=:\!\hat{\mathbf{E}}(\mathbf{r},t)\times\hat{\mathbf{H}}(\mathbf{r},t)\!:=\hat{\mathbf{E}}^{-}(\mathbf{r},t)\times\hat{\mathbf{H}}^{+}(\mathbf{r},t)-\hat{\mathbf{H}}^{-}(\mathbf{r},t)\times\hat{\mathbf{E}}^{+}(\mathbf{r},t)$ Loudon2000 . As detailed in Appendix C and in analogy with the one-dimensional QFED formalism Partanen2014a ; Partanen2016b , we obtain the Poynting-vector-related interference density of states (IFDOS) as

[TABLE]

where $n_{\mathrm{r}}(\mathbf{r},\omega)$ is the real part of the refractive index and we have used the short-hand notation \mathrm{Tr}[\overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{1}(\mathbf{r},\omega,\mathbf{r}^{\prime})\times\overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{2}^{\dagger}(\mathbf{r},\omega,\mathbf{r}^{\prime})]=\sum_{\alpha}[\overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{1}(\mathbf{r},\omega,\mathbf{r}^{\prime})\cdot\hat{\mathbf{e}}_{\alpha}]\times[\overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{2}(\mathbf{r},\omega,\mathbf{r}^{\prime})\cdot\hat{\mathbf{e}}_{\alpha}]^{\dagger}, which is a vector, in contrast to the conventional trace of a matrix. Using the IFDOS, the Poynting vector is given by

[TABLE]

where $v(\mathbf{r},\omega)=c/n_{\mathrm{r}}(\mathbf{r},\omega)$ is the propagation velocity of the field in the direction of the wave vector. The integral of the IFDOS with respect to $\mathbf{r}^{\prime}$ is always zero, i.e., $\int\boldsymbol{\rho}_{\mathrm{IF}}(\mathbf{r},\omega,\mathbf{r}^{\prime})d^{3}r^{\prime}=0$, which is required by the fact that, in a medium in thermal equilibrium, there is no net energy flow. For stratified media, the IFDOS is presented in Appendix D.

II.5 Field-matter interaction operators and thermal balance

A particularly insightful view of the effective photon numbers in the QFED framework is provided by their connection to local thermal balance between the field and matter Partanen2014a . First, we define the normal-ordered emission and absorption operators $\hat{Q}_{\mathrm{em}}(\mathbf{r},t)$ and $\hat{Q}_{\mathrm{abs}}(\mathbf{r},t)$ as

[TABLE]

where the electric and magnetic absorption current operators $\hat{\mathbf{J}}_{\mathrm{e,abs}}(\mathbf{r},t)$ and $\hat{\mathbf{J}}_{\mathrm{m,abs}}(\mathbf{r},t)$ are written in the spectral domain as $\hat{\mathbf{J}}_{\mathrm{e,abs}}^{+}(\mathbf{r},\omega)=-i\omega\varepsilon_{0}\chi_{\mathrm{e}}(\mathbf{r},\omega)\hat{\mathbf{E}}^{+}(\mathbf{r},\omega)$ and $\hat{\mathbf{J}}_{\mathrm{m,abs}}^{+}(\mathbf{r},\omega)=-i\omega\mu_{0}\chi_{\mathrm{m}}(\mathbf{r},\omega)\hat{\mathbf{H}}^{+}(\mathbf{r},\omega)$, where $\chi_{\mathrm{e}}(\mathbf{r},\omega)=\varepsilon(\mathbf{r},\omega)-1$ and $\chi_{\mathrm{m}}(\mathbf{r},\omega)=\mu(\mathbf{r},\omega)-1$ are the electric and magnetic susceptibilities of the medium.

The net emission operator $\hat{Q}(\mathbf{r},t)=\hat{Q}_{\mathrm{em}}(\mathbf{r},t)-\hat{Q}_{\mathrm{abs}}(\mathbf{r},t)$, which describes the energy transfer between the EM field and the local medium, is given by

[TABLE]

where $\hat{\mathbf{J}}_{\mathrm{e,tot}}(\mathbf{r},t)=\hat{\mathbf{J}}_{\mathrm{e}}(\mathbf{r},t)+\hat{\mathbf{J}}_{\mathrm{e,abs}}(\mathbf{r},t)$ and $\hat{\mathbf{J}}_{\mathrm{m,tot}}(\mathbf{r},t)=\hat{\mathbf{J}}_{\mathrm{m}}(\mathbf{r},t)+\hat{\mathbf{J}}_{\mathrm{m,abs}}(\mathbf{r},t)$ correspond to the classical total current densities, which are sums of free and bound current densities. The spectral component of the expectation value of the net emission operator in Eq. (22) can be written in terms of the LDOSs and the electric- and magnetic-field photon numbers as

[TABLE]

This shows that local thermal balance [$\langle\hat{Q}(\mathbf{r},t)\rangle_{\omega}=0$] is generally reached when the source-field photon numbers coincide with the field photon numbers as defined in Eq. (23). In addition, the net emission operator satisfies $\langle\hat{Q}(\mathbf{r},t)\rangle_{\omega}=\nabla\cdot\langle\hat{\mathbf{S}}(\mathbf{r},t)\rangle_{\omega}$. In resonant systems where the energy exchange is dominated by a narrow frequency band, the condition $\langle\hat{Q}(\mathbf{r},t)\rangle_{\omega}=0$ can be used to approximately determine the steady-state temperature of a weakly interacting resonant particle Bohren1998 .

III Results

We apply the QFED formalism presented in Sec. II to the study of an example plasmonic multilayer structure, which has recently been of experimental and theoretical interest Homeyer2013 ; Sadi2013b . In contrast to the previous QFED models, the generalized model can, among other oblique-angle problems, describe the optical properties of plasmonic devices. Here we study the contribution of the evanescent SP modes to the position-dependent LDOSs and the effective-field temperatures in the vicinity of a light-emitting Ag/GaN/In0.15Ga0.85N/GaN/Al2O3 multilayer structure illustrated in Fig. 1. The 2-nm In0.15Ga0.85N QW has a band gap of 2.76 eV ($\lambda=450$ nm), and it acts as the emitter layer. It is deposited 20 nm below the 20-nm silver layer which supports SP modes. The refractive indices of GaN and InN are taken from Refs. 52; 53; 54; 55; 56; 57, and the refractive index of In0.15Ga0.85N is deduced by using Vegard’s law; the refractive index of silver is calculated by using the Drude model plasma frequency $\omega_{\mathrm{p}}=9.04$ eV/$\hbar$ and damping frequency $\omega_{\tau}=0.02125$ eV/$\hbar$ taken from Ref. 58, and the refractive index of sapphire is taken from Ref. 59. For example, in the case of the photon energy $\hbar\omega=2.76$ eV corresponding to the QW band gap, the refractive indices of air, silver, GaN, In0.15Ga0.85N, and sapphire are 1.00, $0.013+3.119i$, $2.51+0.0029i$, $2.51+0.094i$, and 1.78, respectively.

The background temperature of the materials is $T=300$ K. We compare the emission of the structure in two cases: (1) the QW is thermally excited to temperature $T_{\mathrm{ex}}=350$ K, and (2) the QW is electrically or optically excited to a state corresponding to direct excitation by a $U=2.6$ V voltage source. In the first case, the QW source-field photon-number expectation value is modeled using the Bose-Einstein distribution $\langle\hat{\eta}_{\text{\tiny QW}}\rangle=1/(e^{\hbar\omega/(k_{\mathrm{B}}T_{\mathrm{ex}})}-1)$. In other words, we apply the local thermal equilibrium (LTE) approximation. The LTE approximation is justified when the gradients in the temperature are expected to be small compared to a material-dependent current-current correlation length scale, which is of the order of atomic scale or the phonon mean free path Polimeridis2015 . In the second case, the source-field photon number of the QW is modeled using $\langle\hat{\eta}_{\text{\tiny QW}}\rangle=1/(e^{(\hbar\omega-eU)/(k_{\mathrm{B}}T)}-1)$ for photon energies above the band gap $\hbar\omega\geq E_{\mathnormal{g}}$ and the background value $\langle\hat{\eta}_{\text{\tiny BG}}\rangle=1/(e^{\hbar\omega/(k_{\mathrm{B}}T)}-1)$ for photon energies below the band gap $\hbar\omega\leq E_{\mathnormal{g}}$ corresponding to the interactions with the free carriers. For example, in the case of the photon energy $\hbar\omega=E_{\mathnormal{g}}+k_{\mathrm{B}}T=2.786$ eV, the source-field photon number of the electrically or optically excited QW is $\langle\hat{\eta}_{\text{\tiny QW}}\rangle=7.51\times 10^{-4}$, which is very large in comparison with the photon number of the thermal 300 K background $\langle\hat{\eta}_{\text{\tiny BG}}\rangle=1.57\times 10^{-47}$. As the photon-number expectation values are relatively small and depend strongly on the frequency, it is convenient to illustrate the results by using the effective field temperature that is defined in terms of the photon-number expectation value as $T_{\mathrm{eff}}(z,K,\omega)=\hbar\omega/\{(k_{\mathrm{B}}\ln[1+1/\langle\hat{n}(z,K,\omega)\rangle]\}$ Partanen2014c ; Partanen2014b . The corresponding effective-source-field temperature of the electrically or optically excited QW ranges from 5175 K (compare with $\sim$6000 K of solar radiation on earth) to 625 K as the photon energy ranges from 2.76 to 5 eV.

Figure 2(a) shows the base-10 logarithm of the total EM LDOS for photon energy $\hbar\omega=E_{\mathnormal{g}}+k_{\mathrm{B}}T=2.786$ eV as a function of position and the in-plane component of the wave vector. The sapphire substrate lies on the left and air on the right. The light cones for sapphire, GaN, and air are defined by the in-plane wave vector component values $K<nk_{0}$, where $n$ is the real part of the refractive index of the respective material. The light cones of the different material layers are clearly visible in the figure. Due to the evanescent fields, the LDOSs are slightly elevated also beyond the material interfaces. One can also see the very large LDOS associated with the GaN/Ag SP resonance near the position $z=0$ and $K/k_{0}=5.0$. The less visible air/Ag SP resonance is near the position $z=0$ and $K/k_{0}=1.0$. The GaN guided modes and the associated interference patterns can be seen between the GaN light cone and the sapphire light cone with $1.78<K/k_{0}<2.51$.

Figure 2(b) shows the effective-field temperature of the total EM field corresponding to the LDOS in Fig. 2(a) in the case of a thermally excited QW. For the narrow In0.15Ga0.85N layer located slightly left from the position $z=0$ $\mu$m, the source-field temperature is 350 K; for other material layers it is 300 K. It can be seen that the light cones of each material are visible also in the effective-field temperature. The evanescent fields near the material interfaces are even more pronounced when compared to the LDOS in Fig. 2(a). At high values of $K$, the effective-field temperatures approach the source-field temperature in each layer, whereas for K within the light cones of air and sapphire, the effective temperature is reduced due to the strong coupling to the semi-infinite air and sapphire layers. Figure 2(c) shows the corresponding effective-field temperature of the total EM field in the case of an electrically or optically excited QW. The figure clearly resembles the case of thermal excitation in Fig. 2(b), but the values of the effective-field temperature are significantly higher, as expected.

Figure 3(a) presents the base-10 logarithm of the total EM LDOS in the QW as a function of the photon energy and the in-plane component of the wave vector. The figure clearly shows the GaN/Ag SP resonance, as well as the GaN guided modes corresponding to the Fabry-Pérot resonances of the cavity. At photon energy slightly above 3 eV, the GaN becomes absorptive, and therefore, there are no resonances visible above this energy. If, instead of the total EM LDOS, we were to plot the LDOS parts corresponding to the TE and TM polarizations, then the SP modes would be visible only in the TM case, as previously discussed, e.g., in Ref. Sadi2013b . Otherwise, the LDOSs of the TE and TM polarizations are qualitatively very similar. Figure 3(b) shows the corresponding base-10 logarithm of the total EM LDOS in air at 1 nm above the structure. In addition to the resonances visible in Fig. 3(a), in Fig. 3(b), one can also see the Ag/air SP mode just above the light cone of air. Figure 3(c) presents the base-10 logarithm of the EM LDOS in air at 1 $\mu$m above the structure. One can clearly see that there is only a small contribution of the evanescent fields remaining, especially at high frequencies, and the only significant contribution to the EM LDOS arises from the propagating modes in the light cone of air.

Figure 3(d) shows the effective-field temperature of the total EM field in the middle of the QW as a function of energy and $K/k_{0}$ for the case of a thermally excited QW. The effective temperature is essentially above the background temperature of 300 K when the imaginary part of the refractive index of the InGaN QW significantly deviates from zero either due to band-to-band or other absorption and emission mechanisms. At low frequencies, the thin InGaN becomes nearly transparent, and therefore, the effective temperature reaches the background temperature. The emissivity peak near photon energy $\hbar\omega=0.5$ eV follows from the peak in the infrared absorption coefficient of the QW Tansley1986 . The corresponding effective-field temperature in the case of an electrically or optically excited QW is shown in Fig. 3(g). The emission begins at the photon energy corresponding to the band gap, where the resulting effective-field temperature also obtains its highest values, as expected. At high energies well above the band gap, the effective-field temperature again reaches the source-field temperature of 300 K. The effective temperature of the field generally increases as the optical confinement of the mode increases: In the case of an electrically or optically excited QW and photon energy $\hbar\omega=E_{\mathnormal{g}}+k_{\mathrm{B}}T=2.786$ eV, the modes extending into the light cone of air have $T_{\mathrm{eff}}\approx 2200$ K, whereas the modes bound in the light cone of GaN reach $T_{\mathrm{eff}}\approx 2700$ K, while the evanescent InGaN modes reach values as high as $T_{\mathrm{eff}}\approx 3500$ K. For all these cases, however, $T_{\mathrm{eff}}$ remains well below the source-field temperature of the QW due to the losses caused by the surrounding lossy materials.

Figure 3(e) presents the effective-field temperature of the total EM field in air at 1 nm above the structure, corresponding to the LDOS in Fig. 3(b) in the case of a thermally excited QW. The values of the effective-field temperature are somewhat lower than the values of the effective-field temperature in the QW in Fig. 3(d). This is mainly due to the attenuation related to the increased distance to the excited QW. The effective-field temperature in Fig. 3(e), however, resembles the effective-field temperature in the QW. Also, the effective-field temperatures in the case of an electrically or optically excited QW at the two positions presented in Figs. 3(g) and Fig. 3(h) are quite similar. In the case of a thermally excited QW at low frequencies, the infrared emission of the QW is not visible in air as the silver layer between air and the QW becomes very lossy at low frequencies.

Figures 3(f) and 3(i) show the effective-field temperatures of the total EM field in air at 1 $\mu$m above the structure corresponding to the LDOS in Fig. 3(c) in the cases of thermally and electrically or optically excited QWs. The contribution of the evanescent fields is reduced as in the case of the EM LDOS. Due to the longer distance to the structure and reflections at the interfaces, the values of the effective-field temperatures are also consequently lower compared to the values of the effective-field temperatures in the QW in Figs. 3(d) and 3(g).

IV Conclusions

We have developed a three-dimensional QFED method to describe the photon-number quantization and thermal balance in general lossy and lossless geometries. By appropriately defining the photon ladder operators and the densities of states, we were able to present the ladder operators and the photon-number expectation values using formulas that are equivalent to the forms previously obtained by using a one-dimensional formalism. The resulting generalized QFED method allows studying, e.g., plasmonic structures and defining an effective-field temperature that realistically describes the excitation of the optical field.

To demonstrate the applicability and physical implications of the presented QFED method, we have used the model to study the energy and position dependencies of the EM LDOSs and effective-field temperatures in a light-emitting InGaN QW structure, which has recently been of experimental and theoretical interest. The results show that the developed method is well suited for analyzing the emission of electrically, optically, or thermally excited QWs. The effective temperatures were studied both as a function of position and as a function of photon energy. Electrical and optical excitations of the QW produce high effective-field temperatures, whose energy spectrum is quite narrow, whereas the effective-field temperature of a thermally excited QW has a significantly broader emission spectrum, as expected.

In addition to providing further insight into the classical fluctuational electrodynamics theory, the QFED method enables interesting further studies as it bridges the classical propagating wave picture of the EM field and the fluctuational electrodynamics, which is widely used to model near-field effects. Therefore, we expect that using the QFED, method one could, for instance, find a radiative transfer equation that allows describing interference effects, thus widening the applicability of the conventional radiative transfer equation beyond its main limitation in describing interference effects. This would make it possible to use the radiative transfer equation to describe also near-field effects in resonant structures.

Acknowledgements.

This work has been funded in part by the Academy of Finland and the Aalto Energy Efficiency Research Programme.

Appendix A Green’s functions

Here we briefly review the known relations between the electric, magnetic, and exchange Green’s functions. We first define the electric Green’s function \overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{\mathrm{ee}}(\mathbf{r},\omega,\mathbf{r}^{\prime}) that satisfies Paulus2000

[TABLE]

where \overset{\text{\tiny\leftrightarrow}}{\mathbf{I}} is the unit dyadic and $k_{0}=\omega/c$ is the wavenumber in vacuum with the vacuum velocity of light $c$. The subscript $\mathbf{r}$ in $\nabla_{\mathbf{r}}$ highlights that the differentiation is here performed with respect to $\mathbf{r}$ instead of $\mathbf{r}^{\prime}$. The solution of Eq. (5) is then written in terms of the electric Green’s function \overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{\mathrm{ee}}(\mathbf{r},\omega,\mathbf{r}^{\prime}) by integrating the product of the Green’s function and the source terms over all the source points $\mathbf{r}^{\prime}$ as

[TABLE]

where the subscript $\mathbf{r}^{\prime}$ in $\nabla_{\mathbf{r}^{\prime}}$ inside the integral indicates that the differentiation is performed with respect to the source point $\mathbf{r}^{\prime}$. In the case of the second term, we have applied the Stokes’ theorem resulting in the integration by parts formula $\int_{V}\mathbf{G}^{\alpha}\cdot(\nabla_{\mathbf{r}^{\prime}}\times\mathbf{J})d^{3}r^{\prime}=\int_{V}(\nabla_{\mathbf{r}^{\prime}}\times\mathbf{G}^{\alpha})\cdot\mathbf{J}d^{3}r^{\prime}-\int_{\partial V}(\mathbf{G}^{\alpha}\times\mathbf{J})\cdot d\mathbf{S}^{\prime}$ separately for each row vector $\mathbf{G}^{\alpha}$ of the matrix representation of \overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{\mathrm{ee}}(\mathbf{r},\omega,\mathbf{r}^{\prime}) with the boundary condition that the Green’s functions go to zero when the separation between the source point $\mathbf{r}^{\prime}$ and the field point $\mathbf{r}$ tends to infinity. Using the short hand notation \overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{\mathrm{ee}}(\mathbf{r},\omega,\mathbf{r}^{\prime})\times\nabla_{\mathbf{r}^{\prime}}=-[\nabla_{\mathbf{r}^{\prime}}\times\mathbf{G}^{1},\nabla_{\mathbf{r}^{\prime}}\times\mathbf{G}^{2},\nabla_{\mathbf{r}^{\prime}}\times\mathbf{G}^{3}]^{T}, where $T$ denotes transpose, we then define the exchange Green’s function \overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{\mathrm{em}}(\mathbf{r},\omega,\mathbf{r}^{\prime}) as

[TABLE]

Solving for the magnetic field by using Faraday’s law in Eq. (3) and substituting the electric field operator in terms of the Green’s functions in Eq. (25) give

[TABLE]

where we have first substituted the expression for $\hat{\mathbf{E}}^{+}(\mathbf{r},\omega)$ from Eq. (25) and incorporated the separate $\hat{\mathbf{J}}_{\mathrm{m}}^{+}(\mathbf{r},\omega)$ term into the integral using a suitable $\delta$-function presentation and then defined the exchange Green’s function \overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{\mathrm{me}}(\mathbf{r},\omega,\mathbf{r}^{\prime}) and the magnetic Green’s function \overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{\mathrm{mm}}(\mathbf{r},\omega,\mathbf{r}^{\prime}) as

[TABLE]

[TABLE]

By using Eqs. (26) and (29), one also obtains an expression of the magnetic Green’s function \overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{\mathrm{mm}}(\mathbf{r},\omega,\mathbf{r}^{\prime}) directly in terms of the electric Green’s function \overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{\mathrm{ee}}(\mathbf{r},\omega,\mathbf{r}^{\prime}) as

[TABLE]

Appendix B Green’s functions for stratified media

To gain more insight and analytical formulas directly applicable to common planar geometries, and to partly lift the divergences associated with absorbing media, we apply the formalism developed above for stratified media. In the case of stratified media, it is convenient to use the plane wave representation for the dyadic Green’s functions: A point in space is denoted in the Cartesian basis $(\hat{\mathbf{x}},\hat{\mathbf{y}},\hat{\mathbf{z}})$ by $\mathbf{r}=x\hat{\mathbf{x}}+y\hat{\mathbf{y}}+z\hat{\mathbf{z}}=\mathbf{R}+z\hat{\mathbf{z}}$, where $\mathbf{R}=x\hat{\mathbf{x}}+y\hat{\mathbf{y}}$ is the in-plane coordinate and the surface normals are along the $z$ coordinate. Similarly, a wave vector of a plane wave is denoted by $\mathbf{k}=\mathbf{K}+k_{z}\,\mathrm{sgn}(z-z^{\prime})\hat{\mathbf{z}}$ where the component $\mathbf{K}$ is in the $x$-$y$ plane and $k_{z}$ is given by $k_{z}=\sqrt{k_{0}^{2}n^{2}-K^{2}}$, with $\mathrm{Im}(k_{z})\geq 0$. For convenience, we also define the unit vector $\hat{\mathbf{K}}=\mathbf{K}/K$.

The above notation is convenient since at the $x$-$y$ plane, the dyadic Green’s functions of stratified media depend only on the relative in-plane coordinate $\mathbf{R}-\mathbf{R}^{\prime}$. Therefore, in the plane-wave representation, the dyadic Green’s functions \overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{jk}(\mathbf{r},\omega,\mathbf{r}^{\prime}), $j,k\in\{\mathrm{e,m}\}$, can be written as Tomas1995 ; Intravaia2015

[TABLE]

where the terms \overset{\text{\tiny\leftrightarrow}}{\mathbf{R}}^{\mathrm{T}}\overset{\text{\tiny\leftrightarrow}}{\mathbf{g}}_{jk}(z,K,\omega,z^{\prime})\overset{\text{\tiny\leftrightarrow}}{\mathbf{R}} are the Fourier transforms of \overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{jk}(\mathbf{r},\omega,\mathbf{r}^{\prime}) that have been obtained by rotating the dyadic plane-wave Greens functions \overset{\text{\tiny\leftrightarrow}}{\mathbf{g}}_{jk}(z,K,\omega,z^{\prime}) calculated using the standard techniques to evaluate the fields in layered structures as presented below. More specifically, \overset{\text{\tiny\leftrightarrow}}{\mathbf{g}}_{jk}(z,K,\omega,z^{\prime}) have been evaluated in a coordinate system where the in-plane noise components are taken to be perpendicular and parallel to $\mathbf{K}$, and the rotation with the rotation matrix \overset{\text{\tiny\leftrightarrow}}{\mathbf{R}} is used to return this convention of direction back to the coordinate system where the direction of the dipoles does not depend on $\mathbf{K}$.

Due to the symmetry properties of the Green’s functions \overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{jk}(\mathbf{r},\omega,\mathbf{r}^{\prime}), also \overset{\text{\tiny\leftrightarrow}}{\mathbf{g}}_{jk}(z,K,\omega,z^{\prime}) obey symmetry relations. When the field and source positions $z$ and $z^{\prime}$ are interchanged, the values of the Green’s functions are changed according to the reciprocity relations as follows. The spectral dyadic Green’s functions \overset{\text{\tiny\leftrightarrow}}{\mathbf{g}}_{jk}(z,K,\omega,z^{\prime}) obey the reciprocity relation \overset{\text{\tiny\leftrightarrow}}{\mathbf{g}}_{kj}(z^{\prime},K,\omega,z)=\mathrm{diag}(-1,1,-1)\overset{\text{\tiny\leftrightarrow}}{\mathbf{g}}_{jk}(z,K,\omega,z^{\prime})^{T}\mathrm{diag}(-1,1,-1), where $\mathrm{diag}(-1,1,-1)$ is a diagonal matrix with diagonal elements $-1$, $1$, and $-1$. In addition, they obey the complex conjugation relation \overset{\text{\tiny\leftrightarrow}}{\mathbf{g}}_{jk}^{*}(z,K,\omega,z^{\prime})=\overset{\text{\tiny\leftrightarrow}}{\mathbf{g}}_{jk}(z,K,-\omega,z^{\prime}).

B.1 Multi-interface reflection and transmission coefficients

In order to write the spectral dyadic Green’s functions for a multi-interface geometry in a compact form, we first define the multi-interface reflection and transmission coefficients that take into account all the reflections in the geometry. First, the single-interface reflection and transmission coefficients following from the boundary conditions of the tangential and normal polarizations ($\sigma\in\{\parallel,\perp\}$) requiring the tangential components of the electric and magnetic fields to be continuous at interfaces are given for the electric and magnetic fields by

[TABLE]

where $\varepsilon_{l}$, $\mu_{l}$, $n_{l}$, $k_{z,l}$, $l=1,2$, are the relative permittivities, permeabilities, refractive indices, and the $z$ components of the wave vectors in the two materials. The single-interface coefficients in Eq. (32) equal the conventional reflection and transmission coefficients used, e.g., in Ref. Paulus2000 . In the following, with primed reflection and transmission coefficients we denote the reflection and transmission coefficients for the incidence from medium 2 to medium 1, and they are obtained by switching indices 1 and 2 in Eq. (32).

The multi-interface geometry is defined by interface positions $z_{l}$, $l=1,2,...,N$, separating material layers with relative permittivities and permeabilities $\varepsilon_{l}$ and $\mu_{l}$, $l=1,2,...,N+1$. The layer thicknesses are denoted by $d_{l}=z_{l}-z_{l-1}$, where $l=2,...,N$. The multi-interface reflection and transmission coefficients $\mathcal{R}_{l,j}$ and $\mathcal{T}_{l,j}$, which account for the multiple reflections in different medium layers, are recursively given in terms of the single-interface reflection and transmission coefficients as

[TABLE]

where $l=1,2,...,N$, $j\in\{\mathrm{e},\mathrm{m}\}$, $\sigma\in\{\parallel,\perp\}$ $\nu_{l,j,\sigma}=1/(1-\mathcal{R}_{l-1,j,\sigma}^{\prime}\mathcal{R}_{l,j,\sigma}e^{2ik_{z,l}d_{l}})$, and $\mathcal{R}_{0,j,\sigma}^{\prime}=\mathcal{R}_{N+1,j,\sigma}=0$. As in the case of single-interface coefficients in Eq. (32) the primed coefficients denote the coefficients for right incidence. The layers are indexed such that $\mathcal{R}_{l,j,\sigma}^{\prime}$ corresponds to the same interface as $\mathcal{R}_{l,j,\sigma}$. The propagation coefficient for a certain material layer of thickness $d_{l}$ is given as $\mathcal{P}_{l}=e^{ik_{z,l}d_{l}}$ when the transmission coefficient from layer $l^{\prime}$ to layer $l>l^{\prime}+1$ is recursively given by $\mathcal{T}_{l^{\prime},l,j,\sigma}=\mathcal{T}_{l^{\prime},l-1,j,\sigma}\mathcal{T}_{l-1,j,\sigma}e^{ik_{z,l-1}d_{l-1}}$, with $\mathcal{T}_{l^{\prime},l^{\prime}+1,j,\sigma}=\mathcal{T}_{l^{\prime},j,\sigma}$, and that from layer $l^{\prime}$ to layer $l<l^{\prime}-1$ by $\mathcal{T}_{l^{\prime},l,j,\sigma}^{\prime}=\mathcal{T}_{l^{\prime},l+1,j,\sigma}^{\prime}\mathcal{T}_{l,j,\sigma}^{\prime}e^{ik_{z,l+1}d_{l+1}}$, with $\mathcal{T}_{l^{\prime},l^{\prime}-1,j,\sigma}^{\prime}=\mathcal{T}_{l^{\prime}-1,j,\sigma}^{\prime}$.

B.2 Spectral dyadic Green’s functions

Here we give a compact componentwise representation of the spectral dyadic Green’s functions for general stratified media. The presentation adapts the dyadic Green’s functions given, e.g., in Ref. 45 or, in the case of purely dielectric structures, in Appendix A of Ref. 47. However, the chosen presentation has a few differences: (1) We use the orthonormal basis $(\hat{\mathbf{K}}\times\hat{\mathbf{z}},\hat{\mathbf{K}},\hat{\mathbf{z}})$, where the in-plane noise components are taken to be perpendicular and parallel to $\mathbf{K}$. Then, the rotation matrix \overset{\text{\tiny\leftrightarrow}}{\mathbf{R}} in Eq. (31) is used to return this convention of direction back to the coordinate system where the direction of the dipoles does not depend on $\mathbf{K}$. (2) Instead of using the orthonormal system of complex-valued unit dyads of Refs. 45 and 47, we write the dyadic Green’s functions as matrices. (3) We use the scaled forms of the dyadic Green’s functions; for example, \overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{\mathrm{ee}}(\mathbf{r},\omega,\mathbf{r}^{\prime}) is obtained as a solution to the differential equation in Eq. (24) instead of the corresponding equation in Ref. 45, whose right-hand side contains an additional factor $1/\mu(\mathbf{r},\omega)$. Thus, our notation corresponds to the notation used in the case of normal incidence in Ref. 34.

The spectral dyadic Green’s functions \overset{\text{\tiny\leftrightarrow}}{\mathbf{g}}_{\mathrm{ee}}(z,K,\omega,z^{\prime}) and \overset{\text{\tiny\leftrightarrow}}{\mathbf{g}}_{\mathrm{mm}}(z,K,\omega,z^{\prime}), in our notation, are given in terms of the scaled dyadic Green’s functions \overset{\text{\tiny\leftrightarrow}}{\boldsymbol{\xi}}_{\mathrm{e}}(z,K,\omega,z^{\prime}) and \overset{\text{\tiny\leftrightarrow}}{\boldsymbol{\xi}}_{\mathrm{m}}(z,K,\omega,z^{\prime}) as

[TABLE]

The scaled dyadic Green’s functions \overset{\text{\tiny\leftrightarrow}}{\boldsymbol{\xi}}_{\mathrm{e}}(z,\omega,z^{\prime}) and \overset{\text{\tiny\leftrightarrow}}{\boldsymbol{\xi}}_{\mathrm{m}}(z,\omega,z^{\prime}) in Eqs. (35) and (36) are given in the orthonormal basis $(\hat{\mathbf{K}}\times\hat{\mathbf{z}},\hat{\mathbf{K}},\hat{\mathbf{z}})$ by

[TABLE]

where the primed and unprimed quantities correspond to the quantities at positions $z^{\prime}$ and $z$, respectively, and $\xi_{j,\sigma}^{\pm}(z,\omega,z^{\prime})$ are the scaled scalar Green’s functions, which are presented in the case of normal incidence in Ref. 34. For non-normal incidence, the generalization of the scaled scalar Green’s functions of Ref. 34 is obtained by substituting the wave number $k$ with its $z$ component $k_{z}$ and the reflection and transmission coefficients of normal incidence with the corresponding quantities for non-normal incidence, given in Eqs. (33) and (34). Assuming that the source point $z^{\prime}$ is located in layer $l^{\prime}$ ($z_{l^{\prime}-1}<z^{\prime}<z_{l^{\prime}}$) and the field point $z$ is located in layer $l$ ($z_{l-1}<z<z_{l}$), in the three cases $l=l^{\prime}$, $l>l^{\prime}$, and $l<l^{\prime}$, the scaled scalar Green’s functions are compactly given by

[TABLE]

[TABLE]

[TABLE]

The spectral dyadic Green’s functions \overset{\text{\tiny\leftrightarrow}}{\mathbf{g}}_{\mathrm{me}}(z,K,\omega,z^{\prime}) and \overset{\text{\tiny\leftrightarrow}}{\mathbf{g}}_{\mathrm{em}}(z,K,\omega,z^{\prime}) following from Eqs. (26), (28), and (31) are, respectively, presented in terms of the scaled dyadic exchange Green’s functions \overset{\text{\tiny\leftrightarrow}}{\boldsymbol{\xi}}_{\mathrm{ex,e}}(z,K,\omega,z^{\prime}) and \overset{\text{\tiny\leftrightarrow}}{\boldsymbol{\xi}}_{\mathrm{ex,m}}(z,K,\omega,z^{\prime}) as

[TABLE]

The scaled dyadic exchange Green’s functions \overset{\text{\tiny\leftrightarrow}}{\boldsymbol{\xi}}_{\mathrm{ex,e}}(z,K,\omega,z^{\prime}) and \overset{\text{\tiny\leftrightarrow}}{\boldsymbol{\xi}}_{\mathrm{ex,m}}(z,K,\omega,z^{\prime}) are given in terms of the scaled scalar Green’s functions in Eqs. (41)–(43) by

[TABLE]

Appendix C Derivation of the densities of states

C.1 Nonlocal densities of states

The time-domain field operators are obtained from the frequency-domain operators by Fourier transforms. For example, the time-domain electric field operator is given by

[TABLE]

where $\hat{\mathbf{E}}^{-}(\mathbf{r},\omega)$ is the negative-frequency part obtained by a Hermitian conjugate of the positive-frequency part $\hat{\mathbf{E}}^{+}(\mathbf{r},\omega)$ in Eq. (6).

The frequency-space correlation functions are given by

[TABLE]

[TABLE]

In the time domain, we have

[TABLE]

which then becomes

[TABLE]

Using $|j_{\mathrm{0,e}}(\mathbf{r}^{\prime},\omega)|^{2}=4\pi\hbar\omega^{2}\varepsilon_{0}\varepsilon_{\mathrm{i}}(\mathbf{r}^{\prime},\omega)$ and $|j_{\mathrm{0,m}}(\mathbf{r}^{\prime},\omega)|^{2}=4\pi\hbar\omega^{2}\mu_{0}\mu_{\mathrm{i}}(\mathbf{r}^{\prime},\omega)$ gives

[TABLE]

This allows defining the NLDOS for the electric field as

[TABLE]

Corresponding equations can be written for the magnetic field. The NLDOS of the magnetic field is then given by

[TABLE]

C.2 Interference density of states

For an optical mode the quantum optical Poynting vector is defined as a normal ordered operator in terms of the positive- and negative-frequency parts of the electric and magnetic field operators as $\hat{\mathbf{S}}(\mathbf{r},t)=:\!\hat{\mathbf{E}}(\mathbf{r},t)\times\hat{\mathbf{H}}(\mathbf{r},t)\!:=\hat{\mathbf{E}}^{-}(\mathbf{r},t)\times\hat{\mathbf{H}}^{+}(\mathbf{r},t)-\hat{\mathbf{H}}^{-}(\mathbf{r},t)\times\hat{\mathbf{E}}^{+}(\mathbf{r},t)$ Loudon2000 . Substituting the time-space forms of the electric and magnetic field operators in Eqs. (6) and (7) gives

[TABLE]

The frequency-space correlation functions are given by

[TABLE]

[TABLE]

Here \mathrm{Tr}[\overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{jj}^{\dagger}(\mathbf{r},\omega,\mathbf{r}^{\prime})\times\overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{kj}(\mathbf{r},\omega,\mathbf{r}^{\prime})]=\sum_{\sigma}[\hat{\mathbf{e}}_{\sigma}\cdot\overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{jj}^{\dagger}(\mathbf{r},\omega,\mathbf{r}^{\prime})]\times[\overset{\text{\tiny\leftrightarrow}}{\mathbf{G}}_{kj}(\mathbf{r},\omega,\mathbf{r}^{\prime})\cdot\hat{\mathbf{e}}_{\sigma}], which is a vector, in contrast to the conventional trace of a matrix. The Poynting vector then becomes

[TABLE]

Using $|j_{\mathrm{0,e}}(\mathbf{r}^{\prime},\omega)|^{2}=4\pi\hbar\omega^{2}\varepsilon_{0}\varepsilon_{\mathrm{i}}(\mathbf{r}^{\prime},\omega)$ and $|j_{\mathrm{0,m}}(\mathbf{r}^{\prime},\omega)|^{2}=4\pi\hbar\omega^{2}\mu_{0}\mu_{\mathrm{i}}(\mathbf{r}^{\prime},\omega)$ gives

[TABLE]

This allows defining the IFDOS as

[TABLE]

where $n_{\mathrm{r}}(\mathbf{r},\omega)$ is the real part of the refractive index.

Appendix D Densities of states for stratified media

Here we present the densities of states for stratified media by using the components $g_{jk}^{\alpha\beta}$, $\alpha,\beta\in\{1,2,3\}$, of the matrix representations of the spectral dyadic Green’s functions \overset{\text{\tiny\leftrightarrow}}{\mathbf{g}}_{jk}.

D.1 Nonlocal densities of states

Using Eq. (55) with $\varepsilon(\mathbf{r}^{\prime},\omega)=\varepsilon(z^{\prime},\omega)$, $\mu(\mathbf{r}^{\prime},\omega)=\mu(z^{\prime},\omega)$, $\langle\hat{\eta}(\mathbf{r}^{\prime},\omega)\rangle=\langle\hat{\eta}(z^{\prime},\omega)\rangle$, and

[TABLE]

where $j,k\in\{\mathrm{e,m}\}$, gives

[TABLE]

Then, the NLDOS for the electric field can be written as

[TABLE]

where $g_{\mathrm{ee}}^{\alpha\beta}$ and $g_{\mathrm{em}}^{\alpha\beta}$, with $\alpha,\beta\in\{1,2,3\}$, are components of the matrix representations of the spectral dyadic Green’s functions \overset{\text{\tiny\leftrightarrow}}{\mathbf{g}}_{\mathrm{ee}} and \overset{\text{\tiny\leftrightarrow}}{\mathbf{g}}_{\mathrm{em}}. The NLDOS of the magnetic field is given by

[TABLE]

D.2 Local densities of states

As integrals of the electric and magnetic NLDOSs in Eqs. (66) and (67), the electric and magnetic LDOSs are given by

[TABLE]

D.3 Interference density of states

Using Eq. (62) with $\varepsilon(\mathbf{r}^{\prime},\omega)=\varepsilon(z^{\prime},\omega)$, $\mu(\mathbf{r}^{\prime},\omega)=\mu(z^{\prime},\omega)$, $\langle\hat{\eta}(\mathbf{r}^{\prime},\omega)\rangle=\langle\hat{\eta}(z^{\prime},\omega)\rangle$, and

[TABLE]

[TABLE]

gives

[TABLE]

Note that the Poynting vector points purely in the $z$ direction, which is natural due to the symmetry with respect to the $z$ axis. Hence, the IFDOS can be written as $\boldsymbol{\rho}_{\mathrm{IF}}(z,K,\omega,z^{\prime})=\hat{\mathbf{z}}\rho_{\mathrm{IF}}(z,K,\omega,z^{\prime})$, where the scalar IFDOS $\rho_{\mathrm{IF}}(z,K,\omega,z^{\prime})$ is given by

[TABLE]