Few-mode Field Quantization of Arbitrary Electromagnetic Spectral Densities

We develop a framework that provides a few-mode master equation description of the interaction between a single quantum emitter and an arbitrary electromagnetic environment. The field quantization requires only the fitting of the spectral density, obtained through classical electromagnetic simulations, to a model system involving a small number of lossy and interacting modes. We illustrate the power and validity of our approach by describing the population and electric field dynamics in the spontaneous decay of an emitter placed in a complex hybrid plasmonic-photonic structure.

We develop a framework that provides a few-mode master equation description of the interaction between a single quantum emitter and an arbitrary electromagnetic environment. The field quantization requires only the fitting of the spectral density, obtained through classical electromagnetic simulations, to a model system involving a small number of lossy and interacting modes. We illustrate the power and validity of our approach by describing the population and electric field dynamics in the spontaneous decay of an emitter placed in a complex hybrid plasmonic-photonic structure.
Over the last years, there has been large interest in developing strategies for quantizing electromagnetic (EM) modes in open, dispersive and absorbing photonic environments. This has been motivated largely by the desire to use nanophotonic devices for quantum optics and quantum technology applications. This is a theoretical challenge, as standard ways of obtaining quantized modes are not valid in these systems [1,2]. In principle, macroscopic quantum electrodynamics (QED) is the framework for such quantization in material structures described by EM constitutive relations [3][4][5][6][7][8][9]. However, in this formalism, electromagnetic fields are described by a continuum of harmonic oscillators, restricting its applicability to cases where they can be treated perturbatively or eliminated by Laplace transform or similar techniques. Work on specific structures has focused on (possibly approximately) obtaining quantized few-mode descriptions for plasmonic geometries such as a metal surface [10], metallic spheres [11][12][13][14] or sphere dimers [15,16]. It is thus desirable to develop tractable but versatile models using only a small number of EM modes. During the past few decades, there has been extensive work in this direction, with one notable development given by pseudomode theory [17][18][19].
Within nanophotonics increasing attention has recently focused on hybrid metallodielectric setups [20][21][22][23], with the objective of combining the strong field confinement and enhanced light-matter interactions of plasmonic resonances with the long-lived nature (large quality factors) of microcavity or photonic crystal modes. In these systems, EM field quantization is particularly complex due to the inherent coexistence of modes with very different properties and their mutual coupling. Quasinormal modes [22,[24][25][26] can be useful to unveil the EM mode structure, but due to their lossy nature, direct quantization remains challenging, and has only been carried out within very limited spectral windows [23,27]. A complementary technique developed very recently in the context of X-ray quantum optics is based on a partition of the physical space [28].
It is well-known that the interaction of a single emitter with an arbitrary EM environment can be described by means of the spectral density J(ω), which encodes the den- sity of EM modes and their coupling to the emitter. In this Letter, we present a simple and easily implementable framework for obtaining a few-mode quantum description of any given spectral density. We start from a macroscopic QED treatment, which already provides the basis of quantization in terms of a frequency continuum, and then construct a model system consisting of a discrete number of interacting modes coupled to independent flat background baths, see Fig. 1. Making use of Fano diagonalization [29,30], we obtain the general form for the model spectral density, J mod (ω). This can be fitted to any level of accuracy to J(ω), usually obtained by means of classical EM calculations. We illustrate the power and validity of this procedure in a hybrid structure comprising a plasmonic nanocavity and a high-refractive-index microresonator. We show that our approach enables accurate calculations of any far-and near-field observables, proving that replacing the full environment by our few-mode model does not lead to a loss of information.
For a single emitter, the EM mode basis in macroscopic QED can be chosen such that the Hamiltonian becomes ( = 1 here and in the following) [13,[31][32][33] where H e is the bare emitter Hamiltonian andμ e is its dipole operator, while a ω are the bosonic annihilation operators of the EM mode at frequency ω. These fulfill [a ω , a † ω ] = δ(ω − ω ) and correspond to the "true modes" in the nomenclature of Refs. [18,19]. The coupling between the emitter and the EM modes is encoded by where r e is the emitter position, n is the orientation of its dipole moment (for simplicity, we assume that all relevant transitions are oriented identically), and G( r, r , ω) is the classical dyadic Green's function. This is directly related to the spectral density of the environment, J(ω) = µ 2 g(ω) 2 for a given transition dipole moment µ [34]. As discussed above, Eq. (1) requires the treatment of an EM continuum. Without approximations, this is only possible with advanced and costly computational techniques, such as tensor network approaches [35,36]. Our goal is thus to construct an equivalent system that gives rise to dynamical equations that can be solved easily. Our model environment (sketched in Fig. 1) consists of N interacting EM modes with ladder operators a i , a † i , linearly coupled to the quantum emitter. Each of them is also coupled to an independent Markovian (spectrally flat) background bath. The resulting Hamiltonian is (3b) H S is the system (emitter + discrete modes) Hamiltonian, where the real symmetric matrix ω ij describes the mode energies and their interactions, and the real positive vector g i describes their coupling to the emitter. The bath Hamiltonian H B contains both the continuous bath modes, described by the bosonic operators b i,Ω and b † i,Ω at frequency Ω, and the coupling between the baths and the system, characterized by the rates κ i . The power of our approach lies in the fact that the Hamiltonian above can be analytically treated in two different ways: First, since the background baths are completely flat, the Markov approximation is exact and the dynamics described by H is identically reproduced, as proven recently [37], by a Lindblad master equatioṅ where ρ is the system density matrix and ρ} is a standard Lindblad dissipator. Secondly, the linearly coupled system of N interacting modes and continua can be diagonalized by adapting Fano diagonalization for autoionizing states of atomic systems [29], which in this context is related to the theory of quasimodes and pseudomodes [18,19]. This strategy allows us to obtain a simple, closed expression for J mod (ω).
We first formally discretize the EM continua and rewrite Eq. (3) as where A T = (a 1 , . . . , a N , b 1,Ω1 , . . . , b N,Ωn ) collects the annihilation operators of both modes and baths, the matrix H describes their mutual couplings, and the vector M T = (g 1 , . . . , g N , 0, . . . , 0) collects their couplings to the emitter (note that they vanish for the baths). Here, N is the number of discrete modes and associated continua, while n is the number of modes used to discretize each continuum. Diagonalizing H gives the energies E j and eigenmodes ψ j of the environment expressed in the basis of the model system, and determines the spectral Using the Sokhotski-Plemelj formula, we can write the spectral density in terms of the resolvent of Eq. (5) as Applying a formalism developed for the calculation of the absorption spectrum in atomic systems [30] and recovering the continuum limit for the baths (n → ∞), we finally obtain where g T = (g 1 , g 2 , . . . , g N ) is now an N -element vector and the N × N matrixH has entriesH ij = ω ij − i 2 κ i δ ij . Note that we have absorbed Lamb shifts of the modes due to the coupling with the baths into the mode frequencies ω ii (as we have also implicitly done in Eq. (4)). We note that as required for a spectral density, this form is non-negative, i.e., J mod (ω) ≥ 0 for all ω [38].
The last step in our approach consists in using Eq. (6) to fit J(ω) for a given EM environment, which allows us to parameterize Eq. (4) for that system. Although the number of unknowns in J mod (ω) is relatively large (N 2 + 2N real numbers for ω ij , κ i , and g i ), the fit procedure turns out to be stable even for large numbers of modes (up to N = 19 in the example below). A similar procedure based on the fitting of the bath correlation function has been reported as problematic recently [39]. Note that for non-interacting modes, ω ij = ω i δ ij , Eq. (6) simplifies to a sum of Lorentzians of the form J mod (ω) = i . This reproduces the well-known relation between Lorentzian spectral densities and lossy modes [17,40] that has been widely used to quantize simpler EM environments such as plasmonic cavities in the quasi-static approximation [10,12,14,15]. The introduction of interactions allows significantly more freedom in fitting J mod (ω), and in particular allows for the representation of interference effects and the associated Fano-like line shapes. In the supplemental material [38], explicit expressions of Eq. (6) for 2-4 interacting modes are presented, as well as their fitting to J(ω) in two recent studies [22,23]. In the following, we illustrate the power of our approach by considering a hybrid nanophotonic structure, with a significantly more complex spectral density. Fig. 2(a) shows the system under study: a 600 nm radius GaP [41] microsphere (ε sph = 9) embedding two 120 nm long silver nanorods (with permittivity taken from Ref. [42]) separated by a 3 nm gap, substantially displaced from the center of the sphere. The microsphere by itself supports many long-lived and delocalized Mie resonances, while the plasmonic dimer sustains confined surface plasmons with strongly sub-wavelength effective volumes [43]. The interaction between these different modes leads to a complex EM spectrum, shown in Fig. 2(b) through the Purcell factor P (ω) = J(ω)/J 0 (ω) for an emitter located in the center of the nanorods. Here, J 0 (ω) = ω 3 µ 2 6π 2 ε0c 3 is the spectral density in free space. The thick black line plots classical EM simulations performed with the Maxwell's equation solver implemented in COMSOL Multiphysics. This Purcell factor, and the corresponding J(ω), presents a large number of maxima, with several Fano-like profiles that indicate interference effects as typical for hybrid metallodielectric systems [22,23].
In order to obtain a stable fit of J mod (ω) to J(ω) despite the large number of parameters, we consider the noninteracting model (where ω ij = ω i δ ij ) as a starting point. This fit converges rapidly by using the spectral positions, curvature, and amplitudes of the local maxima in J(ω) to obtain initial guesses for the frequencies ω i , loss rates κ i , and coupling strengths g i of the non-interacting model. This gives a good fit for many of the peaks [light blue line in Fig. 2(b)], but strongly overestimates the background at lower frequencies, and fails to reproduce the Fano-like asymmetric profiles that originate from the hybridization between sphere and dimer modes. Using the non-interacting fit as a starting point for Eq. (6) leads to rapid convergence, and as shown by the orange line in Fig. 2(b), the resulting spectrum is in almost perfect agreement with the numerical Purcell factor over the full frequency range. Thus, we have constructed a compact model with a relatively small number of quantized interacting modes, each coupled to a flat background bath, that fully represents the EM environment in the nanophotonic structure in Fig. 2(a). Such an accurate fitting is not possible by means of non-interacting modes, at least not without significantly increasing the amount of modes considered. The thin grey lines in Fig. 2(b) indicate the (real part) of the eigenenergies of the matrixH. These correspond to the complex resonance positions (poles) of J(ω) [19].
We next demonstrate that the model system indeed gives a faithful representation of the EM environment, i.e., that the emitter dynamics with the model and with the original spectral density are equivalent. To do so, we treat the canon- ical spontaneous emission (Wigner-Weisskopf) problem for a two-level emitter initially in its excited state [44]. We thus haveĤ e = ω eg σ + σ − ,μ e = µ(σ + + σ − ), where σ ± are Pauli matrices. The emitter parameters are chosen to represent InAs/InGaAs quantum dots [45], with transition energy ω eg = 1.145 eV, indicated by the dashed red line in Fig. 2(b), and transition dipole moment µ = 0.55 e nm.
Since in the Wigner-Weisskopf problem, there is at most one excitation in the system (either in the emitter or in one of the EM modes), it can be solved easily for arbitrary spectral densities [40]. The exact excited-state population σ + σ − (t) obtained through this approach is shown in Fig. 3 (thick black line). The emitter dynamics in the model Lindblad master equation Eq. (4) is obtained using QuTiP [46]. The excited-state population in the full model (orange line) reproduces the exact results perfectly, while the non-interacting model (light blue) fails to do so and shows significant deviations from the correct result after about 100 fs.
In order to gain further insight into the relevance and meaning of the discrete modes obtained in our fit, we also show the populations of the modesã α obtained by diagonalizing ω ij . Since ω ij is a real symmetric matrix, this corresponds to an orthogonal transformation of the modes a i , withã α = i V αi a i [47]. The dashed lines in Fig. 3 show the mode populations of modesã 6 andã 7 , which we find to be the only significantly populated modes for the chosen emitter parameters. They are exactly the modes close to resonance to the emitter frequency,ω 6 = 1.129 eV andω 7 = 1.149 eV (see Fig. 2). The sum over all other mode populations, shown as a dashed magenta line in Fig. 3, remains small during the whole propagation. This demonstrates that our approach also allows the identification of the relevant modes in the dynamics.
Finally, we show that although the model system is written in terms of discrete lossy modes, it retains the full information about the EM near and far field. The electric field operator for the modes a ω within the formulation we use, Eq. (1), can be written as [48] where the field mode profile is given by Im{G( r, r e , ω)} · n. (8) The calculation then proceeds by solving the Heisenberg equations of motion for the mode operators a ω (t), which can be formally integrated to yield Inserting into Eq. (7) and defining the temporal kernel gives compact expressions for, e.g., the electric field intensity where, for simplicity, we have assumed that the field is initially in the vacuum state. Eq. (11) enables the calculation of the field intensity anywhere in space through the emitter correlation functions, which can be easily obtained from Eq. (4). This is displayed for the spontaneous emission case of Fig. 3 in Fig. 4, with panel (a) showing the temporal dependence of the electric field intensity at various points in space, with locations indicated by numbered white circles in panel (c). The intensities at each point are normalized to their maximum value, which is given in the figure caption. Panels (b) and (c) show snapshots of the field intensity profile in space at t = 10 fs (b) and t = 100 fs (c). A movie showing the full field distribution evolving in time is available in the supplemental material [38]. We note explicitly that in this spontaneous emission problem, there is no coherent field, E (+) + E (−) , and it is thus necessary to calculate the field intensity to observe the emission dynamics. Interestingly, the dynamics at points (1), next to the emitter, and (4) in the far-field (at a distance of 1.5 µm from the emitter), are quite similar, reaching their maximum value within a few tens of femtoseconds and then decaying rapidly. In contrast, points (2) and (3) inside and just outside the dielectric sphere (but at some distance to the nanorod dimer) show a much slower FIG. 4. Electric field intensity E (−) ( r) · E (+) ( r) for the initially excited emitter as in Fig. 3. (a) Field intensity at four points, shown by numbered white circles in panel (c), normalized to their maximum values Imax, given by 2.0 · 10 8 W/m 2 , 1.8 · 10 6 W/m 2 , 1.1 · 10 6 W/m 2 , and 4.0 · 10 3 W/m 2 , respectively. (b) Field intensity distribution in the x-z plane at time t = 10 fs, and (c) at t = 100 fs.
build-up and decay of the field intensity in time. The comparison against Fig. 3 reveals that the largest contribution to the field intensity at positions 2 and 3 is given by the hybrid modes 6 and 7, while the initial fast decay is due to the contribution of other modes and leads to an intense initial pulse radiated from the system, as seen at position 4.
To conclude, we have presented a simple and insightful procedure to quantize the electromagnetic field in arbitrary nanophotonic systems. Our approach works at the level of the spectral density, calculated through the solution of Maxwell's equations. This is fitted to a model spectral density, obtained through Fano diagonalization and involving only a small number of lossy and interacting electromagnetic modes. This makes it possible to construct and parameterize a few-mode master equation accurately describing the interaction of a quantum emitter with the original EM environment. We have illustrated the power and validity of our ideas by calculating the spontaneous emission population dynamics and near-and far-field intensity for an emitter placed within a hybrid structure comprising a dielectric microresonator and a plasmonic cavity. Our findings offer a versatile and easily implementable framework for the theoretical description of quantum nano-optical phenomena with Dyadic Green's function calculations as the single input.
The authors thank Diego Martín-Cano for interesting discussions and sharing their data. This work has been funded by the European Research Council through grant ERC-2016-StG-714870 and by the Spanish Ministry for Science, Innovation, and Universities -Agencia Estatal de Investigación through grants RTI2018-099737-B-I00, PCI2018-093145 (through the QuantERA program of the European Commission), and MDM-2014-0377 (through the María de Maeztu program for Units of Excellence in R&D). It was also supported by a 2019 Leonardo Grant for Researchers and Cultural Creators, BBVA Foundation. I. M. thanks the nanophotonics group for the warm hospitality during his stay at UAM, and acknowledges funding from the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior -Brazil (Capes) through grant 88887.368031/2019-00.