Description of ultrastrong light–matter interaction through coupled harmonic oscillator models and their connection with cavity-QED Hamiltonians

2.1 Derivation of the classical models from the Hamiltonians

In this subsection, we introduce the classical harmonic oscillator models. To this purpose, we first analyze the light–matter interaction using the cavity-QED framework. The cavity modes and the matter excitations are quantized using bosonic operators. The use of bosonic operators is valid for the cavity modes, and for matter excitations such as vibrations or phonons associated with a potential with a harmonic dependence on the degrees of freedom. The correspondence with classical harmonic oscillators (and thus the use of bosonic operators) is also valid to treat the coupling with matter excitations of fermionic nature provided that the number of excitations is much smaller than the number of quantum emitters (molecules, quantum dots…) and that any other effects induced by the saturation of the fermionic states can be discarded. Under these conditions, for example, the Quantum Rabi model (a generalization of the Jaynes–Cummings model to the ultrastrong coupling regime that includes a fermionic excitation [21]) becomes analogous to an appropriate bosonic Hamiltonian with a single matter excitation. Under this prescription based on bosonic operators, we can use a Hopfield-type Hamiltonian [48] in the form

as shown in the Supplementary Material. In this Hamiltonian, the creation operator a ̂ † and the annihilation operator a ̂ act on the cavity mode, while the equivalent operators b ̂ † and b ̂ are associated to the matter excitation, obeying commutation rules [ a ̂ , a ̂ † ] = [ b ̂ , b ̂ † ] = 1 . The first two terms on the right-hand side of Eq. (1) indicate the energy of the uncoupled (or bare) cavity modes and matter excitations at (angular) frequencies ω _cav and ω _mat, respectively, with ℏ the reduced Planck constant. The third term describes the light–matter coupling, which is parameterized by the coupling strength g _QED, and where we include the antiresonant terms a ̂ b ̂ and a ̂ † b ̂ † required to describe the ultrastrong coupling regime correctly. g _QED can in principle depend on ω _cav and ω _mat in specific systems (Section S6 in Supplementary Material). Last, we introduce the diamagnetic term, scaled by a parameter D that is initially considered to have an arbitrary value (including the zero value). This diamagnetic term, which is included in many (but not all) studies of ultrastrong coupling, is negligible in the strong coupling regime, but becomes important under ultrastrong coupling. It typically originates from the |A _⊥|² term of the minimal coupling Hamiltonian, where A _⊥ is the transverse vector potential. In the main text, we work in the Coulomb gauge, where the vector potential is completely transverse (∇⋅A = 0, and thus, A _⊥ = A), so that hereafter we omit the symbol ⊥ in A for brevity.

From the Hopfield Hamiltonian, we can obtain the equations of motion of the displacements (or oscillation amplitudes) of two quantum oscillators. With this aim, we connect the creation and annihilation operators from the Hamiltonian in Eq. (1) with the quantum operators x ̂ cav = ℏ 2 ω cav ( a ̂ + a ̂ † ) , x ̂ mat = ℏ 2 ω mat ( b ̂ + b ̂ † ) , p ̂ cav = − i ℏ ω cav 2 ( a ̂ − a ̂ † ) , and p ̂ mat = − i ℏ ω mat 2 ( b ̂ − b ̂ † ) . These operators correspond to the canonical position and momentum operators of harmonic oscillators (except that no mass has been included in their definitions). They fulfill the canonical commutation relations [ x ̂ mat , x ̂ cav ] = [ p ̂ mat , p ̂ cav ] = [ x ̂ mat , p ̂ cav ] = [ x ̂ cav , p ̂ mat ] = 0 , and [ x ̂ mat , p ̂ mat ] = [ x ̂ cav , p ̂ cav ] = i ℏ . The dynamics of these operators are calculated from the general Heisenberg equation of motion of an operator O ̂ , d d t O ̂ = 1 i ℏ [ O ̂ , H ̂ ] . We convert the four resulting first-order differential equations into two second-order equations by eliminating the momentum operators and obtain the following equations of motion for the expectation values ⟨ x ̂ cav ⟩ and ⟨ x ̂ mat ⟩ :

These are not the only classical equations that could describe the spectra of the coupled system. Any Hamiltonian H ̂ 2 related to the Hopfield Hamiltonian H ̂ 1 by a unitary transformation will have the same eigenfrequencies but will lead to different Heisenberg equations of motion. We perform a unitary transformation to H ̂ 1 with the operator U ̂ = e − i π 2 b ̂ † b ̂ . In the new reference frame, H ̂ 2 = U ̂ H ̂ 1 U ̂ † + i ℏ ∂ U ̂ ∂ t U ̂ † is expressed as:

where the prime ′ denotes that the matter operators are transformed ( b ̂ → i b ̂ ′ and b ̂ † → − i b ̂ † ′ ). In the representation of position and momentum operators, this transformation can be understood as a rotation in phase space so that the canonical variables transform as

In this new reference frame, we can calculate the equations of motion for the expectation values ⟨ x ̂ cav ⟩ and ⟨ x ̂ mat ′ ⟩ :

We find that, in contrast to Eq. (2), the coupling term is now proportional to the time derivative of the oscillation amplitudes.

To obtain the classical harmonic oscillator models, it is just necessary to associate the expectation values of the quantum operators to classical oscillation amplitudes, e.g., ⟨ x ̂ cav ⟩ → x cav , so that Eq. (2) becomes

and Eq. (5) becomes

where we do not make an explicit distinction between x mat ≡ ⟨ x ̂ mat ′ ⟩ used in Eqs. (5) and (7) and x mat ≡ ⟨ x ̂ mat ⟩ in Eqs. (2) and (6). However, the physical interpretation of the expectation values ⟨ x ̂ mat ⟩ and ⟨ x ̂ mat ′ ⟩ (or the oscillation amplitudes x _mat in each set of equations) is different, as discussed in more detail in Section 3 when applying each equation to specific coupled systems. Loss mechanisms are not included in these equations (friction terms proportional to the time derivatives x ̇ cav and x ̇ mat ), because they were derived from Hermitian cavity-QED Hamiltonians. Neglecting losses is usually an excellent approximation for calculating the eigenfrequencies and eigenvectors of the system in the ultrastrong coupling regime, where the coupling strength can be much larger than the system losses (the inclusion of dissipation in cavity-QED descriptions is discussed in Refs. [49], [50]).

Importantly, once the different interpretation of ⟨ x ̂ mat ⟩ and ⟨ x ̂ mat ′ ⟩ is accounted for, the two sets of coupled harmonic oscillator equations can be used to obtain the same result for any physical magnitude of a given system, as they correspond to Hamiltonians related by a unitary transformation. In particular, the eigenfrequencies of Eqs. (6) and (7) are identical.

Thus, it is always possible to obtain the optical response of the coupled system by considering the coupling to be proportional to either the oscillation amplitudes or their time derivatives. An important point to notice is that, in both Eqs. (6b) and (7b), the “matter resonant frequency” (square root of the term proportional to x _mat) is the bare frequency ω _mat. However, the “cavity resonant frequency” (square root of the term proportional to x _cav) is different in the different oscillator models. In the model characterized by Eq. (6a), the shifted cavity frequency is ω cav 2 + 4 D ω cav , while in Eq. (7a), the shifted cavity frequency is ω cav 2 + 4 D ω cav − 4 g QED 2 ω cav ω mat (shifts of the matter excitation are discussed in Section S2 of the Supplementary Material). In the following, we use the term dressed cavity/excitation (or dressed/renormalized frequency) when the shift is not zero, and thus the value of the shifted frequency does not coincide with the original value ω _cav before coupling (notice that ω cav 2 + 4 D ω cav = ω cav when using Eq. (6a) with D = 0 and ω cav 2 + 4 D ω cav − 4 g QED 2 ω cav ω mat = ω cav when using Eq. (7a) with D ω mat = g QED 2 , so that in these two cases we will refer to bare cavity frequencies).

In nanophotonics, coupled harmonic oscillator equations have been widely used to fit data without considering frequency renormalization so we adhere to this procedure, i.e., we consider harmonic oscillator models where the frequency of the cavity and matter excitations are the bare ones. This approach gives preference to the model with coupling constant proportional to the oscillation amplitude (Eq. (6)) or to its derivative (Eq. (7)), depending on the value of D, as discussed next. Thus, throughout the remaining of this paper (including the Supplementary Material unless otherwise stated), we analyze these two preferred models using bare frequencies. We denote these models the Spring Coupling (SpC) model and the Momentum Coupling (MoC) model, respectively. Other models are discussed in Section S2 and summarized in Section S4 of the Supplementary Material. Additionally, Section S1 of the Supplementary Material details how to obtain the classical coupled harmonic oscillator equations directly from the classical electromagnetic Lagrangian.

2.2 Spring Coupling (SpC) model

We consider first a system without diamagnetic term, D = 0. This choice is appropriate, for example, when the interaction between the emitter and cavity excitations is mediated by Coulomb coupling, as discussed in more detail in Section 3.2. Eq. (6) then becomes

where the coupling is proportional to the classical oscillation amplitudes x _cav and x _mat and we have changed the notation g _SpC = g _QED (using a different symbol for the coupling strength in the classical and quantum descriptions becomes useful in Section 2.3). The ω cav ω mat prefactor appears directly from the Hamiltonian and ensures that g have units of frequency. Other choices of prefactor have been used (such as using the arithmetic mean of the bare frequencies instead of the geometric mean [51]), which are equivalent in the strong coupling regime but not in the ultrastrong one. In the frequency (ω) domain, these equations are transformed to

We refer to Eqs. (8) and (9) as the Spring Coupling (SpC) model because they are analogous to the equations describing the movement of two coupled springs (sketch in Figure 2(a)) (we emphasize that we could also describe the same physics of ultrastrongly coupled systems by setting D = 0 in Eq. (7), but, in this case, the dressed frequency ω cav 2 − 4 g SpC 2 ω cav ω mat would appear in the equations instead of the bare one, contrary to our previous choice). The eigenfrequencies ω _±,SpC of the SpC model are obtained by diagonalizing the matrix associated with Eq. (9), which leads to

Figure 2:

Comparison of the Spring Coupling (SpC) and Momentum Coupling (MoC) models. (a) Schematics of the SpC model in analogy to an oscillator model in classical mechanics. The coupling mechanism of strength g _SpC is analogous to a force F _SpC exerted by a spring and proportional to the oscillator displacements x _cav and x _mat. (b) Schematics of the MoC model. The coupling mechanism of strength g _MoC is analogous to a force F _MoC proportional to the time derivatives of the oscillator displacements x ̇ cav and x ̇ mat . We represent the coupling with dashed lines to highlight the different coupling mechanism compared with the SpC model, but we are not aware of any system described by the MoC model in classical mechanics. (c) Eigenfrequencies ω _± of the hybrid states calculated from the bare values ω _cav and ω _mat, with ω _mat fixed and ω _cav/ω _mat changing. ω _± are obtained from the SpC model (blue solid line, corresponding to Eq. (10)) and the MoC model (red dashed line, Eq. (13)) for coupling strength g = g _SpC = g _MoC = 0.1 ω _mat. The thin gray lines correspond to the bare cavity frequency ω _cav and the bare frequency of the matter excitation, ω _mat. (d) Same as panel (c), for coupling strength g = g _SpC = g _MoC = 0.3 ω _mat. (e) Minimum splitting between the hybrid modes Ω^min = ω ₊ − ω ₋, as a function of the coupling strength g for the SpC model (blue solid line) and the MoC model (red solid line). All frequencies in panels (c–e) are normalized with respect to the fixed frequency of the matter excitation ω _mat, so that the results do not depend on the particular value of ω _mat, only on the ω _cav/ω _mat and g/ω _mat ratios.

We note that frequencies given by Eq. (10) correspond to the energy difference between the first excited and ground state, and not to the absolute values of the eigenfrequencies themselves. This distinction is not necessary in classical descriptions that set the energy of the ground state to zero (or a fixed value). However, the cavity-QED model indicates a g _QED-dependent shift of the ground-state energy from zero, which is a fully quantum phenomenon. The information of this shift is lost when we take the expectation value of the operators ⟨ x ̂ cav ⟩ and ⟨ x ̂ mat ⟩ in Eq. (2). The g _QED dependence of this shift can be found, for instance, in Figure 2(f) of Ref. [21].

2.3 Momentum Coupling (MoC) model

For a diamagnetic term with D = g QED 2 ω mat (this value normally appears in atomic physics and in cavity-QED models [23] in the Coulomb Gauge and is discussed in Ref. [21] and Section 3.1), Eq. (7) takes the form

with the coupling term proportional to the time derivative of the oscillation amplitudes (the “velocities”) so that we call this model the Momentum Coupling (MoC) model (sketch in Figure 2(b)). The coupling strength g _MoC in these equations is related to the constant g _QED in the cavity-QED Hamiltonian as g MoC = ω cav ω mat g QED (and thus D = g QED 2 ω mat = g MoC 2 ω cav ). We introduce this new coupling strength because, in this way, (i) Eqs. (11a) and (11b) take the same form as in previous work [52], [53], [54] and (ii) g _MoC becomes independent of the resonant frequencies ω _mat and ω _cav for the systems studied in Section 3. However, it is also possible to write Eqs. (11a) and (11b) in terms of g _QED as long as one is consistent in all the derivation. We further emphasize that g _MoC = g _QED in resonance (ω _cav = ω _mat), and these two parameters only take significantly different values for strong detuning. In the frequency domain, Eq. (11) becomes

and the corresponding eigenfrequencies are

Although the MoC is used to describe the coupling between matter excitations and cavity modes (Figure 2(b)), we are not aware of any equivalent mechanical system in classical mechanics that follows the equations of motion in Eqs. (11a) and (11b) (with coupling terms proportional to the time derivatives of the oscillation amplitude, similarly to friction terms but describing the interaction between two different oscillators). This is in contrast to the SpC model where the equivalent system, composed of masses and springs, is shown in Figure 2(a).

2.4 Comparison of the MoC and SpC models

As mentioned above, the MoC and SpC models are appropriate when D = g QED 2 / ω mat and D = 0, respectively (we emphasize again that the resonant frequencies ω _cav, ω _mat in these models are the bare resonant frequencies). Further, regardless of whether the diamagnetic term should be included in the description or not, these models have been used in the past as phenomenological tools for extracting coupling parameters by fitting the spectra of the coupled system obtained from experimental data or simulations [17], [18], [51], [52], [53], [54], [55], [56], [57], [58], [59]. In this section, we compare the results provided by both models as a function of the coupling strength.

The MoC and SpC models are known to give very different results for g ≫ 0.1ω _mat, as we briefly illustrate in this section (we use g in this subsection to refer to g _SpC or g _MoC in discussions that are valid for both models). Figure 2(c) compares the eigenfrequencies of the SpC (blue solid line) and MoC (red dashed line) models for g = 0.1ω _mat, as given by Eqs. (10) and (13), respectively. For simplicity, we consider that the coupling strength is the same for all values of ω _cav (a different parameter choice is discussed in Section S6 of the Supplementary Material). The eigenfrequencies of the hybrid modes ω _± are calculated as a function of the bare cavity frequency ω _cav, with the bare ω _mat frequency fixed (all frequencies are normalized by ω _mat, so that the figures are independent of the value of this parameter). The eigenfrequencies obtained within the MoC and SpC models follow a nearly identical dependence on ω _cav, and the agreement is even better for g < 0.1 ω _mat. Thus, when analyzing systems not in the ultrastrong coupling regime, the two models can generally be used interchangeably with minimal impact on the results, although exceptions can exist [60].

In contrast, the choice of the model is crucial for even larger coupling strengths, such that the system is well into the ultrastrong coupling regime. The differences between the two models are illustrated in Figure 2(d) for coupling strength g = 0.3 ω _mat. In this case, the two models predict significantly different eigenfrequencies of the coupled system. The difference is smaller for larger cavity frequencies, ω _cav ≫ ω _mat, because the oscillators become uncoupled and the eigenfrequencies approach the bare frequencies ω _cav and ω _mat in the two models. However, even for a relatively large ω cav ω mat = 1.5 , the difference between the values of ω _± according to the two models is around 10 %.

We compare next the splitting Ω = ω ₊ − ω ₋ between the two eigenmodes at zero detuning, ω _cav = ω _mat. In the MoC model, the splitting equals twice the coupling strength, i.e., Ω = 2g, which is the minimum splitting in this model [61], [62]. On the other hand, in the SpC model, the relation between Ω and the coupling strength for zero detuning is

We find Ω_SpC = 2.11g _SpC for the values used in Figure 2(d). Further, according to the SpC model, the minimum splitting between the branches does not happen at zero detuning but at cavity frequencies larger than the matter excitation frequencies. To further emphasize the difference between the models, Figure 2(e) shows the minimum splitting as a function of coupling strength, with a linear dependence for the MoC (red solid line) model, Ω^min = 2g, in contrast with the strong deviation from nonlinearity of the SpC model results (blue line) for g ≫ 0.1ω _mat. As a consequence, close to the so-called deep strong coupling regime g ω mat ≈ 1 , Ω SpC min is approximately twice that of the MoC model.

Last, Figure 2(d) shows important differences at small cavity frequencies, ω _cav ≪ ω _mat. The dispersion of the MoC model shows two hybrid modes for all values of the detuning, with the lower mode frequency ω _−,MoC tending toward ω _cav for decreasing value of ω _cav. In contrast, for the SpC model, the lower mode ceases to exist (ω _−,SpC becomes imaginary) under the condition ω cav ω mat < 2 g SpC ω mat 2 (for fixed g _SpC = 0.3ω _mat; see Section S6 of the Supplementary Material where a different choice is discussed). Further, in the SpC description, the upper branch approaches the bare matter frequency at ω _cav → 0, but this is not the case in the MoC model, where the corresponding asymptotic limit is ω + , MoC = ω mat 2 + 4 g MoC 2 . Thus, in the MoC model, the coupling affects the upper hybrid mode even in this highly detuned situation.

The two models’ different asymptotic limits of the upper branch determine the predicted range of energies where hybrid modes can exist. The MoC results show a frequency band between ω _mat and ω mat 2 + 4 g MoC 2 with no modes available. This forbidden band is not present in the SpC dispersion. In Section 3 and Section S8 of the Supplementary Material, we connect this result with the Reststrahlen band of polar materials and show that we can reproduce the experimental dispersion of these materials by using the MoC [45] and alternative models but not the SpC model.

The connection between classical and quantum models is summarized in Table 1. The classical SpC and MoC models result in the same eigenfrequencies as cavity-QED Hamiltonians without the diamagnetic term (D = 0) and with D = g QED 2 ω mat = g MoC 2 ω cav , respectively. Other classical coupled harmonic oscillator models where dressed frequencies are used instead of the bare ones (with an associated change of the coupling term) are discussed in Section S2 of the SI. For completeness, we also discuss in Section S5 of the Supplementary Material an often-used linearized model that is a good approximation to the MoC and SpC models for low coupling strengths (especially for the anticrossing region of the spectrum corresponding to small detunings). However, this linearized model is not appropriate in the ultrastrong coupling regime.

Table 1:

Summary of the correspondences of the classical SpC and MoC models with the cavity-QED description without diamagnetic term D = 0 (second column) and with diamagnetic term and D = g MoC 2 ω cav (third column). The second row shows the two considered cavity-QED Hamiltonians. The third row indicates the equations of motion of the oscillation amplitudes x _cav and x _mat obtained with the classical SpC (second column) and MoC (third column) harmonic oscillator models. The fourth row provides the frequencies of the two resulting hybrid modes, which are the same for the cavity-QED and classical models for the value of D and choice of coupled harmonic oscillator model indicated in each column. The last row indicates the relationship between the coupling constant g _QED in the cavity-QED Hamiltonian and those in the classical coupled harmonic oscillator models (g _MoC and g _SpC). For the system in Section 3.1, g _MoC is constant and thus g QED ∝ ω mat / ω cav .

At this point, we have discussed the connections between a general quantum description and classical equations of motion. However, we still need to determine how to choose between the MoC and SpC models for a given system (or equivalently, whether the Hamiltonian has D ≠ 0 or D = 0). In the next section, we consider three representative systems to explore this question and highlight the key role played by the nature of the matter–cavity interaction (Coulomb coupling or coupling with transverse electromagnetic modes in dielectric cavities).

Additionally, we have focused thus far on the eigenfrequencies, which can be extracted directly from the equations of coupled harmonic oscillators without needing an exact understanding of what the oscillation amplitudes x _cav and x _mat represent. However, a clear physical interpretation of these parameters is necessary to evaluate magnitudes of interest, such as the electric field at a given location inside or outside the optical cavity. In Section 3, we also address how x _cav and x _mat relate to relevant physical quantities in the representative systems of choice.

3.1 A quantum emitter interacting with a transverse mode of a dielectric cavity

We consider first a dipole interacting with a single transverse mode of a resonant dielectric cavity (Figures 1(a) and 3(a)). The dipole is associated with matter excitations, and it can represent an excitonic transition of a molecule or quantum dot or a transition between vibrational states, for example. For concreteness, we consider the coupling with a molecular emitter in the following. Cavity-QED models of this system have successfully described phenomena such as the modification of the spontaneous emission rate of the emitter [63], [64], of the photon statistics of the emitted light [60], [65], [66], or of the coherence time of the quantum states [67].

Figure 3:

Interaction of a quantum emitter with a transverse cavity mode within the classical MoC model. (a) Schematics of the system. The two oscillators are associated with the vector potential A of the cavity mode and the induced dipole moment d of the excitation in the quantum emitter, which we consider to be a molecule. The oscillators are coupled with each other with strength g _MoC. The bottom sketch indicates the cavity dimensions that we analyze in the rest of the panels. The emitter is placed at the center of the cavity. The green shaded areas in the sketches represent the field distribution of the cavity mode. (b) Spatial distribution of the electric field for the upper (blue) and the lower (red) hybrid modes at frequencies ω _+,MoC and ω _−,MoC, respectively, for coupling strength g _MoC = 2.5 ⋅ 10⁻⁴ ω _cav. The electric field is calculated along the cavity axis (along the x direction in panel (a), with x = y = z = 0 corresponding to the cavity center). The inset is a zoom of the region near the emitter. (c) Contribution to the electric field from the cavity Σ cav ± (dots) and from the emitter Σ mat ± (solid lines), for the hybrid mode at frequency ω _+,MoC (blue) and the hybrid mode at frequency ω _−,MoC (red), as a function of the detuning ω _mat − ω _cav. The fields are real and are evaluated at the position (x, y, z) = (10.5 nm, 0, 0), i.e., at 10.5 nm distance from the center of the cavity where the molecular emitter is located (see sketch in (a) for directions), which corresponds to the position indicated by the dashed line in the inset of panel (b). The coupling strength is g _MoC = 2.5 ⋅ 10⁻⁴ ω _cav. (d) Same as in (c), for g _MoC = 0.2 ω _cav.

The whole derivation of the equations of motion of the classical variables within the Coulomb gauge is discussed in the Supplementary Material (Section S1), but we summarize it in the following. We represent the molecular emitter as two point charges with relative position l (forming a dipole), which couple through Coulomb interactions determined by the potential V _Cou(l) approximated as a harmonic one, V Cou ( l ) = 1 2 m red ω mat 2 l 2 , with l = |l| the distance and m _red the reduced mass of this two-body system. The dipole moment induced in the molecular emitter is d = q l, where q is the absolute value of the charge of the particles in the dipole. On the other hand, the cavity mode is characterized by the vector potential A, which is the canonical position variable of the transverse electromagnetic fields [68]. This description does not include nonlinear effects, being thus valid for harmonic molecular vibrations, and also for anharmonic vibrations or excitonic transitions under weak illumination.

In cavity-QED models, the standard approach to describe light–matter interactions in this system is to use the minimal-coupling classical Hamiltonian in the Coulomb gauge of the form H min - c = q 2 ( p − A ) 2 2 m red , where p is the classical canonical momentum of the dipolar matter excitation. To obtain the quantum Hamiltonian, we use the following quantization relations [48], [68]:

where Π ̂ ( r ) is the canonical momentum associated to the vector potential A ̂ ( r ) (see Sections S1, S2 of the Supplementary Material). The function Ξ(r) accounts for the spatial distribution of the vector potential of the cavity mode and is normalized so that its maximum value is 1. Further, we have introduced the effective mode volume of the cavity field [69], V _eff, and the oscillator strength of the dipolar excitation f mat = q 2 m red . From the minimal-coupling Hamiltonian H _min-c, the light–matter interaction term is H int = − q 2 p ⋅ A m red . Considering that the induced dipole moment and the vector potential form an angle θ, and using Eqs. (15) and (18), the interaction term of the quantized Hamiltonian becomes

Comparing this expression with the third term of the Hopfield Hamiltonian (Eq. (3)), we directly obtain that the coupling strength in the cavity-QED formalism is g QED = 1 2 f mat ε 0 V eff ω mat ω cav Ξ ( r ) cos ⁡ θ . We consider from now on the maximum coupling strength g QED = 1 2 f mat ε 0 V eff ω mat ω cav , which is achieved in the position of the maximum field (Ξ(r) = 1) for optimal orientation (θ = 0). Further, the A ² term in the minimal-coupling Hamiltonian leads to a diamagnetic term (fourth term on the right-hand side of Eq. (3)) with D = g QED 2 ω mat . Following the discussion of Section 2, the presence of the diamagnetic term in the cavity-QED Hamiltonian with this exact value of D indicates that this system can be described by the classical MoC model in Eq. (11).

Next, we use the connection between the classical and cavity-QED approaches to illustrate the procedure to obtain the value of physical observables from the classical oscillation amplitudes of the cavity x _cav and of the molecular excitation x _mat. The classical coupling strength g _MoC is directly obtained from the quantum value as g MoC = ω cav ω mat g QED = 1 2 f mat ε 0 V eff . Further, the quantum position operators ( ∝ a ̂ + a ̂ † and ∝ b ̂ + b ̂ † ) and the classical oscillation amplitudes (x _cav and x _mat) are related by the standard quantization relationship

where, for an appropriate comparison between classical amplitudes and quantum operators, the real part of the oscillator amplitudes must be taken: Re(x _cav) = Re(|x _cav|e^−iωt+ϕ ) ∝|x _cav| cos(ωt + ϕ), with ϕ a phase. Equations (20a) and (15) indicate that the oscillation amplitude x _cav in the MoC model (Eq. (11)) is given by x cav = A ε 0 V eff , where A is the maximum amplitude of the classical vector potential (i.e., in the position where Ξ(r) = 1). Therefore, the oscillation amplitude x _cav can be used to calculate the spatial distribution of this potential as A ( r ) = A Ξ ( r ) = x cav ε 0 V eff Ξ ( r ) ( A ( r ) = ⟨ A ̂ ( r ) ⟩ is the classical counterpart of the quantum operator of the vector potential). Equivalently, from Eqs. (20b) and (17), the amplitude of the oscillator corresponding to the matter excitation is directly connected with the induced classical dipole moment (d = |d|) as x mat = d f mat . These relations are schematically shown in Figure 3(a).

We are finally in conditions to obtain the value of physical observables such as the electric field from the classical harmonic MoC model. We first consider the spatial distribution of the electric fields of each hybrid mode. The transverse cavity mode field (given by A(r, t)) must be added to the longitudinal near field induced by the induced dipole,^[1] which is obtained from the scalar Coulomb potential

with unit vectors n d = d | d | and n r = r | r | . The total electric field is, therefore, given as

and the electric field at frequencies ω _±,MoC of each hybrid mode (given by Eq. (13)) corresponds to

with n A = A ( r ) | A ( r ) | . This equation indicates the contribution of the cavity E _cav(r, ω _±,MoC) ∝ x _cav and of the matter excitation E _mat (r, ω _±,MoC) ∝ x _mat to the electric field. Further, we use the eigenvectors (Eq. (12a)) to obtain the ratio between the amplitudes x _cav and x _mat of the classical harmonic oscillators:

Inserting Eq. (24) into (23), we obtain the ratio between the contributions of the cavity electric field and the matter excitation.

Equations (23) and (24) are a main result of this subsection and can be used to obtain the electric field at any position and for an arbitrary transverse mode with field distribution given by Ξ(r). We consider for illustration the particular case of a molecule (as an example of quantum emitter) introduced in the center of a dielectric cavity consisting in a rectangular vacuum box enclosed in the three dimensions by perfect mirrors, as sketched in Figure 3(a). The cross-section of the box is square, with size L _x = L _y = 292 nm and its height is L _z = 215 nm, which results in a fundamental lowest-order cavity mode at frequency ω _cav = 3 eV and an effective volume V _eff = 4.483 ⋅ 10⁶ nm³ (for an easier comparison between classical frequencies ω and quantum energies ℏω, in this paper, we use eV as a unit for both of them). This value of V _eff is calculated from the general expression of dielectric structures [70]

where ɛ(r) refers to the permittivity of the system at position r, and in this particular case, we consider ɛ(r) = 1 inside the cavity. The molecular excitation is nearly resonant with the cavity, ω _mat ≈ ω _cav = 3 eV, but its exact frequency is changed to study the effects of detuning. The transition dipole moment μ mat = ℏ f mat 2 ω mat (associated with the transition from the ground state to the first excited state) is parallel to the z axis and is relatively strong, μ _mat = 15 Debye, achievable with organic molecules such as nonacene, for example, [71]. This value of the transition dipole moment implies that this molecular emitter has an oscillator strength of f mat = ( 118.74 e ) 2 m p , where e is the electron charge and m _p the mass of the proton. By placing the molecular emitter in the center of the cavity where the electric field of the mode is maximum, this choice of parameters leads to a coupling strength g _MoC ≈ 2.5 ⋅ 10⁻⁴ ω _cav, far from the ultrastrong coupling regime (a larger value of g _MoC is considered at the end of this subsection).

We show in Figure 3(b) the distribution of the z component of the electric field inside this cavity for the upper hybrid mode E _z (x, ω _+,MoC) and for the lower hybrid mode E _z (x, ω _−,MoC), as obtained from Eq. (23). We plot the fields as a function of the position in the x direction with respect to the location of the molecular emitter at the center of the cavity. To highlight the differences between the contributions of the cavity and the induced dipole in the two modes, we choose a slight detuning of ω _cav − ω _mat = 1.5 meV. Since the classical MoC model does not give the absolute value of the eigenmode fields, we choose arbitrary units so that the contribution of the cavity mode to the electric field of the upper hybrid mode (E _cav(r, ω _+,MoC) in Eq. (23)) has a maximum absolute value of 1. This choice fixes all the other values according to Eq. (24).^[2] The fields are dominated by the cavity mode far from the cavity center and by the contribution from the molecular dipole close to x = 0. The field distribution shows a clear difference in the behavior of the two hybrid modes. For the upper mode, the induced dipole points in the same direction as the cavity field ( x cav ( ω + , MoC ) x mat ( ω + , MoC ) > 0 ) , but in the inverse direction for the lower mode ( x cav ( ω − , MoC ) x mat ( ω − , MoC ) < 0 ) . Further, at the detuning considered, the relative contribution of the cavity to the fields is larger for the upper than for the lower mode, as indicated by the values of the electric field far from the molecular emitter at ω _+,MoC and ω _−,MoC. In contrast, as shown in the inset, the relative contribution from the molecular dipole to the field close to the molecule (x = 0) is stronger for the lower mode. Figure 3(b) thus confirms that the classical harmonic oscillator model allows for calculating the relative contribution of cavity and matter for each mode, as desired.

Further, Eqs. (23) and (24) also enable to examine the dependence of the field E(r, ω _±,MoC) inside the cavity with detuning ω _mat − ω _cav. Figure 3(c) shows the contributions to this electric field of the cavity and the molecular emitter for each hybrid mode, normalized with respect to the sum of both contributions, according to Σ cav ± = | E cav ( ω ± , MoC ) | 2 | E cav ( ω ± , MoC ) | 2 + | E mat ( ω ± , MoC ) | 2 (dots) and Σ mat ± = | E mat ( ω ± , MoC ) | 2 | E cav ( ω ± , MoC ) | 2 + | E mat ( ω ± , MoC ) | 2 (solid lines). These ratios play a similar role as the Hopfield coefficients from cavity-QED descriptions. The blue (red) dots and solid lines correspond to the upper (lower) hybrid mode. We obtain E _cav(ω _±,MoC) and E _mat(ω _±,MoC) by replacing Eq. (24) into (23), for a fixed coupling strength g _MoC = 2.5 ⋅ 10⁻⁴ ω _cav and for a distance of 10.5 nm from the molecular emitter in the x direction. This position (indicated by the dashed line in the inset of Figure 3(b)) is chosen because it is where the matter and cavity contributions have the same weight for the two hybrid modes at zero detuning and very small coupling strengths ( Σ cav ± = Σ mat ± = 0.5 ). For ω _cav > ω _mat, the field of the lower mode is predominantly given by the matter excitation ( Σ mat − > Σ cav − as indicated by the red dots and the red solid line). In contrast, for the upper mode, the cavity contribution dominates ( Σ cav + > Σ mat + , blue). Further, already at detunings as small as ω _cav − ω _mat ≳ 15meV = 5 ⋅ 10⁻³ ω _cav, the modes are essentially uncoupled for this small coupling strength ( Σ mat + ≪ Σ cav + and Σ mat − ≫ Σ cav − ).

The coupling strength we have considered in this subsection corresponds to the strong coupling regime (we have neglected losses) but is far from the ultrastrong coupling regime so that the phenomena studied can also be explained with the classical linearized model (Section S5 of the Supplementary Material). On the other hand, we consider again in Figure 3(d) the contributions to the electric field Σ cav ± and Σ mat ± as a function of the detuning, but in this case for a considerably larger coupling strength g _MoC = 0.2ω _cav. This value of g _MoC is not currently achievable with dielectric cavities at the single molecule or single emitter level (it would correspond to a transition dipole moment μ _mat = 1.2 ⋅ 10⁴ Debye), but we choose it to illustrate the analysis of ultrastrongly coupled systems within the classical MoC model. Further, such large g _MoC can be achieved in dielectric cavities fully filled by a material or many molecular emitters, as discussed in Section 3.3. For zero detuning ω _cav = ω _mat, the contributions of the induced dipole and the cavity are no longer identical in the ultrastrong coupling regime, with Σ cav + ≈ 0.6 and Σ mat + ≈ 0.4 for the upper hybrid mode at frequency ω _+,MoC (and the opposite for the lower hybrid mode). More strikingly, the results in Figure 3(d) indicate a very different tendency of the modes at large detunings as compared to strong coupling, especially in the case of the upper hybrid mode at frequency ω _+,MoC. In the ultrastrong coupling regime, in the ω _mat → 0 limit (ω _mat − ω _cav → −3 eV), this mode (blue solid line and dots) has significant contributions from both the cavity and the matter ( Σ cav + ≈ 0.9 and Σ mat + ≈ 0.1 ). Thus, these two excitations do not decouple in this limit. This behavior is consistent with the discussion of the dispersion in Figure 2(d), where at large detunings, the upper mode frequency does not reach the bare frequency ω _mat. The SpC model (not shown) does not reproduce this behavior because the modes become uncoupled ( Σ cav + ≈ 1 and Σ mat + ≈ 0 ).

The described methodology thus enables obtaining results equivalent to those of the cavity-QED description (Hopfield Hamiltonian with the diamagnetic term) by using an intuitive classical model of coupled harmonic oscillators. In summary, we have shown in this section how to use the classical MoC model to characterize the fields in a hybrid system composed of a molecular emitter coupled to a transverse mode of a cavity.

3.2 A quantum emitter interacting with the longitudinal field of a metallic nanoparticle through Coulomb coupling

Next, we consider a quantum emitter placed close to a metallic nanoparticle to analyze how to model an alternative system and obtain physical observables in the strong and ultrastrong coupling regimes. These nanoparticles are attractive in nanophotonics because they support localized surface plasmon modes characterized by very low effective volumes [18], [71], [72], [73], [74]. Since the coupling strength is inversely proportional to the square root of the effective mode volume, very large coupling strengths can be obtained even when the nanoparticle interacts with a single molecule or quantum dot. We consider again a molecule as a representative quantum emitter.

In order to analyze the interaction of the dipolar plasmonic mode of the nanoparticle with a molecular (harmonic) excitation of dipole moment d _mat, we consider that the size of the nanoparticle and the molecule-nanoparticle distance are much smaller than the light wavelength and treat the system within the quasistatic approximation. Under this approximation, the temporal variation of the vector potential A in Eq. (22) is negligible. Therefore, the coupling between the nanoparticle and the molecular emitter is governed by Coulomb interactions expressed by a scalar potential V _Cou. The coupling is then mediated by longitudinal fields, in contrast to the coupling with transverse fields in Section 3.1.

In this context, the emitter-nanoparticle coupling cannot be modeled with the minimal coupling Hamiltonian as in Section 3.1, and it is described instead through the interaction Hamiltonian [75]

E ̂ cav ‖ is the electric field associated with the dipolar mode of the nanocavity, which in the quasistatic approximation is completely longitudinal (we indicate this explicitly with the symbol ‖) and d _cav is the induced plasmonic dipole moment (operator d ̂ cav ). For simplicity, we consider small spherical particles of radius R _cav composed by a Drude metal with plasma frequency ω _p, but this approach could be generalized to other systems. The spherical particles present a dipolar plasmonic resonance of Lorentzian lineshape at frequency ω cav = ω p 3 , and oscillator strength f cav = 4 π ε 0 R cav 3 ω cav 2 [76]. The quasistatic field outside them is directly determined by d _cav according to E ̂ cav ‖ ( r ) = 3 ( d ̂ cav ⋅ n r cav ) n r cav − d ̂ cav 4 π ε 0 | r cav − r | 3 , where r _cav is the center of the nanoparticle, |r − r _cav| > R _cav, and we define the unit vector n r cav = r − r cav | r − r cav | .

We insert the quantized expressions of the induced dipole moments d ̂ cav and d ̂ mat of Eq. (17) into the Hamiltonian in Eq. (26) and the expression of E ̂ cav ‖ ( r ) and obtain

with coupling strength

where we have defined the unit vectors as n d cav = d cav | d cav | , n d mat = d mat | d mat | , and n r rel = r cav − r mat | r cav − r mat | . The total Hamiltonian is thus the sum of H ̂ int 2 and the terms related to the energy of the uncoupled plasmon and molecular excitation, corresponding to the Hopfield Hamiltonian of Eq. (1) without the diamagnetic term (D = 0). Thus, the corresponding classical model to be adopted is the SpC model in Section 2.2, with the equations of motion in Eq. (8). Additional details can be found in Section S1 of the Supplementary Material.

The representation of the plasmon-molecule system with the SpC model is schematically shown in Figure 4(a). To obtain the observables in this system, we use the equivalence of the oscillation amplitudes x _cav and x _mat with the induced dipole moments of the cavity and the molecular (or matter) excitation. This equivalence can be obtained from Eqs. (17) and (20), and it follows x cav = d cav f cav and x mat = d mat f mat . Further, this treatment can be extended to other dipole–dipole interactions beyond the coupling of a molecular emitter with a plasmon (direct dipole–dipole interactions between molecules are considered in Section 3.3).

Figure 4:

Modeling of the coupling between a quantum emitter and a metallic spherical nanoparticle (a plasmonic nanocavity) within the classical SpC model. (a) Schematics of the system. The quantum emitter is considered to be a molecule. The molecular excitation (of induced dipole moment d _mat) and the dipolar mode of the plasmonic nanocavity (of induced dipole moment d _cav) are described as two harmonic oscillators (of oscillation amplitudes x _mat and x _cav) that are coupled with strength g _SpC. The system is excited by a laser of electric field amplitude E _inc. The radius of the spherical nanoparticle is 5 nm, and the molecular emitter is placed at a 1 nm distance from the nanoparticle surface along the x axis (the center of the nanoparticle corresponds to x = y = z = 0). d _cav, d _mat, and E _inc are polarized along x. (b) Electric field distribution along the x axis (y = z = 0) when the system is excited at the frequency of the upper hybrid mode ω _+,SpC (top panel) and of the lower hybrid mode ω _−,SpC (bottom panel). The fields are real and are evaluated only outside the nanocavity, with the positions inside highlighted by the green-shaded area. The position of the molecular emitter is indicated by the vertical brown line. We evaluate the fields for coupling strength g _SpC = 0.1 ω _cav, and ω _cav = ω _mat = 3 eV. For each hybrid mode, the cavity contribution to the field is indicated by dots, the contribution from the emitter by dashed lines, and the total field by blue solid lines. (c) Scattering cross section of the same system, with g _SpC = 0.1 ω _cav, as a function of the detuning of the laser ω − ω _cav. Solid lines: tuned system with frequencies ω _cav = ω _mat = 3 eV. Dashed lines: detuned system with frequencies ω _cav = 3 eV and ω _mat = 3.2 eV. (d) Scattering cross section of the tuned system (ω _cav = ω _mat = 3 eV), comparing the result of the SpC model (blue line) to the results of the MoC model (black line), in the strong coupling regime, g = 10⁻² ω _cav. (e) Same as in (d) for the ultrastrong coupling regime, g = 0.3 ω _cav. In all results f _cav = (4345e)²/m _p (where m _p is the mass of the proton), F cav = f cav | E inc | , f _mat = (118.74e)²/m _p, F mat = f mat | E inc | , κ = 20 meV and γ = 10 meV (except that we modify f _cav in panels (d) and (e) to achieve the desired values of g _SpC).

We consider next that the dipolar mode of the metallic nanoparticle is illuminated by an external field of amplitude E _inc and frequency ω. We introduce this field in the SpC model as a forcing term that acts both onto the nanoparticle and onto the molecular emitter. Specifically, this is done by adding terms F α e − i ω t = f α | E inc | e − i ω t (α = “cav” or α = “mat”) on the right-hand side of Eq. (8), i.e., the amplitude F _α of the time-dependent force is proportional to the induced dipole moments d _α and the electric field of the illumination (see Section S1 in the Supplementary Material for further details). By solving the equations of motion of the SpC model (Eq. (8)) in the frequency domain with this external force included, we can calculate the induced dipole moments of the cavity plasmon and matter excitation:

These expressions are consistent with an alternative classical model that describes the nanocavity and the molecular emitter as dipolar polarizable objects (Section S1 of the Supplementary Material), supporting the validity of the general approach presented here. In the absence of losses [49], the induced dipole moments d _cav and d _mat diverge at the eigenfrequencies ω _±,SpC of the SpC model (Eq. (10)). To avoid these divergences, we add an imaginary part to the bare cavity and matter frequencies in this section. These imaginary parts are related to the decay rates of the cavity, κ, and of the matter excitation, γ, as I m ( ω cav ) = − κ 2 and I m ( ω mat ) = − γ 2 , respectively.

As an example, we consider a metallic spherical nanoparticle of radius R _cav = 5 nm and with a cavity mode of frequency ω _cav = 3 eV. We consider the same molecular emitter of Section 3.1, with a strong transition dipole moment of magnitude μ _mat = 15 Debye. As indicated by Eq. (28), the coupling strength of the system can be adjusted based on the position and orientation of the molecular emitter. We choose that the dipolar molecular transition is polarized perpendicularly to the surface of the nanoparticle and parallel to the amplitude of the incident field E _inc, to maximize the coupling strength (as a consequence d _cav, d _mat, E ̂ cav ‖ ( r mat ) , and n _rrel are all oriented in the same direction in, e.g., Eqs. (26) and (28)). With this choice, and placing the molecular emitter at 1 nm from the surface of the nanoparticle, we obtain a coupling strength g _SpC ≈ 300 meV = 0.1 ω _cav and thus reach the limit of the ultrastrong coupling regime. This large value of g _SpC is possible due to the small size of the nanoparticle (large field confinement) and to the strong transition dipole moment considered for the molecular emitter, which lies slightly beyond the values of μ _mat = 3 − 5 Debyes corresponding to typical molecules used in combination with plasmonic systems. Even larger field confinement may be possible in nonspherical experimental nanostructure configurations that exploit very narrow gaps [18], [77]. To ensure that the system is also in the strong coupling regime when considering lower values of g _SpC below, we choose γ = 10 meV and a damping rate of the plasmonic cavity κ = 20 meV that is small compared to those of usual plasmonic metals.

The induced dipole moments obtained from Eq. (29) can be used, for example, to calculate the near-field distribution under excitation at frequency ω. The total electric field is the sum of the cavity E cav ‖ and molecular or matter contribution E mat ‖ . Under the quasistatic approximation, with d cav ( ω ) = f cav x cav ( ω ) and d mat ( ω ) = f mat x mat ( ω ) , we obtain that the fields at position r outside the metallic nanoparticle, |r − r _cav| > R _cav, depend on the amplitude of the harmonic oscillators as

From this expression, the fields at the frequency of each hybrid mode are calculated by replacing into Eq. (30) the oscillation amplitudes in Eq. (29) induced at the mode frequencies ω _±,SpC.

The electric fields along the x-axis associated with the upper and lower mode frequencies are plotted in the top and bottom panels of Figure 4(b) (blue lines), respectively. These fields are real and polarized along the x direction. We further show the decomposition of the fields into the contribution of the cavity (black dots) and the molecular emitter (black dashed line) as given by the first and second terms on the right-hand side of Eq. (30), respectively. It can be appreciated from Figure 4(b) that, for example, when the upper hybrid mode is excited, the dipoles associated with the cavity and the molecular emitter are oriented in the same direction (same sign). In contrast, for the lower mode, the dipoles point toward the opposite direction.

The near field plotted in Figure 4(b) is useful for analyzing the behavior of the hybrid modes. Still, it is difficult to measure, and most experiments focus on far-field spectral information, such as in the scattering cross section spectral σ _sca. Due to the small emitter-nanocavity distance, we neglect retardation effects so that σ _sca is related to the total induced dipole moment of the system as [78]

We show in Figure 4(c) the scattering cross section for the same nanoparticle-molecular emitter system in the outset of the ultrastrong coupling regime (g _SpC = 0.1 ω _cav). Since the oscillator strength of the cavity is much larger than that of the emitter (f _cav ≫ f _mat), the spectrum is entirely dominated by the cavity contribution, obtained from Eq. (29a) (however, in other systems, where both oscillator strengths are similar, f _cav ≈ f _mat, it is crucial to consider both contributions in Eq. (31)). The scattering cross section spectra are shown for two different detunings between the nanocavity and the molecular emitter. At zero detuning (ω _cav = ω _mat = 3 eV, solid lines in Figure 4(c)) the upper hybrid mode has a (moderately) larger cross section than the lower hybrid mode, mostly due to the ω ⁴ factor in Eq. (31). However, when the molecular excitation is blue detuned with respect to the cavity (ω _cav = 3 eV and ω _mat = 3.2 eV, dashed line), the strength of the peak in the cross section spectra associated with the lower hybrid mode increases and the upper hybrid mode becomes weaker. This behavior occurs because, for this detuning, the lower hybrid mode acquires a larger contribution of the cavity resonance that dominates the scattering spectra, while the predominantly emitter-like behavior of the upper mode results in a smaller cross section due to f _mat ≪ f _cav.

To assess the importance of using the classical SpC model to describe this system, we compare the results of the scattering cross section spectra calculated with this model with those obtained using the MoC model. For this purpose, it is necessary to obtain the expressions of the scattering cross section for the latter model under external illumination. By introducing forcing terms in the equations of motion of the MoC model (Eq. (11)) to account for the external field, we obtain the corresponding oscillation amplitudes