Strong coupling dynamics of driven quantum systems with permanent dipoles

Burgess, Adam; Florescu, Marian; Rouse, Dominic M.

doi:10.1116/5.0157714

Many optically active systems possess spatially asymmetric electron orbitals. These generate permanent dipole moments, which can be stronger than the corresponding transition dipole moments, significantly affecting the system dynamics and creating polarized Fock states of light. We derive a master equation for these systems with an externally applied driving field by employing an optical polaron transformation that captures the photon mode polarization induced by the permanent dipoles. This provides an intuitive framework to explore their influence on the system dynamics and emission spectrum. We find that permanent dipoles introduce multiple-photon processes and a photon sideband, which causes substantial modifications to single-photon transition dipole processes. In the presence of an external drive, permanent dipoles lead to an additional process that we show can be exploited to control the decoherence and transition rates. We derive the emission spectrum of the system, highlighting experimentally detectable signatures of optical polarons, and measurements that can identify the parameters in the system Hamiltonian, the magnitude of the differences in the permanent dipoles, and the steady-state populations of the system.

I. INTRODUCTION

In general, the interaction of atomic systems with light through transition dipoles is well understood, and the optical master equation describing exciton creation and annihilation through photon emission and absorption, respectively, is derived in many introductory texts dedicated to open quantum systems.^1–4 Atomic systems have highly symmetric electron orbitals and so possess negligible permanent dipoles. However, many physical systems do not share this property and can possess permanent dipoles stronger than their transition dipole moments. Such systems include molecules with parity mixing of the molecular state,^5–12 quantum dots with asymmetric confining potentials,^13–23 nanorods with non-centrosymmetric crystallographic lattices,^24–26 and superconducting circuits.²⁷

Permanent dipoles introduce additional pure dephasing interactions into the Hamiltonian. The non-additivity of the pure dephasing and transition dipole interactions yields unique physical effects, including modifications to decoherence,^28,29 steady-state coherence,^28,30,31 laser-driven population inversion,³² multiphoton conversion,^33,34 entanglement generation,^35,36 second-harmonic generation,^37,38 and bathochromic shifts.^39,40 Understanding the role of strong permanent dipoles and their influence on control schemes, such as external driving fields, is highly relevant for the design of new quantum technologies, as decoherence induced by the local electromagnetic field limits the capabilities of quantum computation implementations. Studying the impacts of permanent dipoles may provide additional channels for control over such quantum systems, as well as advancements in quantum chemistry that may uncover novel biochemical processes.

Previous studies have either neglected the transition dipole moments assuming only a pure dephasing interaction,¹² have considered single-mode fields,^{29,34,41–43} or have treated the permanent dipoles in a perturbative manner.^{12,28,29,35,41–43} Furthermore, most studies neglect an additional identity interaction also induced by the permanent dipoles, which, as we show here, modifies the initial state of the environment and can have a significant impact on the system dynamics. Consequently, such treatments do not capture the role of strong permanent dipoles in asymmetric systems under illumination by multimode fields, such as common thermal fields.

In this paper, we utilize a polaron transformation to derive a master equation for driven quantum dipole systems, possessing strong permanent dipoles, interacting with a thermal electromagnetic field. Importantly, we make no assumptions about the dipole matrix beyond perturbative transition dipoles. The polaron transformation is a unitary, state-dependent displacement transformation widely used when dealing with strong, pure dephasing interactions. In the polaron frame, pure dephasing interactions are absorbed into the definition of the basis and treated to all orders in the coupling strength. The basis describes an optical polaron quasiparticle, which is a hybridization of the matter excitations and the displaced harmonic oscillator states of the multimode photonic field. These photonic states, called polarized Fock states, correspond to a non-zero (polarized) vector potential field and were first introduced to explain how permanent dipoles can generate multiphoton conversion in a single-mode cavity³⁴ and are useful in polaritonic chemistry.^44,45

The optical polaron formalism provides us with an essential intuition for the very complex phenomena introduced by permanent dipoles, allowing us to unpick the role of permanent dipoles in the dynamics. We find that the formation of optical polarons results in unique physical phenomena such as modifications to single-photon transition dipole processes, entirely new multiple photon processes, the appearance of photonic sidebands, and a novel interplay with external driving that allows for control over the dynamics. By careful choice of a driving field, we can strongly modulate the decoherence rate of the system induced by the local environment by completely suppressing the transition and dephasing rates, effectively decoupling the system from the field. We also derive the emission spectrum for the system, highlighting experimentally detectable signatures of the optical polarons.

This paper is organized as follows. In Sec. II, we introduce the Hamiltonian and transform it into the polaron frame, and in Sec. III, we discuss realistic parameters. In Sec. IV, we derive the polaron frame master equation (PFME) and discuss the new physical processes using the analytical expressions. Following this, in Sec. V, we compare the PFME to perturbative dynamics and to numerically exact dynamics using the time-evolving matrix product operators (TEMPO)^46–54 algorithm, through the open source code OQuPy.⁵⁵ In Sec. VI, we derive the emission spectrum, and in Sec. VII, we briefly discuss the role of the initial state and identity type interactions. Finally, in Sec. VIII, we present concluding remarks.

II. THE MODEL AND POLARON FRAME

We consider a driven asymmetric emitter with a single quantized dipole coupled to a long wavelength multimode field. After truncating the material subsystem to its two lowest energy levels, the fundamental multipolar-gauge Hamiltonian is

H = \frac{ϵ}{2} σ_{z} + V σ_{+} + V^{*} σ_{-} + \sum_{k} ν_{k} a_{k}^{†} a_{k} + d \cdot Π + E_{dip},

(1)

where ϵ is the transition energy of the emitter and

V = | V | e^{i ϑ_{V}}

is the complex valued amplitude of the external drive.^56–58

Π = i \sum_{k} e_{k} f_{k} (a_{k}^{†} - a_{k})

is the electric displacement field, where

k = {k, λ}

is a four-vector representing both the wavevector k and polarization state λ of the mode with polarization vector

e_{k}

⁠, energy

ν_{k}

⁠, and coupling strength

f_{k} = \sqrt{ν_{k} / (2 V)}

⁠, where

V

is the field volume. a_k and

a_{k}^{†}

are the field mode annihilation and creation operators, and the dipolar self-energy term is

E_{dip} = \sum_{k} \frac{f_{k}^{2}}{ν_{k}} {(d \cdot e_{k})}^{2} .

(2)

In the truncated system Hilbert space, the dipole operator is

d = (\begin{matrix} d_{e e} & d_{e g} \\ d_{g e} & d_{g g} \end{matrix}) = d_{Δ} σ_{z} + d_{D} I + d_{μ} σ_{+} + d_{μ}^{*} σ_{-},

(3)

where the Pauli operators are

σ_{z} = | e ⟩ ⟨ e | - | g ⟩ ⟨ g |, σ_{+} = | e ⟩ ⟨ g |, σ_{-} = | g ⟩ ⟨ e |, I = | g ⟩ ⟨ g | + | e ⟩ ⟨ e |

⁠, and

d_{i j} = ⟨ i | d | j ⟩

for

i, j \in {e, g}

⁠. We have defined the following combinations of dipole matrix elements:

d_{Δ} = \frac{d_{e e} - d_{g g}}{2}, d_{D} = \frac{d_{e e} + d_{g g}}{2}, d_{μ} = d_{e g},

(4)

which play an essential role in our analysis. The

d_{p}

vectors where

p \in {μ, Δ, D}

are not guaranteed to be co-linear and Hermiticity of H require that

d_{p} \in ℝ^{3}

for

p \in {Δ, D}

and

d_{μ} \in ℂ^{3}

⁠.

At this point, it is often assumed that either $| d_{Δ} | \approx | d_{D} | \approx 0$ ⁠, which leads to the standard optical master equation (SOME),^1–4 or that $| d_{D} | \approx | d_{μ} | \approx 0$ leading to a pure dephasing interaction.¹² In both of those limits, the dipolar self-energy term is proportional to the identity. In this study, we make no assumptions about the size of the permanent dipoles.

Substituting Eq. into Eq. and absorbing the drive phase into a basis

| e' ⟩ = e^{i ϑ_{V} / 2} | e ⟩

and

| g' ⟩ = e^{- i ϑ_{V / 2}} | g ⟩

⁠, we find

\begin{matrix} H = \frac{ϵ}{2} σ_{z}^{'} + | V | σ_{x}^{'} + \sum_{k} ν_{k} a_{k}^{†} a_{k} + E_{dip} + π_{D D} I' \\ + π_{Δ Δ} σ_{z}^{'} + π_{μ \bar{μ}} e^{- i ϑ_{V}} σ_{+}^{'} + π_{\bar{μ} μ} e^{i ϑ_{V}} σ_{-}^{'}, \end{matrix}

(5)

where primed operators are in the

{| e' ⟩, | g' ⟩}

basis and

π_{p q} = \sum_{k} (p_{k} a_{k}^{†} + q_{k}^{*} a_{k}),

(6)

p_{k} = i f_{k} (d_{p} \cdot e_{k}),

(7)

with

p, q \in {μ, \bar{μ}, Δ, D}

⁠, and we denote

d_{\bar{μ}} = d_{μ}^{*}

⁠. We refer to Eq. as the lab frame Hamiltonian, which is equivalent to Eq. .

The photon-only part of Eq. can be diagonalized by the displacement transformation

H_{d} = B (D / ν) H B (- D / ν)

⁠, where subscript d denotes the displaced frame, and displacement operators are given by

B (α) = e^{\sum_{k} (α_{k} a_{k}^{†} - α_{k}^{*} a_{k})},

(8)

and

B {(α)}^{†} = B (- α)

⁠. These act on photon operators by

B (\pm α) a_{k} B (\mp α) = a_{k} \mp α_{k},

(9)

and we further analyze the displacement operators in Appendix A. Ignoring constant terms, the resulting Hamiltonian is

\begin{matrix} H_{d} = \frac{ϵ}{2} σ_{z}^{'} + | V | σ_{x}^{'} + \sum_{k} ν_{k} a_{k}^{†} a_{k} + π_{Δ Δ} σ_{z}^{'} + π_{μ \bar{μ}} e^{- i ϑ_{V}} σ_{+}^{'} + π_{\bar{μ} μ} e^{i ϑ_{V}} σ_{-}^{'} . \end{matrix}

(10)

Notice that the dipolar self-energy term has canceled exactly with the terms that result from displacing the light–matter interactions, such that Eq. is independent of D_k. This is a remarkable result, as it tells us that the standard perturbative master equation would be considerably inaccurate for most values of

d_{D}

as this term would have a non-negligible contribution to the dynamics.

Throughout this paper, we make the standard assumption that the transition dipole is small enough to permit an accurate second-order perturbative expansion in its magnitude, $| d_{μ} |$ ⁠. Even for perturbative transition dipole moments, the master equation derived using H_d in Eq. (10), referred to as the displaced frame master equation (DFME), will become inaccurate if $| d_{Δ} |$ is large.

To overcome this challenge, we make a polaron transformation prior to deriving the Redfield master equation. In this frame, the polarizing effect of the pure dephasing interaction on the field is absorbed into the definition of a new basis, called the polarized Fock states.^34,44 The polarization direction of the Fock state is dependent on the matter state, and it hybridizes with the excitation to create an optical polaron quasiparticle. In Fig. 1, we illustrate the various frames introduced and the optical polaron concept. A Redfield master equation derived in the new basis, termed the polaron frame master equation (PFME), will be robust to all magnitudes of the permanent dipoles and will recover the DFME in the limit of $Δ_{k} \to 0$ ⁠.

Fig. 1.

View large Download slide

Illustration of the model and the optical polaron concept. (a)–(c) depicts the molecular energies and interactions in the (a) lab frame, (b) displaced frame, and (c) polaron frame. In the lab frame, the dipolar self-energy term $E_{dip}$ causes the renormalizations $ϵ \to ϵ_{dip}$ and $V \to V_{dip}$ ⁠, which cancel out in the displaced and polaron frames. (d) is an illustration of an optical polaron: a quasiparticle formed of the matter excitation and the polarized Fock states of the displaced photon modes.

The polaron frame Hamiltonian is

H_{p} = U H_{d} U^{†}

⁠, where

U = B (\frac{Δ}{ν}) | e' ⟩ ⟨ e' | + B (- \frac{Δ}{ν}) | g' ⟩ ⟨ g' | .

(11)

Using Eq. , and ignoring constant terms, we obtain

H_{p} = \frac{ϵ}{2} σ_{z}^{'} + κ | V | σ_{x}^{'} + \sum_{k} ν_{k} a_{k}^{†} a_{k} + C σ_{+}^{'} + C^{†} σ_{-}^{'},

(12)

where the coupling operator is

C = C - ⟨ C ⟩

with

C = B (\frac{Δ}{ν}) π_{μ \bar{μ}} B (\frac{Δ}{ν}) e^{- i ϑ_{V}} + | V | B (2 \frac{Δ}{ν}) .

(13)

We denote

⟨ \cdot ⟩ = {Tr}_{E} (\cdot ρ_{E})

⁠, where

ρ_{E} = \frac{1}{Z_{E}} e^{- β \sum_{k} ν_{k} a_{k}^{†} a_{k}},

(14)

and

Z_{E} = Tr [exp (- β \sum_{k} ν_{k} a_{k}^{†} a_{k})]

with

β = 1 / (k_{B} T)

the inverse temperature. In Appendix B, we prove that

⟨ C ⟩ = κ | V |

⁠, where

κ = ⟨ B (\pm 2 Δ / ν) ⟩

⁠. In Eq. , we moved a factor of

⟨ C ⟩ σ_{x}^{'}

from the coupling operator and into the unperturbed part of H_p, which ensures that the perturbation theory yields the Redfield equation.

In the polaron frame Hamiltonian in Eq. (12), we can now interpret $σ_{\pm}^{'}$ as causing transitions between the optical polaron states in Fig. 1(d). Compared to matter excitations, described by the Hamiltonian in the displaced frame in Eq. (10), polarons experience a more complicated interaction, albeit one without a pure dephasing term. The polaron frame Hamiltonian is diagonal when the transition dipoles and driving vanish indicating that permanent dipoles generate trivial dynamics if considered in isolation, described by the independent boson model.

Finally, the unperturbed polaron Hamiltonian can be diagonalized using a unitary rotation,

(ϵ / 2) σ_{z}^{'} + κ | V | σ_{x}^{'} = (η / 2) τ_{z}

⁠, where

η = \sqrt{ϵ^{2} + 4 κ^{2} {| V |}^{2}},

(15)

τ_{z} = | + ⟩ ⟨ + | - | - ⟩ ⟨ - | .

(16)

The eigenbasis relates to the original basis by

(\begin{matrix} | e' ⟩ \\ | g' ⟩ \end{matrix}) = (\begin{matrix} cos (\frac{φ}{2}) & - sin (\frac{φ}{2}) \\ sin (\frac{φ}{2}) & cos (\frac{φ}{2}) \end{matrix}) (\begin{matrix} | + ⟩ \\ | - ⟩ \end{matrix})

(17)

with

cos (φ) = ϵ / η

and

sin (φ) = 2 κ | V | / η

⁠. In the eigenbasis, the polaron frame Hamiltonian is

H_{p} = \frac{η}{2} τ_{z} + \sum_{k} ν_{k} a_{k}^{†} a_{k} + \sum_{α \in {z, +, -}} g_{α} τ_{α},

(18)

where we have defined

τ_{+} = | + ⟩ ⟨ - |, τ_{-} = | - ⟩ ⟨ + |

⁠, and the coupling operators

(19)

g_{z} = \frac{1}{2} sin (φ) [C + C^{†}],

(19a)

g_{+} = [{cos}^{2} (\frac{φ}{2}) C - {sin}^{2} (\frac{φ}{2}) C^{†}],

(19b)

and

g_{-} = g_{+}^{†}

⁠.

To derive the master equation, we will take the continuum limit of the photon modes, in which summations over an arbitrary function

F (ν)

transform according to

\begin{matrix} \sum_{k} f_{k}^{2} F (ν_{k}) (d_{p} \cdot e_{k}) (d_{q} \cdot e_{k}) \to h_{p q} \int_{0}^{\infty} d ν J (ν) F (ν), \end{matrix}

(20)

where

J (ν)

is the spectral density,

p, q \in {μ, \bar{μ}, Δ}

⁠, and

h_{p q} = \int_{Ω_{k}} d Ω_{k} \sum_{λ} ({\tilde{d}}_{p} \cdot e_{k}) ({\tilde{d}}_{q} \cdot e_{k}),

(21)

where

\int_{Ω_{k}} d Ω_{k} = \int_{0}^{π} d θ_{k} sin (θ_{k}) \int_{0}^{2 π} d ϕ_{k}

⁠.² We have defined dimensionless dipole vectors

{\tilde{d}}_{p} = d_{p} / d_{ref}

measured against an arbitrary reference value which for convenience we take to be

d_{ref} = 1 D \approx 9 \times 10^{- 6} {eV}^{- 1}

⁠. The purpose of introducing

d_{ref}

is to ensure that the spectral density has dimensions of energy. For free space and unpolarized light,

(22)

h_{μ Δ} = Ω_{μ Δ} e^{i ϑ_{μ}} cos (θ_{μ Δ}),

(22a)

h_{μ \bar{μ}} = Ω_{μ μ},

(22b)

h_{μ μ} = Ω_{μ μ} e^{2 i ϑ_{μ}},

(22c)

where

Ω_{p q} = Ω_{q p} = (8 π / 3) | {\tilde{d}}_{p} | | {\tilde{d}}_{q} |, ϑ_{μ}

is the complex phase of

d_{μ}, θ_{μ Δ}

is the angle between

d_{Δ}

and

d_{μ}

⁠, and one can obtain, e.g.,

h_{\bar{μ} Δ} = h_{μ Δ}^{*}

by suitable complex conjugation. Notice that the dynamics only depend on the relative phase

ϑ_{μ V} = ϑ_{μ} - ϑ_{V}

⁠. Also, in the continuum limit

κ = exp [- ϕ (0) / 2]

where the photon propagator is

\begin{matrix} ϕ (s) = 4 Ω_{Δ Δ} \int_{0}^{\infty} d ν \frac{J (ν)}{ν^{2}} [coth (\frac{β ν}{2}) cos (ν s) \\ - i sin (ν s)] . \end{matrix}

(23)

III. REALISTIC PARAMETERS

Emitters in free space have weaker permanent dipole interactions than, for example, those in solvents or cavities. However, solvent and cavity models introduce substantial complexity to polaron type master equations due to infrared divergences and non-Markovianity originating from their Drude–Lorentz type spectral densities.¹² Conversely, the multipolar-gauge free-space spectral density

J (ν) = S \frac{ν^{3}}{ν_{c}^{2}} e^{- \frac{ν}{ν_{c}}} Θ (ν),

(24)

where

Θ (ν)

is the Heaviside step function and does not suffer from these problems provided that the high-frequency cutoff ν_c is large. For simplicity, in this paper, we use a free space model with permanent dipole strengths inspired by realistic molecules embedded in solvents and cavities.

As we derive in Appendix C, $S = \int_{0}^{\infty} d ν J (ν) / ν^{2} = s_{0} {(ν_{c} d_{ref})}^{2} / [2 {(2 π)}^{3}]$ is a dimensionless Huang–Rhys parameter where s₀ is a Purcell enhancement factor. $s_{0} = 1$ corresponds to free space, and $s_{0} > 1$ to restricted environmental geometries and solvents. Multiplying the spectral density by s₀ is accurate within the weak field regime.⁵⁹ The cutoff frequency in Eq. (24) is required to enforce the electric dipole approximation implicit in Eq. (1), such that $ν_{c} \sim 1 / r_{d}$ ⁠, where $r_{d} ≳ 50^{2} a_{0}$ is the size of the dipole, giving $ν_{c} \sim [1, 10] eV$ ⁠.⁵⁷ Unless otherwise stated, we use $ν_{c} = 1 eV$ ⁠.

In all our numerical results—except when comparing to TEMPO which we discuss in that section—we take parameter values using the gamma-globulin protein as an example: $| d_{Δ} | \sim 100 D, | d_{μ} | = 0.1 D$ ⁠, and emitter energy $ϵ = 3.16 eV$ ⁠.^32,60 Taking the lifetime of an individual dipole—i.e., when $| d_{Δ} | = 0 = V$ —to be $1 ns$ ⁠, requires that $s_{0} = 5.705 \times 10^{6}$ for the values introduced so far.⁶¹ These values fix $S = 9.341 \times 10^{- 7}$ and so the field remains weak, thus justifying our introduction of the Purcell enhancement factor. Moreover, these values predict that polaron formation (⁠ $κ ≪ 1$ ⁠) will occur at $| d_{Δ} | \sim 100 D$ ⁠. Finally, we choose the inverse temperature as $β = 2 {eV}^{- 1}$ corresponding to $T \approx 5800 K$ ⁠.

To summarize, the important parameters in this model are the dipole strengths $| d_{μ} |$ and $| d_{Δ} |$ ⁠, their relative angle $θ_{μ Δ}$ ⁠, the eigenenergy η, the driving strength $| V |$ ⁠, the relative complex phase $ϑ_{μ V}$ ⁠, and the cutoff frequency ν_c.

IV. EFFECTS OF STRONG PERMANENT DIPOLES

In this section, we derive the secularized PFME, provide a brief analytical review of it, and analyze the population transfer rates, decoherence rate, and Lamb shift appearing in the master equation, which are given by Fourier transforms of the environment correlation functions (ECFs). This will allow us to analytically explore the role of permanent dipoles.

Following the analytical discussion of the PFME, we will compare the PFME to the DFME, and to the numerically exact TEMPO,^46–54 which both serve as benchmarks. In the numerical approach, we use the full non-secular Redfield master equations derived in Appendix D for both frames. We emphasize that the non-secular PFME depends only on the physical processes discussed in the main text. Moreover, the non-secular terms in the PFME are perturbative in $| V | / ϵ$ and so are negligible when the polaron transformation is valid.

In the eigenbasis, the secularized, time-local, Redfield master equation is

(25)

\frac{\partial ρ_{+ +} (t)}{\partial t} = - γ_{↓} ρ_{+ +} (t) + γ_{↑} ρ_{- -} (t),

(25a)

\frac{\partial ρ_{+ -} (t)}{\partial t} = - [γ_{d} + i \bar{η}] ρ_{+ -} (t),

(25b)

and

(\partial / \partial t) ρ_{- -} (t) = - (\partial / \partial t) ρ_{+ +} (t), (\partial / \partial t) ρ_{- +} (t) = (\partial / \partial t) ρ_{+ -} {(t)}^{†}

⁠, where

ρ_{i j} (t) = ⟨ i | ρ_{S} (t) | j ⟩

for

i, j \in {+, -}

and

ρ_{S} (t) = {Tr}_{E} [ρ (t)]

⁠. Equations describe population transfer from

| + ⟩

to

| - ⟩

at decay rate

γ_{↓}

⁠, vice versa at an excitation rate

γ_{↑}

⁠, decoherence at a rate γ_d, and oscillations in the coherence at frequency

\bar{η}

⁠. The transition and decoherence rates are

γ_{\overset{↓}{↑}} = 2 ℜ [Γ_{\mp \mp} (\pm η)],

(26a)

γ_{d} = \frac{1}{2} [γ_{↑} + γ_{↓}] + 4 ℜ [Γ_{z z} (0)],

(26b)

where

Γ_{α β} (ω) = \int_{0}^{\infty} d s e^{i ω s} ⟨ g_{α}^{†} (s) g_{β} (0) ⟩,

(27)

for

α, β \in {z, +, -}

and

g_{α} (s)

is the interaction picture form of

g_{α}

in Eqs. . Finally, the Lamb shifted eigenenergy is

\bar{η} = η + ℑ [Γ_{- -} (η) - Γ_{+ +} (- η)] .

(28)

Equation includes an assumption that the environment is Markovian, which is accurate for free space fields—described by the spectral density in Eq. —if the high-frequency cutoff is large.¹

In the PFME, we choose the initial state to be $ρ_{p} (0) = | g ⟩ ⟨ g | \otimes ρ_{E}$ ⁠, where ρ_E is given in Eq. (14) and assume that the environment state does not change throughout the dynamics. We discuss the implications of this initial state in Sec. VII.

Evaluating the ECFs in Eq. ,

⟨ g_{α}^{†} (s) g_{β} (0) ⟩

⁠, is a laborious process, and we present the full derivations in Appendix E. Here, we focus the discussion on new physical processes introduced by the presence of the strong permanent dipoles. The ECFs,

⟨ g_{α}^{†} (s) g_{β} (0) ⟩

⁠, depend on linear combinations of

⟨ C^{†} (s) C (0) ⟩, ⟨ C (s) C^{†} (0) ⟩, ⟨ C^{†} (s) C^{†} (0) ⟩

⁠, and

⟨ C (s) C (0) ⟩

⁠, with the relative weighting of each dependent on the eigenbasis angle

φ

⁠. Importantly, as we show in Appendix E, each of these correlation functions describes the same physics. Therefore, for the analytics in the main text, we focus the discussion on

Γ (ω) = \int_{0}^{\infty} d s e^{i ω s} ⟨ C^{†} (s) C (0) ⟩ .

(29)

This function has contributions from four distinct processes

Γ (ω) = Γ_{1} (ω) + Γ_{2} (ω) + Γ_{V, 1} (ω) + Γ_{V, 0} (ω) .

(30)

The number in the subscripts of each term on the right-hand side of Eq. denotes how many photons are involved in the process (in the absence of the sideband), and terms with a subscript V are induced by the driving. The other three Fourier transforms of the

C (s)

two-time correlation functions also depend on these four contributions, except that

⟨ C^{†} (s) C^{†} (0) ⟩

and

⟨ C (s) C (0) ⟩

do not have a

Γ_{V, 1} (ω)

-type contribution. The prefactor of each term is dependent on the relative dipole angle

θ_{μ Δ}

⁠, the driving amplitude

| V |

⁠, the relative phase

ϑ_{μ V}

⁠, and dipole magnitudes Ω_pq, as

(31)

Γ_{1} (ω) \propto Ω_{μ μ},

(31a)

Γ_{2} (ω) \propto Ω_{μ Δ}^{2} {cos}^{2} (θ_{μ Δ}),

(31b)

Γ_{V, 1} (ω) \propto Ω_{μ Δ} | V | cos (ϑ_{μ V}) cos (θ_{μ Δ}),

(31c)

Γ_{V, 0} (ω) \propto {| V |}^{2} .

(31d)

Equations show that when the system is driven,

Γ (ω)

is not an even function of

θ_{μ Δ}

or

ϑ_{μ V}

⁠. We will later show that this can be utilized to control the system, for example, to prevent decoherence and transitions. Moreover, for perpendicular dipoles

cos (θ_{μ Δ}) = 0

⁠, so

Γ (ω)

becomes equivalent to the analogous function in the driven spin boson model where the transition dipole and permanent dipole interactions are with independent baths. Therefore, processes dependent on

θ_{μ Δ}

arise because the interactions are non-commutative.

For additional context, when $| V | ≪ ϵ$ —the parameter regime in which polaron theory is valid—the eigenbasis angle is $φ \approx 0$ and so $g_{z} \approx 0$ and $g_{+} \approx C$ ⁠. In this limit, the secular approximation used in the master equation in Eqs. (25) becomes exact, and $γ_{↑} \approx 2 ℜ [Γ (- η)]$ with $Γ (ω)$ in Eq. (30) and $γ_{↓} \approx \int_{0}^{\infty} d s e^{i η s} ⟨ C (s) C^{†} (0) ⟩$ which, as we show in Appendix E, evaluates to $γ_{↓} \approx 2 ℜ [\bar{Γ} (η)]$ where $\bar{Γ} (ω)$ is equal to Eq. (30) but with $Γ_{V, 1} (ω) \to - Γ_{V, 1} (ω)$ ⁠.

We will now discuss the four contributions to Eq. (30) in turn. Each is derived in Appendix E.

A. One photon processes, Γ₁(ω)

The function describing driving-independent one photon processes is

\begin{matrix} Γ_{1} (ω) = π Ω_{μ μ} \int_{0}^{\infty} d ν [J_{A} (ν) K (ω + ν) + J_{E} (ν) K (ω - ν)], \end{matrix}

(32)

where

K (ε) = \frac{1}{π} \int_{0}^{\infty} d s e^{i ε s} e^{ϕ (s) - ϕ (0)}

(33)

contains the influence of the permanent dipoles within this term, and we have introduced the absorption and emission spectral densities

(34)

J_{A} (ν) = J (ν) N (ν),

(34a)

J_{E} (ν) = J (ν) \tilde{N} (ν),

(34b)

where

N (ν) = 1 / [exp (β ν) - 1]

is the Bose–Einstein distribution and

\tilde{N} (ν) = N (ν) + 1

will be a useful notation throughout this paper.

To understand the physical processes associated with this term, it is convenient to temporarily ignore the factor of

exp [ϕ (s) - ϕ (0)]

in Eq. . The resulting term is the typical function appearing in the standard optical master equation (SOME)¹—our model in the absence of permanent dipoles and driving—which is

\begin{matrix} Γ_{SOME} (ω) = Ω_{μ μ} [π (J_{A} (- ω) + J_{E} (ω)) + i P \int_{0}^{\infty} d ν (\frac{J_{A} (ν)}{ω + ν} + \frac{J_{E} (ν)}{ω - ν})], \end{matrix}

(35)

where

P

denotes the principal value. Clearly, when

ω = - η

⁠, twice the real part of

Γ_{SOME} (ω)

describes excitation-by-absorption at a rate

2 π Ω_{μ μ} J_{A} (η)

⁠, and when

ω = η

⁠, it describes decay-by-emission at a rate

2 π Ω_{μ μ} J_{E} (η)

⁠. The imaginary part of

Γ_{SOME} (ω)

will determine the Lamb shift.

We now return to Eq. (32). In Ref. 62, an analytic solution to Eq. (33) was found by exploiting the fact that $K (ε)$ only depends on $J (ν)$ through its moments $μ_{m} = 4 Ω_{Δ Δ} \int_{0}^{\infty} d ν [J (ν) / ν^{2}] ν^{m}$ for $m = 1, 2, \dots, \infty$ ⁠, and that moments of lower order contribute relatively more. We can then evaluate $K (ε)$ by replacing $J (ν)$ with a truncated spectral density $J' (ν) = \sum_{k = 1}^{N_{*}} | f'_{k} |^{2} δ (ν - ν_{k}^{'})$ ⁠, as long as we choose the coupling strengths ${f_{k}^{'}}$ and energies ${ν_{k}^{'}}$ of the modes such that $J' (ν) / ν^{2}$ has the same lowest moments as μ_m for $m = 1, 2, \dots, 2 N_{*}$ ⁠. The solution converges rapidly for increasing $N_{*}$ ⁠, and a single-mode truncation is often very accurate, with truncation mode parameters: $ν'_{1} \equiv ν_{s} = μ_{2} / μ_{1}$ and $f'_{1} \equiv f_{s} = μ_{2}^{1 / 2} \propto Ω_{Δ Δ}^{1 / 2}$ ⁠. A single-mode truncation captures all of the essential physics described by $K (ω)$ ⁠. Thus, for the analytical analysis in this work, we use a single-mode truncation to evaluate Eq. (33); however, in all of the simulations in this paper, we use converged values.

From Ref. 62, the single-mode solution of Eq. is

K (ε) = \sum_{ℓ = - \infty}^{\infty} A_{ℓ} [δ (ε - ℓ ν_{s}) + \frac{i}{π} \frac{P}{ε - ℓ ν_{s}}],

(36)

where

A_{ℓ} = \sum_{n = | ℓ |}^{\infty'} \sum_{m = \frac{n - ℓ}{2}}^{n} (\begin{matrix} n \\ m \end{matrix}) (\begin{matrix} m \\ m - \frac{n - ℓ}{2} \end{matrix}) W_{n} V_{m},

(37)

and the prime on the first summation indicates that only every other term is included, i.e.,

n = | ℓ |, | ℓ | + 2, \dots,

and

W_{n} = S_{s}^{n} exp [- S_{s}] / n!

is the Franck–Condon factor of the mode in the truncation,

S_{s} = | f_{s} |^{2} / ν_{s}^{2} = μ_{1}^{2} / μ_{2}

is its Huang–Rhys factor, and

V_{m} = N {(ν_{s})}^{m} exp [- 2 S_{s} N (ν_{s})]

⁠.

A_{ℓ}

has the normalization property

\sum_{ℓ = - \infty}^{\infty} A_{ℓ} = 1

⁠, is maximized for

ℓ = Round (S')

⁠, and becomes

A_{ℓ} = δ_{ℓ 0}

if

| d_{Δ} | = 0

⁠, and at zero temperature

V_{m} = δ_{m 0}

⁠, which leads to

A_{ℓ < 0} = 0

and

A_{ℓ \geq 0} = W_{ℓ}

⁠.

Substituting Eq. into Eq. yields

Γ_{1} (ω) = \frac{1}{2} γ_{1} (ω) + i S_{1} (ω),

(38)

where

γ_{1} (ω) = 2 π Ω_{μ μ} \sum_{ℓ = - \infty}^{\infty} A_{ℓ} [J_{A} (ℓ ν_{s} - ω) + J_{E} (ω - ℓ ν_{s})]

(39)

and

\begin{matrix} S_{1} (ω) = i Ω_{μ μ} \sum_{ℓ = - \infty}^{\infty} A_{ℓ} P \int_{0}^{\infty} d ν (\frac{J_{A} (ν)}{ω + ν - ℓ ν_{s}} + \frac{J_{E} (ν)}{ω - ν - ℓ ν_{s}}) . \end{matrix}

(40)

Equations and show that the effect of the permanent dipoles within

Γ_{1} (ω)

is to introduce manifolds of harmonic levels [in the single-mode truncation of

K (ε)

⁠, there is only one manifold] to/from which transitions can occur, and which influence the Lamb shift value. This is commonly observed when vibrational displacement interactions occur simultaneously with transition dipole interactions and manifest in spectra as vibrational side bands and Stokes's shift. However, in our model, this is a purely photonic effect and it results in a photon sideband in the spectrum of the system.

Γ_{SOME} (ω)

in Eq. is recovered from

Γ_{1} (ω)

in the limit of no permanent dipoles because

A_{ℓ} \to δ_{ℓ 0}

⁠.

In Fig. 2, we illustrate decay-by-emission transitions corresponding to $ℓ = 1$ ⁠, and decay-by-absorption transitions for $ℓ = ℓ_{*}$ where $Round (η / ν_{s}) < ℓ_{*} < Round (η / ν_{s}) + 1$ ⁠. The $A_{ℓ}$ values can be interpreted as the probabilities for a decay between levels with energy difference $η - ℓ ν_{s}$ to occur.

Fig. 2.

View large Download slide

Illustration of decay processes (⁠ $ω = + η$ ⁠) captured by the $ℓ = 1$ and $ℓ = ℓ_{*}$ terms in Eq. (39), within the single-mode truncation of $K (ε)$ ⁠. Decay-by-emission processes have dotted–dashed arrows, and decay-by-absorption processes have solid arrows. The type of decay that occurs for a given value of $ℓ$ depends on the sign of $η - ℓ ν_{s}$ ⁠. In the absence of permanent dipoles, $A_{ℓ} = δ_{ℓ 0}$ and so only the $ℓ = 0$ decay processes are possible.

B. Two-photon processes, Γ₂(ω)

We now consider the two-photon processes in Eq. that occur independent of the driving,

Γ_{2} (ω)

⁠. These processes require that the σ_z and

σ_{\pm}

interactions do not commute, and so are unique to permanent dipole interactions. The rate function is

\begin{matrix} Γ_{2} (ω) = 4 κ^{2} Ω_{μ Δ}^{2} {cos}^{2} (θ_{μ Δ}) \sum_{n, m \in {- 1, 1}} \int_{0}^{\infty} d s e^{i ω s + ϕ (s)} {\tilde{χ}}_{n, m} (s), \end{matrix}

(41)

where

{\tilde{χ}}_{n, m} (s) = \int_{0}^{\infty} d ν \int_{0}^{\infty} d ν' χ_{n, m} (ν, ν') e^{in ν s} e^{i m ν' s}

(42)

is the inverse Fourier transform of

χ_{n, m} (ν, ν')

which we write in matrix notation as

χ (ν, ν') = \bar{J} (ν) \bar{J} (ν') (\begin{matrix} \tilde{N} (ν) \tilde{N} (ν') & - \tilde{N} (ν) N (ν') \\ - N (ν) \tilde{N} (ν') & N (ν) N (ν') \end{matrix}),

(43)

where, for example,

χ_{- 1, 1} (ν, ν') = - \bar{J} (ν) \bar{J} (ν') \tilde{N} (ν) N (ν')

⁠. We have introduced the polarized spectral density

\bar{J} (ν) = \frac{J (ν)}{ν},

(44)

which will be relevant to all transitions resulting from permanent dipole interactions.

The two-frequency Fourier transform in Eq. (42) indicates that $Γ_{2} (ω)$ describes two simultaneous processes. From the phase factors in Eq. (42), when n or m indices are equal to −1, we deal with emission processes, and when they are equal to +1, we deal with absorption, into the frequency channel ν and $ν'$ for n and m, respectively. This interpretation is reinforced by the positions of the factors of $N (ν), N (ν'), \tilde{N} (ν)$ ⁠, and $\tilde{N} (ν')$ in the matrix in Eq. (43). Similar to the case of the one photon processes, the factor of $exp [ϕ (s) - ϕ (0)]$ in Eq. (41) will introduce a photon sideband to the overall process, again enabling processes with more photons.

Before we move onto the processes induced by driving, in Fig. 3, we plot the decay rate $γ_{↓}$ as a function of ν_c and $| d_{Δ} |$ for V = 0. Since V = 0, the eigenbasis coincides with the ${| e ⟩, | g ⟩}$ basis and the decay rate is $γ_{↓} = 2 ℜ [Γ (ϵ)]$ ⁠, where $Γ (ω)$ is in Eq. (30).

Fig. 3.

View large Download slide

Polaron frame decay rate without driving. The total decay rate (dashed black), and the contributions of the one photon processes (solid green) and two photon processes (dot–dashed purple) vs ν_c. In (a)–(d), the value of $| d_{Δ} |$ is $0 D, 100 D, 200 D$ ⁠, and $400 D$ ⁠, respectively. Additionally, V = 0, $ϵ = 3.16 \sim eV$ ⁠, and $θ_{μ Δ} = 0$ (parallel dipole moments). For the other parameters, see the discussion in Sec. III. When varying ν_c, we keep S fixed to its value that ensures a bare emitter lifetime of $1 ns$ at $ν_{c} = 1 eV$ ⁠, and the gray shaded region corresponds to realistic values of ν_c.

Figure 3 shows that strong permanent dipoles suppress one photon transition rates and, for large enough $| d_{Δ} |$ ⁠, cause two photon processes to dominate the rate. However, for the parameter regime shown, corresponding to gamma-globulin, two-photon processes only dominate at ν_c values smaller than are permitted by the electric dipole approximation, indicated by the gray region. As we will show later, if the transition energy is reduced to $ϵ ≲ 1 eV$ ⁠, two-photon processes can contribute significantly and negatively for realistic values of ν_c and $| d_{Δ} | \sim 100 D$ ⁠.

C. Driving-induced one-photon processes, Γ_V,1(ω)

We now move onto the processes in Eq. induced by the driving. The function for such one-photon processes is

\begin{matrix} Γ_{V, 1} (ω) = 4 π Ω_{μ Δ} | V | cos (ϑ_{μ V}) cos (θ_{μ Δ}) \\ \times \int_{0}^{\infty} d ν (- {\bar{J}}_{A} (ν) K (ω + ν) + {\bar{J}}_{E} (ν) K (ω - ν)), \end{matrix}

(45)

where

K (ε)

is given in Eq. and, analogously to Eqs. , we have introduced the polarized absorption and emission spectral densities

{\bar{J}}_{A} (ν) = \bar{J} (ν) N (ν),

(46)

{\bar{J}}_{E} (ν) = \bar{J} (ν) \tilde{N} (ν),

(47)

where

\bar{J} (ν)

is in Eq. . Similarly to the driving-independent two-photon processes,

Γ_{2} (ω)

in Eq. , these processes also arise due to the non-commutativity of the permanent and transition dipole interactions and vanish when

d_{μ} \cdot d_{Δ} = 0

⁠.

Using the single-mode truncation solution of

K (ε)

in Eq. , we can decompose

Γ_{V, 1} (ω)

into real and imaginary parts

Γ_{V, 1} (ω) = \frac{1}{2} γ_{V, 1} (ω) + i S_{V, 1} (ω),

(48)

which we do not write explicitly. Therefore, driving the system generates additional one photon (plus sideband photons) transitions, similar in nature to the transitions described by

Γ_{1} (ω)

but here scaling with

Ω_{μ Δ} | V |

⁠, and dependent on the phase

ϑ_{μ V}

and polarized spectral density.

Importantly, the prefactor in Eq. (45) is proportional to $cos (ϑ_{μ V}) cos (θ_{μ Δ})$ ⁠, which allows one to control whether this function suppresses or enhances the rates and Lamb shifts. As we discuss later, this control could be easy to achieve in practice by tuning the driving phase.

D. Driving-induced zero-photon processes, Γ_V,0(ω)

Finally, the zero-photon processes induced by the driving in Eq. are described by

Γ_{V, 0} (ω) = κ^{2} | V |^{2} \int_{0}^{\infty} d s e^{i ω s} (e^{ϕ (s)} - 1) .

(49)

This term vanishes when

ϕ (s) = 0

and so is generated entirely through the photon sideband. Equation also appears in the driven spin boson model where it describes phonon mediated transfer between eigenstates and so is an effect of pure dephasing.

E. Numerical analysis of the rates

In Fig. 4, panel (a) shows the total decay rate for different driving strengths V as a function of $| d_{Δ} |$ ⁠, and panel (b) shows the four contributions to the rate with $V = 10^{- 4} eV$ ⁠. Due to detailed balance, the excitation rate can be obtained from Fig. 4 by $γ_{↑} = e^{- β η} γ_{↓}$ and, provided pure dephasing is negligible (which is true for the parameters in Fig. 4 because $2 κ | V | ⪡ ϵ$ ⁠), the decoherence rate can be obtained from Fig. 4 by $γ_{d} = (1 + e^{- β η}) γ_{↓} / 2$ ⁠. The key points resulting from the analysis of Fig. 4 are as follows: (1) one can both increase and decrease rates through the physical mechanisms enabled by the presence of permanent dipoles and driving, particularly by controlling $Γ_{V, 1} (ω)$ processes, and (2) all physical processes found in the polaron ECFs contribute non-negligibly to the total rates.

Fig. 4.

View large Download slide

Polaron frame decay rate for different permanent dipole and driving strengths. (a) shows the total decay rate for different driving strengths. In all curves, $ϑ_{μ} = 0$ and $ϑ_{V} \in {0, π}$ control the sign of V given in the legend. (b) shows the four contributions (colored curves) to the total rate (dashed black), given in Eq. (30), for the $V = 10^{- 4} eV$ parameter set. The dashed black curves in both panels (a) and (b) are the same. The markers in panel (b) refer to $| d_{Δ} |$ values used in Fig. 5. The left-hand side of each panel has anti-parallel $d_{μ}$ and $d_{Δ}$ ⁠, while the right-hand side has parallel dipole moments, which is indicated by the arrows at the top of the figure. $ϵ = 0.1 \sim eV$ and other parameters are discussed in Sec. III.

Through control of $| d_{Δ} |$ and V, one can completely suppress the transition and decoherence rates. The decay rate at each of the minima in the curves for $V \in {- 10^{- 4}, 10^{- 4}, 2 \times 10^{- 4}} eV$ in Fig. 4(a) is zero to within numerical error, and one can show that the excitation and decoherence rates are also zero at these minima. As well as tuning $| d_{Δ} |$ ⁠, similar control can be achieved through the phase of the driving. In Fig. 5, we demonstrate this by varying $ϑ_{μ V}$ with $V = 10^{- 4} eV$ and the $| d_{Δ} |$ values marked in Fig. 4(b) and given in the legend of Fig. 5. Maximum suppression is achieved when the driving is in, or out of, phase with the transition dipole moment, depending on the alignment of the dipole moments. This is because of the factor of $cos (θ_{μ Δ}) cos (ϑ_{μ V})$ in Eq. (45).

Fig. 5.

View large Download slide

Phase control of the polaron frame decay rate. The relative phase $ϑ_{μ V}$ is varied with $| d_{μ} | = 0.1 D$ and $| V | = 10^{- 4} eV$ fixed. Each curve has a different value of $| d_{Δ} |$ indicated in the legend, which correspond to the markers in Fig. 4(b). Other parameters are discussed in Sec. III.

Figures 4 and 5 exemplify the impact of permanent dipoles on the ability to control quantum systems. Beyond reducing the decoherence in the quantum system, another critical challenge in quantum computing is the ability to relax systems faster to be reset for computation. The inverse optimization of maximizing the decay rate can also be achieved, allowing for a scheme for faster relaxation of such systems. For example, let us take the $| d_{Δ} | = 50 D$ case in Fig. 5 with aligned dipoles. During computation, one could drive the system with phase $ϑ_{μ V} = 0$ so that decoherence from optical sources is minimized. Once computation is complete, one could then increase the drive phase to $ϑ_{μ V} = π$ so that the decay rate is maximized. Moreover, this same control could be used to operate an energy storage device: by changing the phase of the drive, one can change the operation of the device from storage to discharge.

V. COMPARISON TO EXACT DYNAMICS

In this section, we compare the predictions for the system density operator found using the PFME to those with the DFME and a numerically exact approach, TEMPO. (See Refs. 46–54 for discussions on the TEMPO algorithm.) In Fig. 6, we plot the population of the ground state and the coherence between the excited and ground state against time calculated using the three methods. For this comparison, we do not use the dipole moments and emitter energies discussed in Sec. III, and instead choose these to emphasize the accuracy of the PFME in a strong coupling regime.

Fig. 6.

View large Download slide

Comparison of the PFME predictions to exact numerical approach and the DFME. (a) and (c) are a comparison of $ρ_{g g} (t)$ calculated using the PFME (green), TEMPO (dot-dashed black), and the DFME (dashed orange). (b) and (d) are a similar comparison but for $| ρ_{e g} (t) |$ and note the discontinuity in the time axis. In (d), the cyan dotted curve is the polaron prediction for the coherence when artificially accounting for the initial non-Markovian slip. In (a) and (b), $ϵ = 0.363 eV, V = - 0.006 37 eV$ ⁠, and $| d_{Δ} | = 292 D$ ⁠, and in (c) and (d), $ϵ = 1 eV, V = 0 eV$ ⁠, and $| d_{Δ} | = 584 D$ ⁠. In all plots, $ϑ_{μ V} = θ_{μ Δ} = 0$ and $| d_{μ} | = 5.84 D$ ⁠, and other parameters are discussed in Sec. III. To make the comparison to TEMPO easier, we have assumed that the field aligns with the dipole moments such that all dipoles are parallel and we have ignored the factor of $8 π / 3$ in the Ω_pq in Eqs. (22). We have taken this into account when evaluating the inequality in Eq. (50) in the text for these parameter sets. We use the non-secular PFME and DFME in our results, although the non-secular terms merely introduce oscillations to the population dynamics that cannot be resolved on the scales shown.

Since the DFME does not capture many of the unique processes attributed to the permanent dipoles, we expect it to break down when these are strong. Motivated by the Hamiltonian in Eq. , an approximate boundary of this regime is

\sum_{k} f_{k}^{2} {(d_{Δ} \cdot e_{k})}^{2} = Ω_{Δ Δ} \int_{0}^{\infty} d ν J (ν) ≳ η^{2} .

(50)

For the parameters used in Figs. 6(a) and 6(b) (see the caption), the inequality in Eq. reads

0.69 > 0.36

and

1.4 > 1.0

⁠, respectively. Thus, both parameter sets constitute strong permanent dipole interactions. Despite this, the PFME predicts qualitatively correct dynamics and agrees well with TEMPO, and, as expected, the DFME does not. This emphasizes that each new process in Fig. 4 is required to correctly describe strong permanent dipoles.

A phenomenon not captured by the PFME or any Markovian master equation is the so-called “slip”⁶³ at short times, which is most evident in the coherence evolution in panel (d). However, as shown by the cyan dotted curve, if this is artificially accounted for in the PFME by starting in a coherent initial state, the PFME again agrees well over the long time duration, indicating it captures the main decoherence mechanisms.

VI. EMISSION SPECTRUM

Due to the many new physical processes generated by the presence of permanent dipoles discussed in Sec. IV, we expect alterations to the emission spectrum. In the absence of permanent dipoles, the spectrum consists of a Mollow triplet with peaks at frequencies $- η, 0, η$ ⁠. When V = 0, the triplet becomes a single peak at frequency ϵ with a width determined by the decoherence rate in the standard optical master equation: $π J (ϵ) [1 + 2 N (ϵ)]$ ⁠. In the presence of strong permanent dipoles, we expect that there will be a photon sideband extending to negative frequencies. Moreover, the positions of the Mollow triplet peaks will shift from $\pm η$ to $\pm \bar{η}$ given in Eq. (28), and the widths of the peaks will be determined by γ_d in Eq. (26b).

As we prove in Appendix G, the emission spectrum is given by

I (ω) = α_{prop} (r, R, ω) I_{0} (ω),

(51)

where

α_{prop} (r, R, ω) = | d_{μ} \cdot G (r, R, ω) |^{2}

and

G (r, R, ω)

is the Green's function of the medium.^64–68

α_{prop} (r, R, ω)

accounts for propagation and filtering of the light from the dipole to the detector at positions r and R, respectively.⁶⁴ The polarization spectrum is

I_{0} (ω) = lim_{t \to \infty} ℜ [\int_{0}^{\infty} d τ ⟨ σ_{+} (t + τ) σ_{-} (t) ⟩ e^{- i ω τ}] .

(52)

In the following, we focus on evaluating the polarization spectrum, which captures entirely the effects associated with the permanent dipoles. The specifics of the experimental setup, for example, whether the dipole is coupled to a cavity or a waveguide, are described by the

α_{prop} (r, R, ω)

term, which can be calculated separately.

The expectation value in Eq. is taken with respect to the lab frame density operator ρ_l. Using

ρ_{l} = B (- D / ν) U^{†} ρ_{p} U B (D / ν)

⁠, where U is the polaron transformation in Eq. ,

B (\pm D / ν)

are displacement operators, and ρ_p is the polaron frame density operator, we have

\begin{matrix} ⟨ σ_{+} (t + τ) σ_{-} (t) ⟩ = κ^{2} e^{ϕ (τ)} Tr [σ_{+} (t + τ) σ_{-} (t) ρ_{p} (0)] \\ \equiv κ^{2} e^{ϕ (τ)} {⟨ σ_{+} (t + τ) σ_{-} (t) ⟩}_{p}, \end{matrix}

(53)

where

ϕ (τ)

is defined in Eq. . Substituting Eq. into Eq. yields

\begin{matrix} I_{p, 0} (ω) = κ^{2} lim_{t \to \infty} ℜ [\int_{0}^{\infty} d τ e^{ϕ (τ)} \\ {⟨ σ_{+} (t + τ) σ_{-} (t) ⟩}_{p} e^{- i ω τ}] . \end{matrix}

(54)

The two-time correlation function in Eq. can be calculated using the quantum regression theorem (QRT).^1,69 The QRT utilizes the cyclicity of the trace to rewrite the two-time expectation value as

lim_{t \to \infty} {⟨ σ_{+} (t + τ) σ_{-} (t) ⟩}_{p} = Tr [σ^{+} Λ_{p} (τ)],

(55)

where

Λ_{p} (τ) = U_{0}^{†} (τ) Λ_{p} (0) U_{0} (τ)

is a modified density operator with the initial state

Λ_{p} (0) = σ^{-} U_{0} (\infty) ρ_{p} (0) U_{0}^{†} (\infty),

(56)

and, because

U_{0} (τ)

is the time evolution operator defined by the polaron frame Hamiltonian [see Eq. ], the PFME we have derived holds identically for the

Λ_{p} (τ)

operator but with the initial condition given in Eq. . The QRT allows one to convert two-time correlation functions into one-time expectation values of density operators with modified initial conditions.

Due to the polaron transformation, Eq. (54) captures the photonic sideband. However, the QRT contains an implicit Born approximation⁶⁹ as well as any approximations used in the derivation of the master equation. Moreover, the QRT produces spectral lines at the Lambshifted eigenenergies of H_S because the ground-excited coherence ρ_eg oscillates at these frequencies [see Eq. (25b)]. This means that the QRT is unable to capture spectral lines associated with multiphoton transitions. However, multiphoton transitions scale with the square of the spectral density, and so for weak light–matter coupling, we expect the associated spectral lines to be much smaller than the photonic sideband.

Owing to the modified initial state in Eq. , the power of the polarization spectrum contains information about the steady state of the system

P = \int_{- \infty}^{\infty} d ω I_{p, 0} (ω) = π ρ_{e e} (\infty) .

(57)

Additionally, by noting that the polarization spectrum without the photon sideband, denoted by

I_{p, \times} (ω)

⁠, is given by Eq. with the replacement

exp [ϕ (τ)] \to 1

⁠, one can show that

\frac{\int_{- \infty}^{\infty} d ω I_{p, \times} (ω)}{P} = κ^{2} .

(58)

Hence,

κ^{2}

can be interpreted as the fractional emission into the sideband. Remarkably, an integrated spectrum measurement can be used to determine κ, which in turn provides an effective measurement of

Ω_{Δ Δ}

and, therefore, of the strength of the

| d_{Δ} |

dipole moment. In measurements,

I_{p, \times} (ω)

cannot be separated from

I_{p, 0} (ω)

⁠, but as we show in the following,

\int_{- \infty}^{\infty} d ω I_{p, \times} (ω)

can be accurately approximated as the integral of the measured spectrum over the Lorentzian parts of the dominant peak(s).⁷⁰

In Fig. 7, we plot the polarization spectrum calculated in the polaron frame [ $I_{p, 0} (ω)$ ], the spectrum for the same parameters but without the photonic sideband [ $I_{p, \times} (ω)]$ which serves to highlight the sideband, and the polarization spectrum in the absence of permanent dipoles [ $I_{p, npd} (ω)$ ]. Notably, the photon sideband substantially changes the emission spectrum. As seen in the inset, the widths of the peaks are increased by the permanent dipoles; however, the effect of the permanent dipoles on the Lamb shifts is small.

Fig. 7.

View large Download slide

Polarization spectrum of the system. The total polarization spectrum $I_{p, 0} (ω)$ is shown in solid blue, the polarization spectrum without the photon sideband $I_{p, \times} (ω)$ in dotted purple, and the polarization spectrum with $| d_{Δ} | = 0$ (and other parameters unchanged) $I_{p, npd} (ω)$ is shown in dot-dashed green. For clarity, $I_{p, npd} (ω)$ is only shown in the inset, which is a zoom-in of the peaks plotted on a logarithmic scale. We use $| d_{Δ} | = 100 D, V = 0.01 eV$ ⁠, and $ϵ = 0.1 eV$ ⁠, and other parameters discussed in Sec. III. The vertical dashed lines in the inset indicate the boundaries of R in the approximation $\int_{- \infty}^{\infty} d ω I_{p, \times} (ω) \approx \int_{R} d ω I_{p, 0} (ω)$ used to estimate κ. We use the secularized master equation to produce the spectra, which eliminates a known inconsistency near to the negative frequency peak and has no other impact (see Appendix H).

We now introduce a procedure of directly obtaining essential information about the system from the experimentally measurable spectrum, $I_{p, 0} (ω)$ ⁠, and compare the values obtained for the data shown in Fig. 7 to analytic values in Table I. We show that it is possible to evaluate $ρ_{e e} (\infty)$ and κ from $I_{p, 0} (ω)$ using Eqs. (57) and (58) and, provided that the temperature of the experiment is known, these values can be used to obtain ϵ and $| V |$ ⁠. If the spectral density is also known, it is also possible to obtain $| d_{Δ} |$ from κ.

Table I.

Values of parameters obtained from the measurable spectrum. Each column shows analytic values, obtained directly from the equations in this paper, and “measured” values obtained from $I_{p, 0} (ω)$ in Fig. 7.

Param.	Analytic	Meas.
$κ^{2}$	0.632	0.621
$ρ_{e e} (\infty)$	0.450	0.450
$ϵ / eV$	0.100	0.100
$\| V \| / eV$	0.010	0.010
$\| d_{Δ} \| / D$	100	101

As shown in Eq. , the steady-state population of the excited state can straightforwardly be obtained by integrating

I_{p, 0} (ω)

over all frequencies. κ can be approximately found using Eq. with

\int_{- \infty}^{\infty} d ω I_{p, \times} (ω) \approx \int_{R} d ω I_{p, 0} (ω)

⁠, where R are the frequency regions for which the dominant peaks in the spectrum are approximately Lorentzian before the photon sideband begins. In the example in Fig. 7, this region is between the vertical dashed lines in the inset. One could increase the accuracy of the

κ^{2}

estimation in Fig. 7 by integrating over the other two peaks, but these contributions are small. If the spectral density of the photon bath is known, the relation

κ = exp [- ϕ (0) / 2]

can be inverted to obtain

| d_{Δ} |

⁠. The remaining parameters of the system, ϵ and

| V |

⁠, can be estimated by assuming that the system thermalizes with respect to

H_{S} = η τ_{z} / 2

⁠, yielding the steady state

ρ_{e e} (\infty) = \frac{1}{2} [1 - \frac{ϵ}{η} tanh (\frac{β η}{2})] .

(59)

Using Eq. , the measured value of

ρ_{e e} (\infty)

⁠, and assuming negligible Lamb shifts such that the frequency of the peak in the spectrum

\bar{η} \approx η

⁠, one can estimate ϵ. Then, using

\bar{η} \approx η = {(ϵ^{2} + 4 κ^{2} | V |^{2})}^{1 / 2}

⁠, one can estimate

| V |

⁠.

VII. IMPORTANCE OF INITIAL CONDITIONS

It is typical in theoretical quantum optics to consider the system and environment initially in an uncorrelated state,

ρ (0) = ρ_{S} (0) \otimes ρ_{E},

(60)

where the environment is assumed to be in a thermal Gibbs state

ρ_{E} = exp (- β H_{E}) / Z_{E}

with

H_{E} = \sum_{k} ν_{k} a_{k}^{†} a_{k}

⁠. This assumption is valid in weakly coupled systems where the system and the environment are only weakly correlated even after thermalization. However, for strongly coupled systems, this assumption is invalidated, and the total Gibbs state

ρ_{β} = exp (- β H) / Z

should be used as the initial state, where H now refers to the total Hamiltonian. Due to the interaction term, this will no longer be a separable state and so is difficult to model.

A benefit of the polaron transformation is that the separable initial state $ρ_{S} (0) \otimes ρ_{E}$ in the polaron frame models an initial state in the lab frame that is more similar to $ρ_{β}$ ⁠. For example, in our calculations in Fig. 6, $ρ_{S} (0) = | g ⟩ ⟨ g |$ and so $ρ_{S} (0) \otimes ρ_{E}$ in the polaron frame becomes $ρ_{l} (0) = | g ⟩ ⟨ g | \otimes B ([D + Δ] / ν) ρ_{E} B (- [D + Δ] / ν)$ in the lab frame, which is equal to $ρ_{β}$ in the limit of negligible transition dipoles and projected onto the ground state. This is the state of the system after decaying to the ground state from a thermalized state, and so is a good initial state in which to model excitation.

Therefore, to compare the PFME to TEMPO (which operates in the lab frame), we must use the Hamiltonian in Eq. (5) with the initial state $ρ_{l} (0)$ ⁠. However, most numerical techniques, including TEMPO, assume that the initial state is $ρ_{S} (0) \otimes ρ_{E}$ ⁠. This discrepancy can be overcome by deriving an effective Hamiltonian to use in TEMPO such that the environment appears to be in the correct state. For $ρ_{S} (0) = | g ⟩ ⟨ g |$ in the polaron frame, we show in Appendix I that the required Hamiltonian is $\tilde{H} = {\tilde{H}}_{0} + {\tilde{H}}_{I}$ ⁠, where

(61)

{\tilde{H}}_{0} = \frac{\tilde{ϵ}}{2} σ_{z} + \tilde{V} σ_{+} + {\tilde{V}}^{*} σ_{-} + \sum_{k} ν_{k} a_{k}^{†} a_{k},

(61a)

{\tilde{H}}_{I} = π_{Δ Δ} (I + σ^{z}) + π_{μ \bar{μ}} σ_{+} + π_{\bar{μ} μ} σ_{-},

(61b)

and π_pq is given in Eq. ,

\tilde{ϵ} = ϵ + 2 G_{Δ Δ}, \tilde{V} = V + G_{μ \bar{μ}}

⁠, where

G_{p q} = \sum_{k} (p_{k} Δ_{k}^{*} + q_{k}^{*} Δ_{k}) / ν_{k}

⁠. We note that because the displacement direction of the polaron transformation is state dependent, the necessary effective Hamiltonian depends on

ρ_{S} (0)

⁠.

By comparing the Hamiltonians in Eqs. (61) and (5), we see stark differences if both are assumed to have the initial state $ρ_{S} (0) \otimes ρ_{E}$ ⁠. This includes a renormalization of the transition energy $ϵ \to \tilde{ϵ}$ and the driving term $V \to \tilde{V}$ ⁠. We can understand these differences by noting that the environment is far from equilibrium for a strongly coupled system, and initially, no optical polarons exist in our system. Dynamically created optical polarons effectively introduce a strong restoring force in the system that scales with the strength of the permanent dipoles. The ensuing dynamics are thus considerably different, exacerbating the need for careful consideration of the initial conditions of the physical models we deploy.

Identity type interactions, such as $π_{D D} I'$ in Eq. (5), have a similar effect. As we showed in Eq. (10), this type of interaction can be removed through a displacement transformation and so, ultimately, the role of an identity interaction is to change the initial environment state from a thermal Gibbs state to a displaced one.

VIII. CONCLUSION

We have studied a driven quantum optical system with strong permanent dipole moments associated with molecular orbital asymmetry. The optical polaron transformation, which captures the polarization of photonic modes caused by the permanent dipoles, allows for the construction of an optical master equation perturbative only in the transition dipole moment and driving strength and provides an intuitive formalism to understand the effects associated with the presence of the permanent dipoles on the system dynamics and emission spectrum.

We have shown three key results. (1) In Sec. IV, we showed that transition and decoherence rates can be engineered for practical application by exploiting permanent dipoles, for example, to eliminate optical decoherence and transitions. By using the novel physical processes explicit in the polaron rate equations we derive, one can design systems to exploit these effects. The novel physical processes arising from the permanent dipoles are as follows: a photon sideband with a relative contribution to the emission spectrum scaling as $κ^{2}$ ⁠, two-photon processes scaling as $Ω_{μ Δ}^{2}$ ⁠, and a term linear in the driving amplitude scaling as $| V | Ω_{μ Δ} cos (ϑ_{μ Ω})$ ⁠. (2) In Sec. V, we proved that the optical polaron description provides a much more accurate master equation by comparing it to the DFME and TEMPO. (3) In Sec. VI, we indicated distinguishable features of permanent dipoles in emission spectra and described possible measurements to obtain $ρ_{e e} (\infty)$ ⁠, κ, $\bar{η}$ and by extension the bare energy splitting ϵ, the driving strength $| V |$ ⁠, and the permanent dipole magnitude $| d_{Δ} |$ ⁠. The distinguishable features are the photonic sideband and the altered widths and positions of peaks.

Furthermore, two of the processes we identified, namely, the $Γ_{2} (ω)$ and $Γ_{V, 1} (ω)$ processes, originate from the non-commutativity of the bosonic operators on the transition and permanent dipole interactions. These vanish if $d_{μ}$ and $d_{Δ}$ are perpendicular and are unique to pure dephasing interactions that are induced by permanent dipoles, as opposed to vibrational environments. The linear driving term, $Γ_{V, 1} (ω)$ ⁠, may be particularly useful for quantum technologies. As we have shown, one could engineer devices in which the optical decoherence and transition rates can be varied from essentially zero to faster than those in the same system but without permanent dipoles, by simply changing the driving phase by π. This has potential applications in quantum computing and energy storage.

There are many interesting features of optical polarons that warrant future exploration. For example, how the transition rates are affected by embedding the dipole in a structured,⁷¹ or anisotropic,⁷² dielectric medium, common to many biological systems. As discussed in the main text, such systems usually possess stronger permanent dipoles but have more complicated spectral densities, requiring more sophisticated theories—such as variational polaron theory⁶⁹—to model. Moreover, many asymmetric systems couple strongly to vibrational baths and the interplay between photonic and vibrational physics leads to non-additive and non-equilibrium phenomena such as population inversion.^50,73 How permanent dipoles affect these phenomena are still an open question. Additionally, it has been shown in Refs. 28, 30, and 31 that the interplay between the pure dephasing and dissipative interactions leads to non-zero coherences in the steady state. It would be of great interest to explore the nature of the steady-state coherences in the optical polaron formalism.

ACKNOWLEDGMENTS

We would like to thank the Tempo Collaboration for use of the open-source code Oqupy.⁵⁵ In particular, we thank Gerald Fux for extremely insightful conversations on the use of Oqupy. D.M.R. also thanks Ahsan Nazir and Owen Diba for helpful discussions.

The work by A.B. was supported by the Leverhulme Quantum Biology Doctoral Training Centre at the University of Surrey funded by a Leverhulme Trust training center (Grant No. DS-2017-079), the EPSRC Strategic Equipment [Grant No. EP/L02263X/1 (EP/M008576/1)], and EPSRC (Grant No. EP/M027791/1) awards to M.F. D.M.R. was supported by the EPSRC (Grant No. EP/T517896/1).

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Adam Burgess: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Marian Florescu: Conceptualization (equal); Data curation (supporting); Formal analysis (supporting); Funding acquisition (equal); Investigation (supporting); Methodology (supporting); Project administration (supporting); Resources (equal); Software (equal); Supervision (supporting); Validation (supporting); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (equal). Dominic Michael Rouse: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

DATA AVAILABILITY

The data that support the findings of this study are openly available at http://dx.doi.org/10.5525/gla.researchdata.1437, Ref. 74.

APPENDIX A: DISPLACEMENT OPERATORS

Within the appendixes, we will regularly use many identities involving displacement operators, which are proven in Ref. 75. This appendix is dedicated to listing the necessary displacement operator identities. For the purpose of this appendix, we use displacement operators with a single-photon mode

B (α) = exp [α a^{†} - α^{*} a],

(A1)

and note that

B {(α)}^{†} = B (- α)

and

B (- α) B (α) = I

⁠. The first identity is the action of a displacement operator on harmonic operators

(B2)

B (\pm α) a^{†} B (\mp α) = a^{†} \mp α^{*},

(A2a)

B (\pm α) a B (\mp α) = a \mp α .

(A2b)

The second is that the product of two displacement operators is

B (α_{1}) B (α_{2}) = B (α_{1} + α_{2}) e^{\frac{1}{2} (α_{1} α_{2}^{*} - α_{1}^{*} α_{2})} .

(A3)

The third is the action of a displacement operator on the vacuum state to generate a coherent state

B (α) | 0 ⟩ = | α ⟩ .

(A4)

Fourth, that the expectation value of an operator with respect to the thermal state

ρ_{E} = exp [- β ν a^{†} a] / Z_{E}

can be written as an integral over coherent states as

{Tr}_{E} [O ρ_{E}] = \frac{1}{π N (ν)} \int_{- \infty}^{\infty} d^{2} α e^{- \frac{{| α |}^{2}}{N (ν)}} ⟨ α | O | α ⟩,

(A5)

where

\int_{- \infty}^{\infty} d^{2} α = \int_{- \infty}^{\infty} d ℑ [α] \int_{- \infty}^{\infty} d ℜ [α]

and

N (ν) = {(e^{β ν} - 1)}^{- 1}

is the Bose–Einstein distribution. The final identity is the expectation value of a displacement operator with respect to the vacuum state

⟨ 0 | B (α) | 0 ⟩ = e^{- \frac{1}{2} {| α |}^{2}} .

(A6)

APPENDIX B: CALCULATION OF $⟨ C ⟩$

Recall from Eq. that

C = B (δ) π_{μ \bar{μ}} B (δ) e^{- i ϑ_{V}} + | V | B (2 δ)

⁠, where

δ_{k} = Δ_{k} / ν_{k}

⁠. Using the properties of displacement operators listed in Appendix A, this can be written as

C = B (2 δ) (\sum_{k} [μ_{k} (a_{k}^{†} + δ_{k}^{*}) + {\bar{μ}}_{k}^{*} (a_{k} + δ_{k})] e^{- i ϑ_{V}} + | V |) .

(B1)

To calculate

⟨ C ⟩ = {Tr}_{E} [C ρ_{E}]

⁠, we require

⟨ B (2 δ) ⟩, ⟨ B (2 δ) a_{k}^{†} ⟩

⁠, and

⟨ B (2 δ) a_{k} ⟩

⁠. We will perform these calculations explicitly, because the polaron frame environment correlation functions require analogous, but more algebraically involved, mathematics. We do the calculations with a single photonic mode, where

B (α) = exp [α a^{†} - α a]

⁠, and reinstate the multimode summations at the end. To calculate the expectation values, we aim to use Eq. .

Starting with

⟨ B (2 δ) ⟩

⁠, we first evaluate the integrand of Eq. as

\begin{matrix} ⟨ α | B (2 δ) | α ⟩ = ⟨ 0 | B (- α) B (2 δ) B (α) | 0 ⟩ \\ = ⟨ 0 | B (2 δ) | 0 ⟩ e^{2 δ α^{*} - 2 δ^{*} α} \\ = e^{- 2 {| δ |}^{2}} e^{2 δ α^{*} - 2 δ^{*} α}, \end{matrix}

(B2)

where in the first equality, we have used Eq. twice, in the second equality Eq. twice, and in the final equality Eq. . Substituting Eq. into Eq. then yields

⟨ B (2 δ) ⟩ = e^{- 2 {| δ |}^{2} coth (β ν / 2)} \equiv κ .

(B3)

Evaluating

⟨ B (2 δ) a^{†} ⟩

follows a similar procedure; however, one must now use Eqs. to move the creation operator such that it annihilates with the vacuum bra-state

⟨ 0 |

⁠. That is,

\begin{matrix} ⟨ α | B (2 δ) a^{†} | α ⟩ = ⟨ 0 | B (- α) B (2 δ) a^{†} | α ⟩ \\ = ⟨ 0 | B (- α) [a^{†} - 2 δ^{*}] B (2 δ) | α ⟩ \\ = [α^{*} - 2 δ^{*}] ⟨ α | B (2 δ) | α ⟩ \\ = [α^{*} - 2 δ^{*}] e^{- 2 {| δ |}^{2}} e^{2 δ α^{*} - 2 δ^{*} α} . \end{matrix}

(B4)

Substituting Eq. into Eq. and performing the integrations yields

⟨ B (2 δ) a^{†} ⟩ = - 2 κ δ^{*} [1 + N (ν)] .

(B5)

The final expectation value,

⟨ B (2 δ) a ⟩

⁠, is easier to calculate because a annihilates with

| 0 ⟩

⁠. One finds that

⟨ α | B (2 δ) a | α ⟩ = α e^{- 2 {| δ |}^{2}} e^{2 δ α^{*} - 2 δ^{*} α},

(B6)

and so,

⟨ B (2 δ) a ⟩ = 2 κ δ N (ν) .

(B7)

Collecting the expectation values in Eqs. , , and and substituting these into

⟨ C ⟩

with C given in Eq. yield

⟨ C ⟩ = κ [coth (\frac{β ν}{2}) ({\bar{μ}}^{*} δ - μ δ^{*}) e^{- i ϑ_{V}} + | V |] .

(B8)

Since

d_{Δ} \in ℜ

⁠, we find that

{\bar{μ}}^{*} δ - μ δ^{*} = 0

⁠, and so,

⟨ C ⟩ = κ | V |

as stated in the main text.

APPENDIX C: THE HUANG–RHYS PARAMETER

In the main text, we write the dipole vectors in terms of a reference value, as $d_{p} = d_{ref} {\tilde{d}}_{p}$ ⁠, which ensures that the spectral density in Eq. (24) has dimensions of energy. Consequently, the value of S is a function of $d_{ref}$ ⁠, and in this appendix, we derive this relationship.

The proof starts from the free space continuum limit

\sum_{k} \to \frac{V}{{(2 π)}^{3}} \sum_{λ} \int_{0}^{\infty} d ν \int_{Ω_{k}} d Ω_{k},

(C1)

which can be found, for example, in Eq. (2.30) of Ref. 2. Consequently, the continuum limit of the left-hand side of Eq. is

\sum_{k} f_{k}^{2} F (ν_{k}) (d_{p} \cdot e_{k}) (d_{q} \cdot e_{k}) \to h_{p q} \frac{d_{ref}^{2}}{2 {(2 π)}^{3}} \int_{0}^{\infty} d ν ν^{3} e^{- (ν / ν_{c})} F (ν) = h_{p q} A_{0} .

(C2)

We have used that

f_{k} = \sqrt{ν_{k} / (2 V)} exp [- ν / (2 ν_{c})]

where the exponential factor accounts for the phenomenological high-frequency cutoff, and defined the spectral area

A_{0} = \int_{0}^{\infty} d ν J (ν) = 6 S ν_{c}^{2} .

(C3)

The other parameters are defined in the main text. By comparing the right-hand sides of Eqs. and , one finds that the Huang–Rhys parameter in the free space spectral density must be related to the reference dipole magnitude by

S = \frac{{(ν_{c} d_{ref})}^{2}}{2 {(2 π)}^{3}} .

(C4)

In the main text, we multiply the Huang–Rhys parameter by the dimensionless Purcell enhancement factor s_0, which is justified for weak field interactions.

APPENDIX D: NON-SECULAR MASTER EQUATIONS

The non-secular master equations in both the polaron and displaced frames have the same forms, only differing by which operators enter the two-time correlation functions, and therefore, the rates and energies in the master equation below will be different. We find that

(E1)

\partial_{t} ρ_{+ +} (t) = - γ_{↓} ρ_{+ +} (t) + γ_{↑} ρ_{- -} (t) + \bar{γ} ρ_{- +} (t) + {\bar{γ}}^{*} ρ_{+ -} (t),

(D1a)

\partial_{t} ρ_{+ -} (t) = - [γ_{d} + i \bar{η}] ρ_{+ -} (t) + k_{1} ρ_{- +} (t) + k_{-} ρ_{- -} (t) + k_{+}^{*} ρ_{+ +} (t),

(D1b)

and

\partial_{t} ρ_{- -} (t) = - \partial_{t} ρ_{+ +} (t)

and

\partial_{t} ρ_{- +} (t) = \partial_{t} ρ_{- +} {(t)}^{†}

⁠. As defined in the main text, the secular terms are

γ_{\overset{↓}{↑}} = 2 ℜ [Γ_{\mp \mp} (\pm η)],

(D2)

γ_{d} = \frac{1}{2} [γ_{↑} + γ_{↓}] + 4 ℜ [γ_{z z} (0)],

(D3)

\bar{η} = η + ℑ [Γ_{- -} (η) - Γ_{+ +} (- η)] .

(D4)

These are written in terms of the Markovian ECFs as

Γ_{α β} (ω) = \int_{0}^{\infty} d s e^{i ω s} ⟨ g_{α}^{†} (s) g_{β} (0) ⟩,

(D5)

for

α \in {z, +, -}

⁠, and it is the expressions for

g_{α}

that vary between the PMFE and DFME, and these are given in the main text in Eqs. and Appendix F, respectively. The non-secular rates are

\bar{γ} = Γ_{- z} (0) + Γ_{+ z} {(0)}^{*},

(D6a)

k_{1} = Γ_{- +} (- η) + Γ_{+ -} {(η)}^{*},

(D6b)

k_{\pm} = \mp Γ_{\pm z} (0) \pm Γ_{\mp z} {(0)}^{*} \pm 2 Γ_{z \mp} (\pm η) .

(D6c)

We evaluate

Γ_{α β} (ω)

for the polaron frame in Appendix E and for the displaced frame in Appendix F.

Notice that if the coherences are initially zero (⁠ $ρ_{+ -} (0) = ρ_{- +} (0) = 0$ ⁠), then if $k_{-} = k_{+} = 0$ as well, the coherences will be zero at all times t. This occurs in the polaron frame if $V \to 0$ since $g_{z} \propto sin (φ) \to 0$ ⁠, i.e., the eigenstates fully localize, but this does not happen in the displaced frame in the same limit because $g_{z}^{d}$ does not become zero.

APPENDIX E: CALCULATION OF THE POLARON FRAME TWO-TIME CORRELATION FUNCTIONS

The environment correlation functions, which determine the second-order Born–Markov rates, depend on the two-time correlation functions

⟨ g_{α}^{†} (s) g_{β} (0) ⟩

for

α, β \in {z, +, -}

where recalling from Eq. ,

g_{z} = \frac{1}{2} sin (φ) [C + C^{†}],

(E1a)

g_{+} = [{cos}^{2} (\frac{φ}{2}) C - {sin}^{2} (\frac{φ}{2}) C^{†}],

(E1b)

and

g_{-} = g_{+}^{†}

where

C = C - κ | V |

and

C = B_{+} (\sum_{k} [μ_{k} (a_{k}^{†} + δ_{k}^{*}) + {\bar{μ}}_{k}^{*} (a_{k} + δ_{k})] e^{- i ϑ_{V}} + | V |),

(E2)

and

B_{\pm} \equiv B (\pm 2 δ)

with

δ_{k} = Δ_{k} / ν_{k}

⁠. The interaction picture forms of the relevant operators are

a_{k} (s) = a_{k} e^{- i ν_{k} s},

(E3)

B_{\pm} (s) = B (\pm 2 δ e^{i ν s}) .

(E4)

1. Preliminary calculations

All $⟨ g_{α}^{†} (s) g_{β} (0) ⟩$ depend on four two-time correlation functions: $⟨ C^{†} (s) C (0) ⟩, ⟨ C (s) C^{†} (0) ⟩, ⟨ C (s) C (0) ⟩$ ⁠, and $⟨ C^{†} (s) C^{†} (0) ⟩$ ⁠. Each of these in turn depends on either seven or nine unique two-time correlation functions involving a_k, $a_{k}^{†}$ ⁠, and $B_{\pm}$ ⁠. In this appendix, we list the results of all necessary two-time correlation functions, which are each derived using the same mathematics as in the explicit examples in Appendix B.

a. $⟨ C^{†} (s) C (0) ⟩$

To calculate this two-time correlation function, we require

(F5)

⟨ B_{-} (s) B_{+} (0) ⟩ = κ^{2} e^{ϕ},

(E5a)

⟨ a_{k}^{†} B_{-} (s) B_{+} (0) ⟩ = - N_{k} x_{k}^{*} κ^{2} e^{ϕ},

(E5b)

⟨ a_{k} B_{-} (s) B_{+} (0) ⟩ = {\tilde{N}}_{k} x_{k} κ^{2} e^{ϕ},

(E5c)

⟨ B_{-} (s) B_{+} (0) a_{k}^{†} ⟩ = - {\tilde{N}}_{k} x_{k}^{*} κ^{2} e^{ϕ},

(E5d)

⟨ B_{-} (s) B_{+} (0) a_{k} ⟩ = N_{k} x_{k} κ^{2} e^{ϕ},

(E5e)

⟨ a_{k}^{†} B_{-} (s) B_{+} (0) a_{q}^{†} ⟩ = N_{k} {\tilde{N}}_{q} x_{k}^{*} x_{q}^{*} κ^{2} e^{ϕ},

(E5f)

⟨ a_{k} B_{-} (s) B_{+} (0) a_{q} ⟩ = {\tilde{N}}_{k} N_{q} x_{k} x_{q} κ^{2} e^{ϕ},

(E5g)

⟨ a_{k}^{†} B_{-} (s) B_{+} (0) a_{q} ⟩ = (- N_{k} N_{q} x_{k}^{*} x_{q} + N_{k} δ_{k q}) κ^{2} e^{ϕ},

(E5h)

⟨ a_{k} B_{-} (s) B_{+} (0) a_{q}^{†} ⟩ = (- {\tilde{N}}_{k} {\tilde{N}}_{q} x_{k} x_{q}^{*} + {\tilde{N}}_{k} δ_{k q}) κ^{2} e^{ϕ},

(E5i)

where, for brevity,

N_{k} \equiv N (ν_{k})

is the Bose–Einstein distribution,

{\tilde{N}}_{k} \equiv \tilde{N} (ν_{k}) = 1 + N (ν_{k}), ϕ \equiv ϕ (s)

given in Eq. ,

κ^{2} = exp [- ϕ (0)]

⁠, and finally,

x_{k} \equiv x_{k} (s) = 2 δ_{k} (1 - e^{i ν_{k} s}) .

(E6)

Substituting Eqs. into the two-time correlation function, we find the form

⟨ C^{†} (s) C (0) ⟩ = C_{1}^{(†, \cdot)} (s) + C_{2}^{(†, \cdot)} (s) + C_{V, 1}^{(†, \cdot)} (s) + C_{V, 0}^{(†, \cdot)} (s),

(E7)

where the subscripts denote the same processes as in Eq. , and

C_{1}^{(†, \cdot)} (s) = κ^{2} e^{ϕ (s)} \sum_{k} f_{k}^{2} μ_{k}^{0} (N_{k} e^{i ν_{k} s} + {\tilde{N}}_{k} e^{- i ν_{k} s}),

(E8a)

C_{2}^{(†, \cdot)} (s) = 4 κ^{2} e^{ϕ (s)} \sum_{n, m \in {- 1, 1}} \sum_{k q} μ_{k}^{0 *} μ_{q}^{0} Δ_{k}^{0} Δ_{q}^{0} χ_{n m} (k, q) e^{i n ν_{k} s} e^{i m ν_{q} s},

(E8b)

C_{V, 1}^{(†, \cdot)} (s) = 2 | V | κ^{2} e^{ϕ (s)} \sum_{n \in {- 1, 1}} \sum_{k} Δ_{k}^{0} (e^{- i ϑ_{V}} μ_{k}^{0} + e^{i ϑ_{V}} μ_{k}^{0 *}) ξ_{n} (k) e^{i n ν_{k} s},

(E8c)

C_{V, 0}^{(†, \cdot)} (s) = κ^{2} | V |^{2} (e^{ϕ (s)} - 1),

(E8d)

where we have written

μ_{k} = i μ_{k}^{0} f_{k}, {\bar{μ}}_{k} = i μ_{k}^{0 *} f_{k}

⁠, and

Δ_{k} = i Δ_{k}^{0} f_{k}

with

μ_{k}^{0} = d_{μ} \cdot e_{k}

and

Δ_{k}^{0} = d_{Δ} \cdot e_{k}

⁠. The values of

χ_{n m} (k, q)

and

ξ_{n} (k)

are written in matrix notation for brevity, using the format

χ (k, q) = (\begin{matrix} χ_{- 1, - 1} (k, q) & χ_{- 1, 1} (k, q) \\ χ_{1, - 1} (k, q) & χ_{1, 1} (k, q) \end{matrix}),

(E9)

ξ (k) = {ξ_{- 1} (k), ξ_{1} (k)}

(E10)

with the elements

χ (k, q) = \frac{f_{k}^{2}}{ν_{k}} \frac{f_{q}^{2}}{ν_{q}} (\begin{matrix} {\tilde{N}}_{k} {\tilde{N}}_{q} & - {\tilde{N}}_{k} N_{q} \\ - N_{k} {\tilde{N}}_{q} & N_{k} N_{q} \end{matrix}),

(E11a)

ξ (k) = \frac{f_{k}^{2}}{ν_{k}} {{\tilde{N}}_{k}, - N_{k}} .

(E11b)

We now take the continuum limit of the photon wavenumber [see Eq. ] to obtain

C_{1}^{(†, \cdot)} (s) = κ^{2} e^{ϕ (s)} Ω_{μ μ} \int_{0}^{\infty} d ν J (ν) (N (ν) e^{i ν s} + \tilde{N} (ν) e^{- i ν s}),

(E12a)

C_{2}^{(†, \cdot)} (s) = 4 κ^{2} e^{ϕ (s)} Ω_{μ Δ}^{2} {cos}^{2} (θ_{μ Δ}) \sum_{n, m \in {- 1, 1}} \int_{0}^{\infty} d ν \int_{0}^{\infty} d ν' χ_{n m} (ν, ν') e^{in ν s} e^{i m ν' s},

(E12b)

C_{V, 1}^{(†, \cdot)} (s) = 4 | V | κ^{2} e^{ϕ (s)} Ω_{μ Δ} cos (θ_{μ Δ}) cos (ϑ_{μ V}) \int_{0}^{\infty} d ν \bar{J} (ν) (\tilde{N} (ν) e^{- i ν s} - N (ν) e^{i ν s}),

(E12c)

C_{V, 0}^{(†, \cdot)} (s) = κ^{2} | V |^{2} (e^{ϕ (s)} - 1),

(E12d)

where in the continuum limit,

χ_{n m} (ν, ν')

are given in Eq. . The Fourier transform of Eq. leads to

Γ_{1} (ω)

in the main text, of Eq. to

Γ_{2} (ω)

⁠, of Eq. to

Γ_{V, 1} (ω)

⁠, and of Eq. to

Γ_{V, 0} (ω)

⁠. Lastly, we note that upon expansion of the summation in Eq. , one obtains

C_{2}^{(†, \cdot)} (s) = {[2 κ e^{\frac{1}{2} ϕ (s)} Ω_{μ Δ} cos (θ_{μ Δ}) \int_{0}^{\infty} d ν \bar{J} (ν) (\tilde{N} (ν) e^{- i ν s} - N (ν) e^{i ν s})]}^{2} .

(E13)

b. $⟨ C (s) C^{†} (0) ⟩$

To calculate this two-time correlation function, we require

(F14)

⟨ B_{+} (s) B_{-} (0) ⟩ = κ^{2} e^{ϕ},

(E14a)

⟨ B_{+} (s) a_{k}^{†} B_{-} (0) ⟩ = (- 2 δ_{k}^{*} - {\tilde{N}}_{k} y_{k}^{*}) κ^{2} e^{ϕ},

(E14b)

⟨ B_{+} (s) a_{k} B_{-} (0) ⟩ = (- 2 δ_{k} + N_{k} y_{k}) κ^{2} e^{ϕ},

(E14c)

⟨ B_{+} (s) a_{k}^{†} a_{q}^{†} B_{-} (0) ⟩ = (- 2 δ_{k}^{*} - {\tilde{N}}_{k} y_{k}^{*}) (- 2 δ_{q}^{*} - {\tilde{N}}_{q} y_{q}^{*}) κ^{2} e^{ϕ},

(E14d)

⟨ B_{+} (s) a_{k} a_{q} B_{-} (0) ⟩ = (- 2 δ_{k} + N_{k} y_{k}) (- 2 δ_{q} + N_{q} y_{q}) κ^{2} e^{ϕ},

(E14e)

⟨ B_{+} (s) a_{k}^{†} a_{q} B_{-} (0) ⟩ = [N_{k} δ_{k q} + (- 2 δ_{k}^{*} - {\tilde{N}}_{k} y_{k}^{*}) (- 2 δ_{q} + N_{q} y_{q})] κ^{2} e^{ϕ},

(E14f)

⟨ B_{+} (s) a_{k} a_{q}^{†} B_{-} (0) ⟩ = [{\tilde{N}}_{k} δ_{k q} + (- 2 δ_{k} + N_{k} y_{k}) (- 2 δ_{q}^{*} - {\tilde{N}}_{q} y_{q}^{*})] κ^{2} e^{ϕ},

(E14g)

where

y_{k} \equiv y_{k} (s) = 2 δ_{k} (e^{i ν_{k} s} - 1) .

(E15)

Substituting Eqs. into the two-time correlation function, one finds the general form

⟨ C (s) C^{†} (0) ⟩ = C_{1}^{(\cdot, †)} (s) + C_{2}^{(\cdot, †)} (s) + C_{V, 1}^{(\cdot, †)} (s) + C_{V, 0}^{(\cdot, †)} (s),

(E16)

where,

C_{1}^{(\cdot, †)} (s) = C_{1}^{(†, \cdot)} (s), C_{2}^{(\cdot, †)} (s) = C_{2}^{(†, \cdot)} (s), C_{V, 1}^{(\cdot, †)} (s) = - C_{V, 1}^{(†, \cdot)} (s)

⁠, and

C_{V, 0}^{(\cdot, †)} (s) = C_{V, 0}^{(†, \cdot)} (s)

⁠.

c. $⟨ C^{†} (s) C^{†} (0) ⟩$

To calculate this two-time correlation function, we require

(F17)

⟨ B_{-} (s) B_{-} (0) ⟩ = κ^{2} e^{- ϕ},

(E17a)

⟨ a_{k}^{†} B_{-} (s) B_{-} (0) ⟩ = - N_{k} z_{k}^{*} κ^{2} e^{- ϕ},

(E17b)

⟨ a_{k} B_{-} (s) B_{-} (0) ⟩ = {\tilde{N}}_{k} z_{k} κ^{2} e^{- ϕ},

(E17c)

⟨ B_{-} (s) a_{k}^{†} B_{-} (0) ⟩ = (- {\tilde{N}}_{k} z_{k}^{*} - 2 δ_{k}^{*}) κ^{2} e^{- ϕ},

(E17d)

⟨ B_{-} (s) a_{k} B_{-} (0) ⟩ = (N_{k} z_{k} - 2 δ_{k}) κ^{2} e^{- ϕ},

(E17e)

⟨ a_{k}^{†} B_{-} (s) a_{q}^{†} B_{-} (0) ⟩ = - N_{k} z_{k}^{*} (- {\tilde{N}}_{q} z_{q}^{*} - 2 δ_{q}^{*}) κ^{2} e^{- ϕ},

(E17f)

⟨ a_{k} B_{-} (s) a_{q} B_{-} (0) ⟩ = {\tilde{N}}_{k} z_{k} (N_{q} z_{q} - 2 δ_{q}) κ^{2} e^{- ϕ},

(E17g)

⟨ a_{k}^{†} B_{-} (s) a_{q} B_{-} (0) ⟩ = [N_{k} δ_{k q} - N_{k} z_{k} (N_{q} z_{q} - 2 δ_{q})] κ^{2} e^{- ϕ},

(E17h)

⟨ a_{k} B_{-} (s) a_{q}^{†} B_{-} (0) ⟩ = [{\tilde{N}}_{k} δ_{k q} + {\tilde{N}}_{k} z_{k} (N_{q} z_{q} - 2 δ_{q})] κ^{2} e^{- ϕ},

(E17i)

where

z_{k} \equiv z_{k} (s) = - 2 δ_{k} (e^{i ν_{k} s} + 1) .

(E18)

Substituting Eqs. into the two-time correlation function and taking the continuum limit, one finds the general form

⟨ C^{†} (s) C^{†} (0) ⟩ = C_{1}^{(†, †)} (s) + C_{2}^{(†, †)} (s) + C_{V, 0}^{(†, †)} (s),

(E19)

which does not have a driving-induced one-photon process. The components of the correlation functions are the same as for

⟨ C^{†} (s) C (0) ⟩

except for overall phase and the replacement

ϕ (s) \to - ϕ (s)

⁠. Explicitly,

C_{1}^{(†, †)} (s) = κ^{2} e^{- ϕ (s)} Ω_{μ μ} e^{2 i ϑ_{μ V}} \int_{0}^{\infty} d ν J (ν) (N (ν) e^{i ν s} + \tilde{N} (ν) e^{- i ν s}),

(E20a)

C_{2}^{(†, †)} (s) = 4 κ^{2} e^{- ϕ (s)} Ω_{μ Δ}^{2} {cos}^{2} (θ_{μ Δ}) e^{2 i ϑ_{μ V}} \sum_{n, m \in {- 1, 1}} \int_{0}^{\infty} d ν \int_{0}^{\infty} d ν' χ_{n m} (ν, ν') e^{in ν s} e^{i m ν' s},

(E20b)

C_{V, 0}^{(†, †)} (s) = κ^{2} | V |^{2} (e^{- ϕ (s)} - 1) .

(E20c)

The Fourier transform of Eq. will depend on the function,

K_{-} (ε) = \frac{1}{π} \int_{0}^{\infty} d s e^{i ε s} e^{- ϕ (s) - ϕ (0)},

(E21)

which is equivalent to Eq. except that

ϕ (s) \to - ϕ (s)

⁠. The finite-mode truncation method solution to Eq. is the same as for Eq. but with the replacement

W_{n} = S^{n} exp (- S) / n! \to {(- 1)}^{n} S^{n} exp (- S) / n!

⁠.

d. $⟨ C (s) C (0) ⟩$

To calculate this two-time correlation function, we require

(E22)

⟨ B_{+} (s) B_{+} (0) ⟩ = κ^{2} e^{- ϕ},

(E22a)

⟨ B_{+} (s) a_{k}^{†} B_{+} (0) ⟩ = (z_{k}^{*} {\tilde{N}}_{k} + 2 δ_{k}^{*}) κ^{2} e^{- ϕ},

(E22b)

⟨ B_{+} (s) a_{k} B_{+} (0) ⟩ = (- z_{k} N_{k} + 2 δ_{k}) κ^{2} e^{- ϕ},

(E22c)

⟨ B_{+} (s) B_{+} (0) a_{k}^{†} ⟩ = {\tilde{N}}_{k} z_{k}^{*} κ^{2} e^{- ϕ},

(E22d)

⟨ B_{+} (s) B_{+} (0) a_{k} ⟩ = - N_{k} z_{k} κ^{2} e^{- ϕ},

(E22e)

⟨ B_{+} (s) a_{k}^{†} B_{+} (0) a_{q}^{†} ⟩ = {\tilde{N}}_{k} z_{k}^{*} ({\tilde{N}}_{q} z_{q}^{*} + 2 δ_{q}^{*}) κ^{2} e^{- ϕ},

(E22f)

⟨ B_{+} (s) a_{k} B_{+} (0) a_{q} ⟩ = N_{k} z_{k} (N_{q} z_{q} - 2 δ_{q}) κ^{2} e^{- ϕ},

(E22g)

⟨ B_{+} (s) a_{k}^{†} B_{+} (0) a_{q} ⟩ = [N_{k} δ_{k q} - (z_{k}^{*} {\tilde{N}}_{k} + 2 δ_{k}^{*}) N_{q} z_{q}] κ^{2} e^{- ϕ},

(E22h)

⟨ B_{+} (s) a_{k} B_{+} (0) a_{q}^{†} ⟩ = [{\tilde{N}}_{k} δ_{k q} + (- z_{k} N_{k} + 2 δ_{k}) {\tilde{N}}_{q} z_{q}^{*}] κ^{2} e^{- ϕ},

(E22i)

where

z_{k} \equiv z_{k} (s)

is given in Eq. . Substituting Eqs. into the two-time correlation function and taking the continuum limit, one finds that

⟨ C (s) C (0) ⟩

is equal to

⟨ C^{†} (s) C^{†} (0) ⟩

with the replacement

ϑ_{μ V} \to - ϑ_{μ V}

⁠.

2. Two-time correlation functions $⟨ g_{α}^{†} (s) g_{β} (0) ⟩$

With the expressions for

⟨ C^{†} (s) C (0) ⟩, ⟨ C (s) C^{†} (0) ⟩, ⟨ C^{†} (s) C^{†} (0) ⟩

⁠, and

⟨ C (s) C (0) ⟩

⁠, we can now write down the two-time correlation functions

⟨ g_{α}^{†} (s) g_{β} (0) ⟩

for

α, β \in {z, +, -}

⁠. We will write these expressions in the continuum limit using the generic operators

g_{α} = a_{α} C + b_{α} C^{†},

(E23)

and the equivalent for

α \to β

⁠. One can recover the desired two-time correlation functions of

g_{-}, g_{+}

⁠, and g_z by using the coefficients written in Table II. The Fourier transforms of the two-time correlation functions are

Γ_{α β} (ω) \equiv \int_{0}^{\infty} d s e^{i ω s} ⟨ g_{α}^{†} (s) g_{β} (0) ⟩ = a_{α} a_{β} Γ^{(†, \cdot)} (ω) + b_{α} b_{β} Γ^{(\cdot, †)} (ω) + a_{α} b_{β} Γ^{(†, †)} (ω) + b_{α} a_{β} Γ^{(\cdot, \cdot)} (ω),

(E24)

where, for example,

Γ^{(†, \cdot)} (ω) = \int_{0}^{\infty} d s e^{i ω s} ⟨ C^{†} (s) C (0) ⟩ .

(E25)

In the main text, we demonstrate how to evaluate

Γ^{(†, \cdot)} (ω)

[there denoted

Γ (ω)

] and the remaining Fourier transforms follow similarly, but with the addition of the various minus signs we have already mentioned.

Table II.

The coefficients used in $g_{α}$ in Eq. (E23) to obtain the polaron frame coupling operators $g_{-}, g_{+}$ ⁠, and g_z. Recall that ${cos}^{2} (φ / 2) = (1 + ϵ / η) / 2$ and ${sin}^{2} (φ / 2) = (1 - ϵ / η) / 2$ ⁠.

α	−	+	z
$a_{α}$	$- {sin}^{2} (\frac{φ}{2})$	${cos}^{2} (\frac{φ}{2})$	${sin}^{2} (\frac{φ}{2}) {cos}^{2} (\frac{φ}{2})$
$b_{α}$	${cos}^{2} (\frac{φ}{2})$	$- {sin}^{2} (\frac{φ}{2})$	${sin}^{2} (\frac{φ}{2}) {cos}^{2} (\frac{φ}{2})$

APPENDIX F: CORRELATION FUNCTIONS IN THE DISPLACED FRAME MASTER EQUATION

If one does not make the polaron transformation of the Hamiltonian in Eq. and instead moves straight to the eigenbasis, one finds the displaced frame Hamiltonian

H = \frac{η^{d}}{2} τ_{z}^{d} + \sum_{k} ν_{k} a_{k}^{†} a_{k} + \sum_{μ \in {z, +, -}} g_{μ}^{d} τ_{μ}^{d},

(F1)

where

η^{d} = \sqrt{ϵ^{2} + 4 | V |^{2}}

[there is no κ renormalization of V] and

g_{z}^{d} = [{cos}^{2} (\frac{φ^{d}}{2}) - {sin}^{2} (\frac{φ^{d}}{2})] π_{Δ Δ} + cos (\frac{φ^{d}}{2}) sin (\frac{φ^{d}}{2}) [π_{μ \bar{μ}} e^{- i ϑ_{V}} + π_{\bar{μ} μ} e^{i ϑ_{V}}],

(F2a)

g_{+}^{d} = - 2 cos (\frac{φ^{d}}{2}) sin (\frac{φ^{d}}{2}) π_{Δ Δ} + {cos}^{2} (\frac{φ^{d}}{2}) π_{μ \bar{μ}} e^{- i ϑ_{V}} - {sin}^{2} (\frac{φ^{d}}{2}) π_{\bar{μ} μ} e^{i ϑ_{V}},

(F2b)

and

g_{-}^{d} = g_{+}^{d †}

⁠, where

π_{p q} = \sum_{k} (p_{k} a_{k}^{†} + q_{k}^{*} a_{k}) .

(F3)

The Pauli matrices in the eigenbasis are

τ_{\pm}^{d} = | \pm^{d} ⟩ ⟨ \mp^{d} |

and

τ_{z}^{d} = | +^{d} ⟩ ⟨ +^{d} | - | -^{d} ⟩ ⟨ -^{d} |

where

(\begin{matrix} | e ⟩ \\ | g ⟩ \end{matrix}) = (\begin{matrix} cos (\frac{φ^{d}}{2}) & - sin (\frac{φ^{d}}{2}) \\ sin (\frac{φ^{d}}{2}) & cos (\frac{φ^{d}}{2}) \end{matrix}) (\begin{matrix} | +^{d} ⟩ \\ | -^{d} ⟩ \end{matrix})

(F4)

with

cos (φ^{d}) = ϵ / η^{d}

and

sin (φ^{d}) = 2 | V | / η^{d}

⁠.

Since Eq. has the same structure as Eq. , the master equation in the displaced frame has the same algebraic form as in the polaron frame, i.e., given within the secular approximation by Eqs. and in full in Appendix C, but the environment correlation functions (ECFs) are different. The displaced frame ECFs depend on linear combinations of Fourier transforms of the form

Γ_{p q, r s}^{d} (ω) = \int_{0}^{\infty} d s e^{i ω s} ⟨ π_{p q}^{†} (s) π_{r s} (0) ⟩ \equiv \frac{1}{2} γ_{p q, r s}^{d} (ω) + i S_{p q, r s}^{d} (ω) .

(F5)

Substituting Eq. into the ECF and using that

⟨ a_{k}^{†} a_{k'} ⟩ = δ_{k k'} N (ν_{k}), ⟨ a_{k} a_{k'}^{†} ⟩ = δ_{k k'} \tilde{N} (ν_{k})

⁠, where

\tilde{N} (ν) = 1 + N (ν)

and that other combinations equal zero, leads to

(G6)

γ_{p q, r s}^{d} (ω) = 2 π [cos (θ_{p r}) e^{i (ϑ_{r} - ϑ_{p})} Ω_{p r} J (ω) \tilde{N} (ω) + cos (θ_{q s}) e^{i (ϑ_{q} - ϑ_{s})} Ω_{q s} J (- ω) N (- ω)],

(F6a)

S_{p q, r s}^{d} (ω) = P \int_{0}^{\infty} d ν J (ν) [cos (θ_{p r}) e^{i (ϑ_{r} - ϑ_{p})} Ω_{p r} \frac{\tilde{N} (ν)}{ω - ν} + cos (θ_{q s}) e^{i (ϑ_{q} - ϑ_{s})} Ω_{q s} \frac{N (ν)}{ω + ν}] .

(F6b)

To exemplify how these functions relate to the rates in Eqs. , we will derive

γ_{- -}^{d} (η^{d})

explicitly. This function is given by

γ_{- -}^{d} (η^{d}) = 2 ℜ \int_{0}^{\infty} d s e^{i η^{d} s} ⟨ g_{-}^{d †} (s) g_{-}^{d} (0) ⟩ .

(F7)

Substituting in

g_{-}^{d}

and using Eqs. , we find that [ignoring arguments

(η^{d})

on the right-hand side and temporarily denoting

cos (φ^{d} / 2) = c

and

sin (φ^{d} / 2) = s

],

\begin{matrix} γ_{- -}^{d} (η^{d}) = 4 c^{2} s^{2} γ_{Δ Δ, Δ Δ}^{d} + 2 c s^{3} e^{- i ϑ_{μ V}} γ_{Δ Δ, μ μ}^{d} - 2 c^{3} s e^{i ϑ_{μ V}} γ_{Δ Δ, μ μ}^{d} \\ + 2 c s^{3} e^{i ϑ_{μ V}} γ_{μ μ, Δ Δ}^{d} + s^{4} γ_{μ μ, μ μ}^{d} - s^{2} c^{2} e^{2 i ϑ_{μ V}} γ_{μ μ, μ μ}^{d} \\ - 2 c^{3} s e^{- i ϑ_{μ V}} γ_{μ μ, Δ Δ}^{d} - c^{2} s^{2} e^{- 2 i ϑ_{μ V}} γ_{μ μ, μ μ}^{d} + c^{4} γ_{μ μ, μ μ}^{d} . \end{matrix}

(F8)

Using the symmetry that

Ω_{a b} = Ω_{b a}

and that

d_{Δ} \in ℜ

⁠, we can read off from Eq. that

γ_{Δ Δ, Δ Δ}^{d} (η^{d}) = Ω_{Δ Δ} γ_{0} (η^{d}),

(F9)

γ_{μ μ, μ μ}^{d} (η^{d}) = Ω_{μ μ} γ_{0} (η^{d}),

(F10)

γ_{μ μ, Δ Δ}^{d} (η^{d}) = γ_{Δ Δ, μ μ}^{d} (η^{d}) = cos (θ_{μ Δ}) Ω_{μ Δ} γ_{0} (η^{d}),

(F11)

where

γ_{0} (ω) = 2 π [J (ω) \tilde{N} (ω) + J (- ω) N (- ω)

]. Using these, we find that

γ_{- -}^{d} (η^{d}) = k_{- -}^{0} γ_{0} (η^{d}),

(F12)

where the rate coefficient is

k_{- -}^{0} = [c^{4} + s^{4}] Ω_{μ μ} + 4 c^{2} s^{2} Ω_{Δ Δ} - s^{2} c^{2} [Ω_{μ μ} e^{- 2 i ϑ_{μ V}} + Ω_{μ μ} e^{2 i ϑ_{μ V}}] - 2 s c [c^{2} - s^{2}] [Ω_{μ Δ} e^{- i ϑ_{μ V}} + Ω_{μ Δ} e^{i ϑ_{μ V}}] cos (θ_{μ Δ}) .

(F13)

This is rather complicated but only because of the transition dipole and driving phases. If we assume that

ϑ_{μ} = ϑ_{V}

such that

ϑ_{μ V} = 0

⁠, then

k_{- -}^{0} = (8 π / 3) (d_{-} \cdot d_{-})

where we have defined a new dipole vector

d_{-} = [c^{2} - s^{2}] d_{μ} + 2 c s d_{Δ} .

(F14)

A final limit worth checking is when the eigenstates fully localize, $| +^{d} ⟩ \to | e ⟩$ and $| -^{d} ⟩ \to | g ⟩$ ⁠, i.e., $c \to 1$ and $s \to 0$ ⁠. In this case, $k_{- -}^{0} \to Ω_{μ μ}$ ⁠, and so the decay rate is $γ_{- -}^{d} (η^{d}) = γ_{- -}^{d} (ϵ) = Ω_{μ μ} γ_{0} (ϵ)$ ⁠, which is the decay rate in the standard optical master equation.

APPENDIX G: DERIVATION OF THE EMISSION SPECTRUM

In this appendix, we derive Eq. (51), which gives the emission spectrum for the single emitter system. Our derivation closely follows those provided for single emitter systems without permanent dipoles in Refs. 2, 64–66, and 68. Due to the permanent dipoles, there are additional light–matter interaction terms in the Hamiltonian in Eq. (5), and so, the derivation is slightly more cumbersome. However, as we show here, these additional terms do not affect the expression for the emission spectrum of the emitter within the standard approximations.

The emission spectrum is given exactly by

I (ω) = lim_{t \to \infty} ℜ \int_{0}^{\infty} d τ e^{- i ω τ} ⟨ E_{-} (R, t + τ) \cdot E_{+} (R, t) ⟩,

(G1)

where R is the position of the detector and the positive frequency component of the electric field is

E_{+} (R, t) = i \sum_{k} e_{k} f_{k} a_{k} e^{- i k \cdot R},

(G2)

and

E_{-} (R, t) = E_{+} {(R, t)}^{†}

⁠. Our aim is to express the expectation value

⟨ a_{k}^{†} a_{q} ⟩

in terms of dipole operators

σ_{\pm}

⁠, which is achieved through the Heisenberg equation of motion

\frac{\partial}{\partial t} a_{k} (t) = - i [H, a_{k} (t)],

(G3)

where H is the lab frame Hamiltonian in Eq. . By obtaining

⟨ a_{k}^{†} a_{q} ⟩

in terms of

σ_{\pm}

⁠, we will have arrived at the form of the emission spectrum in Eq. .

Using

[a_{k}, a_{q}^{†}] = δ_{k q}

⁠, one can show that

\frac{\partial}{\partial t} a_{k} (t) = - i ν_{k} a_{k} (t) - Δ_{k} σ_{z} (t) - μ_{k} σ_{+} (t) - μ_{k}^{*} σ_{-} (t),

(G4)

where

σ_{α} (t) = U {(t)}^{†} σ_{α} U (t)

and

U (t) = exp (- iHt)

⁠, which can be solved to yield

a_{k} (t) = a_{k} (0) e^{- i ν_{k} t} - \int_{0}^{t} d τ e^{i ν_{k} (τ - t)} [Δ_{k} σ_{z} (τ) + μ_{k} σ_{+} (τ) - μ_{k}^{*} σ_{-} (τ)] .

(G5)

The first term in Eq. is the free field term, which does not contribute to the spectrum of the emitter and is henceforth ignored. After substituting the second term of Eq. into Eq. and taking the continuum limit, we obtain

E_{+} (R, t) = \int_{0}^{\infty} d ν \sqrt{J (ν)} \int_{0}^{t} d τ e^{i ν_{k} (τ - t)} [O_{z} σ_{z} (τ) + O_{μ} σ_{+} (τ) + O_{\bar{μ}} σ_{-} (τ)],

(G6)

where

O_{p} = \int_{d Ω_{k}} d Ω_{k} \sum_{λ} e_{k} (d_{p} \cdot e_{k}) e^{i k \cdot (r - R)},

(G7)

where we briefly use the notation

d_{z} = d_{Δ}

for convenience.

To progress analytically, we need to know how the emitter operators

σ_{α} (t)

evolve in time. However, this is very complicated. Instead, we make the so-called harmonic decomposition (see Chap. 2.2 of Ref. 2) in which we assume that the timescale over which the emitter evolves unitarily is much faster than the timescale over which spontaneous emission occurs. Within this approximation, we write that

σ_{+} (τ) \approx σ_{+} (t) e^{i ϵ (τ - t)}, σ_{-} (τ) \approx σ_{-} (t) e^{- i ϵ (τ - t)}, and σ_{z} (τ) \approx σ_{z} (t) .

(G8)

Making the harmonic decomposition in Eq. yields

E_{+} (R, t) = \int_{0}^{\infty} d ν \sqrt{J (ν)} [O_{z} σ_{z} (t) j_{0} (ν, t) + O_{μ} σ_{+} (t) j_{+} (ν, t) + O_{\bar{μ}} σ_{-} (t) j_{-} (ν, t)],

(G9)

where

j_{α} (ν, t) = \int_{0}^{t} d τ e^{i (τ - t) (ν + α ϵ)} = - i \frac{1 - e^{- i t (ν + α ϵ)}}{ν + α ϵ} .

(G10)

The function

j_{α} (ν, t)

is dominated by the contribution near to

ν = - α ϵ

and so we approximate it as a delta function,

j_{α} (ν, t) \approx 2 π δ (ν + α ϵ)

⁠.² After performing this approximation, one finds that the delta function corresponding to

j_{+} (ν, t)

lies out with the integration domain

ν \in [0, \infty]

⁠, and so, the term going as

σ_{+} (t)

in Eq. vanishes. Moreover, the term proportional to

j_{0} (ν, t)

provides a delta function at zero frequency, leading to the term going as

σ_{z} (t)

in Eq. to be proportional to

\sqrt{J (0)}

⁠, which is equal to zero for any well-defined spectral density. Therefore, within the standard approximations outlined in this derivation, the presence of permanent dipoles does not change the expression of the emission spectrum from Eq. . Thus, the only surviving term in Eq. is proportional to

σ_{-} (t)

⁠, and so, we obtain

E_{+} (R, t) \approx 2 π \sqrt{J (ϵ)} O_{\bar{μ}} (r, R) σ_{-} (t),

(G11)

where we have made the dependence of

O_{\bar{μ}}

on r and R explicit. After substituting Eq. and its Hermitian conjugate into Eq. , one obtains Eq. , and an explicit expression for the Green's function

α_{prop} (r, R, ω)

⁠.

APPENDIX H: NON-SECULAR MASTER EQUATION PRODUCES A NONPHYSICAL POLARIZATION SPECTRUM

In this appendix, we show that if a non-secular master equation is used, then the standard optical master equation can produce a nonphysical (negative) polarization spectrum. We will show that, for an isolated transition dipole without external drive, the negative frequency peak occurs with a negative magnitude for the non-secular master equation, proportional to the square of the Lamb shift.

The Hamiltonian for the standard optical master equation is

H = \frac{ϵ}{2} σ_{z} + \sum_{k} ν_{k} a_{k}^{†} a_{k} + σ_{x} \sum_{k} (g_{k} a_{k}^{†} + g_{k}^{*} a_{k}),

(H1)

where

g_{k} = i f_{k} (d_{μ} \cdot e_{k})

and other symbols are defined in the main text. The non-secular master equation (the standard optical master equation) is

{\dot{ρ}}_{e e} (t) = - γ_{↓} ρ_{e e} (t) + γ_{↑} ρ_{g g} (t),

(H2)

{\dot{ρ}}_{g e} (t) = (i ω_{0} - γ_{d}) ρ_{g e} + γ_{χ} ρ_{e g} (t),

(H3)

and

{\dot{ρ}}_{g g} (t) = - {\dot{ρ}}_{e e} (t)

and

{\dot{ρ}}_{e g} (t) = {\dot{ρ}}_{g e} {(t)}^{†}

⁠. The transition and decoherence rates are

γ_{↑} = 2 π Ω_{μ μ} J (ϵ) N (ϵ),

(H4)

γ_{↓} = 2 π Ω_{μ μ} J (ϵ) (1 + N (ϵ)),

(H5)

γ_{d} = \frac{1}{2} (γ_{↑} + γ_{↓}),

(H6)

and the Lamb shifted transition energy is

ω_{0} = ϵ + S_{+} - S_{-},

(H7)

where

S_{\pm} = S (\pm ϵ)

and

S (ω) = P \int_{0}^{\infty} d ν J (ν) (\frac{N (ν)}{ω + ν} + \frac{1 + N (ν)}{ω - ν}) .

(H8)

Finally, the non-secular rate is

γ_{χ} = γ_{d} + i (S_{+} - S_{-}) .

(H9)

To recover the secular theory, one just sets

γ_{χ} = 0

⁠. Note that for this system,

γ_{χ} = 0

is also obtained if one makes the rotating wave approximation in the interaction Hamiltonian.

Using the quantum regression theorem, we find that the polarization spectrum is the real part of the Fourier transform of

ρ_{g e} (t)

with the initial conditions

ρ_{g e} (0) = 0

and

ρ_{g e} (0) = ρ_{e e} (\infty) = N (ϵ) / (1 + 2 N (ϵ))

⁠. The solution to the non-secular master equation is

ρ_{g e} (t) = e^{- t γ_{d}} ρ_{e e} (\infty) (cosh (t \sqrt{ξ}) + \frac{i ω_{0} sinh (t \sqrt{ξ})}{\sqrt{ξ}}),

(H10)

where

ξ^{2} = | γ_{χ} |^{2} - ω_{0}^{2}

⁠. Notice that, as expected, if

γ_{χ} = 0

we recover the secular result that

ρ_{g e} (t) \propto exp [(i ω_{0} - γ_{d}) t]

⁠. The polarization spectrum is then

I_{p, SOME} (ω) = ℜ \int_{0}^{\infty} d t e^{- i ω t} ρ_{g e} (t) = ρ_{e e} (\infty) ℜ [\frac{γ_{d} + i (ω + ω_{0})}{{(γ_{d} + i ω)}^{2} - ξ^{2}}],

(H11)

which has the secular limit

lim_{γ_{χ} \to 0} I_{p, SOME} (ω) = ρ_{e e} (\infty) \frac{γ_{d}}{γ_{d}^{2} + (ω - ω_{0}^{2})}

(H12)

that is always positive valued. Conversely, the value of the non-secular spectrum at the negative frequency peak,

ω = - ω_{0}

⁠, is

I_{p, SOME} (- ω_{0}) = - ρ_{e e} (\infty) \frac{γ_{d} L^{2}}{L^{4} + 4 γ_{d}^{2} ω_{0}^{2}} \sim - \frac{ρ_{e e} (\infty)}{4 γ_{d} ϵ^{2}} L^{2} + O (L^{3}),

(H13)

where

L = S_{+} - S_{-}

is the Lamb shift. Equation shows that the non-secular spectrum is nonphysical at the negative peak if the Lamb shift is finite, while the secular spectrum in Eq. is always physical. This effect is small for weak light–matter coupling (⁠

L / ϵ ≪ 1

⁠) and low temperatures (⁠

β ϵ ≪ 1

⁠) but is nevertheless present. However, within this regime, the secular and non-secular spectra are otherwise very similar and only start to deviate very close to

ω = - ω_{0}

⁠. Therefore, in the main text where our results are well within this regime, we use the secular master equation to produce the spectra. By doing this, we avoid a narrow and negative peak at

ω = - ω_{0}

when the driving is small, but otherwise leave the spectra unchanged.

APPENDIX I: EFFECTIVE HAMILTONIAN

In many numerical schemes to solve the open quantum dynamics of systems coupled to thermal environments, it is assumed that the environment is in a free Gibbs state at the temperature $T = 1 / β$ ⁠. In the polaron framework, the lab frame environment is not in such a convenient form. As such, we derive here an effective Hamiltonian that encapsulates the polaron thermalized state but has the effective environment in a free Gibbs state.

We start with the lab frame Hamiltonian,

H = \frac{ϵ}{2} σ_{z} + V σ_{+} + V^{*} σ_{-} + E_{dip} + \sum_{k} ν_{k} a_{k}^{†} a_{k} + π_{D D} I + π_{Δ Δ} σ_{z} + π_{μ \bar{μ}} σ_{+} + π_{\bar{μ} μ} σ_{-},

(I1)

where

π_{p q} = \sum_{k} (p_{k} a_{k}^{†} + q_{k}^{*} a_{k}),

(I2)

p_{k} = i (d_{p} \cdot e_{k}),

(I3)

with

p, q \in {μ, \bar{μ}, D, Δ}

⁠. In our calculations, the polaron frame initial state is

ρ_{p}^{0} = | g ⟩ ⟨ g | \otimes ρ_{E},

(I4)

where

ρ_{E} = exp [- β \sum_{k} ν_{k} a_{k}^{†} a_{k}] / Z_{E}

is a thermal state. Performing the inversion of the unitary transformations to go from the polaron frame to the displaced frame, the initial state is

ρ_{d}^{0} = U^{†} ρ_{p}^{0} U

where the polaron transformation is

U = B (δ) | e ⟩ ⟨ e | + B (- δ) | g ⟩ ⟨ g |

and

δ_{k} = Δ_{k} / ν_{k}

and

B (α) = exp [\sum_{k} (α_{k} a_{k}^{†} - α_{k}^{*} a_{k})]

⁠. Therefore,

ρ_{d}^{0} = | g ⟩ ⟨ g | \otimes η (δ),

(I5)

where we have defined a displaced thermal state as

η (α) = B (α) ρ_{E} B (- α)

⁠. We can then obtain the lab frame initial state via

ρ_{l}^{0} = B (- d) ρ_{d}^{0} B (d)

⁠, where

d_{k} = D_{k} / ν_{k}

⁠, leading to

ρ_{l}^{0} = | g ⟩ ⟨ g | \otimes η (δ - d) .

(I6)

The expression for

η (α)

can be rewritten by making use of the identity

exp (e^{S} X e^{- S}) = e^{S} e^{X} e^{- S},

(I7)

which holds if

e^{S} e^{- S} = I

⁠. Using this identity yields

η (α) = exp [- β \sum_{k} ν_{k} (a_{k}^{†} - α_{k}^{*}) (a_{k} - α_{k})] .

(I8)

We now rewrite the lab frame Hamiltonian using new ladder operators b_k, which we will relate to the a_k to ensure we model our desired initial state,

η (δ - d)

⁠. The initial environment state in the effective lab frame will be

{\bar{ρ}}_{E} = exp [- β \sum_{k} ν_{k} b_{k}^{†} b_{k}] / Z_{E}

⁠, and so by comparison with Eq. , we know that

b_{k} = a_{k} - (δ_{k} - d_{k}) .

(I9)

Substituting this into our actual lab frame Hamiltonian in Eq. and ignoring any terms that are identity operators in both Hilbert spaces, we arrive at the effective Hamiltonian

\bar{H} = {\bar{H}}_{0} + {\bar{H}}_{I},

(I10)

where

{\bar{H}}_{0} = (\frac{ϵ}{2} + G_{Δ Δ}) σ_{z} + V σ_{+} + V^{*} σ_{-} + G_{μ \bar{μ}} σ_{+} + G_{\bar{μ} μ} σ_{-} + \sum_{k} ν_{k} b_{k}^{†} b_{k},

(I11)

{\bar{H}}_{I} = π'_{Δ Δ} (I + σ^{z}) + π'_{μ \bar{μ}} σ_{+} + π'_{\bar{μ} μ} σ_{-},

(I12)

and

π'_{p q} = \sum_{k} (p_{k} b_{k}^{†} + q_{k}^{*} b_{k}),

(I13)

G_{p q} = \sum_{k} (p_{k} δ_{k}^{*} + q_{k}^{*} δ_{k}) .

(I14)

In order to get numerical agreement between TEMPO and the PFME, both of which assume separable initial states but crucially in different frames, we, therefore, must use

\bar{H}

in Eq. for TEMPO calculations.

REFERENCES

1.

H.-P.

Breuer

and

F.

Petruccione

et al,

The Theory of Open Quantum Systems

(

Oxford University Press on Demand

,

2002

).

Google Scholar

2.

Z.

Ficek

and

S.

Swain

,

Quantum Interference and Coherence: Theory and Experiments

(

Springer Science & Business Media

,

2005

), Vol.

100

.

Google Scholar

3.

A. J.

Leggett

,

S.

Chakravarty

,

A. T.

Dorsey

,

M. P.

Fisher

,

A.

Garg

, and

W.

Zwerger

,

Rev. Mod. Phys.

59

,

1

(

1987

).

https://doi.org/10.1103/RevModPhys.59.1

Google Scholar

Crossref

4.

G. S.

Agarwal

, “

Quantum statistical theories of spontaneous emission and their relation to other approaches

,” in

Quantum Optics

(

Springer

,

1974

).

Google Scholar

Crossref

5.

P.-H.

Chung

,

C.

Tregidgo

, and

K.

Suhling

,

Methods Appl. Fluoresc.

4

,

045001

(

2016

).

https://doi.org/10.1088/2050-6120/4/4/045001

Google Scholar

Crossref

PubMed

6.

C.

Filippi

,

F.

Buda

,

L.

Guidoni

, and

A.

Sinicropi

,

J. Chem. Theory Comput.

8

,

112

(

2012

).

https://doi.org/10.1021/ct200704k

Google Scholar

Crossref

PubMed

7.

V.

Kovarskiĭ

and

O.

Prepelitsa

,

Opt. Spectrosc.

90

,

351

(

2001

).

https://doi.org/10.1134/1.1358440

Google Scholar

Crossref

8.

J.

Deiglmayr

,

A.

Grochola

,

M.

Repp

,

O.

Dulieu

,

R.

Wester

, and

M.

Weidemüller

,

Phys. Rev. A

82

,

032503

(

2010

).

https://doi.org/10.1103/PhysRevA.82.032503

Google Scholar

Crossref

9.

R.

Guérout

,

M.

Aymar

, and

O.

Dulieu

,

Phys. Rev. A

82

,

042508

(

2010

).

https://doi.org/10.1103/PhysRevA.82.042508

Google Scholar

Crossref

10.

C.-Y.

Lin

and

S. G.

Boxer

,

J. Am. Chem. Soc.

142

,

11032

(

2020

).

https://doi.org/10.1021/jacs.0c02796

Google Scholar

Crossref

PubMed

11.

B.

Jagatap

and

W. J.

Meath

,

J. Opt. Soc. Am. B

19

,

2673

(

2002

).

https://doi.org/10.1364/JOSAB.19.002673

Google Scholar

Crossref

12.

J.

Gilmore

and

R. H.

McKenzie

,

J. Phys.: Condens. Matter

17

,

1735

(

2005

).

https://doi.org/10.1088/0953-8984/17/10/028

Google Scholar

Crossref

13.

L.

Garziano

,

V.

Macrì

,

R.

Stassi

,

O.

Di Stefano

,

F.

Nori

, and

S.

Savasta

,

Phys. Rev. Lett.

117

,

043601

(

2016

).

https://doi.org/10.1103/PhysRevLett.117.043601

Google Scholar

Crossref

PubMed

14.

I. Y.

Chestnov

,

V.

Shakhnazaryan

,

I. A.

Shelykh

, and

A. P.

Alodjants

,

JETP Lett.

104

,

169

(

2016

).

https://doi.org/10.1134/S0021364016150066

Strong coupling dynamics of driven quantum systems with permanent dipoles

I. INTRODUCTION

II. THE MODEL AND POLARON FRAME

III. REALISTIC PARAMETERS

IV. EFFECTS OF STRONG PERMANENT DIPOLES

A. One photon processes, Γ1(ω)

B. Two-photon processes, Γ2(ω)

C. Driving-induced one-photon processes, ΓV,1(ω)

D. Driving-induced zero-photon processes, ΓV,0(ω)

E. Numerical analysis of the rates

V. COMPARISON TO EXACT DYNAMICS

VI. EMISSION SPECTRUM

VII. IMPORTANCE OF INITIAL CONDITIONS

VIII. CONCLUSION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: DISPLACEMENT OPERATORS

APPENDIX B: CALCULATION OF ⟨ C ⟩

APPENDIX C: THE HUANG–RHYS PARAMETER

APPENDIX D: NON-SECULAR MASTER EQUATIONS

APPENDIX E: CALCULATION OF THE POLARON FRAME TWO-TIME CORRELATION FUNCTIONS

1. Preliminary calculations

a. ⟨ C † ( s ) C ( 0 ) ⟩

b. ⟨ C ( s ) C † ( 0 ) ⟩

c. ⟨ C † ( s ) C † ( 0 ) ⟩

d. ⟨ C ( s ) C ( 0 ) ⟩

2. Two-time correlation functions ⟨ g α † ( s ) g β ( 0 ) ⟩

APPENDIX F: CORRELATION FUNCTIONS IN THE DISPLACED FRAME MASTER EQUATION

APPENDIX G: DERIVATION OF THE EMISSION SPECTRUM

APPENDIX H: NON-SECULAR MASTER EQUATION PRODUCES A NONPHYSICAL POLARIZATION SPECTRUM

APPENDIX I: EFFECTIVE HAMILTONIAN

REFERENCES

A. One photon processes, Γ₁(ω)

B. Two-photon processes, Γ₂(ω)

C. Driving-induced one-photon processes, Γ_V,1(ω)

D. Driving-induced zero-photon processes, Γ_V,0(ω)

APPENDIX B: CALCULATION OF $⟨ C ⟩$

a. $⟨ C^{†} (s) C (0) ⟩$

b. $⟨ C (s) C^{†} (0) ⟩$

c. $⟨ C^{†} (s) C^{†} (0) ⟩$

d. $⟨ C (s) C (0) ⟩$

2. Two-time correlation functions $⟨ g_{α}^{†} (s) g_{β} (0) ⟩$