Exact Quantum Electrodynamics of Radiative Photonic Environments

Abstract

We present a comprehensive second quantization scheme for radiative photonic devices. We canonically quantize the continuum of photonic eigenmodes by transforming them into a discrete set of pseudomodes that provide a complete and exact description of quantum emitters interacting with electromagnetic environments. This method avoids all reservoir approximations and offers new insights into quantum correlations, accurately capturing all non-Markovian dynamics. This method overcomes challenges in quantizing non-Hermitian systems and is applicable to diverse nanophotonic geometries.

The problem

The geometry of the environment defines a photon's interaction with matter, bringing complexity to its radiative behavior, from the Casimir force between an atom and a surface, or the enhanced fluorescence of a molecule by a nanoparticle, to the tapestry of colors scattered by stained glass. Lately, quantum emitters — such as atoms, fluorescent molecules, and quantum dots — have been coupled to ever more geometrically complex photonic devices, such as optical microcavities, nanobeams, plasmonic nanostructures and hybrid nanophotonic devices. The interplay with multiple photonic modes of varying radiative behavior leads to non-Markovian quantum dynamics. Such dynamics significantly impact the future development of quantum information processing, quantum transport, photochemistry and biological processes such as light harvesting.

The interaction of quantum emitters in photonic environments is usually described by an open quantum system coupled to a reservoir that characterises radiative and material loss. Such descriptions, typical of cavity quantum electrodynamics, simplify the continuous electromagnetic spectrum via a simple phenomenological model of the local density of states, and were well suited to early high-finesse micro resonators and Fabry-Perot cavities that radiate weakly. Lately however, focus has shifted toward ever more intricate photonic devices where extreme subwavelength field confinement leads to extreme light-matter interactions. Such systems typically radiate efficiently to the far field, exhibit broadband overlapping modes, and often have significant material losses. Classical models for the local density of states have been phenomenologically quantised to form an open quantum system by assuming Lorentzian photonic modes. Alternatively, classical quasi-normal modes have been transformed and quantized "canonically". However, such methods lack the generality of established quantum optical theory, and their disconnect between the near and far field precludes straightforward descriptions of photon statistics, input and output field, squeezed light, dressed states, and so on. Furthermore, reservoir correlations — from which complex non-Markovian dynamics arises — are unaccounted for. Hence, there is a need for an a priori method that provides a global picture for open quantum nanophotonic systems.

In this Letter, we develop an a priori complete quantum electrodynamic description of light-matter interactions for radiative photonic devices by transforming the quantized fields into pseudomodes, originally used in an elegant treatment of isolated Lorentzian resonances. We show that our generalized pseudomode expansion gives a complete and exact description of both near and far field without the need for a reservoir. The pseudomodes arise naturally from the quantum dynamical equations of motion, yet maintain a direct relationship to the resonant modes of the photonic system. We derive the pseudomode's quantum equations of motion to obtain all the correlations and non-Markovian dynamics of the field, and accurately describe the propagation of light. We give an example application of a quantum emitter coupled to a microresonator that supports many spectrally overlapping Mie resonances.

Quantizing the field

We start with the second quantization of the electromagnetic field, expanded in terms of its eigenmodes, labelled by a mode index and a wave number. For any nanophotonic system these are the solutions of the Helmholtz equation, Equation 1 in the source: the Laplacian plus the product of the relative permeability, the relative permittivity and the squared wave number, acting on the eigenmode, equals zero. From these one constructs the corresponding electric field — the imaginary unit times the speed of light times the wave number times the eigenmode — and magnetic field, the curl of the eigenmode. The nanostructure geometry is specified by the relative electric permittivity as a function of position, assuming non-magnetic materials.

To canonically quantize the field via this mode decomposition, the eigenmodes must satisfy an orthogonality condition — the volume integral of the permittivity times the product of two eigenmodes equals the product of Kronecker deltas in the mode index and the wave number — and the energy of each mode must be conserved in time. For free space this is achieved using periodic or zero Dirichlet boundary conditions within a box of some volume, which is subsequently taken to infinity. A similar approach is adopted here, but the shape of the bounding volume is chosen to match the geometry of the nanophotonic device, and the boundary condition that the Poynting flux through its surface vanishes is applied. The solutions cover all of space, and the discrete spectrum of wave numbers becomes continuous when we take the infinite-volume limit, which for convenience is performed after quantization.

The quantum field operators are then expanded over the eigenmodes. Equation 2 in the source gives the electric field operator as the imaginary unit times a sum over mode index and wave number, of the square root of Planck's reduced constant times the speed of light times the wave number, divided by twice the vacuum permittivity times the local relative permittivity, multiplied by the eigenmode and by the sum of the creation and annihilation operators.

We consider a two-level quantum emitter that interacts via the dipole interaction in the length gauge, where the system is described by the Hamiltonian of Equation 3 in the source: the Hamiltonian divided by Planck's reduced constant is the emitter transition frequency times the raising-lowering product, plus a sum over modes of the speed of light times the wave number times the photon number operator, plus a sum over modes of the coupling strength multiplying the annihilation operator times the raising operator plus the creation operator times the lowering operator. Each respective term corresponds to the emitter with its transition frequency, the field, and their interaction under the dipole and rotating wave approximations.

Equation 4 in the source gives the coupling strength for an emitter at a given position: the square root of Planck's reduced constant times the speed of light times the wave number, divided by twice the vacuum permittivity times the local relative permittivity, multiplied by the dipole moment dotted into the eigenmode. While here we only consider two-level emitters, our approach can be generalized to complex multilevel emitters starting from the minimal coupling Hamiltonian, keeping terms linear in the field.

The pseudomode transformation

We transform the system into the "pseudomode picture" to reveal the full range of quantum dynamics. The mode functions are transformed from the continuous set of Helmholtz solutions to a discrete set of pseudomodes, determined by the integral transformation of Equation 5 in the source: the integral from zero to infinity, over the wave number weighted by the mode density, of the product of the eigenmode at two positions times a decaying complex exponential in the delay, equals a sum over resonances of a residue factor times the product of the pseudomodes at those two positions, times a Heaviside step function of the delay minus the propagation time between the points, times a complex exponential in the pseudomode frequency. This is evaluated by the residue theorem over the poles in the lower half plane for positive delay. The complex pseudomode frequencies correspond to the successive resonances of each mode. The integrand's asymptotic behavior determines the pole location and which half plane the contour integral encloses, and leads to the Heaviside step function. Physically, this accounts for the finite time delay for light to propagate from one point to the other via the photonic structure.

The pseudomode transformation arises naturally from the Schrödinger equations of motion with the emitter initially excited; we consider general initial conditions using the Heisenberg equations elsewhere. Here the state vector is described by the amplitude of the emitter excited state and the amplitudes of the singly excited field modes. These obey a pair of coupled equations derived from the Schrödinger equation in the interaction picture. Formally integrating the second with zero initial field amplitude, inserting into the first, and taking the continuum limit in which the sum over wave numbers becomes an integral weighted by the mode density, gives Equation 6 in the source: the time derivative of the emitter amplitude equals minus a time integral, over a sum of mode indices and an integral over wave number, of the mode density times the squared coupling strength times a complex exponential in the detuning times the elapsed interval, multiplying the emitter amplitude at the earlier time.

The first Markov approximation — taking the coupling strength to be a constant square root of a decay rate over two pi — is often applied here to simplify the integral over wave number. Instead we evaluate the wave-number integral exactly, using the coupling strength of Equation 4 and the transformation of Equation 5, to obtain Equation 7 in the source: the same time derivative equals minus a time integral of a sum over discrete pseudomodes of the squared pseudomode coupling times a complex exponential in the detuning from the complex pseudomode frequency, multiplying the earlier emitter amplitude.

This is equivalent to the discrete set of coupled pseudomode equations, Equations 8a and 8b in the source: the imaginary unit times the time derivative of the emitter amplitude equals the sum over pseudomodes of the pseudomode coupling times the phase factor times the pseudomode amplitude; and the imaginary unit times the time derivative of each pseudomode amplitude equals its coupling times the conjugate phase factor times the emitter amplitude. Equation 9 in the source gives the pseudomode-emitter interaction strength: the square root of Planck's reduced constant times the speed of light times the complex pseudomode frequency, divided by twice the vacuum permittivity times the local relative permittivity, multiplied by the dipole moment dotted into the pseudomode.

By solving these equations for an initially excited emitter one obtains the full quantum dynamical evolution of the system. These equations are non-Hermitian, since the pseudomode frequencies and couplings are complex valued, which causes the pseudomode amplitudes to decay in time as energy is radiated to the far field.

The electromagnetic fields are expressed exactly by a time-dependent superposition of pseudomodes. The intensity is found from the field quadrature acting on the state. Equation 10 in the source gives the positive-frequency field acting on the state as the imaginary unit times a sum over pseudomodes of the same square-root amplitude factor, multiplied by the pseudomode, by its amplitude taken at the retarded time, and by a complex exponential in its frequency. Expanding that exponential into a retarded part and a delay part highlights the decay by a factor set by the imaginary part of the pseudomode frequency times the propagation delay. This decay ensures the field remains regular even though the pseudomodes themselves diverge at large radius, as expected for a non-Hermitian theory. Hence our pseudomode approach overcomes the mode divergence problem for radiating photonic systems to give a global description of the field.

Worked example: a silicon microsphere

To comprehensively demonstrate this new approach, we apply it to a quantum emitter coupled to a spherical silicon resonator of radius 1 micrometre and refractive index 3.446, surrounded by vacuum. We use dimensionless units in which Planck's reduced constant and the speed of light are one. This geometry supports Mie resonances with the near field enhanced due to the confinement of modes by the dielectric interface at the sphere surface. For microspheres much larger than the wavelength the Mie modes interfere constructively, forming high-finesse whispering gallery modes. For radii around a micrometre or less, individual Mie modes form a spectrum of distinct but overlapping resonances.

The eigenmodes take the form of vector spherical harmonics: one family built as the curl of the radial unit vector times a scalar function, the other as the curl of the first divided by the wave number, where the scalar function is an azimuthal phase times an associated Legendre polynomial in the polar angle times a radial function. The angular distributions for different integer orders are orthogonal.

Equation 11 in the source gives the radial functions as piecewise continuous solutions of the spherical Bessel equation, normalised by a mode normalisation factor: inside the sphere, a coefficient times the spherical Bessel function of the first kind with the refractive index folded into its argument; outside, a combination of the spherical Bessel and spherical Neumann functions with their own coefficients. These solutions are regular at the origin, zero on the surface of the bounding volume of radius much greater than the sphere, and produce a continuous tangential field at the dielectric interface. The interface conditions determine the three coefficients.

These produce the resonances of the spectrum, determined by the normalisation factor of Equation 12 in the source: in the infinite-volume limit it is the bounding radius over twice the squared wave number, multiplied by the product of the sum and the difference of the two exterior coefficients weighted by the imaginary unit — a quantity analogous to the mode volume of the photonic system. Resonances occur for real wave numbers adjacent to the complex roots of Equation 12, which are the poles of Equation 11. Similarly, the scattering cross-section is resonant at the same poles, which have long been identified as the natural modes. Furthermore, the radial functions are orthogonal due to the conditions at the origin, at the sphere surface and at the bounding radius, and therefore the vector spherical harmonics satisfy the orthonormality condition required of the eigenmodes, with the mode index now a triple of angular and radial labels. Finally, for a finite bounding radius, the allowed wave numbers form a countably infinite set, with asymptotic mode density equal to the bounding radius over pi; in the infinite limit the wave number becomes continuous.

We now quantize these eigenmodes and subsequently perform the pseudomode transformation. The second quantized electric field operator of Equation 2 is now obtained using these orthonormal modes, and the Hamiltonian is given by Equation 3 with coupling strengths given by Equation 4. For simplicity we choose a radially oriented emitter dipole, which only couples to the transverse magnetic modes. We now transform the system into a discrete set of pseudomodes that interact with the emitter, by evaluating the integral of Equation 5 over those modes. Separating out the wave-number dependent terms, this becomes an outer product of differential operators multiplying a radial integral, Equation 13 in the source: the integral from zero to infinity over the wave number, of the mode density divided by the wave number, times the product of the two radial functions, times a decaying complex exponential in the delay.

To evaluate that integral, initially taking the lower limit to minus infinity, we extend the integrand into the complex plane. We expand the radial function into its outgoing and incoming components, proportional to the spherical Hankel functions of the first and second kind. This splits the integrand into incoming, outgoing, and mixed components that converge to zero at large complex argument in the lower half plane under stated conditions on the delay, and in the upper half plane otherwise. After careful consideration of the poles and contours, Equation 14 in the source gives the integral as a step function times a lower-half-plane contour integral of the incoming component, plus a step function times an upper-half-plane contour integral of the outgoing component, evaluated using the residue theorem. To evaluate the integral with the physical lower limit of zero, we split the integrand further into symmetric and antisymmetric components, integrate the symmetric part, or its product with the sign function, as before, then divide the result by two. Consequently, for positive times, Equation 15 in the source gives the result as two pi times the imaginary unit times a sum of residues of the outgoing component over poles in the fourth quadrant of the complex plane, multiplied by a step function.

Equation 16 in the source gives the resulting pseudomodes: the differential operator acting on pi times the square root of the difference of the exterior coefficients divided by the derivative of their sum, evaluated at the pole, all multiplying the outgoing spherical Hankel function at that complex frequency. The term in brackets defines the pseudomodes' radial behavior. These pseudomodes differ from the original eigenmodes in their radial behavior, their complex nature, and their discrete spectrum over the roots of Equation 12 in the fourth quadrant. Two of them, labelled by angular and radial indices five, zero, four and eight, zero, three, correspond to the two narrow overlapping resonances near 1.72 micrometres in the spectrum.

The quantum dynamics of the system are found using Equation 8. Choosing the emitter to sit on the polar axis simplifies further calculations, since the spherical harmonics vanish there for non-zero azimuthal index. We set the emitter's shifted transition frequency to be resonant with the eight-three pseudomode, where Equation 17 in the source gives the Lamb shift due to off-resonant modes: a sum over those modes of the squared coupling multiplied by the detuning minus the imaginary unit times the linewidth, divided by the sum of the squared detuning and the squared linewidth. The dynamical equations for the 613 pseudomodes with angular index up to 30 and real part of the frequency below 20 inverse micrometres are solved rapidly via a similarity transformation that diagonalizes Equation 8.

The time evolution of the emitter excited state and the pseudomode populations shows the following. For a dipole moment of 10 debye, Rabi oscillations — inherently non-Markovian by nature — occur between the emitter and the strongly coupled eight-three mode at a frequency set by twice its coupling, and their decay is primarily due to the radiation of that mode. By contrast, the weakly coupled five-four mode's dynamics are Markovian, oscillating in phase with the emitter six orders of magnitude lower. All other modes have small populations and are initially rapidly excited to values set by the square of the coupling over the detuning. For 100 debye the Rabi oscillations are ten times faster, while the decay rate remains largely unchanged. For ten thousand debye the dynamics are more complex: the five-four mode, on the cusp of strong coupling, now oscillates out of phase with the emitter and radiates energy more efficiently than the eight-three mode, significantly accelerating the emitter's decay. The remaining non-resonant modes' behavior depends on the ratio of their linewidth to the emitter decay rate: strongly radiating pseudomodes follow the emitter amplitude, demonstrating Markovian behavior, while high-quality pseudomodes are non-Markovian, evolving with their own complex exponential at long times.

We can approximate our system, dominated by the eight-three and five-four modes, by a two-mode model given by Equation 18 in the source: the imaginary unit times the time derivative of the emitter amplitude equals the Lamb shift times that amplitude, plus a sum over just those two modes of their couplings times their phase factors times their amplitudes, with Equation 9 used for those two modes only, while the Lamb shift is calculated from Equation 17 over all the other modes. These approximate solutions show remarkable accuracy provided the emitter population is larger than the off-resonant mode populations.

Finally, we calculate the expected field intensity using Equation 10 for a 10 debye emitter, plotted against distance from the sphere and time, which clearly shows the light cone. The contribution from the high-quality eight-three mode, confined to the near field within about a micrometre of the surface, shows decaying Rabi oscillations. In all other modes, a short pulse originates at the sphere surface at early times due to transient excitation of the field. Radiation from the weakly coupled five-four mode then dominates over an intermediate range and interferes with more weakly excited modes. Later still, the five-four mode has radiatively decayed, revealing non-Markovian dynamics in the near field due to high-finesse off-resonant modes. The total intensity is dominated at short times by the off-resonant modes and at longer times by the coherent oscillations of the eight-three mode.

What the transformation shows

By transforming the continuum into a discrete set of pseudomodes, we solve the dynamics without the need for a reservoir or its accompanying approximations, and therefore we retain all the information about the continuum. Our theory demonstrates that emitter decay arises from the continuous nature of the quantized field: the system never decays to the ground state; its energy is merely dispersed over the continuous spectrum of eigenmodes, for which there is no inherent decay. The "decay" of the pseudomode amplitudes describes the tail of the photon wave packet as it propagates, when taken correctly at the retarded time. This arises naturally in our infinite yet closed quantum system, which avoids outgoing wave boundary conditions on the normal modes that necessitate nonstandard quantization of divergent modes for non-Hermitian systems.

Conclusion

In conclusion, we present a general theory that gives a complete and exact description for the quantum electrodynamics of a quantum emitter strongly interacting with a radiating photonic device. We quantize the continuous Helmholtz eigenmodes, which we then transform with one-to-one correspondence into a discrete set of non-Hermitian pseudomodes. Thus, we solve common problems met when quantizing non-Hermitian systems, such as mode divergence, defining mode volumes, and identifying canonical field variables. Furthermore, our approach precisely captures all quantum correlations of the field and emitter, avoiding common Markovian approximations, and, unlike other methods, accurately captures the light propagation to the far field. This new method can be further extended for arbitrary photonic geometries through the analytic continuation of the local density of states and can reveal the non-Markovian behavior exhibited in experimentally realizable systems at the nanoscale.

Acknowledgments

A. D. gratefully acknowledges support from the Royal Society University Research Fellowship URF/R1/180097 and URF/R/231024, and funding from EPSRC EP/X012689/1, EP/Y008774/1 and the CDT in Topological Design EP/S02297X/1. A. D. and B. Y. gratefully acknowledge support funding from Research Fellows Enhancement Award RGF/EA/181038.

---

Ben Yuen and Angela Demetriadou, School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham, United Kingdom. Published as Physical Review Letters 133, 203604 (2024), doi.org/10.1103/PhysRevLett.133.203604, published by the American Physical Society under CC BY 4.0.

(Running heads, page numbers and reference-number markers have been dropped; display equations are rendered in words keyed to their source equation numbers and inequalities written out; the five figures, the Supplemental Material and the references are at the source.)

(On this site: Milonni’s account of why an excited atom emits at all — the same balance of radiation reaction against vacuum fluctuations, read from the emitter’s side — is at /library/stm-c9f20044be and /library/stm-a55cf0aad0, with his two books at /library/stm-d0a2779af6 and /library/stm-9b9f932b26. Puthoff’s reading of the ground state as an equilibrium with that same field is at /library/stm-c7c1082f9b. For what happens when the geometry itself is engineered: the Casimir effect in microstructured geometries is at /library/stm-d41136300f, non-monotonic Casimir forces between silicon nanostructures at /library/stm-5d1ea02d6d, the dynamical Casimir effect in a superconducting waveguide at /library/stm-d6682d44c9, and White’s worldline computation of a custom Casimir cavity at /library/stm-1775ebeff1. Cavity- and vacuum-enhanced superconductivity, where reshaping the photonic environment changes a material’s state, is at /library/stm-812175a230, /library/stm-94b2666369 and /library/stm-b7a1a66f71; the direct sampling of vacuum field fluctuations is at /library/stm-6095b3deaf.)