Interference-exact radiative transfer equation

The Purcell effect, i.e., the modification of the spontaneous emission rate by optical interference, profoundly affects the light-matter coupling in optical resonators. Fully describing the optical absorption, emission, and interference of light hence conventionally requires combining the full Maxwell’s equations with stochastic or quantum optical source terms accounting for the quantum nature of light. We show that both the nonlocal wave and local particle features associated with interference and emission of propagating fields in stratified geometries can be fully captured by local damping and scattering coefficients derived from the recently introduced quantized fluctuational electrodynamics (QFED) framework. In addition to describing the nonlocal optical interference processes as local directionally resolved effects, this allows reformulating the well known and widely used radiative transfer equation (RTE) as a physically transparent interference-exact model that extends the useful range of computationally efficient and quantum optically accurate interference-aware optical models from simple structures to full optical devices.


Results
Homogeneous medium. In general, the damping and scattering coefficients in Eq. (2) can be position dependent. This is naturally not the case in a homogeneous space, where the damping and scattering coefficients are constant and separately equal for fields propagating in different directions, i.e., α +,σ (z, K, ω) = α −,σ (z, K, ω) = α ±,σ and β +,σ (z, K, ω) = β −,σ (z, K, ω) = β ±,σ . As shown in the Supplemental Material, substituting the densities of states corresponding to the homogeneous space Green's functions into the damping and scattering coefficients in Eq. (2) leads to damping and scattering coefficients where k z,i is the imaginary part of the wave vector z component = − k k K z 2 2 , the wavenumber is k = nω/c, n is the refractive index, c is the speed of light in vacuum, and the parameter ψ σ is given for the TE and TM polarizations by Here k r is the real part of the wavenumber and ε and μ are the relative permittivity and permeability of the medium, which are related to the refractive index as εµ = n . In a lossless uniform medium, the damping and scattering coefficients are all zero for propagating modes as the imaginary part of the z-component of the wave vector is zero. In homogeneous lossy media, on the other hand, the damping and scattering coefficients are both positive and the damping coefficients α ±,σ are larger than the scattering coefficients β ±,σ . For normal incidence with K = 0, the coefficients for the TE and TM polarizations are equal as expected. In a purely dielectric medium with ε = n , the damping and scattering coefficients in Eq. (3) simplify for K = 0 to α = | | + | | ± k n n n n ( / / ) i 2 r 2 r 2 2 and β = | | − | | ± k n n n n ( / / ) i 2 r 2 r 2 2 where k i is the imaginary part of the wavenumber and n r is the real part of the refractive index. In the limit of small losses, one can approximate |n| ≈ n r and the damping and scattering coefficients become approximately equal to the classical expressions α ± ≈ 2k i and β ± ≈ 0. For larger losses, however, the scattering coefficients become nonzero. For example, in the case of a dielectric material with refractive index n = 2 + i, the coefficients are given by α ± = 2.05k 0 and β ± = 0.45k 0 , where k 0 = ω/c is the wavenumber in vacuum. The damping coefficient is still close to the classical value, but the scattering coefficient clearly deviates from the classical result of zero, indicating that a part of the photons is scattered backwards due to the induced electric and magnetic dipoles in the medium.
Single-interface geometry. To illustrate the general position dependence of the damping and scattering coefficients, we next study the damping and scattering coefficients for photon energy ħω = 1 eV (λ = 1.24 μm) in the vicinity of an interface between two lossy media with refractive indices ε = = . + . . The normal incidence damping and scattering coefficients for the single interface structure are given in Supplemental Material and, in Fig. 1, they are plotted as a function of the position. In Fig. 1(a), the damping coefficients α + and α − reach homogeneous field values 1.00077k 0 and 0.60046k 0 far from the interface, whereas, near the interface, they oscillate. The above homogeneous field values of the damping coefficients are seen to be very close to the classical values 2k 1,i = k 0 and 2k 2,i = 0.6k 0 . The oscillations of the damping coefficients near the interface originate from the interference and the modified position-dependent emission and absorption rates in analogy with the Purcell effect 6 . The oscillations in both the scattering and damping coefficients are substantially larger in the propagation directions away from the interface, suggesting that the photons propagating away from the interface experience stronger interference effects. For certain material combinations, the oscillations in the damping coefficients can become much larger and it is also possible for the damping coefficients to obtain negative values suggesting that the field may experience local amplification.
The scattering coefficients β + and β − in Fig. 1(b) have significantly smaller values than the damping coefficients α + and α − in Fig. 1(a). This is expected as the change of the field propagating in one direction generally depends more on the field itself than on the field propagating in the other direction. In addition, also the scattering coefficients β + and β − can obtain negative values near the interfaces due to interference. On the left and right, the oscillations in the scattering coefficients die out and saturate to the homogeneous space values 0.03923k 0 and 0.02354k 0 , which are nonzero, thus slightly deviating from the classical results. Two-interface resonator. Next we study the damping and scattering coefficients in a two-interface resonator formed by a dielectric slab with a refractive index ε = = + . n i 2 0 1 placed in vacuum. We also compare the results of our interference-exact RTE model and the classical field-based methods directly solving Maxwell's equations with appropriate boundary conditions. For a concise comparison, we use the negative divergence of the Poynting vector −∇ · S (here S is calculated using Eq. (6) in ref. 14), which describes the net absorption rate that is well-known to oscillate inside lossy resonant structures due to interference and the related Purcell effect 6 . These oscillations cannot be described correctly by using models such as the conventional RTE model which neglects essentially all interference effects, e.g., coherent backscattering 5 . Figure 2(a) shows the damping and scattering coefficients as a function of position for photon energy ħω = 0.46 eV (λ = 2.68 μm) and for normal incidence in the vicinity of the dielectric slab. The used photon energy corresponds to the second constructive interference of the field reflected from the slab, i.e., the intensity of the reflected field obtains its second maximum when it is plotted as a function of photon energy. One can clearly see that the damping and scattering coefficients are oscillating in the slab. Outside the slab, the coefficients are zero as there are no losses in vacuum.
When comparing the results of the derived interference-exact RTE model in Fig. 2(b) with the results of the classical solution of Maxwell's equations, we assume normal incidence, set the source-field temperature of the resonator to zero as T = 0 K, and use the initial condition that, on the right, there is only a right propagating field with a fixed average photon number, i.e., In the classical method, this corresponds to fixing the electric and magnetic fields on the right such that the resulting Poynting vector is the same as that in our interference-exact RTE model. The net absorption rates calculated by using the interference-exact RTE model and the classical method in Fig. 2(b) are equal within the precision of the numerical accuracy of the computations. This clearly demonstrates that by using the derived position-dependent damping and scattering coefficients, the conventional RTE model can be extended to account for all interference effects. More extensive comparisons between the interference-exact RTE model and full FED based calculations have additionally indicated that the results are also generally valid for other angles of incidence and resonator source fields T > 0 K.

Conclusions
In conclusion, we have used the newly developed QFED framework to derive the interference modified local field-matter coupling strengths of propagating fields and to extend the widely deployed RTE model so that it receives the ability to also fully capture interference effects in stratified geometries. The approach involves deriving the quantum optically exact damping and scattering coefficients from the position dependent expectation values of the propagating photon-number operators provided by the QFED framework. This approach allows providing an accurate and transparent local physical picture of interference as a mechanism that modulates the strength of the light-matter interactions. In addition to the physical transparency, the approach is expected to be very useful from the computational point of view when the Green's functions are known analytically or can be easily solved for, as the photon numbers no longer exhibit the strong oscillations throughout the simulation space, as is the case for the electric fields in conventional Maxwell's equation based models. This allows substantially relaxing the requirements set on the problem discretization away from interfaces. Furthermore, the model is also expected to enable new possibilities for modeling quantum effects (like coherence and collective effects) of the fields in macroscopic structures or devices. Overall, the presented interference exact RTE model therefore allows substantial widening of the use of RTE-based models to a wide variety of new geometries involving e.g. thin-films and resonators which have not been previously accessible to the simple RTE-based methods.

Number of photons and densities of states.
In contrast to other approaches used to describe fields in lossy resonant media, the fundamental requirement of the QFED is the preservation of the canonical commutation relation , NL , Here K is the wave vector component in the x−y plane, ρ σ (z, K, ω) is the local density of states (LDOS) 9, 14, 15 and we refer to the weighting coefficients ρ NL±,σ (z, K, ω, z′) as the nonlocal densities of states (NLDOSs) of the left and right propagating fields. These NLDOSs are given as sums and differences of the NLDOSs of the total electromagnetic field ρ NL,σ (z, K, ω, z′) and the interference densities of states (IFDOSs) ρ IF,σ (z, K, ω, z′) as ρ NL±,σ (z, K, ω, z′) = ρ NL,σ (z, K, ω, z′) ± ρ IF,σ (z, K, ω, z′) 14 , which are all related to the electromagnetic Green's functions of the system. These densities of states have been originally derived in ref. 14 and 9 and, in Supplemental Material, they are explicitly given in terms of the spectral dyadic Green's function components for stratified media. In Eq. (5), η ω 〈 〉 σ z K ( , , ) is the source-field photon-number expectation value given for thermal fields by the Bose-Einstein distribution Derivation of the RTE coefficients. Substituting the integral expressions for the photon-number expectation values of the QFED framework as given by Eq. (5) into the RTE model in Eq. (1) and omitting the arguments z, z′, K, and ω for brevity, we obtain where the source-field terms η 〈 〉 σ of Eq. (1) have been imported within the integrals by using suitable δ-function representations. To determine the RTE-coefficients α ±,σ and β ±,σ , we require that the integrands on the left and right side of Eq. (6) must be equal at all positions z and z′, which corresponds to requiring that the RTE model is valid for arbitrary material temperature distributions represented by the position-dependent source-field photon-number expectation value η σ ⟨ ⟩ . Although it is not immediately apparent, Eq. (6) only provides two linearly independent components, one for the case z = z′ and one for any z ≠ z′. This enables us to separate Eq. (6) into an equation group of two linearly independent equations. For the case z = z′ we obtain an equation