Ultrastrong coupling effects in molecular cavity QED

2.1 The model

We consider one molecule inside a cavity (see Figure 1A for a schematic diagram). This system can be described as the single-molecule version of the Holstein-Tavis-Cummings Hamiltonian, which has been analyzed in depth in the past for collections of molecules [16], [17], [18], [19] (throughout this article, we denote this single-molecule case as the Holstein-Jaynes-Cummings (HJC) case). However, this Hamiltonian is obtained after neglecting the counter-rotating (CR) terms that arise when quantizing the dipole Hamiltonian [1], [13]. This is correct in the usual case where the molecule-photon interaction is weak, but as the coupling increases when the photon modal volume decreases, the CR terms may be relevant when considering ultrasmall plasmonic cavities. We thus retain the CR terms and propose the Holstein-Quantum-Rabi (HQR) Hamiltonian:

Figure 1:

Energy levels in molecular cavity QED.

(A) Schematic representation of the energy scales in the problem of a molecule in a cavity. (B) Lower eigen energies as a function of the Huang-Rhys factor. Red lines correspond to the vibrational ladder of the photon-exciton ground-state, while black lines correspond to the 1 excitation light-matter (polariton) sector. The parameters used are ω_c=Δ=1, ω_v=0.075, and g=0.05. (C) Zoom of panel (B) in the polariton sector. Solid lines were obtained with the HQR model, data points marked with crosses were obtained with the HJC Hamiltonian (neglecting the CR terms), and the discontinuous lines are obtained with the effective QR model described by Eq. (3).

where the operators a⁺, σ⁺, and b⁺ create one cavity photon (with energy ω_c), one exciton in the molecule (with energy Δ), and one molecular vibrational quantum (with energy ω_v), respectively, while their adjoint operators (a, σ⁻, and b) annihilate the corresponding excitations. The operators a and b are bosonic, while the σ’s are Pauli matrices operating in the molecular ground state-exciton two-level manifold.

The Holstein exciton-phonon interaction takes into account that the molecule vibrates differently in the ground and excited states, and it is characterized by the Huang-Rhys factor λ². The coefficient g sets the exciton-photon interaction strength and depends on both the molecular transition dipole moment and the photon modal volume. When g is small enough compared to both ω_c and Δ, the CR term H_CR=g(σ⁺a⁺+σ⁻a⁻) can be safely neglected, arriving at the HJC model. On the contrary, for large enough g, the CR term is relevant to the dynamics of the system (situation termed as “ultrastrong coupling regime” or USC). In CQED, the rule of thumb is that reaching the USC requires g≳0.1 Δ [12], [13], [20]; here we will show that this condition is modified in molecular CQED.

Notice that the diamagnetic term (A²~(a+a⁺)²)) has not been included in the Hamiltonian (1), as we assume that its effect has already been taken into account in the values of ω_c and g [17].

Hamiltonian (1) is expressed in the base of vibrational levels of the electronic ground state {n}. It is possible to go into a representation where vibrations in the electronic ground state are expressed in the base {n} while the vibrations in the exciton sector are expressed in their own eigenfunctions: the displaced oscillators {n˜}. In this case, the base vectors are {|↓, i, n〉,|↑, i, n˜〉}, where ↓, ↑ refers to the electronic degree of freedom and i is the number of photons. This basis change is implemented by a polaron transformation, which takes b→b–λσ⁺σ⁻ via the unitary transformation H→UP†HUP, with U_P=exp(–λσ⁺σ⁻(b⁺–b)). After standard manipulations, we obtain a Hamiltonian that is exactly equivalent to Hamiltonian (1):

In this representation, the vibrations “dress” the exciton-photon coupling through the Frank-Condon factors 〈n|D(λ)|m˜ 〉  ≡  〈n|eλ(b+−b)|m˜〉 [21].

Although the motivation behind the presentation of the HQR Hamiltonian is the application to molecular CQED, note that it could more generally apply to cases where a 2LS is coupled both to a cavity photon and to another bosonic degree of freedom, a situation that may occur in circuit QED [22]. With this in mind, in what follows we present results over a wide range of Huang-Rhys factors (which in molecules typically range from 0 to ~2 [23], [24]).

2.2 The spectrum

In this article we focus on the modification of the dynamics of a 2LS in a cavity due to the presence of vibrational modes. We thus consider the cavity to be in resonance with the zero-phonon excitonic transition (ω_c=Δ, which is taken as the energy unit). Figure 1B renders the numerically computed spectrum for the case g=0.05 (the chosen values for ω_c and g representative for vibrations in organic molecules [25] and ultrasmall plasmonic cavities [10], [11], respectively). Results for both HQR and HJC models are shown, demonstrating that the CR terms have only a minimal impact in the eigen energies for that value of g. The spectrum shows a series of vibrational modes associated with the exciton-photon ground state (with energies virtually independent of λ) and another set associated with the vibrational dressed polaritonic states. At λ=0, when |m〉= |m˜〉, this set comprises m–phonon replicas of the polaritonic states |Pm±〉= |±〉 ⊗|m〉=2−1/2{|↓, 1〉 ± |↑, 0〉}⊗|m〉 (in the Rabi model, polaritons may have a more complex structure, but for the considered values of g, this “Jaynes-Cummings” expression is an excellent approximation). For finite λ, these vibrational states couple and the eigenstates do not have a well-defined number of phonons. Notably, as shown in Figure 1C, some eigenstates trend toward degeneracy at large λ (e.g. the two lowest polaritonic states in Figure 1C, arising from |P₀₋〉 and |P₁₋〉). In order to understand this feature, which as we will show has consequences on the Rabi oscillation, we assume that ω_c=Δ, and (i) neglect the CR terms in the HQR Hamiltonian in Eq. (2), (ii) project the Hamiltonian on the polariton basis on the Jaynes-Cummings light-matter interaction |±〉, and (iii) perform a polaron transformation with U_P=exp(–λ(b⁺–b)/2). After this, we obtain that the polaritonic sector is governed by an effective Quantum-Rabi (QR) Hamiltonian

where the σ_P operators work in the two-level subspace spanned by the exciton-photon polaritons, Δ^=2g,g^=λωv/2, and ε= ω_c+ω_vλ²/4–g is just an energy shift of all eigen energies. We emphasize that Eq. (3) has been derived assuming the resonant condition ω_c=Δ, as we wanted to stress that in the Jaynes-Cummings polariton sector the system behaves according to the QR Hamiltonian (or, in other words, it could be used as a “quantum simulator” for this Hamiltonian). The validity of this effective Hamiltonian can be seen in Figure 1C: the approximate spectrum is very close to those obtained with both HQR and HJC models. Thus, the dynamics in the polaritonic subspace mimic that of CQED, with the polaritons playing the part of the 2LS and the phonons the part of the cavity photons! The effective description also unveils another feature that was not evident in the HJC one: at photon-exciton resonance, the number of excitations in Hamiltonian (3) is not conserved, but they uncouple in sectors with odd and even numbers of excitations. The only exception is near degeneracies, for instance, the crossing between P₀₊ and P₂₋ at λ≈0.74. While both HJC and effective Rabi models predict a crossing of levels (notice that those levels have different parity), the full HQR model produces an anti-crossing with a small gap (not shown). This gap originates from the CR terms, which in the effective Rabi model induce additional parity-breaking terms. This ultrastrong coupling effect is small in the spectrum but, as will be shown later, affects the dynamics.

Thus, away from these degeneracy points, the dynamics in the polariton subspace is governed by an effective QR Hamiltonian, which can even reach the deep-ultrastrong coupling regime g^≳Δ^ (a parameter regime that has proven difficult to access experimentally), even when g is small enough for the CR terms in the original Hamiltonian to be negligible. The analogy also explains why the evolution of the energy levels with λ strongly resembles that in the QR model with coupling strength g. It must be stressed nonetheless that the effective QR Hamiltonian applies to the dynamics of excited states but, additionally, the system presents a manifold of energy states associated with the molecular ground state. Obviously, this manifold should be taken into account in the presence of decay channels for the polaritons.

2.3 Dynamics: ultrastrong effects

We consider the situation where one photon enters the system, in resonance with the zero-phonon exciton, and study the subsequent dynamics. In this work we assume that the decay rates are small enough to be safely neglected in the timescales we examine; the effect of losses will be analyzed in a subsequent publication. It can be anticipated, nevertheless, that the high losses present in today’s room-temperature plasmonic cavities would have to be drastically reduced in order to observe any ultrastrong coupling effects. This can perhaps be achieved by lowering the temperature, considering metallodielectric cavities with high dielectric index, or using quantum circuits [22].

Figure 2 renders, for different values of λ and g, the time evolution of the photon number P(t)≡〈a⁺a〉(t) (the exciton number E(t)=〈σ⁺σ⁻〉(t) is complementary to P(t), as their sum is 1). Each panel shows the comparison between the calculations using the full HQR and the HJC models. In the λ=0 case, the vibrational degrees of freedom decouple and, for the considered initial condition, the system is always in the zero-vibration state. Thus, the molecule behaves as a 2LS and the system maps into traditional CQED, where ultrastrong coupling effects are negligible for g=0.05 (Figure 2B) and very small even for g=0.2 (Figure 2C). As shown in the figure, the frequency of the Rabi oscillations Ω_R strongly decreases with λ. This occurs because the oscillations mainly involve the two lowest polaritonic states, whose energy decreases with λ (as shown in Figure 1C). But, notably, the influence of the CR terms on the dynamics is strongly enhanced for larger values of λ, as shown by the incompleteness of Rabi oscillations in the lower panels of Figure 2. This is highlighted in Figure 3, which renders the comparison between the time-averaged values for P(t), E(t), and V(t)≡〈b⁺b〉(t) when the CR terms have been either considered or neglected, as a function of λ. In the last (“Jaynes-Cummings”) case, P_JC=E_JC=1/2 for all λ. The presence of CR terms change the occupations in two ways. First, they “dress” the bare energies of the states (“Bloch-Siegert” effect). This can be taken into account considering H_CR as a perturbation to the HJC Hamiltonian. Within second order, the Bloch-Siegert corrections to the bare eigen energies are

Figure 2:

Time evolution of the average number of photons, P(t), for different values of λ and g.

The photon and the zero-phonon exciton are considered to be in resonance, ω_c=Δ, and ω_v=0.075. Each panel shows the comparison between the full HQR model (red curves) and the one without the CR terms (Holstein-Jaynes-Cummings model, black curves). Top panels (A–C) are for λ=0, and bottom ones (D–F) are for λ=2.5. The case g=0.05 is rendered for both values of λ, while the other panels are representative of the values g needed for ultrastrong coupling effects to appear for each λ.

Figure 3:

Time-averaged number of excitons (E), photons (P), and molecular vibrations (V) as a function of λ.

Continuous lines are computed using the Holstein-Quantum-Rabi model, while dashed lines were computed by neglecting the counter-rotating terms (HJC model). The inset shows the Bloch-Siegert correction to the photon energy, ΔE_BS(|↓, 1, 0〉), and the effective value of g (i.e. half their energy difference between the two lowest polaritons). The parameters used are ω_c=Δ=1, ω_v=0.075, and g=0.05.

where, in the approximation to ΔE_BS(|↓, 1, 0〉), we have used the following properties of the Frank-Condon factors: (i) 〈n˜|0〉 is peaked at n˜=λ2 and (ii) ∑n˜|〈n˜|m〉|2=1. The important point is that the Bloch-Siegert corrections dress the exciton and photon states differently. This “ultrastrong” mechanism can be incorporated into an effective Jaynes-Cummings model by “renormalizing” the photon frequency ω_c→ ω_c+ΔE_BS(|↓, 1, 0〉), which clearly affects whether the photon is in resonance with the exciton or not. This shift, combined with the strong renormalization of the effective coupling that occurs at large λ, brings the exciton and photon out of resonance. This is illustrated in the inset to Figure 3. Assuming the bare resonant condition ω_c=Δ, the system essentially remains at resonance for values of λ such that |ΔE_BS(|↓, 1, 0〉)|<g_eff, thus developing complete Rabi oscillations. But these cease to happen when |ΔE_BS(|↓, 1, 0〉)| and g_eff are comparable (for λ≳2.5 in the inset to Figure 3, computed for g=0.05). As g_eff~g, while ΔE_BS~g², the BS corrections are more relevant for larger exciton-photon interactions but, admittedly, this effect plays a role only for large values of λ. It is worth noting that the CR parity-breaking terms we mentioned when discussing the effective Hamiltonian (3) cancel when the “renormalized” cavity is in resonance with the 2LS, i.e. when ω_c+ΔE_BS=Δ. In that case, Hamiltonian (3) is virtually exact in the polariton subspace, for all values of g, ω_v, and λ.

The second way in which the CR terms modify the occupations occurs at smaller values of λ. It works via mixing states which would be orthogonal within the HJC Hamiltonian but anti-cross when the CR terms are considered (which, as mentioned before and shown in Figure 1, occurs for the states P₀₊ and P₂₋ at λ≈0.74). This mixing allows P₀₊ to couple to P₂₋, thus enhancing the average number of phonons P (see the peak in P in Figure 3, at λ≈0.74).