# Optical polaron formation in quantum systems with permanent dipoles

Adam Burgess<EMAIL_ADDRESS>Leverhulme Quantum Biology Doctoral
Training Centre, University of Surrey, Guildford, GU2 7XH, United Kingdom
Department of Physics and Advanced Technology Institute, University of Surrey,
Guildford, GU2 7XH, United Kingdom  Marian Florescu<EMAIL_ADDRESS>Department of Physics and Advanced Technology Institute, University of Surrey,
Guildford, GU2 7XH, United Kingdom  Dominic M. Rouse
<EMAIL_ADDRESS>School of Physics and Astronomy, University of
Glasgow, Glasgow G12 8QQ, Scotland, United Kingdom  Department of Physics and
Astronomy, University of Manchester, Oxford Road, Manchester M13 9PL, United
Kingdom

###### Abstract

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 polarised Fock states of light. We derive a master
equation for these systems by employing an optical polaron transformation that
captures the photon mode polarisation 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 optimise 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, 2, 3, 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, 6, 7, 8, 9, 10, 11, 12], quantum dots with
asymmetric confining potentials [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 22,
23], nanorods with non-centrosymmetric crystallographic lattices [24, 25, 26],
and superconducting circuits [27].

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 [32], multiphoton conversion [33, 34], entanglement
generation [35, 36], and second-harmonic generation [37, 38]. Understanding
the role of strong permanent dipoles is highly relevant for the design of new
quantum technologies, as well as in exploring novel biochemical processes.
Previous studies have either neglected the transition dipole moments assuming
only a pure dephasing interaction [12], have considered single-mode fields
[29, 39, 40, 41, 34], or have treated the permanent dipoles in a perturbative
manner [12, 28, 35, 29, 39, 40, 41]. 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 utilise a polaron transformation to derive a master equation
for systems with strong permanent dipoles interacting with a thermal 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 hybridisation of the matter excitations and the
displaced harmonic oscillator states of the multi-mode photonic field. These
photonic states, called polarised Fock states, correspond to a non-zero
(polarised) vector potential field and were first introduced to explain how
permanent dipoles can generate multiphoton conversion in a single-mode cavity
[34] and are useful in polaritonic chemistry [42, 43].

The optical polaron formalism provides us with an essential intuition for the
very complex phenomena introduced by the permanent dipoles, allowing us to
unpick the role of permanent dipoles in open quantum system 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 control over the dynamics.
We also derive the emission spectrum for the system, highlighting
experimentally detectable signatures of the optical polarons.

This paper is organised as follows. In Sec. II, we introduce the Hamiltonian
and transform it into the polaron frame. In Sec. III, we derive the polaron
frame master equation and discuss the new physical processes using the
analytical expressions. Following this, in Sec. IV, we compare the polaron
frame master equation to perturbative dynamics and to numerically exact
dynamics using the time-evolving matrix product operators (TEMPO) [44, 45, 46,
47, 48, 49, 50, 51, 52] algorithm, through the open source code OQuPy [53]. In
Sec. V we derive the emission spectrum, and in Sec. VI we briefly discuss the
role of the initial state and identity type interactions. Finally in Sec. VII
we present concluding remarks.

##  II The model and polaron frame

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

$H=\frac{\epsilon}{2}\sigma_{z}+V\sigma_{+}+V^{*}\sigma_{-}+\sum_{k}\nu_{\mathbf{k}}a_{k}^{\dagger}a_{k}+\mathbf{d}\cdot\mathbf{\Pi}+E_{\text{dip}},$
(1)

where $\epsilon$ is the transition energy of the quantum emitter and
$V=|V|\mathrm{e}^{i\vartheta_{V}}$ is a complex drive.
$\mathbf{\Pi}=i\sum_{k}\mathbf{e}_{k}f_{\mathbf{k}}(a_{k}^{\dagger}-a_{k})$ is
the electric displacement field where $k=\\{\mathbf{k},\lambda\\}$ is a four-
vector representing both the wavevector $\mathbf{k}$ and polarisation state
$\lambda$ of the mode with polarisation vector $\mathbf{e}_{k}$, energy
$\nu_{\mathbf{k}}$, and coupling strength
$f_{\mathbf{k}}=\sqrt{\nu_{\mathbf{k}}/2\mathcal{V}}$ where $\mathcal{V}$ is
the field volume. $a_{k}$ and $a^{\dagger}_{k}$ are the field mode
annihilation and creation operators [54, 55, 56]. The dipolar self energy term
is

$E_{\text{dip}}=\sum_{k}\frac{f_{\mathbf{k}}^{2}}{\nu_{\mathbf{k}}}\left(\mathbf{d}\cdot\mathbf{e}_{k}\right)^{2}.$
(2)

In the truncated system Hilbert space, the dipole operator is

$\displaystyle\mathbf{d}$
$\displaystyle=\begin{pmatrix}\mathbf{d}_{ee}&\mathbf{d}_{eg}\\\
\mathbf{d}_{ge}&\mathbf{d}_{gg}\end{pmatrix}=\mathbf{d}_{\Delta}\sigma_{z}+\mathbf{d}_{D}\mathcal{I}+\mathbf{d}_{\mu}\sigma_{+}+\mathbf{d}_{\mu}^{*}\sigma_{-},$
(3)

where the Pauli operators are $\sigma_{z}=\ket{e}\bra{e}-\ket{g}\bra{g}$,
$\sigma_{+}=\ket{e}\bra{g}$, $\sigma_{-}=\ket{g}\bra{e}$,
$\mathcal{I}=\ket{g}\bra{g}+\ket{e}\bra{e}$, and $\mathbf{d}_{ij}=\langle
i|\mathbf{d}|j\rangle$ for $i,j\in\\{e,g\\}$. We have defined the following
combinations of dipole matrix elements,

$\displaystyle\mathbf{d}_{\Delta}=\frac{\mathbf{d}_{ee}-\mathbf{d}_{gg}}{2},\quad\mathbf{d}_{D}=\frac{\mathbf{d}_{ee}+\mathbf{d}_{gg}}{2},\quad\mathbf{d}_{\mu}=\mathbf{d}_{eg},$
(4)

which play an essential role in our analysis. The $\mathbf{d}_{p}$ vectors are
not guaranteed to be co-linear and Hermiticity of $H$ requires that
$\mathbf{d}_{p}\in\mathbb{R}^{3}$ for $p\in\\{\Delta,D\\}$ and
$\mathbf{d}_{\mu}\in\mathbb{C}^{3}$.

At this point it is often assumed that either
$|\mathbf{d}_{\Delta}|\approx|\mathbf{d}_{D}|\approx 0$ which leads to the
standard optical master equation [1, 2, 3, 4], or that
$|\mathbf{d}_{D}|\approx|\mathbf{d}_{\mu}|\approx 0$ leading to a pure
dephasing interaction [12]. 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. (3) into Eq. (1) and absorbing the drive phase into a basis
$\ket{e^{\prime}}=\mathrm{e}^{i\vartheta_{V}/2}\ket{e}$ and
$\ket{g^{\prime}}=\mathrm{e}^{-i\vartheta_{V}}\ket{g}$ we find

$\displaystyle H=$ $\displaystyle\
\frac{\epsilon}{2}\sigma_{z}^{\prime}+|V|\sigma_{x}^{\prime}+\sum_{k}\nu_{\mathbf{k}}a_{k}^{\dagger}a_{k}+E_{\text{dip}}+\pi_{DD}\mathcal{I}^{\prime}$
$\displaystyle+\pi_{\Delta\Delta}\sigma_{z}^{\prime}+\pi_{\mu\bar{\mu}}\mathrm{e}^{-i\vartheta_{V}}\sigma_{+}^{\prime}+\pi_{\bar{\mu}\mu}\mathrm{e}^{i\vartheta_{V}}\sigma_{-}^{\prime},$
(5)

where primed operators are in the $\\{\ket{e^{\prime}},\ket{g^{\prime}}\\}$
basis and

$\displaystyle\pi_{pq}$
$\displaystyle=\sum_{k}\left(p_{k}a_{k}^{\dagger}+q_{k}^{*}a_{k}\right),$ (6)
$\displaystyle p_{k}$
$\displaystyle=if_{\mathbf{k}}\left(\mathbf{d}_{p}\cdot\mathbf{e}_{k}\right),$
(7)

with $p,q\in\\{\mu,\bar{\mu},\Delta,D\\}$ and we denote
$\mathbf{d}_{\bar{\mu}}=\mathbf{d}_{\mu}^{*}$. We refer to Eq. (II) as the lab
frame Hamiltonian, which is equivalent to Eq. (1).

The photon-only part of Eq. (II) can be diagonalised by the displacement
transformation $H_{d}=B(D/\nu)HB(-D/\nu)$, where subscript $d$ denotes the
displaced frame, and displacement operators are given by

$B(\alpha)=\mathrm{e}^{\sum_{k}\left(\alpha_{k}a_{k}^{\dagger}-\alpha_{k}^{*}a_{k}\right)},$
(8)

and $B(\alpha)^{\dagger}=B(-\alpha)$. These act on photon operators by

$B(\pm\alpha)a_{k}B(\mp\alpha)=a_{k}\mp\alpha_{k},$ (9)

and we further analyse the displacement operators in Appendix A. Ignoring
terms proportional to the identity, the resulting Hamiltonian is

$\displaystyle H_{d}=$
$\displaystyle\frac{\epsilon}{2}\sigma_{z}^{\prime}+|V|\sigma_{x}^{\prime}+\sum_{k}\nu_{\mathbf{k}}a_{k}^{\dagger}a_{k}$
$\displaystyle+\pi_{\Delta\Delta}\sigma_{z}^{\prime}+\pi_{\mu\bar{\mu}}\mathrm{e}^{-i\vartheta_{V}}\sigma_{+}^{\prime}+\pi_{\bar{\mu}\mu}\mathrm{e}^{i\vartheta_{V}}\sigma_{-}^{\prime}.$
(10)

Notice that the dipolar self energy term has cancelled exactly with the terms
that result from displacing the light-matter interactions, such that Eq. (II)
is independent of $D_{k}$.

Throughout this paper, we make the standard assumption that the transition
dipole moment $\mathbf{d}_{\mu}$ is small enough to permit an accurate second-
order perturbative expansion in its magnitude $|\mathbf{d}_{\mu}|$. Even for
perturbative transition dipole moments, the master equation derived using
$H_{d}$ in Eq. (II), referred to as the displaced frame master equation
(DFME), will become inaccurate if $|\mathbf{d}_{\Delta}|$ is large.

To overcome this challenge we make a polaron transformation prior to deriving
the Redfield master equation. In this frame, the polarising effect of the pure
dephasing interaction on the field is absorbed into the definition of a new
basis, called the polarised Fock states [34, 42]. The polarisation direction
of the Fock state is dependent on the matter state, and it hybridises 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 $\Delta_{k}\to 0$.

Figure 1: 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_{\text{dip}}$ causes the renormalisations
$\epsilon\to\epsilon_{\text{dip}}$ and $V\to V_{\text{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 polarised
Fock states of the displaced photon modes.

The polaron frame Hamiltonian is $H_{p}=UH_{d}U^{\dagger}$, where

$U=B\left(\frac{\Delta}{\nu}\right)\ket{e^{\prime}}\bra{e^{\prime}}+B\left(-\frac{\Delta}{\nu}\right)\ket{g^{\prime}}\bra{g^{\prime}}.$
(11)

Using Eq. (9), and ignoring constant terms, we obtain

$\displaystyle
H_{p}=\frac{\epsilon}{2}\sigma_{z}^{\prime}+\kappa|V|\sigma_{x}^{\prime}+\sum_{k}\nu_{\mathbf{k}}a_{k}^{\dagger}a_{k}+\mathcal{C}\sigma_{+}^{\prime}+\mathcal{C}^{\dagger}\sigma_{-}^{\prime},$
(12)

where the coupling operator is $\mathcal{C}=C-\langle C\rangle$ with

$\displaystyle
C=B\left(\frac{\Delta}{\nu}\right)\pi_{\mu\bar{\mu}}B\left(\frac{\Delta}{\nu}\right)\mathrm{e}^{-i\vartheta_{V}}+|V|B\left(2\frac{\Delta}{\nu}\right).$
(13)

We denote $\langle\cdot\rangle=\text{Tr}_{E}(\cdot\rho_{E})$ where

$\rho_{E}=\frac{1}{\mathcal{Z}_{E}}\mathrm{e}^{-\beta\sum_{k}\nu_{\mathbf{k}}a_{k}^{\dagger}a_{k}},$
(14)

and
$\mathcal{Z}_{E}=\text{Tr}[\exp(-\beta\sum_{k}\nu_{\mathbf{k}}a_{k}^{\dagger}a_{k})]$
with $\beta=1/(k_{B}T)$ the inverse temperature. In Appendix B we prove that
$\langle C\rangle=\kappa|V|$ where $\kappa=\langle B(\pm 2\Delta/\nu)\rangle$.
In Eq. (12) we moved a factor of $\langle C\rangle\sigma_{x}^{\prime}$ from
the coupling operator and into the unperturbed part of $H_{p}$. This is to
ensure that the perturbation theory yields the Redfield equation.

In the polaron frame Hamiltonian in Eqs. (12) we can now interpret
$\sigma_{\pm}^{\prime}$ 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 Eqs. (II), polarons experience a more
complicated interaction, albeit one without a pure dephasing term. We also
note that the polaron frame Hamiltonian is diagonal when the transition
dipoles and driving vanish. Hence, the permanent dipoles generate trivial
dynamics if considered in isolation, described by the independent boson model.

Finally, the unperturbed polaron Hamiltonian can be diagonalised using a
unitary rotation,
$(\epsilon/2)\sigma_{z}^{\prime}+\kappa|V|\sigma_{x}^{\prime}=(\eta/2)\tau_{z}$,
where

$\displaystyle\eta=$ $\displaystyle\
\sqrt{\epsilon^{2}+4\kappa^{2}\left|V\right|^{2}},$ (15)
$\displaystyle\tau_{z}=$ $\displaystyle\ \ket{+}\bra{+}-\ket{-}\bra{-}.$ (16)

The eigenbasis relates to the original basis by

$\begin{pmatrix}\ket{e^{\prime}}\\\
\ket{g^{\prime}}\end{pmatrix}=\begin{pmatrix}\cos\left(\frac{\varphi}{2}\right)&-\sin\left(\frac{\varphi}{2}\right)\\\
\sin\left(\frac{\varphi}{2}\right)&\cos\left(\frac{\varphi}{2}\right)\end{pmatrix}\begin{pmatrix}\ket{+}\\\
\ket{-}\end{pmatrix},$ (17)

with $\cos(\varphi)=\epsilon/\eta$ and $\sin(\varphi)=2\kappa|V|/\eta$. In the
eigenbasis, the polaron frame Hamiltonian is

$H_{p}=\frac{\eta}{2}\tau_{z}+\sum_{k}\nu_{\mathbf{k}}a_{k}^{\dagger}a_{k}+\sum_{\alpha\in\\{z,+,-\\}}g_{\alpha}\tau_{\alpha},$
(18)

where we have defined $\tau_{+}=\ket{+}\bra{-}$, $\tau_{-}=\ket{-}\bra{+}$,
and the coupling operators

$\displaystyle g_{z}=$ $\displaystyle\
\frac{1}{2}\sin\left(\varphi\right)\left[\mathcal{C}+\mathcal{C}^{\dagger}\right],$
(19a) $\displaystyle g_{+}=$ $\displaystyle\
\left[\cos^{2}\left(\frac{\varphi}{2}\right)\mathcal{C}-\sin^{2}\left(\frac{\varphi}{2}\right)\mathcal{C}^{\dagger}\right],$
(19b)

and $g_{-}=g_{+}^{\dagger}$.

To derive the master equation we will take the continuum limit of the photon
modes, in which summations over an arbitrary function $F(\nu)$ transform
according to

$\displaystyle\sum_{k}f_{\mathbf{k}}^{2}F(\nu_{\mathbf{k}})$
$\displaystyle\left(\mathbf{d}_{p}\cdot\mathbf{e}_{k}\right)\left(\mathbf{d}_{q}\cdot\mathbf{e}_{k}\right)$
$\displaystyle\to h_{pq}\int_{0}^{\infty}\text{d}\nu\ J(\nu)F(\nu),$ (20)

where $p,q\in\\{\mu,\bar{\mu},\Delta\\}$ and

$h_{pq}=\int_{\Omega_{\mathbf{k}}}\text{d}\Omega_{\mathbf{k}}\
\sum_{\lambda}\left(\tilde{\mathbf{d}}_{p}\cdot\mathbf{e}_{k}\right)\left(\tilde{\mathbf{d}}_{q}\cdot\mathbf{e}_{k}\right).$
(21)

We have defined a spectral density $J(\nu)$, a dimensionless dipole vector
$\tilde{\mathbf{d}}_{p}=\mathbf{d}_{p}/|\mathbf{d}_{\text{ref}}|$ measured
against an arbitrary reference value $|\mathbf{d}_{\text{ref}}|$, and the
solid angle integral
$\int_{\Omega_{\mathbf{k}}}\text{d}\Omega_{\mathbf{k}}=\int_{0}^{\pi}\text{d}\theta_{\mathbf{k}}\
\sin(\theta_{\mathbf{k}})\int_{0}^{2\pi}\text{d}\phi_{\mathbf{k}}\ $ [2]. For
free space and unpolarised light,

$\displaystyle h_{\mu\Delta}$
$\displaystyle=\Omega_{\mu\Delta}\mathrm{e}^{i\vartheta_{\mu}}\cos\left(\theta_{\mu\Delta}\right),$
(22a) $\displaystyle h_{\mu\bar{\mu}}$ $\displaystyle=\Omega_{\mu\mu},$ (22b)
$\displaystyle h_{\mu\mu}$
$\displaystyle=\Omega_{\mu\mu}\mathrm{e}^{2i\vartheta_{\mu}},$ (22c)

where
$\Omega_{pq}=\Omega_{qp}=(8\pi/3)|\tilde{\mathbf{d}}_{p}||\tilde{\mathbf{d}}_{q}|\in\mathbb{R}$,
$\vartheta_{\mu}$ is the complex phase of $\mathbf{d}_{\mu}$,
$\theta_{\mu\Delta}$ is the angle between $\mathbf{d}_{\Delta}$ and
$\mathbf{d}_{\mu}$, and one can obtain e.g.
$h_{\bar{\mu}\Delta}=h_{\mu\Delta}^{*}$ by suitable complex conjugation.
Notice that the dynamics only depend on the relative phase $\vartheta_{\mu
V}=\vartheta_{\mu}-\vartheta_{V}$.

We use the free space multipolar-gauge spectral density

$J(\nu)=S\frac{\nu^{3}}{\nu_{c}^{2}}\mathrm{e}^{-\frac{\nu}{\nu_{c}}}\Theta(\nu),$
(23)

which has dimensions of energy. $\Theta(\nu)$ is the Heaviside step function
and $S=\int_{0}^{\infty}\text{d}\nu\ J(\nu)/\nu^{2}$ is a dimensionless
parameter called the Huang-Rhys parameter. We have introduced a
phenomenological cut-off to the spectral density $\nu_{c}$, which is justified
for finite systems [55]. In the continuum limit, $\kappa=\exp[-\phi(0)/2]$
where the photon propagator is

$\displaystyle\phi(s)=4\Omega_{\Delta\Delta}\int_{0}^{\infty}\text{d}\nu\
\frac{J(\nu)}{\nu^{2}}$
$\displaystyle\big{[}\coth\left(\frac{\beta\nu}{2}\right)\cos\left(\nu
s\right)$ $\displaystyle\qquad\qquad-i\sin\left(\nu s\right)\big{]}.$ (24)

To summarise, the important parameters in this model are the dipole strengths
$|\tilde{\mathbf{d}}_{\mu}|$ and $|\tilde{\mathbf{d}}_{\Delta}|$, their
relative angle $\theta_{\mu\Delta}$, the eigenenergy $\eta$, the driving
strength $|V|$, and the relative complex phase $\vartheta_{\mu V}$. We denote
by $\varphi$ the angle determining the eigenbasis, by $\theta_{\mu\Delta}$ the
angle between the dipole vectors, and by $\vartheta_{\mu V}$ the relative
complex phases. Whether or not the transition or permanent dipole interactions
are deemed strong is approximately determined by
$|\tilde{\mathbf{d}}_{p}|\lambda>\epsilon$ for $p\in\\{\mu,\Delta\\}$ where
the reorganisation energy is

$\lambda=\int_{0}^{\infty}\text{d}\nu\ \frac{J(\nu)}{\nu}=2\nu_{c}S.$ (25)

##  III Effects of strong permanent dipoles

In this section we derive the secularised PFME, provide a brief analytical
review of it, and analyse 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 [44, 45, 46, 47, 48, 49, 50, 51,
52], which both serve as benchmarks. In the numerical approach, we use the
full non-secular Redfield master equations derived in Appendix C for both
frames. We emphasise that the non-secular PFME depends only on the physical
processes discussed in the main text.

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

$\displaystyle\frac{\partial\rho_{++}(t)}{\partial t}$
$\displaystyle=-\gamma_{\downarrow}\rho_{++}(t)+\gamma_{\uparrow}\rho_{--}(t),$
(26a) $\displaystyle\frac{\partial\rho_{+-}(t)}{\partial t}$
$\displaystyle=-[\gamma_{d}+i\bar{\eta}]\rho_{+-}(t),$ (26b)

and $(\partial/\partial t)\rho_{--}(t)=-(\partial/\partial t)\rho_{++}(t)$,
$(\partial/\partial t)\rho_{-+}(t)=(\partial/\partial
t)\rho_{+-}(t)^{\dagger}$, where $\rho_{ij}(t)=\braket{i}{\rho_{S}(t)}{j}$ for
$i,j\in\\{+,-\\}$ and $\rho_{S}(t)=\text{Tr}_{E}[\rho(t)]$. The transition and
decoherence rates are

$\displaystyle\gamma_{\begin{subarray}{c}\downarrow\\\
\uparrow\end{subarray}}$
$\displaystyle=2\Re\left[\Gamma_{\mp\mp}(\pm\eta)\right],$ (27a)
$\displaystyle\gamma_{d}$
$\displaystyle=\tfrac{1}{2}\left[\gamma_{\uparrow}+\gamma_{\downarrow}\right]+4\Re\left[\Gamma_{zz}(0)\right],$
(27b)

where

$\displaystyle\Gamma_{\alpha\beta}(\omega)=\int_{0}^{\infty}\text{d}s\
\mathrm{e}^{i\omega s}\left\langle
g_{\alpha}^{\dagger}(s)g_{\beta}(0)\right\rangle,$ (28)

for $\alpha,\beta\in\\{z,+,-\\}$ and $g_{\alpha}(s)$ is the interaction
picture form of $g_{\alpha}$ in Eqs. (19). Finally, the Lamb shifted
eigenenergy is

$\displaystyle\bar{\eta}=\eta+\Im\left[\Gamma_{--}(\eta)-\Gamma_{++}(-\eta)\right].$
(29)

The secular master equation in Eqs. (26) describes population transfer from
$\ket{+}$ to $\ket{-}$ at decay rate $\gamma_{\downarrow}$, vice-versa at an
excitation rate $\gamma_{\uparrow}$, decoherence at a rate $\gamma_{d}$, and
oscillations in the coherence at frequency $\bar{\eta}$.

In the PFME we choose the initial state to be
$\rho_{p}(0)=\left|g\right>\left<g\right|\otimes\rho_{E}$ where $\rho_{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. VI.

Evaluating the ECFs in Eq. (28), $\langle
g_{\alpha}^{\dagger}(s)g_{\beta}(0)\rangle$, is a laborious process and we
present the full derivations in Appendix D. Here, we focus the discussion on
new physical processes introduced by the presence of the strong permanent
dipoles. The ECFs, $\langle g_{\alpha}^{\dagger}(s)g_{\beta}(0)\rangle$,
depend on linear combinations of
$\langle\mathcal{C}^{\dagger}(s)\mathcal{C}(0)\rangle$,
$\langle\mathcal{C}(s)\mathcal{C}^{\dagger}(0)\rangle$,
$\langle\mathcal{C}^{\dagger}(s)\mathcal{C}^{\dagger}(0)\rangle$ and
$\langle\mathcal{C}(s)\mathcal{C}(0)\rangle$, with the relative weighting of
each dependent on the eigenbasis angle $\varphi$. Importantly, as we show in
Appendix D, each of these correlation functions describes the same physics.
Therefore, for the analytics in the main text, we focus the discussion on

$\Gamma(\omega)=\int_{0}^{\infty}\text{d}s\ \mathrm{e}^{i\omega
s}\left\langle\mathcal{C}^{\dagger}(s)\mathcal{C}(0)\right\rangle.$ (30)

This function has contributions from four distinct processes,

$\Gamma(\omega)=\Gamma_{1}(\omega)+\Gamma_{2}(\omega)+\Gamma_{V,1}(\omega)+\Gamma_{V,0}(\omega).$
(31)

The number in the subscripts of each term on the right-hand-side of Eq. (31)
denote 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 $\mathcal{C}(s)$ two-time correlation
functions also depend on these four contributions, except that
$\langle\mathcal{C}^{\dagger}(s)\mathcal{C}^{\dagger}(0)\rangle$ and
$\langle\mathcal{C}(s)\mathcal{C}(0)\rangle$ do not have a
$\Gamma_{V,1}(\omega)$ type contribution. The prefactor of each term is
dependent on the relative dipole angle $\theta_{\mu\Delta}$, the driving
$|V|$, the relative phase $\vartheta_{\mu V}$, and dipole magnitudes
$\Omega_{pq}$, as

$\displaystyle\Gamma_{1}(\omega)$ $\displaystyle\propto\Omega_{\mu\mu},$ (32a)
$\displaystyle\Gamma_{2}(\omega)$
$\displaystyle\propto\Omega_{\mu\Delta}^{2}\cos^{2}(\theta_{\mu\Delta}),$
(32b) $\displaystyle\Gamma_{V,1}(\omega)$
$\displaystyle\propto\Omega_{\mu\Delta}\left|V\right|\cos\left(\vartheta_{\mu
V}\right)\cos(\theta_{\mu\Delta}),$ (32c) $\displaystyle\Gamma_{V,0}(\omega)$
$\displaystyle\propto\left|V\right|^{2}.$ (32d)

From Eqs. (32) we see that when the system is driven, $\Gamma(\omega)$ is not
an even function of $\theta_{\mu\Delta}$ or $\vartheta_{\mu V}$. We will later
show that this can be utilised to control the system, for example to minimise
decoherence. Moreover, for perpendicular dipoles $\cos(\theta_{\mu\Delta})=0$,
so $\Gamma(\omega)$ becomes equivalent to the analogous function in the driven
spin boson model where the transition dipole and permanent dipole interactions
are with independent baths. Interactions proportional to
$\cos(\theta_{\mu\Delta})$ are therefore arise because the interactions are
non-commutative.

The $\Gamma_{V,1}(\omega)$ contribution in Eq. (47) varies cosinusoidally with
the relative complex phase $\vartheta_{\mu V}$. This means that through
controlling the phase of $V$ one can determine the phase $\vartheta_{\mu}$ of
$\mathbf{d}_{\mu}$, for example by minimising the spectral linewidth, allowing
an estimation of the relative phase between the excited and ground state. This
will work provided that $V$ is not generated via a coupling to the system’s
transition dipole moments, in which case
$\vartheta_{V}=\vartheta_{\mu}+\vartheta_{\text{external}}$ and so
$\vartheta_{\mu V}$ would only depend on the phase of the external drive.

For additional context, when $|V|\ll\epsilon$ – the parameter regime in which
polaron theory is valid – the eigenbasis angle is $\varphi\approx 0$ and so
$g_{z}\approx 0$ and $g_{+}\approx\mathcal{C}$. In this limit, the secular
approximation used in the master equation in Eqs. (26) becomes exact, and
$\gamma_{\uparrow}\approx 2\Re[\Gamma(-\eta)]$ with $\Gamma(\omega)$ in Eq.
(31) and $\gamma_{\downarrow}\approx\int_{0}^{\infty}\text{d}s\
\mathrm{e}^{i\eta s}\langle\mathcal{C}(s)\mathcal{C}^{\dagger}(0)\rangle$
which, as we show in Appendix D, evaluates to $\gamma_{\downarrow}\approx
2\Re[\bar{\Gamma}(\eta)]$ where $\bar{\Gamma}(\omega)$ is equal to Eq. (31)
but with $\Gamma_{V,1}(\omega)\to-\Gamma_{V,1}(\omega)$. (Note that this
asymmetry cannot lead to population inversion due to minus signs within
$\Gamma_{V,1}(\omega)$.)

We will now discuss the four contributions to Eq. (31) in turn. Each are
derived in Appendix D.

###  III.1 One photon processes, $\Gamma_{1}(\omega)$

The function describing driving-independent one photon processes is

$\displaystyle\Gamma_{1}(\omega)=\pi\Omega_{\mu\mu}\int_{0}^{\infty}\text{d}\nu\
\big{[}$ $\displaystyle J_{A}(\nu)\mathcal{K}(\omega+\nu)$
$\displaystyle\qquad+J_{E}(\nu)\mathcal{K}(\omega-\nu)\big{]},$ (33)

where

$\mathcal{K}(\varepsilon)=\frac{1}{\pi}\int_{0}^{\infty}\text{d}s\
\mathrm{e}^{i\varepsilon s}\mathrm{e}^{\phi(s)-\phi(0)},$ (34)

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

$\displaystyle J_{A}(\nu)$ $\displaystyle=J(\nu)N(\nu),$ (35) $\displaystyle
J_{E}(\nu)$ $\displaystyle=J(\nu)\tilde{N}(\nu),$ (36)

where $N(\nu)=1/[\exp(\beta\nu)-1]$ is the Bose-Einstein distribution and
$\tilde{N}(\nu)=N(\nu)+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[\phi(s)-\phi(0)]$ in Eq.
(34). The resulting term is the typical function appearing in the standard
optical master equation (SOME) [1] \- our model in the absence of permanent
dipoles and driving - which is

$\displaystyle\Gamma_{\text{SOME}}(\omega)=\Omega_{\mu\mu}\bigg{[}\pi\left(J_{A}(-\omega)+J_{E}(\omega)\right)$
$\displaystyle\qquad\qquad\qquad+i\mathcal{P}\int_{0}^{\infty}\text{d}\nu\
\left(\frac{J_{A}(\nu)}{\omega+\nu}+\frac{J_{E}(\nu)}{\omega-\nu}\right)\bigg{]},$
(37)

where $\mathcal{P}$ denotes the principal value. Clearly, when $\omega=-\eta$,
twice the real part of $\Gamma_{\text{SOME}}(\omega)$ describes excitation-by-
absorption at a rate $2\pi\Omega_{\mu\mu}J_{A}(\eta)$, and when $\omega=\eta$
it describes decay-by-emission at a rate $2\pi\Omega_{\mu\mu}J_{E}(\eta)$. The
imaginary part of $\Gamma_{\text{SOME}}(\omega)$ will determine the Lamb
shift.

We now return to Eq. (III.1). In Ref. [57] an analytic solution to Eq. (34)
was found by exploiting the fact that $\mathcal{K}(\varepsilon)$ only depends
on $J(\nu)$ through its moments
$\mu_{m}=4\Omega_{\Delta\Delta}\int_{0}^{\infty}\text{d}\nu\
[J(\nu)/\nu^{2}]\nu^{m}$ for $m=1,2,\ldots,\infty$, and that moments of lower
order contribute relatively more. We can then evaluate
$\mathcal{K}(\varepsilon)$ by replacing $J(\nu)$ with a truncated spectral
density
$J^{\prime}(\nu)=\sum_{\mathbf{k}=1}^{N_{*}}|f^{\prime}_{\mathbf{k}}|^{2}\delta(\nu-\nu_{\mathbf{k}}^{\prime})$,
as long as we choose the coupling strengths $\\{f_{\mathbf{k}}^{\prime}\\}$
and energies $\\{\nu_{\mathbf{k}}^{\prime}\\}$ of the modes such that
$J^{\prime}(\nu)/\nu^{2}$ has the same lowest moments as $\mu_{m}$ for
$m=1,2,\ldots,2N_{*}$ [57]. The solution converges rapidly for increasing
$N_{*}$, and a single-mode truncation is often very accurate, with truncation
mode parameters: $\nu^{\prime}_{1}\equiv\nu_{s}=\mu_{2}/\mu_{1}$ and
$f^{\prime}_{1}\equiv f_{s}=\mu_{2}^{1/2}\propto\Omega_{\Delta\Delta}^{1/2}$
[57]. A single-mode truncation captures all of the essential physics described
by $\mathcal{K}(\omega)$. Thus, for the analytical analysis in this work, we
use a single-mode truncation to evaluate Eq. (34); however, in all of the
simulations in this paper, we increase $N_{*}$ until convergence.

From Ref. [57], the single mode solution of Eq. (34) is

$\mathcal{K}(\varepsilon)=\sum_{\ell=-\infty}^{\infty}A_{\ell}\left[\delta(\varepsilon-\ell\nu_{s})+\frac{i}{\pi}\frac{\mathcal{P}}{\varepsilon-\ell\nu_{s}}\right],$
(38)

where

$A_{\ell}=\sum_{n=|\ell|}^{\infty^{\prime}}\sum_{m=\frac{n-\ell}{2}}^{n}\binom{n}{m}\binom{m}{m-\frac{n-\ell}{2}}W_{n}V_{m},$
(39)

and the prime on the first summation indicates that only every other term is
included, i.e. $n=|\ell|,|\ell|+2,\ldots,$ and
$W_{n}=S^{n}_{s}\exp[-S_{s}]/n!$ is the Franck-Condon factor of the mode in
the truncation, $S_{s}=|f_{s}|^{2}/\nu_{s}^{2}=\mu_{1}^{2}/\mu_{2}$ is its
Huang-Rhys factor and $V_{m}=N(\nu_{s})^{m}\exp[-2S_{s}N(\nu_{s})]$ [57].
$A_{\ell}$ has the normalisation property
$\sum_{\ell=-\infty}^{\infty}A_{\ell}=1$, is maximised for
$\ell=\text{Round}(S^{\prime})$, becomes $A_{\ell}=\delta_{\ell 0}$ if
$|\tilde{\mathbf{d}}_{\Delta}|=0$, and the $\ell<0$ terms are only non-zero at
finite temperature (see Ref. [57] for a full discussion on the mode truncation
solution).

Substituting Eq. (38) into Eq. (III.1) yields

$\Gamma_{1}(\omega)=\frac{1}{2}\gamma_{1}(\omega)+iS_{1}(\omega),$ (40)

where

$\gamma_{1}(\omega)=2\pi\Omega_{\mu\mu}\sum_{\ell=-\infty}^{\infty}A_{\ell}\Big{[}J_{A}(\ell\nu_{s}-\omega)+J_{E}(\omega-\ell\nu_{s})\Big{]},$
(41)

and

$\displaystyle S_{1}(\omega)=i\Omega_{\mu\mu}$
$\displaystyle\sum_{\ell=-\infty}^{\infty}A_{\ell}\mathcal{P}\int_{0}^{\infty}\text{d}\nu\
$
$\displaystyle\times\bigg{(}\frac{J_{A}(\nu)}{\omega+\nu-\ell\nu_{s}}+\frac{J_{E}(\nu)}{\omega-\nu-\ell\nu_{s}}\bigg{)}.$
(42)

Eqs. (41) and (III.1) show that the effect of the permanent dipoles within
$\Gamma_{1}(\omega)$ is to introduce manifolds of harmonic levels (in the
single mode truncation of $\mathcal{K}(\varepsilon)$ 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. $\Gamma_{\text{SOME}}(\omega)$ in Eq.
(III.1) is recovered from $\Gamma_{1}(\omega)$ in the limit of no permanent
dipoles because $A_{\ell}\to\delta_{\ell 0}$.

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

Figure 2: Illustration of decay processes ($\omega=+\eta$) captured by the
$\ell=1$ and $\ell=\ell_{*}$ terms in Eq. (41), within the single mode
truncation of $\mathcal{K}(\varepsilon)$. Decay-by-emission processes have
dot-dashed arrows, and decay-by-absorption processes have solid arrows. The
type of decay that occurs for a given value of $\ell$ depends on the sign of
$\eta-\ell\nu_{s}$. In the absence of permanent dipoles $A_{\ell}=\delta_{\ell
0}$ and so only the $\ell=0$ decay processes are possible.

###  III.2 Two photon processes, $\Gamma_{2}(\omega)$

We now consider the two photon processes in Eq. (31) that occur independent of
the driving, $\Gamma_{2}(\omega)$. These processes require that the
$\sigma_{z}$ and $\sigma_{\pm}$ do not commute, and so are unique to permanent
dipole interactions. As shown in Appendix D, the rate function is

$\displaystyle\Gamma_{2}(\omega)=$ $\displaystyle
4\kappa^{2}\Omega_{\mu\Delta}^{2}\cos^{2}\left(\theta_{\mu\Delta}\right)$
$\displaystyle\times\sum_{n,m\in\\{-1,1\\}}\int_{0}^{\infty}\text{d}s\
\mathrm{e}^{i\omega s+\phi(s)}\tilde{\chi}_{n,m}(s),$ (43)

where

$\tilde{\chi}_{n,m}(s)=\int_{0}^{\infty}\text{d}\nu\
\int_{0}^{\infty}\text{d}\nu^{\prime}\
\chi_{n,m}(\nu,\nu^{\prime})\mathrm{e}^{in\nu s}\mathrm{e}^{im\nu^{\prime}s},$
(44)

is the inverse Fourier transform of $\chi_{n,m}(\nu,\nu^{\prime})$ which we
write in matrix notation as

$\displaystyle\chi(\nu,\nu^{\prime})=\bar{J}(\nu)\bar{J}(\nu^{\prime})\begin{pmatrix}\tilde{N}(\nu)\tilde{N}(\nu^{\prime})&-\tilde{N}(\nu)N(\nu^{\prime})\\\
-N(\nu)\tilde{N}(\nu^{\prime})&N(\nu)N(\nu^{\prime})\end{pmatrix},$ (45)

where, for example,
$\chi_{1,1}(\nu,\nu^{\prime})=4\Omega_{\mu\Delta}^{2}\bar{J}(\nu)\bar{J}(\nu^{\prime})N(\nu)N(\nu^{\prime})$.
We have introduced the polarised spectral density,

$\bar{J}(\nu)=\frac{J(\nu)}{\nu},$ (46)

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

The two-frequency Fourier transform in Eq. (44) indicates that
$\Gamma_{2}(\omega)$ describes two simultaneous processes. From the phase
factors in Eq. (44), 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 $\nu$ and $\nu^{\prime}$ for $n$ and $m$,
respectively. This interpretation is reinforced by the positions of the
factors of $N(\nu)$, $N(\nu^{\prime})$, $\tilde{N}(\nu)$ and
$\tilde{N}(\nu^{\prime})$ in the matrix in Eq. (45). Similar to the case of
the one photon processes, the factor of $\exp[\phi(s)-\phi(0)]$ in Eq. (III.2)
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 $\gamma_{\downarrow}$ as a function of
$|\tilde{\mathbf{d}}_{\Delta}|$ and $\nu_{c}$ for $V=0$. In Fig. 3, the
eigenbasis coincides with the $\\{\ket{e},\ket{g}\\}$ basis ($\varphi=0$) and
the decay rate is $\gamma_{\downarrow}=2\text{Re}[\Gamma(\epsilon)]$ where
$\Gamma(\omega)$ is in Eq. (31). One can see that $\Gamma(\omega)$ is
independent of the phase $\vartheta_{\mu}$ when $V=0$.

In Fig. 3(a) and (b), it is clear that at small $\nu_{c}$ it is possible for
permanent dipoles to enhance the decay rate, whilst at large $\nu_{c}$ they
always lead to suppression. This is because of how the functions $J(\nu)$ and
$\bar{J}(\nu)$ vary with $\nu$ at the sideband-modified frequencies compared
with the changes in the mode populations governed by $N(\nu)$. Fig. 3(c) shows
that two photon processes become negligible for large $\nu_{c}$. This is
because two photon processes scale as
$J(\nu)J(\nu^{\prime})\propto\nu_{c}^{-4}$ (see Eq. (23)) whilst one photon
processes scale as $J(\nu)\propto\nu_{c}^{-2}$. For the same reason, this will
also occur in the limit of small Huang-Rhys parameter, $S\to 0$.

Figure 3: Polaron frame decay rate for $V=0$. (a) and (b) shows the total
decay rate (dashed black), and the contributions of the one photon processes
(solid green) and two photon processes (dot-dashed purple) verses
$|\tilde{\mathbf{d}}_{\Delta}|$. In (a) the cutoff frequency is
$\nu_{c}=0.2~{}\text{eV}$ and in (b) $\nu_{c}=1~{}\text{eV}$. (c) shows the
total rate and its contributions verses $\nu_{c}$ for
$|\tilde{\mathbf{d}}_{\Delta}|=0.26$ which corresponds to the maximum in (a).
Other parameters are $|\tilde{\mathbf{d}}_{\mu}|=0.01$,
$\epsilon=1~{}\text{eV}$, $\beta=2~{}\text{eV}^{-1}$ ($T\approx
5800~{}\text{K}$), $S=1/\pi$, and $\theta_{\mu\Delta}=0$ which, unless
otherwise stated, are used in all figures.

###  III.3 Driving induced one photon processes, $\Gamma_{V,1}(\omega)$

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

$\displaystyle\Gamma_{V,1}$
$\displaystyle(\omega)=4\pi\Omega_{\mu\Delta}|V|\cos\left(\vartheta_{\mu
V}\right)\cos\left(\theta_{\mu\Delta}\right)$ (47)
$\displaystyle\times\int_{0}^{\infty}\text{d}\nu\
\bigg{(}-\bar{J}_{A}(\nu)\mathcal{K}(\omega+\nu)+\bar{J}_{E}(\nu)\mathcal{K}(\omega-\nu)\bigg{)},$

where $\mathcal{K}(\varepsilon)$ is given in Eq. (34) and, analogously to Eqs.
(35), we have introduced the polarised absorption and emission spectral
densities,

$\displaystyle\bar{J}_{A}(\nu)$ $\displaystyle=\bar{J}(\nu)N(\nu),$ (48)
$\displaystyle\bar{J}_{E}(\nu)$ $\displaystyle=\bar{J}(\nu)\tilde{N}(\nu),$
(49)

where $\bar{J}(\nu)$ is in Eq. (46). Similarly to the driving-independent two
photon processes, $\Gamma_{2}(\omega)$ in Eq. (III.2), these processes also
arise due to the non-commutativity of the permanent and transition dipole
interactions and vanish when $\mathbf{d}_{\mu}\cdot\mathbf{d}_{\Delta}=0$.

Using the single mode truncation solution in Eq. (38), we can decompose
$\Gamma_{V,1}(\omega)$ into real and imaginary parts,

$\Gamma_{V,1}(\omega)=\frac{1}{2}\gamma_{V,1}(\omega)+iS_{V,1}(\omega),$ (50)

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 $\Gamma_{1}(\omega)$ but here scaling with
$\Omega_{\mu\Delta}$ and dependent on the polarised spectral density.

Importantly, the prefactor in Eq. (47) is proportional to $\cos(\vartheta_{\mu
V})\cos(\theta_{\mu\Delta})$ which allows one to control whether this function
suppresses or enhances the rates and Lamb shifts. In practice, this control
could be easy to achieve by tuning the driving phase.

###  III.4 Driving induced zero photon processes, $\Gamma_{V,0}(\omega)$

Finally, the zero photon processes induced by the driving in Eq. (31) is

$\Gamma_{V,0}(\omega)=\kappa^{2}|V|^{2}\int_{0}^{\infty}\text{d}s\
\mathrm{e}^{i\omega s}\left(\mathrm{e}^{\phi(s)}-1\right).$ (51)

This term vanishes when $\phi(s)=0$ and so is generated entirely through the
photon sideband. Eq. (51) also appears in the driven spin boson model
describing phonon mediated transfer between eigenstates and so is an effect of
pure dephasing.

In Fig. 4, panel (a) shows the total decay rate for different driving
strengths $V$ as a function of $|\tilde{\mathbf{d}}_{\Delta}|$, and panel (b)
shows the four contributions to the rate with $V=0.02~{}\text{eV}$. The dashed
black curves in both panels (a) and (b) are the same. The excitation
$\gamma_{\uparrow}$ and decoherence $\gamma_{d}$ rates can be obtained from
Fig. 3 by scaling the curves by $0.135$ and $0.568$, respectively; all three
rates have the same dependence on the permanent dipoles unless
$|V|\sim\epsilon$, but in that limit polaron theory becomes inaccurate. The
key points resulting from the analysis of Fig. 3 are: (1) one can both
increase and decrease rates through the physical mechanisms enabled by the
presence of permanent dipoles and driving, particularly by controlling
$\Gamma_{V,1}(\omega)$ processes and (2) all physical processes found in the
polaron ECFs contribute non-negligibly to the total rates. Recall from Fig. 3
that the contribution of the two photon processes $\Gamma_{2}(\omega)$ is
larger for smaller $\nu_{c}$.

Figure 4: Polaron frame decay rate for different permanent dipole and driving
strengths. (a) shows the total decay rate for different driving strengths. In
all curves $\vartheta_{\mu}=0$ and $\vartheta_{V}\in\\{0,\pi\\}$ controls the
sign of $V$ given in the legend. (b) shows the four contributions (coloured
curves) to the total rate (dashed black), given in Eq. (31), for the
$V=0.02~{}\text{eV}$ parameter set. The markers in panel (b) refer to
$|\tilde{\mathbf{d}}_{\Delta}|$ values used in Fig. 5. The left-hand-side of
each panel has anti-parallel $\mathbf{d}_{\mu}$ and $\mathbf{d}_{\Delta}$,
whilst the right-hand-side has parallel dipole moments, which is indicated by
the arrows at the top of the figure. $\nu_{c}=1~{}\text{eV}$, and other
parameters are given at the end of the caption for Fig. 3.

Through control of $|\tilde{\mathbf{d}}_{\Delta}|$ and $V$ one can suppress
the transition and decoherence rates. For example, the decay rate at the
minimum when $V=0.05$ in Fig. 3(a), located at
$|\tilde{\mathbf{d}}_{\Delta}|\approx 0.1$, is suppressed by a factor of
$\approx 4000$ compared to its maximum value. Rather than tuning
$|\tilde{\mathbf{d}}_{\Delta}|$, similar control can be achieved through the
phase of the driving $V$. In Fig. 5 we demonstrate this by varying
$\vartheta_{\mu V}$ with $V=0.02$ and the $|\tilde{\mathbf{d}}_{\Delta}|$
values marked in Fig. 4(b). 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(\theta_{\mu\Delta})\cos(\vartheta_{\mu V})$ in Eq. (47).

Figure 5: Phase control of the polaron frame decay rate. The relative phase
$\vartheta_{\mu V}$ is varied whilst keeping $|\tilde{\mathbf{d}}_{\mu}|=0.01$
and $|V|=0.02~{}\text{eV}$ fixed. Each curve has a different value of
$|\tilde{\mathbf{d}}_{\Delta}|$ indicated in the legend which correspond to
the markers in Fig. 4(b). $\nu_{c}=1~{}\text{eV}$, and other parameters are
given at the end of the caption for Fig. 3.

##  IV 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. The DFME does not capture many of the unique processes attributed to
the permanent dipoles as seen in the polaron frame, and we expect these
perturbative approaches to become inaccurate when
$|\tilde{\mathbf{d}}_{\Delta}|\lambda$ approaches $\epsilon$ where
$\lambda=2\nu_{c}S$ is given in Eq. (25). This is indeed manifest in Fig. 6
where we plot the population of the ground state and the coherence between the
excited and ground state against time.

Figure 6: Comparison of the PFME predictions to exact numerical approach and
the DFME. (a) and (c) are a comparison of $\rho_{gg}(t)$ calculated using the
PFME (green) to TEMPO (dot-dashed black), and the DFME (dashed orange) over
short and long times. (b) and (d) are a similar comparison but for
$|\rho_{eg}(t)|$ (note the discontinuity in the time axis). In (d) the cyan
dotted curve is the polaron prediction for the coherence artificially
accounting for the initial non-Markovian slip. In (a) and (b)
$\epsilon=0.36~{}\text{eV}$, $V=-0.0064~{}\text{eV}$, and
$|\tilde{\mathbf{d}}_{\Delta}|=0.5$, and in (c) and (d)
$\epsilon=1~{}\text{eV}$, $V=0~{}\text{eV}$, and
$|\tilde{\mathbf{d}}_{\Delta}|=1$. In both plots $\nu_{c}=1~{}\text{eV}$, and
other parameters are given at the end of the caption for Fig. 3. To make
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\pi/3$ in the $\Omega_{pq}$ in Eqs. (22).

The PFME predicts qualitatively correct dynamics and agrees well with TEMPO,
even for strong permanent dipoles. In Fig. 6(a) and (b),
$|\tilde{\mathbf{d}}_{\Delta}|=0.5$ and $\epsilon=0.36~{}\text{eV}$, and in
(c) and (d) $|\tilde{\mathbf{d}}_{\Delta}|=1$ and $\epsilon=1~{}\text{eV}$,
with $\lambda=2/\pi~{}\text{eV}$ in all panels. Hence,
$|\tilde{\mathbf{d}}_{\Delta}|\lambda$ is $0.32~{}\text{eV}$ and
$0.64~{}\text{eV}$ which is comparable to $\epsilon$ in both panels. This
emphasises that each new process in Fig. 3 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’ [58] 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.

##  V Emission spectrum

Due to the many new physical processes generated by the presence of permanent
dipoles discussed in Sec. III we expect alterations to the emission spectrum.
In the absence of permanent dipoles, the spectrum consists of a Mollow triplet
with peaks at frequencies $-\eta,0,\eta$. When $V=0$, the triplet becomes a
single peak at frequency $\epsilon$ with a width determined by the decoherence
rate in the standard optical master equation: $\pi
J(\epsilon)[1+2N(\epsilon)]$. In the presence of strong permanent dipoles, we
expect that there will be a photon sideband extending to negative frequencies.
If $|\tilde{\mathbf{d}}_{\Delta}|$ is large enough, the positions of the
Mollow triplet peaks will shift noticeably from $\pm\eta$ to $\pm\bar{\eta}$
given in Eq. (29), and the widths of the peaks will be determined by
$\gamma_{d}$ in Eq. (27b).

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

$I(\omega)=\alpha_{\text{prop}}(\mathbf{r},\mathbf{R},\omega)I_{0}(\omega),$
(52)

where
$\alpha_{\text{prop}}(\mathbf{r},\mathbf{R},\omega)=|\mathbf{d}_{\mu}\cdot\mathbf{G}(\mathbf{r},\mathbf{R},\omega)|^{2}$
and $\mathbf{G}(\mathbf{r},\mathbf{R},\omega)$ is the Green’s function of the
medium [59, 60, 61, 62, 63].
$\alpha_{\text{prop}}(\mathbf{r},\mathbf{R},\omega)$ accounts for propagation
and filtering of the light from the dipole to the detector at positions
$\mathbf{r}$ and $\mathbf{R}$, respectively [59]. The polarisation spectrum is

$I_{0}(\omega)=\left\langle\sigma^{+}(\omega)\sigma^{-}(\omega)\right\rangle,$
(53)

or, equivalently,

$I_{0}(\omega)=\lim_{t\to\infty}\Re\left[\int_{0}^{\infty}\text{d}\tau\
\left\langle\sigma_{+}(t+\tau)\sigma_{-}(t)\right\rangle\mathrm{e}^{-i\omega\tau}\right].$
(54)

In the following, we focus on evaluating the polarisation spectrum, which
captures entirely the effects associated with the permanent dipoles. The
specifics of the experimental set-up, for example, whether the dipole is
coupled to a cavity or a waveguide, are described by the
$\alpha_{\text{prop}}(\mathbf{r},\mathbf{R},\omega)$ term, which can be
calculated separately.

The expectation value in Eq. (54) is taken with respect to the lab frame
density operator $\rho_{l}$. Using $\rho_{l}=U^{\dagger}\rho_{p}U$ where $U$
is the polaron transformation in Eq. (11) and $\rho_{p}$ is the polaron frame
density operator, we have

$\displaystyle\left\langle\sigma_{+}(t+\tau)\sigma_{-}(t)\right\rangle$
$\displaystyle=\kappa^{2}\mathrm{e}^{\phi(\tau)}\text{Tr}\left[\sigma_{+}(t+\tau)\sigma_{-}(t)\rho_{p}(0)\right]$
$\displaystyle\equiv\kappa^{2}\mathrm{e}^{\phi(\tau)}\left\langle\sigma_{+}(t+\tau)\sigma_{-}(t)\right\rangle_{p},$
(55)

where $\kappa$ is defined above Eq. (II). Therefore, the final expression for
the polaron frame polarisation spectrum is

$\displaystyle I_{p,0}(\omega)=\kappa^{2}\lim_{t\to\infty}\Re\Big{[}$
$\displaystyle\int_{0}^{\infty}\text{d}\tau\ \mathrm{e}^{\phi(\tau)}$
$\displaystyle\times\left\langle\sigma_{+}(t+\tau)\sigma_{-}(t)\right\rangle_{p}\mathrm{e}^{-i\omega\tau}\Big{]}.$
(56)

The two-time correlation function in Eq. (V) can be calculated using the
quantum regression theorem (QRT) [1, 64]. The QRT utilises the cyclicity of
the trace to rewrite the two-time expectation value as

$\lim_{t\to\infty}\left\langle\sigma_{+}(t+\tau)\sigma_{-}(t)\right\rangle_{p}=\text{Tr}\left[\sigma^{+}\Lambda_{p}(\tau)\right],$
(57)

where $\Lambda_{p}(\tau)=U_{0}^{\dagger}(\tau)\Lambda_{p}(0)U_{0}(\tau)$ is a
modified density operator with the initial state

$\Lambda_{p}(0)=\sigma^{-}U_{0}(\infty)\rho_{p}(0)U_{0}^{\dagger}(\infty),$
(58)

and, because $U_{0}(\tau)$ is the time evolution operator defined by the
polaron frame Hamiltonian (see Eq. (12)), the master equation we have derived
holds identically for the $\Lambda_{p}(\tau)$ operator but with the initial
condition defined by Eq. (58). 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. (V) captures the photonic sideband.
However, the QRT contains an implicit Born approximation [64] as well as any
approximations used in the derivation of the master equation.

Owing to the modified initial state in Eq. (58), the power of the polarisation
spectrum contains information about the steady state of the system,

$P=\int_{-\infty}^{\infty}\text{d}\omega\
I_{p,0}(\omega)=\pi\rho_{ee}(\infty).$ (59)

Additionally, by noting that the polarisation spectrum without the photon
sideband, denoted by $I_{p,\times}(\omega)$, is given by Eq. (V) with the
replacement $\exp[\phi(\tau)]\to 1$, one can show that

$\frac{\int_{-\infty}^{\infty}\text{d}\omega\
I_{p,\times}(\omega)}{P}=\kappa^{2}.$ (60)

Hence, $\kappa^{2}$ can be interpreted as the fractional emission into the
sideband. Remarkably, an integrated spectrum measurement can be used to
determine $\kappa$, which in turn provides an effective measurement of
$\Omega_{\Delta\Delta}$ and, therefore, of the strength of the
$|\tilde{\mathbf{d}}_{\Delta}|$ dipole moment. In measurements,
$I_{p,\times}(\omega)$ cannot be separated from $I_{p,0}(\omega)$, but as we
show in the following, $\int_{-\infty}^{\infty}\text{d}\omega\
I_{p,\times}(\omega)$ can be accurately approximated as the integral of the
measured spectrum over the Lorentzian parts of the dominant peaks [65].

We note that using Eq. (59) to determine $\rho_{ee}(\infty)$ requires the
spectrum to be measured in units of inverse energy such that $P$ is
dimensionless. This may be difficult for experiments that measure the spectrum
in, for example, Watts per area or counts of photodetection events. However,
the measurement of $\kappa^{2}$ in Eq. (60) is independent of the units used
to measure the spectrum.

In Fig. 7, we plot the polarisation spectrum calculated in the polaron frame
[$I_{p,0}(\omega)$] in dot-dashed green, in blue we plot the spectrum for the
same parameters but without the photonic sideband [$I_{p,\times}(\omega)]$
which serves to highlight the sideband, and in dashed black we plot the
polarisation spectrum in the absence of permanent dipoles
[$I_{p,\text{npd}}(\omega)$]. In panel (a)
$|\tilde{\mathbf{d}}_{\Delta}|=0.05$ and $V=0.05~{}\text{eV}$, and in (b)
$|\tilde{\mathbf{d}}_{\Delta}|=0.1$ and $V=0.25~{}\text{eV}$. Notably, the
photon sideband substantially changes the emission spectrum even for small
$|\tilde{\mathbf{d}}_{\Delta}|$. The widths of the peaks are noticeably larger
in panel (b) due to the permanent dipoles; however, the effect of the
permanent dipoles on the Lamb shifts is negligible.

Figure 7: Polarisation spectrum of the system. The total polarisation
$I_{p,0}(\omega)$ is shown in dot-dashed green, the polarisation spectrum
without the photon sideband $I_{p,\times}(\omega)$ in solid blue, and the
polarisation spectrum with $|\mathbf{d}_{\Delta}|=0$ (and other parameters
unchanged) $I_{p,\text{npd}}(\omega)$ is shown in dashed black. In (a)
$|\tilde{\mathbf{d}}_{\Delta}|=0.05$ and $V=0.05~{}\text{eV}$ and in (b)
$|\tilde{\mathbf{d}}_{\Delta}|=0.1$ and $|V|=0.25~{}\text{eV}$. Other
parameters are given in the caption at the end of Fig. 3. The insets show
shaded regions for which $\int_{-\infty}^{\infty}\text{d}\omega\
I_{p,\times}(\omega)\approx\int_{R}\text{d}\omega\ I_{p,0}(\omega)$ where $R$
is the domain of the shaded regions, which can be used in Eq. (60) to estimate
$\kappa$ from a measured spectrum. We use the secularised master equation to
produce the spectra, which eliminates a known inconsistency near to the
negative frequency peak and has no other impact.

We now introduce a procedure of directly obtaining essential information about
the system from the experimentally measurable spectrum, $I_{p,0}(\omega)$, and
compare the values obtained for the data shown in Fig. 7 to analytic values.
We show that it is possible to evaluate $\rho_{ee}(\infty)$ and $\kappa$ from
$I_{p,0}(\omega)$ using Eqs. (59) and (60) and, provided that the temperature
of the experiment is known, these values can be used to obtain $\epsilon$ and
$|V|$. If the spectral density is also known, it is also possible to obtain
$|\tilde{\mathbf{d}}_{\Delta}|$ from $\kappa$. The values of the parameters
obtained from the spectrum (methods explained afterwards) and the analytic
values are given in Table 1.

Param. | Analytic (a) | Meas. (a) | Analytic (b) | Meas. (b)
---|---|---|---|---
$\kappa^{2}$ | 0.962 | 0.953 | 0.855 | 0.832
$\rho_{ee}(\infty)$ | 0.120 | 0.119 | 0.136 | 0.113
$\epsilon/\text{eV}$ | 1.000 | 1.004 | 1.000 | 1.071
$|V|/\text{eV}$ | 0.050 | 0.027 | 0.250 | 0.164
$|\tilde{\mathbf{d}}_{\Delta}|$ | 0.050 | 0.050 | 0.100 | 0.108
Table 1: Values of parameters obtainable 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}(\omega)$ in Fig 7. The
left-hand column is for Fig. 7(a) and the right-hand column for Fig. 7(b).

As shown in Eq. (59), the steady state population of the excited state can
straightforwardly be obtained by integrating $I_{p,0}(\omega)$ over all
frequencies. $\kappa$ can be approximately found using Eq. (60) with
$\int_{-\infty}^{\infty}\text{d}\omega\
I_{p,\times}(\omega)\approx\int_{R}\text{d}\omega\ I_{p,0}(\omega)$ where $R$
are the frequency regions for which the dominant peaks in the spectrum are
approximately Lorentzian before the photon sideband begins. These regions are
shown by the shaded regions in the insets of Fig. 7. In principle, one could
increase the accuracy of the $\kappa^{2}$ estimation in Fig. 7(a) by
integrating over the other two peaks, but these contributions are small. If
the spectral density of the photon bath is known, the relation
$\kappa=\exp[-\phi(0)/2]$ can be inverted to obtain
$|\tilde{\mathbf{d}}_{\Delta}|$. The remaining parameters of the system,
$\epsilon$ and $|V|$, can be estimated by assuming that the system thermalises
with respect to $H_{S}=\eta\tau_{z}/2$, yielding the steady state

$\rho_{ee}(\infty)=\frac{1}{2}\left[1-\frac{\epsilon}{\eta}\tanh\left(\frac{\beta\eta}{2}\right)\right].$
(61)

Using Eq. (61), the measured value of $\rho_{ee}(\infty)$, and assuming
negligible Lamb shifts such that the frequency of the peak in the spectrum
$\bar{\eta}\approx\eta$, one can estimate $\epsilon$. Then, using
$\bar{\eta}\approx\eta=(\epsilon^{2}+4\kappa^{2}|V|^{2})^{1/2}$, one can
estimate $|V|$. The estimations of $|V|$ in Table 1 are the least accurate
owing to the sensitivity of both Eq. (61) and $\eta$ to inaccuracies in
$\rho_{ee}(\infty)$ and $\epsilon$, respectively.

##  VI Importance of initial conditions

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

$\rho(0)=\rho_{S}(0)\otimes\rho_{E},$ (62)

where the environment is assumed to be in a thermal Gibbs state
$\rho_{E}=\exp(-\beta H_{E})/\mathcal{Z}_{E}$ with
$H_{E}=\sum_{k}\nu_{\mathbf{k}}a_{k}^{\dagger}a_{k}$. This assumption is valid
in weakly coupled systems where the system and the environment are only weakly
correlated even after thermalisation. However, for strongly coupled systems,
this assumption is invalidated, and the total Gibbs state
$\rho_{\beta}=\exp(-\beta H)/\mathcal{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
$\rho_{S}(0)\otimes\rho_{E}$ in the polaron frame, models an initial state in
the lab frame that is more similar to $\rho_{\beta}$. In our calculations,
$\rho_{S}(0)=\left|g\right>\left<g\right|$, and so
$\rho_{S}(0)\otimes\rho_{E}$ in the polaron frame becomes
$\rho_{l}(0)=\left|g\right>\left<g\right|\otimes
B([D+\Delta]/\nu)\rho_{E}B(-[D+\Delta]/\nu)$ in the lab frame, which is equal
to $\rho_{\beta}$ 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 thermalised 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. (II) with the initial state $\rho_{l}(0)$.
However, most numerical techniques, including TEMPO, assume that the initial
state is $\rho_{S}(0)\otimes\rho_{E}$. In Appendix G, we show that 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
$\rho_{S}(0)=\left|g\right>\left<g\right|$ in the polaron frame, the required
Hamiltonian is $\bar{H}=\bar{H}_{0}+\bar{H}_{I}$, where

$\displaystyle\bar{H}_{0}$
$\displaystyle=\frac{\tilde{\epsilon}}{2}\sigma_{z}+\tilde{V}\sigma_{+}+\tilde{V}^{*}\sigma_{-}+\sum_{k}\nu_{\mathbf{k}}a_{k}^{\dagger}a_{k},$
(63a) $\displaystyle\bar{H}_{I}$
$\displaystyle=A_{\Delta\Delta}\left(\mathcal{I}+\sigma^{z}\right)+\pi_{\mu\bar{\mu}}\sigma_{+}+\pi_{\bar{\mu}\mu}\sigma_{-},$
(63b)

and $\pi_{pq}$ is given in Eq. (6),
$\tilde{\epsilon}=\epsilon+2G_{\Delta\Delta}$, $\tilde{V}=V+G_{\mu\bar{\mu}}$,
where
$G_{pq}=\sum_{k}(p_{k}\Delta_{k}^{*}+q_{k}^{*}\Delta_{k})/\nu_{\mathbf{k}}$.
We note that because the displacement direction of the polaron transformation
is state dependent, the necessary effective Hamiltonian depends on
$\rho_{S}(0)$.

By comparing the Hamiltonians in Eqs. (63) and Eq. (II), we see stark
differences if both are assumed to have the initial state
$\rho_{S}(0)\otimes\rho_{E}$. This includes a renormalisation of the
transition energy $\epsilon\to\tilde{\epsilon}$ 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 $\pi_{DD}\mathcal{I}^{\prime}$ in Eq.
(II), have a similar effect. As we showed in Eq. (II), this 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.

##  VII 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 polarisation of photonic modes caused by
the permanent dipoles, allows 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. III, we showed that transition
and decoherence rates can be engineered for practical application by
exploiting permanent dipoles, for example, to reduce decoherence. 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 dipole are: a photon sideband with a relative
contribution to the emission spectrum scaling as $\kappa^{2}$, two-photon
processes scaling as $\Omega_{\mu\Delta}^{2}$, and a term linear in the
driving amplitude. The linear driving term is particularly useful for
engineering enhanced or suppressed rates by varying the driving phase. (2) In
Sec. IV, 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. V,
we indicated distinguishable features of permanent dipoles in emission spectra
and described possible measurements to obtain $\rho_{ee}(\infty)$, $\kappa$,
$\bar{\eta}$ and by extension the bare energy splitting $\epsilon$ and the
permanent dipole magnitude $|\tilde{\mathbf{d}}_{\Delta}|$. The
distinguishable features are the photonic sideband and the altered width and
position of peaks.

Furthermore, two of the processes we identified, namely the
$\Gamma_{2}(\omega)$ and $\Gamma_{V,1}(\omega)$ processes, originate from the
non-commutativity of the bosonic operators on the transition and permanent
dipole interactions. These vanish if $\mathbf{d}_{\mu}$ and
$\mathbf{d}_{\Delta}$ are perpendicular in which case our model becomes the
driven spin boson model. We have also shown that in the absence of driving,
this model maps onto the independent boson model in the limit of large cut off
frequency ($\nu_{c}\to\infty$) or small Huang-Rhys parameter ($S\to 0$) where
two photon processes are suppressed by the prefactor of the spectral density
in Eq. (23).

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 [66], or anisotropic [67], dielectric medium,
common to many biological systems. 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 [48, 68]. How permanent dipoles affect these phenomena
is still an open question. Additionally, it has been shown in Refs. [28, 30,
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.

## Acknowledgements

We would like to thank the Tempo Collaboration for use of the open-source code
Oqupy [53]. 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 centre grant number DS-2017-079, and the EPSRC (United Kingdom)
Strategic Equipment Grant No. EP/L02263X/1 (EP/M008576/1) and EPSRC (United
Kingdom) Grant EP/M027791/1 awards to M.F. D.M.R. is supported by EPSRC
(United Kingdom) Grant EP/T517896/1.

## References

  * Breuer  _et al._ [2002] H.-P. Breuer, F. Petruccione, _et al._ ,  _The Theory of Open Quantum Systems_ (Oxford University Press on Demand, 2002).
  * Ficek and Swain [2005] Z. Ficek and S. Swain,  _Quantum Interference and Coherence: Theory and Experiments_ , Vol. 100 (Springer Science & Business Media, 2005).
  * Leggett  _et al._ [1987] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. Fisher, A. Garg, and W. Zwerger,  _Dynamics of the Dissipative Two-State System,_ Reviews of Modern Physics 59, 1 (1987).
  * Agarwal [1974] G. S. Agarwal,  _Quantum Statistical Theories of Spontaneous Emission and Their Relation to Other Approaches,_ Quantum Optics , 1 (1974).
  * Chung  _et al._ [2016] P.-H. Chung, C. Tregidgo, and K. Suhling,  _Determining a Fluorophore’s Transition Dipole Moment from Fluorescence Lifetime Measurements in Solvents of Varying Refractive Index,_ Methods and Applications in Fluorescence 4, 045001 (2016).
  * Filippi  _et al._ [2012] C. Filippi, F. Buda, L. Guidoni, and A. Sinicropi,  _Bathochromic Shift in Green Fluorescent Protein: A Puzzle for QM/MM Approaches,_ Journal of Chemical Theory and Computation 8, 112 (2012).
  * Kovarskiĭ and Prepelitsa [2001] V. Kovarskiĭ and O. Prepelitsa,  _Effect of a Polar Environment on the Resonant Generation of Higher Optical Harmonics by Dipole Molecules,_ Optics and Spectroscopy 90, 351 (2001).
  * Deiglmayr  _et al._ [2010] J. Deiglmayr, A. Grochola, M. Repp, O. Dulieu, R. Wester, and M. Weidemüller,  _Permanent Dipole Moment of LiCs in the Ground State,_ Physical Review A 82, 032503 (2010).
  * Guérout  _et al._ [2010] R. Guérout, M. Aymar, and O. Dulieu,  _Ground State of the Polar Alkali-Metal-Atom–Strontium Molecules: Potential Energy Curve and Permanent Dipole Moment,_ Physical Review A 82, 042508 (2010).
  * Lin and Boxer [2020] C.-Y. Lin and S. G. Boxer,  _Mechanism of Color and Photoacidity Tuning for the Protonated Green Fluorescent Protein Chromophore,_ Journal of the American Chemical Society 142, 11032 (2020).
  * Jagatap and Meath [2002] B. Jagatap and W. J. Meath,  _Contributions of Permanent Dipole Moments to Molecular Multiphoton Excitation Cross Sections,_ JOSA B 19, 2673 (2002).
  * Gilmore and McKenzie [2005] J. Gilmore and R. H. McKenzie,  _Spin Boson Models for Quantum Decoherence of Electronic Excitations of Biomolecules and Quantum Dots in a Solvent,_ Journal of Physics: Condensed Matter 17, 1735 (2005).
  * Garziano  _et al._ [2016] L. Garziano, V. Macrì, R. Stassi, O. Di Stefano, F. Nori, and S. Savasta,  _One Photon can Simultaneously Excite Two or More Atoms,_ Physical Review Letters 117, 043601 (2016).
  * Chestnov  _et al._ [2016] I. Y. Chestnov, V. Shakhnazaryan, I. A. Shelykh, and A. P. Alodjants,  _Ensemble of Asymmetric Quantum Dots in a Cavity as a Terahertz Laser Source,_ JETP letters 104, 169 (2016).
  * Shim and Guyot-Sionnest [1999] M. Shim and P. Guyot-Sionnest,  _Permanent Dipole Moment and Charges in Colloidal Semiconductor Quantum Dots,_ The Journal of Chemical Physics 111, 6955 (1999).
  * Antón  _et al._ [2016] M. Antón, F. Carreño, O. Calderón, S. Melle, and E. Cabrera,  _Radiation Emission from an Asymmetric Quantum Dot Coupled to a Plasmonic Nanostructure,_ Journal of Optics 18, 025001 (2016).
  * Fry  _et al._ [2000a] P. Fry, I. Itskevich, D. Mowbray, M. Skolnick, J. Barker, E. O’Reilly, L. Wilson, P. Maksym, M. Hopkinson, M. Al-Khafaji, _et al._ ,  _Photocurrent Spectroscopy of InAs/GaAs Self-Assembled Quantum Dots: Observation of a Permanent Dipole Moment,_ Physica E: Low-dimensional Systems and Nanostructures 7, 408 (2000a).
  * Chestnov  _et al._ [2017] I. Y. Chestnov, V. A. Shahnazaryan, A. P. Alodjants, and I. A. Shelykh,  _Terahertz Lasing in Ensemble of Asymmetric Quantum Dots,_ ACS Photonics 4, 2726 (2017).
  * Fry  _et al._ [2000b] P. W. Fry, I. E. Itskevich, D. J. Mowbray, M. S. Skolnick, J. J. Finley, J. A. Barker, E. P. O’Reilly, L. R. Wilson, I. A. Larkin, P. A. Maksym, M. Hopkinson, M. Al-Khafaji, J. P. R. David, A. G. Cullis, G. Hill, and J. C. Clark,  _Inverted Electron-Hole Alignment in InAs-GaAs Self-Assembled Quantum Dots,_ Physical Review Letters 84, 733 (2000b).
  * Antón  _et al._ [2017] M. A. Antón, S. Maede-Razavi, F. Carreño, I. Thanopulos, and E. Paspalakis,  _Optical and Microwave Control of Resonance Fluorescence and Squeezing Spectra in a Polar Molecule,_ Physical Review A 96, 063812 (2017).
  * Patane  _et al._ [2000] A. Patane, A. Levin, A. Polimeni, F. Schindler, P. Main, L. Eaves, and M. Henini,  _Piezoelectric Effects in In 0.5 Ga 0.5 As Self-Assembled Quantum Dots Grown on (311) B GaAs Substrates,_ Applied Physics Letters 77, 2979 (2000).
  * Warburton  _et al._ [2002] R. J. Warburton, C. Schulhauser, D. Haft, C. Schäflein, K. Karrai, J. M. Garcia, W. Schoenfeld, and P. M. Petroff,  _Giant Permanent Dipole Moments of Excitons in Semiconductor Nanostructures,_ Physical Review B 65, 113303 (2002).
  * Ostapenko  _et al._ [2010] I. A. Ostapenko, G. Hönig, C. Kindel, S. Rodt, A. Strittmatter, A. Hoffmann, and D. Bimberg,  _Large Internal Dipole Moment in InGaN/GaN Quantum Dots,_ Applied Physics Letters 97, 063103 (2010).
  * Li and Alivisatos [2003] L.-s. Li and A. P. Alivisatos,  _Origin and Scaling of the Permanent Dipole Moment in CdSe Nanorods,_ Physical Review Letters 90, 097402 (2003).
  * Gupta  _et al._ [2006] S. Gupta, Q. Zhang, T. Emrick, and T. P. Russell,  _“Self-Corralling” Nanorods Under an Applied Electric Field,_ Nano Letters 6, 2066 (2006).
  * Mohammadimasoudi  _et al._ [2016] M. Mohammadimasoudi, Z. Hens, and K. Neyts,  _Full Alignment of Dispersed Colloidal Nanorods by Alternating Electric Fields,_ RSC Advances 6, 55736 (2016).
  * Yoshihara  _et al._ [2017] F. Yoshihara, T. Fuse, S. Ashhab, K. Kakuyanagi, S. Saito, and K. Semba,  _Superconducting Qubit–Oscillator Circuit Beyond the Ultrastrong-Coupling Regime,_ Nature Physics 13, 44 (2017).
  * Guarnieri  _et al._ [2018] G. Guarnieri, M. Kolář, and R. Filip,  _Steady-State Coherences by Composite System-Bath Interactions,_ Physical Review Letters 121, 070401 (2018).
  * Greenberg [2007] Y. S. Greenberg,  _Low-Frequency Rabi Spectroscopy of Dissipative Two-Level Systems: Dressed-State Approach,_ Physical Review B 76, 104520 (2007).
  * Román-Ancheyta  _et al._ [2021] R. Román-Ancheyta, M. Kolář, G. Guarnieri, and R. Filip,  _Enhanced Steady-State Coherence Via Repeated System-Bath Interactions,_ Physical Review A 104, 062209 (2021).
  * Purkayastha  _et al._ [2020] A. Purkayastha, G. Guarnieri, M. T. Mitchison, R. Filip, and J. Goold,  _Tunable Phonon-Induced Steady-State Coherence in a Double-Quantum-Dot Charge Qubit,_ NPJ Quantum Information 6, 1 (2020).
  * Macovei  _et al._ [2015] M. Macovei, M. Mishra, and C. H. Keitel,  _Population Inversion in Two-Level Systems Possessing Permanent Dipoles,_ Physical Review A 92, 013846 (2015).
  * Mirzac  _et al._ [2021] A. Mirzac, S. Carlig, and M. A. Macovei,  _Microwave Multiphoton Conversion Via Coherently Driven Permanent Dipole Systems,_ Physical Review A 103, 043719 (2021).
  * Mandal  _et al._ [2020a] A. Mandal, S. Montillo Vega, and P. Huo,  _Polarized Fock states and the dynamical Casimir effect in molecular cavity quantum electrodynamics,_ The Journal of Physical Chemistry Letters 11, 9215 (2020a).
  * Antón  _et al._ [2020] M. Antón, I. Gonzalo, and F. Carreño,  _Bichromatically-Controlled Entanglement Between Asymmetric Quantum Dots in a Photonic Cavity,_ Journal of Physics B: Atomic, Molecular and Optical Physics 54, 015504 (2020).
  * Oster  _et al._ [2012] F. Oster, C. H. Keitel, and M. Macovei,  _Generation of Correlated Photon Pairs in Different Frequency Ranges,_ Physical Review A 85, 063814 (2012).
  * JJuzeliunas  _et al._ [2003] G. JJuzeliunas, L. D. Romero, and D. L. Andrews,  _Eliminating Ground-State Dipole Moments in Quantum Optics Via Canonical Transformation,_ Physical Review A 68, 043811 (2003).
  * Paspalakis  _et al._ [2013] E. Paspalakis, J. Boviatsis, and S. Baskoutas,  _Effects of Probe Field Intensity in Nonlinear Optical Processes in Asymmetric Semiconductor Quantum Dots,_ Journal of Applied Physics 114, 153107 (2013).
  * Zhao  _et al._ [2016] J.-Y. Zhao, L.-G. Qin, X.-M. Cai, Q. Lin, and Z.-Y. Wang,  _The Effect of a Permanent Dipole Moment on the Polar Molecule Cavity Quantum Electrodynamics,_ Chinese Physics B 25, 044202 (2016).
  * Hattori and Kobayashi [1987] T. Hattori and T. Kobayashi,  _Bloch-Siegert Shift in Giant-Dipole Molecules,_ Physical Review A 35, 2733 (1987).
  * Scala  _et al._ [2021] G. Scala, K. Słowik, P. Facchi, S. Pascazio, and F. V. Pepe,  _Beyond the Rabi Model: Light Interactions with Polar Atomic Systems in a Cavity,_ Physical Review A 104, 013722 (2021).
  * Mandal  _et al._ [2020b] A. Mandal, S. M. Vega, and P. Huo,  _Polarized Fock States for Polariton Photochemistry,_ arXiv preprint arXiv:2005.00201 (2020b).
  * Mandal  _et al._ [2023] A. Mandal, M. Taylor, and P. Huo,  _A Theory for Cavity Modified Ground-State Reactivities via Electron-Photon Interactions,_ (2023).
  * Strathearn  _et al._ [2017] A. Strathearn, B. W. Lovett, and P. Kirton,  _Efficient Real-Time Path Integrals for Non-Markovian Spin-Boson Models,_ New Journal of Physics 19, 093009 (2017).
  * Strathearn  _et al._ [2018] A. Strathearn, P. Kirton, D. Kilda, J. Keeling, and B. W. Lovett,  _Efficient Non-Markovian Quantum Dynamics Using Time-Evolving Matrix Product Operators,_ Nature Communications 9, 3322 (2018).
  * Pollock  _et al._ [2018] F. A. Pollock, C. Rodríguez-Rosario, T. Frauenheim, M. Paternostro, and K. Modi,  _Non-Markovian Quantum Processes: Complete Framework and Efficient Characterization,_ Physical Review A 97, 012127 (2018).
  * Jørgensen and Pollock [2019] M. R. Jørgensen and F. A. Pollock,  _Exploiting the Causal Tensor Network Structure of Quantum Processes to Efficiently Simulate Non-Markovian Path Integrals,_ Physical Review Letters 123, 240602 (2019).
  * Gribben  _et al._ [2022a] D. Gribben, D. M. Rouse, J. Iles-Smith, A. Strathearn, H. Maguire, P. Kirton, A. Nazir, E. M. Gauger, and B. W. Lovett,  _Exact Dynamics of Nonadditive Environments in Non-Markovian Open Quantum Systems,_ PRX Quantum 3, 010321 (2022a).
  * Fux  _et al._ [2021] G. E. Fux, E. P. Butler, P. R. Eastham, B. W. Lovett, and J. Keeling,  _Efficient Exploration of Hamiltonian Parameter Space for Optimal Control of Non-Markovian Open Quantum Systems,_ Physical Review Letters 126, 200401 (2021).
  * Gribben  _et al._ [2022b] D. Gribben, A. Strathearn, G. E. Fux, P. Kirton, and B. W. Lovett,  _Using the Environment to Understand Non-Markovian Open Quantum Systems,_ Quantum 6, 847 (2022b).
  * Fux  _et al._ [2022] G. E. Fux, D. Kilda, B. W. Lovett, and J. Keeling,  _Thermalization of a Spin Chain Strongly Coupled to its Environment,_ arXiv preprint, arXiv:2201.05529 (2022).
  * Fowler-Wright  _et al._ [2022] P. Fowler-Wright, B. W. Lovett, and J. Keeling,  _Efficient Many-Body Non-Markovian Dynamics of Organic Polaritons,_ Physical Review Letters 129, 173001 (2022).
  * oqu [2020] _The TEMPO collaboration, OQuPy: A Python 3 Package to Efficiently Compute Non-Markovian Open Quantum Systems,_ http://dx.doi.org/10.5281/zenodo.4428316. (2020).
  * Babiker and Loudon [1983] M. Babiker and R. Loudon,  _Derivation of the Power-Zienau-Woolley Hamiltonian in Quantum Electrodynamics by Gauge Transformation,_ Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 385, 439 (1983).
  * Stokes and Nazir [2018] A. Stokes and A. Nazir,  _A Master Equation for Strongly Interacting Dipoles,_ New Journal of Physics 20, 043022 (2018).
  * Mahan [2013] G. D. Mahan,  _Many-Particle Physics_ (Springer Science & Business Media, 2013).
  * Rouse  _et al._ [2022] D. M. Rouse, E. M. Gauger, and B. W. Lovett,  _Analytic Expression for the Optical Exciton Transition Rates in the Polaron Frame,_ Physical Review B 105, 014302 (2022).
  * Gaspard and Nagaoka [1999] P. Gaspard and M. Nagaoka,  _Slippage of Initial Conditions for the Redfield Master Equation,_ The Journal of Chemical Physics 111, 5668 (1999).
  * Roy-Choudhury and Hughes [2015] K. Roy-Choudhury and S. Hughes,  _Quantum Theory of the Emission Spectrum from Quantum Dots Coupled to Structured Photonic Reservoirs and Acoustic Phonons,_ Physical Review B 92, 205406 (2015).
  * Iles-Smith  _et al._ [2017] J. Iles-Smith, D. P. McCutcheon, A. Nazir, and J. Mørk,  _Phonon Limit to Simultaneous Near-Unity Efficiency and Indistinguishability in Semiconductor Single Photon Sources,_ in  _CLEO: QELS_Fundamental Science_ (Optical Society of America, 2017) pp. FTu3E–3.
  * Hughes and Yao [2009] S. Hughes and P. Yao,  _Theory of Quantum Light Emission from a Strongly-Coupled Single Quantum Dot Photonic-Crystal Cavity system,_ Optics Express 17, 3322 (2009).
  * Bundgaard-Nielsen  _et al._ [2021] M. Bundgaard-Nielsen, J. Mørk, and E. V. Denning,  _Non-Markovian Perturbation Theories for Phonon Effects in Strong-Coupling Cavity Quantum Electrodynamics,_ Physical Review B 103, 235309 (2021).
  * Walls and Milburn [1994] D. F. Walls and G. J. Milburn,  _Quantum Optics (Springer, Berlin,_ (1994).
  * McCutcheon [2016] D. P. S. McCutcheon,  _Optical Signatures of Non-Markovian Behavior in Open Quantum Systems,_ Physical Review A 93, 022119 (2016).
  * Trebbia  _et al._ [2022] J.-B. Trebbia, Q. Deplano, P. Tamarat, and B. Lounis,  _Tailoring the Superradiant and Subradiant Nature of Two Coherently Coupled Quantum Emitters,_ Nature Communications 13, 1 (2022).
  * Burgess and Florescu [2022] A. Burgess and M. Florescu,  _Dynamical Decoherence and Memory Effects in Green Fluorescent Proteins by Dielectric Relaxation,_ arXiv preprint, arXiv:2211.09408 (2022).
  * Messinger  _et al._ [2020] A. Messinger, N. Westerberg, and S. M. Barnett,  _Spontaneous Emission in Anisotropic Dielectrics,_ Physical Review A 102, 013721 (2020).
  * Maguire  _et al._ [2019] H. Maguire, J. Iles-Smith, and A. Nazir,  _Environmental Nonadditivity and Franck-Condon Physics in Nonequilibrium Quantum Systems,_ Physical Review Letters 123, 093601 (2019).
  * Nazir and McCutcheon [2016] A. Nazir and D. P. McCutcheon,  _Modelling Exciton–Phonon Interactions in Optically Driven Quantum Dots,_ Journal of Physics: Condensed Matter 28, 103002 (2016).

##  Appendix A Displacement operators

Within the appendices, we will regularly use many identities involving
displacement operators which are proven in Ref. [69]. This first appendix is
dedicated to listing the necessary displacement operator identities. For the
purposes of this appendix, we use displacement operators with a single photon
mode

$B(\alpha)=\exp[\alpha a^{\dagger}-\alpha^{*}a],$ (64)

and note that $B(\alpha)^{\dagger}=B(-\alpha)$ and
$B(-\alpha)B(\alpha)=\mathcal{I}$. The first identity is the action of a
displacement operator on harmonic operators,

$\displaystyle B(\pm\alpha)a^{\dagger}B(\mp\alpha)$
$\displaystyle=a^{\dagger}\mp\alpha^{*},$ (65a) $\displaystyle
B(\pm\alpha)aB(\mp\alpha)$ $\displaystyle=a\mp\alpha.$ (65b)

The second is that the product of two displacement operators is

$B(\alpha_{1})B(\alpha_{2})=B(\alpha_{1}+\alpha_{2})\mathrm{e}^{\frac{1}{2}\left(\alpha_{1}\alpha_{2}^{*}-\alpha_{1}^{*}\alpha_{2}\right)}.$
(66)

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

$B(\alpha)\ket{0}=\ket{\alpha}.$ (67)

Fourth, that the expectation value of an operator with respect to the thermal
state $\rho_{E}=\exp[-\beta\nu a^{\dagger}a]/Z_{E}$ can be written as an
integral over coherent states as

$\text{Tr}_{E}[O\rho_{E}]=\frac{1}{\pi
N(\nu)}\int_{-\infty}^{\infty}d^{2}\alpha\
\mathrm{e}^{-\frac{\left|\alpha\right|^{2}}{N(\nu)}}\braket{\alpha\left
missing}{O}{\alpha},$ (68)

where
$\int_{-\infty}^{\infty}d^{2}\alpha=\int_{-\infty}^{\infty}d\Im[\alpha]\int_{-\infty}^{\infty}d\Re[\alpha]$
and $N(\nu)=(\mathrm{e}^{\beta\nu}-1)^{-1}$ is the Bose-Einstein distribution.
The final identity is the expectation value of a displacement operator with
respect to the vacuum state,

$\braket{0\left
missing}{B(\alpha)}{0}=\mathrm{e}^{-\frac{1}{2}\left|\alpha\right|^{2}}.$ (69)

##  Appendix B Calculation of $\left\langle C\right\rangle$

Recall from Eq. (13) that $C=B(\delta)A_{\mu\bar{\mu}}B(\delta)+VB(2\delta)$
where $\delta_{k}=\Delta_{k}/\nu_{\mathbf{k}}$. Using properties of
displacement operators this can be written as

$C=B(2\delta)\left(\sum_{k}\left[\mu_{k}\left(a_{k}^{\dagger}+\delta_{k}^{*}\right)+\bar{\mu}_{k}^{*}\left(a_{k}+\delta_{k}\right)\right]+|V|\right).$
(70)

In order to calculate $\langle C\rangle=\text{Tr}_{E}[C\rho_{E}]$ we require
$\langle B(2\delta)\rangle$, $\langle B(2\delta)a^{\dagger}_{k}\rangle$ and
$\langle B(2\delta)a_{k}\rangle$. 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(\alpha)=\exp[\alpha
a^{\dagger}-\alpha a]$, and reinstate the multimode summations at the end. To
calculate the expectation values we aim to use Eq. (68).

Starting with $\langle B(2\delta)\rangle$, we first evaluate the integrand of
Eq. (68) as

$\displaystyle\braket{\alpha\left missing}{B(2\delta)}{\alpha}$
$\displaystyle=\braket{0\left missing}{B(-\alpha)B(2\delta)B(\alpha)}{0}$
$\displaystyle=\braket{0\left
missing}{B(2\delta)}{0}\mathrm{e}^{2\delta\alpha^{*}-2\delta^{*}\alpha}$
$\displaystyle=\mathrm{e}^{-2\left|\delta\right|^{2}}\mathrm{e}^{2\delta\alpha^{*}-2\delta^{*}\alpha},$
(71)

where in the first equality we have used Eq. (67) twice, in the second
equality Eq. (66) twice and in the final equality Eq. (69). Substituting Eq.
(B) into Eq. (67) then yields

$\left\langle
B(2\delta)\right\rangle=\mathrm{e}^{-2\left|\delta\right|^{2}\coth\left(\frac{\beta\nu}{2}\right)}\equiv\kappa.$
(72)

Evaluating $\langle B(2\delta)a^{\dagger}\rangle$ follows a similar procedure.
However one must now use Eqs. (65) to move the creation operator such that it
annihilates with the vacuum bra-state $\bra{0}$. That is,

$\displaystyle\braket{\alpha\left missing}{B(2\delta)a^{\dagger}}{\alpha}$
$\displaystyle=\braket{0\left
missing}{B(-\alpha)B(2\delta)a^{\dagger}}{\alpha}$
$\displaystyle=\braket{0\left
missing}{B(-\alpha)\left[a^{\dagger}-2\delta^{*}\right]B(2\delta)}{\alpha}$
$\displaystyle=\left[\alpha^{*}-2\delta^{*}\right]\braket{\alpha\left
missing}{B(2\delta)}{\alpha}$
$\displaystyle=\left[\alpha^{*}-2\delta^{*}\right]\mathrm{e}^{-2\left|\delta\right|^{2}}\mathrm{e}^{2\delta\alpha^{*}-2\delta^{*}\alpha}.$
(73)

Substituting Eq. (B) into Eq. (68) and performing the integrations yields

$\langle
B(2\delta)a^{\dagger}\rangle=-2\kappa\delta^{*}\left[1+N(\nu)\right].$ (74)

The final expectation value, $\langle B(2\delta)a\rangle$, is easier to
calculate because $a$ annihilates with $\ket{0}$. One finds that

$\braket{\alpha\left
missing}{B(2\delta)a}{\alpha}=\alpha\mathrm{e}^{-2\left|\delta\right|^{2}}\mathrm{e}^{2\delta\alpha^{*}-2\delta^{*}\alpha},$
(75)

and so

$\left\langle B(2\delta)a\right\rangle=2\kappa\delta N(\nu).$ (76)

Collecting the expectation values in Eqs. (72), (74) and (76) and substituting
these into $\langle C\rangle$ with $C$ given in Eq. (70) yields

$\langle
C\rangle=\kappa\left[\coth\left(\frac{\beta\nu}{2}\right)\left(\bar{\mu}^{*}\delta-\mu\delta^{*}\right)+|V|\right].$
(77)

Since $\mathbf{d}_{\Delta}\in\Re$, we find that
$\bar{\mu}^{*}\delta-\mu\delta^{*}=0$, and so $\langle C\rangle=\kappa|V|$ as
stated in the main text.

##  Appendix C Non-secular master equations

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

$\displaystyle\partial_{t}\rho_{++}(t)$
$\displaystyle=-\gamma_{\downarrow}\rho_{++}(t)+\gamma_{\uparrow}\rho_{--}(t)+\bar{\gamma}\rho_{-+}(t)+\bar{\gamma}^{*}\rho_{+-}(t),$
(78a) $\displaystyle\partial_{t}\rho_{+-}(t)$
$\displaystyle=-[\gamma_{d}+i\bar{\eta}]\rho_{+-}(t)+k_{1}\rho_{-+}(t)+k_{-}\rho_{--}(t)+k_{+}^{*}\rho_{++}(t),$
(78b)

and $\partial_{t}\rho_{--}(t)=-\partial_{t}\rho_{++}(t)$ and
$\partial_{t}\rho_{-+}(t)=\partial_{t}\rho_{-+}(t)^{\dagger}$. As defined in
the main text, the secular rates are

$\displaystyle\gamma_{\begin{subarray}{c}\downarrow\\\
\uparrow\end{subarray}}$ $\displaystyle=\gamma_{\mp\mp}(\pm\eta),$ (79)
$\displaystyle\gamma_{d}$
$\displaystyle=\tfrac{1}{2}\left[\gamma_{\uparrow}+\gamma_{\downarrow}\right]+2\gamma_{zz}(0),$
(80) $\displaystyle\bar{\eta}$
$\displaystyle=\eta+\Im\left[\Gamma_{--}(\eta)-\Gamma_{++}(-\eta)\right].$
(81)

The rates are written in terms of the Markovian ECFs as

$\displaystyle\Gamma_{\alpha\beta}(\omega)=\int_{0}^{\infty}\text{d}s\
\mathrm{e}^{i\omega s}\left\langle
g_{\alpha}^{\dagger}(s)g_{\beta}(0)\right\rangle,$ (82)

for $\alpha\in\\{z,+,-\\}$, and it is the expressions for $g_{\alpha}$ that
vary between the PMFE and DFME, and these are given in the main text in Eqs.
(19) and Appendix E, respectively. The non-secular rates are

$\displaystyle\bar{\gamma}$ $\displaystyle=\Gamma_{-z}(0)+\Gamma_{+z}(0)^{*}$
(83a) $\displaystyle k_{1}$
$\displaystyle=\Gamma_{-+}(-\eta)+\Gamma_{+-}(\eta)^{*},$ (83b) $\displaystyle
k_{\pm}$ $\displaystyle=\mp\Gamma_{\pm z}(0)\pm\Gamma_{\mp z}(0)^{*}\pm
2\Gamma_{z\mp}(\pm\eta)$ (83c)

We evaluate $\Gamma_{\alpha\beta}(\omega)$ for the polaron frame in Appendix D
and for the displaced frame in Appendix E.

Notice that if the coherences are initially zero
($\rho_{+-}(0)=\rho_{-+}(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(\varphi)\to 0$, i.e. the eigenstates fully
localise, but this does not happen in the displaced frame in the same limit
because $g_{z}^{d}$ does not become zero.

##  Appendix D 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 $\langle
g^{\dagger}_{\alpha}(s)g_{\beta}(0)\rangle$ for $\alpha,\beta\in\\{z,+,-\\}$
where recalling from Eq. (19),

$\displaystyle g_{z}=$ $\displaystyle\
\frac{1}{2}\sin\left(\varphi\right)\left[\mathcal{C}+\mathcal{C}^{\dagger}\right],$
(84a) $\displaystyle g_{+}=$ $\displaystyle\
\left[\cos^{2}\left(\frac{\varphi}{2}\right)\mathcal{C}-\sin^{2}\left(\frac{\varphi}{2}\right)\mathcal{C}^{\dagger}\right],$
(84b)

and $g_{-}=g_{+}^{\dagger}$ where $\mathcal{C}=C-\kappa|V|$ and

$C=B_{+}\left(\sum_{k}\left[\mu\left(a^{\dagger}_{k}+\frac{\Delta^{*}_{k}}{\nu_{\mathbf{k}}}\right)+\bar{\mu}_{k}^{*}\left(a_{k}+\frac{\Delta_{k}}{\nu_{\mathbf{k}}}\right)\right]+|V|\right),$
(85)

and $B_{\pm}\equiv B(\pm 2\delta)$ with
$\delta_{k}=\Delta_{k}/\nu_{\mathbf{k}}$. The interaction picture form of the
relevant operators are,

$\displaystyle a_{k}(s)$
$\displaystyle=a_{k}\mathrm{e}^{-i\nu_{\mathbf{k}}s},$ (86) $\displaystyle
B_{\pm}(s)$ $\displaystyle=B\left(\pm 2\delta\mathrm{e}^{i\nu s}\right).$ (87)

###  D.1 Preliminary calculations

All $\langle g^{\dagger}_{\alpha}(s)g_{\beta}(0)\rangle$ depend on four two-
time correlation functions: $\langle C^{\dagger}(s)C(0)\rangle$, $\langle
C(s)C^{\dagger}(0)\rangle$, $\langle C(s)C(0)\rangle$, and $\langle
C^{\dagger}(s)C^{\dagger}(0)\rangle$. Each of these in turn depends on either
six or nine unique two-time correlation functions involving $a_{k}$,
$a_{k}^{\dagger}$ 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.

####  D.1.1 $\left\langle\mathcal{C}^{\dagger}(s)\mathcal{C}(0)\right\rangle$

To calculate this two-time correlation function, we require:

$\displaystyle\left\langle B_{-}(s)B_{+}(0)\right\rangle$
$\displaystyle=\kappa^{2}\mathrm{e}^{\phi},$ (88a) $\displaystyle\left\langle
a_{k}^{\dagger}B_{-}(s)B_{+}(0)\right\rangle$
$\displaystyle=-N_{\mathbf{k}}x_{k}^{*}\kappa^{2}\mathrm{e}^{\phi},$ (88b)
$\displaystyle\left\langle a_{k}B_{-}(s)B_{+}(0)\right\rangle$
$\displaystyle=\tilde{N}_{\mathbf{k}}x_{k}\kappa^{2}\mathrm{e}^{\phi},$ (88c)
$\displaystyle\left\langle B_{-}(s)B_{+}(0)a_{k}^{\dagger}\right\rangle$
$\displaystyle=-\tilde{N}_{\mathbf{k}}x_{k}^{*}\kappa^{2}\mathrm{e}^{\phi},$
(88d) $\displaystyle\left\langle B_{-}(s)B_{+}(0)a_{k}\right\rangle$
$\displaystyle=N_{\mathbf{k}}x_{k}\kappa^{2}\mathrm{e}^{\phi},$ (88e)
$\displaystyle\left\langle
a_{k}^{\dagger}B_{-}(s)B_{+}(0)a_{q}^{\dagger}\right\rangle$
$\displaystyle=N_{\mathbf{k}}\tilde{N}_{\mathbf{q}}x_{k}^{*}x_{q}^{*}\kappa^{2}\mathrm{e}^{\phi},$
(88f) $\displaystyle\left\langle a_{k}B_{-}(s)B_{+}(0)a_{q}\right\rangle$
$\displaystyle=\tilde{N}_{\mathbf{k}}N_{\mathbf{q}}x_{k}x_{q}\kappa^{2}\mathrm{e}^{\phi},$
(88g) $\displaystyle\left\langle
a_{k}^{\dagger}B_{-}(s)B_{+}(0)a_{q}\right\rangle$
$\displaystyle=\left(-N_{\mathbf{k}}N_{\mathbf{q}}x_{k}^{*}x_{q}+N_{\mathbf{k}}\delta_{kq}\right)\kappa^{2}\mathrm{e}^{\phi},$
(88h) $\displaystyle\left\langle
a_{k}B_{-}(s)B_{+}(0)a_{q}^{\dagger}\right\rangle$
$\displaystyle=\left(-\tilde{N}_{\mathbf{k}}\tilde{N}_{\mathbf{q}}x_{k}x_{q}^{*}+\tilde{N}_{\mathbf{k}}\delta_{kq}\right)\kappa^{2}\mathrm{e}^{\phi},$
(88i)

where, for brevity, $N_{\mathbf{k}}\equiv N(\nu_{\mathbf{k}})$ is the Bose-
Einstein distribution,
$\tilde{N}_{\mathbf{k}}\equiv\tilde{N}(\nu_{\mathbf{k}})=1+N(\nu_{\mathbf{k}})$,
$\phi\equiv\phi(s)$ given in Eq. (II), $\kappa^{2}=\exp[-\phi(0)]$, and
finally

$x_{k}\equiv
x_{k}(s)=2\delta_{k}\left(1-\mathrm{e}^{i\nu_{\mathbf{k}}s}\right).$ (89)

Substituting Eqs. (88) into the two-time correlation function we find that,

$\displaystyle\left\langle\mathcal{C}^{\dagger}(s)\mathcal{C}(0)\right\rangle=\kappa^{2}\mathrm{e}^{\phi(s)}\bigg{[}$
$\displaystyle\sum_{k}\left(\left|\bar{\mu}_{k}\right|^{2}N_{\mathbf{k}}\mathrm{e}^{i\nu_{\mathbf{k}}s}+\left|\mu_{\mathbf{k}}\right|^{2}\tilde{N}_{\mathbf{k}}\mathrm{e}^{-i\nu_{\mathbf{k}}s}\right)+\sum_{n,m\in\\{-1,1\\}}\sum_{kq}\chi_{nm}(k,q)\mathrm{e}^{in\nu_{\mathbf{k}}s}\mathrm{e}^{im\nu_{\mathbf{q}}s}$
$\displaystyle+\sum_{n\in\\{-1,1\\}}\sum_{k}\xi_{n}(k)\mathrm{e}^{in\nu_{\mathbf{k}}s}\bigg{]}+\kappa^{2}|V|^{2}\left(\mathrm{e}^{\phi(s)}-1\right).$
(90)

$\chi_{nm}(k,q)$ and $\xi_{n}(k)$ are written in matrix notation for brevity,
using the format

$\displaystyle\chi(k,q)$
$\displaystyle=\begin{pmatrix}\chi_{-1,-1}(k,q)&\chi_{-1,1}(k,q)\\\
\chi_{1,-1}(k,q)&\chi_{1,1}(k,q)\end{pmatrix},$ (91) $\displaystyle\xi(k)$
$\displaystyle=\\{\xi_{-1}(k),\xi_{1}(k)\\}.$ (92)

Written in this notation, one finds that

$\displaystyle\chi(k,q)$
$\displaystyle=4\frac{f_{\mathbf{k}}^{2}}{\nu_{\mathbf{k}}}\frac{f_{\mathbf{q}}^{2}}{\nu_{\mathbf{q}}}\mu_{k}^{0*}\mu_{q}^{0}\Delta_{k}^{0}\Delta_{q}^{0}\begin{pmatrix}\tilde{N}_{\mathbf{k}}\tilde{N}_{\mathbf{q}}&-\tilde{N}_{\mathbf{k}}N_{q}\\\
-N_{\mathbf{k}}\tilde{N}_{\mathbf{q}}&N_{\mathbf{k}}N_{\mathbf{q}}\end{pmatrix},$
(93a) $\displaystyle\xi(k)$
$\displaystyle=2|V|\frac{f_{k}^{2}}{\nu_{\mathbf{k}}}\Delta_{k}^{0}\left(\mu_{k}^{0}+\mu_{k}^{0*}\right)\\{\tilde{N}_{\mathbf{k}},-N_{\mathbf{k}}\\},$
(93b)

where we have written $\mu_{k}=i\mu_{k}^{0}f_{\mathbf{k}}$,
$\bar{\mu}_{k}=i\mu_{k}^{0*}f_{\mathbf{k}}$, and=
$\Delta_{k}=i\Delta_{k}^{0}f_{\mathbf{k}}$ with
$\mu_{k}^{0}=\mathbf{d}_{\mu}\cdot\mathbf{e}_{k}$ and
$\Delta_{k}^{0}=\mathbf{d}_{\Delta}\cdot\mathbf{e}_{k}$.

If one ignores the factor of $\kappa^{2}\exp[\phi(s)]$, the first term in Eq.
(D.1.1) will lead to the standard optical master equation rates describing
photon emission and absorption scaling with the magnitude squared of the
transition dipole strength. The second term in Eq. (D.1.1) describes
corrections to the standard optical master equation rates that are second
order in both the permanent dipole moments and the transition dipole moments.
The factor of $\kappa^{2}\exp[\phi(s)]$ accounts for higher order permanent
dipole contributions. The second and third terms in Eq. (D.1.1) along with
Eqs. (93) reveals the physical processes induced by the permanent dipoles. In
the term dependent on $\chi_{nm}(k,q)$, each term $n$ and $m$ describes two
simultaneous processes into the $k$ and $q$ modes, respectively. When $n$ or
$m$ equal $-1$ the processes are photon emission, when $n$ or $m$ equal $0$
the processes are non-radiative, and when $n$ or $m$ equal $+1$ the processes
are photon absorption. This can be seen by both the frequency dependence of
the exponential phase factors in Eq. (D.1.1) and by the dependence of the
matrix elements on $\tilde{N}_{\mathbf{k}}$, $N_{\mathbf{k}}$, or neither.
Similarly, the term dependent on $\xi_{n}(k)$ describes new one photon
emission ($n=-1$) and absorption ($n=+1$) processes induced by the combination
of the driving, permanent dipoles and transition dipoles.

Finally, we take the continuum limit of the photon wavenumber to obtain

$\displaystyle\left\langle\mathcal{C}^{\dagger}(s)\mathcal{C}(0)\right\rangle=$
$\displaystyle\kappa^{2}\mathrm{e}^{\phi(s)}\bigg{[}\Omega_{\mu\bar{\mu}}\int_{0}^{\infty}\text{d}\nu\
J(\nu)\left(N(\nu)\mathrm{e}^{i\nu s}+\tilde{N}(\nu)\mathrm{e}^{-i\nu
s}\right)+\cos^{2}\left(\theta_{\mu\Delta}\right)\sum_{n,m\in\\{-1,1\\}}\tilde{\chi}^{\mu\bar{\mu}}_{nm}(s)$
$\displaystyle+2|V|\cos(\theta_{\mu\Delta})\left(\Omega_{\mu\Delta}+\Omega_{\bar{\mu}\Delta}\right)\int_{0}^{\infty}\text{d}\nu\
\frac{J(\nu)}{\nu}\left(-N(\nu)\mathrm{e}^{i\nu
s}+\tilde{N}(\nu)\mathrm{e}^{-i\nu
s}\right)\bigg{]}+\kappa^{2}|V|^{2}\left(\mathrm{e}^{\phi(s)}-1\right),$ (94)

where we have defined the two-dimensional Fourier transform

$\tilde{\chi}^{ab}_{nm}(s)=\int_{0}^{\infty}\text{d}\nu\
\int_{0}^{\infty}\text{d}\nu^{\prime}\
\chi^{ab}_{nm}(\nu,\nu^{\prime})\mathrm{e}^{in\nu
s}\mathrm{e}^{im\nu^{\prime}s}.$ (95)

The continuum limit of the coefficients can be written in the matrix
representation as

$\chi^{ab}(\nu,\nu^{\prime})=4\frac{J(\nu)}{\nu}\frac{J(\nu^{\prime})}{\nu^{\prime}}\Omega_{a\Delta}\Omega_{b\Delta}\begin{pmatrix}\tilde{N}(\nu)\tilde{N}(\nu^{\prime})&-\tilde{N}(\nu)N(\nu^{\prime})\\\
-N(\nu)\tilde{N}(\nu^{\prime})&N(\nu)N(\nu^{\prime})\end{pmatrix},$ (96)

where $a,b\in\\{\mu,\bar{\mu}\\}$. We have introduced the $a,b$ superscripts
so that we may write all four two-time correlation functions of
$\mathcal{C}(s)$ in a unified notation. Note that
$\chi^{\mu\bar{\mu}}(\nu,\nu^{\prime})$ is the same as Eq. (45) in the main
text. The Fourier transform of the first term ($\propto\Omega_{\mu\bar{\mu}}$)
in Eq. (D.1.1) leads to $\Gamma_{1}(\omega)$ in the main text, of the second
term ($\propto\tilde{\chi}_{nm}^{\mu\bar{\mu}}(\nu,\nu^{\prime})$) leads to
$\Gamma_{2}(\omega)$, of the third term ($\propto|V|$) leads to
$\Gamma_{V,1}(\omega)$, and of the final term ($\propto|V|^{2}$) leads to
$\Gamma_{V,0}(\omega)$.

We now also introduce a more complicated notation for the Fourier transforms
of the two-time correlation functions so that we can write all four down in a
unified notation. Note that this notation differs again from the simpler
notation in the main text. In the new notation, the Fourier transform of
$\langle\mathcal{C}^{\dagger}(s)\mathcal{C}(0)\rangle$ is

$\Gamma^{(\dagger,\cdot)}(\omega)\equiv\int_{0}^{\infty}\text{d}\nu\
\mathrm{e}^{i\omega
s}\left\langle\mathcal{C}^{\dagger}(s)\mathcal{C}(0)\right\rangle=\Gamma_{1}^{\mu\bar{\mu},+}(\omega)+\Gamma_{2}^{\mu\bar{\mu},+}(\omega)+\Gamma_{V,1}^{\mu\bar{\mu},+}(\omega)+\Gamma_{V,0}^{+}(\omega),$
(97)

where we have defined

$\displaystyle\Gamma^{ab,\pm}_{1}(\omega)$
$\displaystyle=\pi\Omega_{ab}\int_{0}^{\infty}\text{d}\nu\
J(\nu)\left[N(\nu)\mathcal{K}^{\pm}(\omega+\nu)+\tilde{N}(\nu)\mathcal{K}^{\pm}(\omega-\nu)\right],$
(98a) $\displaystyle\Gamma^{ab,\pm}_{2}(\omega)$
$\displaystyle=\kappa^{2}\cos^{2}\left(\theta_{\mu\Delta}\right)\sum_{n,m\in\\{-1,1\\}}\int_{0}^{\infty}\text{d}s\
\mathrm{e}^{i\omega
s}\mathrm{e}^{\pm\phi(s)}\tilde{\chi}_{n,m}^{ab}(\nu,\nu^{\prime}),$ (98b)
$\displaystyle\Gamma^{ab,\pm}_{V,1}(\omega)$
$\displaystyle=4\pi\left(\Omega_{a\Delta}+\Omega_{b\Delta}\right)|V|\cos\left(\theta_{\mu\Delta}\right)\int_{0}^{\infty}\text{d}\nu\
\frac{J(\nu)}{\nu}\left[\tilde{N}(\nu)\mathcal{K}^{\pm}(\omega-\nu)-N(\nu)\mathcal{K}^{\pm}(\omega+\nu)\right],$
(98c) $\displaystyle\Gamma^{\pm}_{V,2}$
$\displaystyle=\kappa^{2}|V|^{2}\left(\mathrm{e}^{\pm\phi(s)}-1\right),$ (98d)

and

$\mathcal{K}^{\pm}(\varepsilon)=\frac{1}{\pi}\int_{0}^{\infty}\text{d}s\
\mathrm{e}^{i\varepsilon s}\mathrm{e}^{\pm\phi(s)-\phi(0)},$ (99)

is a generalisation of $\mathcal{K}(\varepsilon)$ in Eq. (34). Note that the
single mode truncation solution of $\mathcal{K}^{-}(\varepsilon)$ is obtained
from the solution for $\mathcal{K}^{+}(\varepsilon)$ in the main text by
replacing $W_{n}=S_{s}^{n}\exp[-S_{s}]/n!$ in Eq. (39) with
$W_{n}^{\prime}=(-S_{s})^{n}\exp[-S_{s}]/n!$ [57].

####  D.1.2 $\left\langle\mathcal{C}(s)\mathcal{C}^{\dagger}(0)\right\rangle$

To calculate this two-time correlation function, we require:

$\displaystyle\left\langle B_{+}(s)B_{-}(0)\right\rangle$
$\displaystyle=\kappa^{2}\mathrm{e}^{\phi},$ (100a) $\displaystyle\left\langle
B_{+}(s)a_{k}^{\dagger}B_{-}(0)\right\rangle$
$\displaystyle=\left(-2\delta_{k}^{*}-\tilde{N}_{\mathbf{k}}y_{k}^{*}\right)\kappa^{2}\mathrm{e}^{\phi},$
(100b) $\displaystyle\left\langle B_{+}(s)a_{k}B_{-}(0)\right\rangle$
$\displaystyle=\left(-2\delta_{k}+N_{\mathbf{k}}y_{k}\right)\kappa^{2}\mathrm{e}^{\phi},$
(100c) $\displaystyle\left\langle
B_{+}(s)a_{k}^{\dagger}a_{q}^{\dagger}B_{-}(0)\right\rangle$
$\displaystyle=\left(-2\delta_{k}^{*}-\tilde{N}_{\mathbf{k}}y_{k}^{*}\right)\left(-2\delta_{q}^{*}-\tilde{N}_{\mathbf{q}}y_{q}^{*}\right)\kappa^{2}\mathrm{e}^{\phi},$
(100d) $\displaystyle\left\langle B_{+}(s)a_{k}a_{q}B_{-}(0)\right\rangle$
$\displaystyle=\left(-2\delta_{k}+N_{\mathbf{k}}y_{k}\right)\left(-2\delta_{q}+N_{\mathbf{q}}y_{q}\right)\kappa^{2}\mathrm{e}^{\phi},$
(100e) $\displaystyle\left\langle
B_{+}(s)a_{k}^{\dagger}a_{q}B_{-}(0)\right\rangle$
$\displaystyle=\left[N_{\mathbf{k}}\delta_{kq}+\left(-2\delta_{k}^{*}-\tilde{N}_{\mathbf{k}}y_{k}^{*}\right)\left(-2\delta_{q}+N_{\mathbf{q}}y_{q}\right)\right]\kappa^{2}\mathrm{e}^{\phi},$
(100f) $\displaystyle\left\langle
B_{+}(s)a_{k}a_{q}^{\dagger}B_{-}(0)\right\rangle$
$\displaystyle=\left[\tilde{N}_{\mathbf{k}}\delta_{kq}+\left(-2\delta_{k}+N_{\mathbf{k}}y_{k}\right)\left(-2\delta_{q}^{*}-\tilde{N}_{\mathbf{q}}y_{q}^{*}\right)\right]\kappa^{2}\mathrm{e}^{\phi},$
(100g)

where

$y_{k}\equiv
y_{k}(s)=2\delta_{k}\left(\mathrm{e}^{i\nu_{\mathbf{k}}s}-1\right).$ (101)

Substituting Eqs. (100) into the two-time correlation function and taking the
continuum limit one obtains

$\displaystyle\left\langle\mathcal{C}(s)\mathcal{C}^{\dagger}(0)\right\rangle=\kappa^{2}$
$\displaystyle\mathrm{e}^{\phi(s)}\bigg{[}\Omega_{\mu\bar{\mu}}\int_{0}^{\infty}\text{d}\nu\
J(\nu)\left(N(\nu)\mathrm{e}^{i\nu s}+\tilde{N}(\nu)\mathrm{e}^{-i\nu
s}\right)+\cos^{2}(\theta_{\mu\Delta})\sum_{n,m\in\\{-1,1\\}}\tilde{\chi}_{nm}(s)$
$\displaystyle-2|V|\cos(\theta_{\mu\Delta})\left(\Omega_{\mu\Delta}+\Omega_{\bar{\mu}\Delta}\right)\int_{0}^{\infty}\text{d}\nu\
\frac{J(\nu)}{\nu}\left(-N(\nu)\mathrm{e}^{i\nu
s}+\tilde{N}(\nu)\mathrm{e}^{i\nu
s}\right)\bigg{]}+\kappa^{2}|V|^{2}\left(\mathrm{e}^{\phi(s)}-1\right),$ (102)

which is the same as Eq. (D.1.1) except for the overall minus sign on the
third term ($\propto V$). Written in terms of the rate functions defined in
Eqs. (98) the Fourier transform of this two-time correlation function is,

$\Gamma^{(\cdot,\dagger)}(\omega)\equiv\int_{0}^{\infty}\text{d}\nu\
\mathrm{e}^{i\omega
s}\left\langle\mathcal{C}(s)\mathcal{C}^{\dagger}(0)\right\rangle=\Gamma_{1}^{\mu\bar{\mu},+}(\omega)+\Gamma_{2}^{\mu\bar{\mu},+}(\omega)-\Gamma_{V,1}^{\mu\bar{\mu},+}(\omega)+\Gamma_{V,0}^{+}(\omega).$
(103)

####  D.1.3
$\left\langle\mathcal{C}^{\dagger}(s)\mathcal{C}^{\dagger}(0)\right\rangle$

To calculate this two-time correlation function, we require:

$\displaystyle\left\langle B_{-}(s)B_{-}(0)\right\rangle$
$\displaystyle=\kappa^{2}\mathrm{e}^{-\phi},$ (104a)
$\displaystyle\left\langle a_{k}^{\dagger}B_{-}(s)B_{-}(0)\right\rangle$
$\displaystyle=-N_{\mathbf{k}}z_{k}^{*}\kappa^{2}\mathrm{e}^{-\phi},$ (104b)
$\displaystyle\left\langle a_{k}B_{-}(s)B_{-}(0)\right\rangle$
$\displaystyle=\tilde{N}_{\mathbf{k}}z_{k}\kappa^{2}\mathrm{e}^{-\phi},$
(104c) $\displaystyle\left\langle
B_{-}(s)a_{k}^{\dagger}B_{-}(0)\right\rangle$
$\displaystyle=\left(-\tilde{N}_{\mathbf{k}}z_{k}^{*}-2\delta_{k}^{*}\right)\kappa^{2}\mathrm{e}^{-\phi},$
(104d) $\displaystyle\left\langle B_{-}(s)a_{k}B_{-}(0)\right\rangle$
$\displaystyle=\left(N_{\mathbf{k}}z_{k}-2\delta_{k}\right)\kappa^{2}\mathrm{e}^{-\phi},$
(104e) $\displaystyle\left\langle
a_{k}^{\dagger}B_{-}(s)a_{q}^{\dagger}B_{-}(0)\right\rangle$
$\displaystyle=-N_{\mathbf{k}}z_{k}^{*}\left(-\tilde{N}_{\mathbf{q}}z_{q}^{*}-2\delta_{q}^{*}\right)\kappa^{2}\mathrm{e}^{-\phi},$
(104f) $\displaystyle\left\langle a_{k}B_{-}(s)a_{q}B_{-}(0)\right\rangle$
$\displaystyle=\tilde{N}_{\mathbf{k}}z_{k}\left(N_{\mathbf{q}}z_{q}-2\delta_{q}\right)\kappa^{2}\mathrm{e}^{-\phi},$
(104g) $\displaystyle\left\langle
a_{k}^{\dagger}B_{-}(s)a_{q}B_{-}(0)\right\rangle$
$\displaystyle=\left[N_{\mathbf{k}}\delta_{kq}-N_{\mathbf{k}}z_{k}\left(N_{\mathbf{q}}z_{q}-2\delta_{q}\right)\right]\kappa^{2}\mathrm{e}^{-\phi},$
(104h) $\displaystyle\left\langle
a_{k}B_{-}(s)a_{q}^{\dagger}B_{-}(0)\right\rangle$
$\displaystyle=\left[\tilde{N}_{\mathbf{k}}\delta_{kq}+\tilde{N}_{\mathbf{k}}z_{k}\left(N_{\mathbf{q}}z_{q}-2\delta_{q}\right)\right]\kappa^{2}\mathrm{e}^{-\phi},$
(104i)

where

$z_{k}\equiv
z_{k}(s)=-2\delta_{k}\left(\mathrm{e}^{i\nu_{\mathbf{k}}s}+1\right).$ (105)

Substituting Eqs. (104) into the two-time correlation function and taking the
continuum limit one obtains

$\displaystyle\left\langle\mathcal{C}^{\dagger}(s)\mathcal{C}^{\dagger}(0)\right\rangle=\kappa^{2}\mathrm{e}^{-\phi(s)}\left[\Omega_{\bar{\mu}\bar{\mu}}\int_{0}^{\infty}\text{d}\nu\
J(\nu)\left(N(\nu)\mathrm{e}^{i\nu s}+\tilde{N}(\nu)\mathrm{e}^{-i\nu
s}\right)-\sum_{n,m\in\\{-1,1\\}}\tilde{\chi}^{\bar{\mu}\bar{\mu}}_{nm}(s)\right]+\kappa^{2}|V|^{2}\left(\mathrm{e}^{-\phi(s)}-1\right),$
(106)

where $\tilde{\chi}_{nm}^{ab}(s)$ is defined in Eqs. (95) and (96). Eq. (106)
has only a few differences to the two-time correlation functions
$\langle\mathcal{C}^{\dagger}(s)\mathcal{C}(0)\rangle$ and
$\langle\mathcal{C}(s)\mathcal{C}^{\dagger}(0)\rangle$ derived earlier. Most
significantly, Eq. (106) has $\phi(s)\to-\phi(s)$, it does not have a term
proportional to $|V|\cos(\theta_{\mu\Delta})$, and there is a minus sign on
the second term. The final difference is that Eq. (106) only depends on
$\mathbf{d}_{\bar{\mu}}=\mathbf{d}_{\mu}^{*}$.

Using the functions defined in Eqs. (98) we can write the Fourier transform of
Eq. (106) as

$\Gamma^{(\dagger,\dagger)}(\omega)\equiv\int_{0}^{\infty}\text{d}\nu\
\mathrm{e}^{i\omega
s}\left\langle\mathcal{C}^{\dagger}(s)\mathcal{C}^{\dagger}(0)\right\rangle=\Gamma_{1}^{\bar{\mu}\bar{\mu},-}(\omega)-\Gamma_{2}^{\bar{\mu}\bar{\mu},-}(\omega)+\Gamma_{V,0}^{-}(\omega).$
(107)

####  D.1.4 $\left\langle\mathcal{C}(s)\mathcal{C}(0)\right\rangle$

To calculate this two-time correlation function we require:

$\displaystyle\left\langle B_{+}(s)B_{+}(0)\right\rangle$
$\displaystyle=\kappa^{2}\mathrm{e}^{-\phi},$ (108a)
$\displaystyle\left\langle B_{+}(s)a_{k}^{\dagger}B_{+}(0)\right\rangle$
$\displaystyle=\left(z_{k}^{*}\tilde{N}_{\mathbf{k}}+2\delta_{k}^{*}\right)\kappa^{2}\mathrm{e}^{-\phi},$
(108b) $\displaystyle\left\langle B_{+}(s)a_{k}B_{+}(0)\right\rangle$
$\displaystyle=\left(-z_{k}N_{\mathbf{k}}+2\delta_{k}\right)\kappa^{2}\mathrm{e}^{-\phi},$
(108c) $\displaystyle\left\langle
B_{+}(s)B_{+}(0)a_{k}^{\dagger}\right\rangle$
$\displaystyle=\tilde{N}_{\mathbf{k}}z_{k}^{*}\kappa^{2}\mathrm{e}^{-\phi},$
(108d) $\displaystyle\left\langle B_{+}(s)B_{+}(0)a_{k}\right\rangle$
$\displaystyle=-N_{\mathbf{k}}z_{k}\kappa^{2}\mathrm{e}^{-\phi},$ (108e)
$\displaystyle\left\langle
B_{+}(s)a_{k}^{\dagger}B_{+}(0)a_{q}^{\dagger}\right\rangle$
$\displaystyle=\tilde{N}_{\mathbf{k}}z_{k}^{*}\left(\tilde{N}_{\mathbf{q}}z_{q}^{*}+2\delta_{q}^{*}\right)\kappa^{2}\mathrm{e}^{-\phi},$
(108f) $\displaystyle\left\langle B_{+}(s)a_{k}B_{+}(0)a_{q}\right\rangle$
$\displaystyle=N_{\mathbf{k}}z_{k}\left(N_{\mathbf{q}}z_{q}-2\delta_{q}\right)\kappa^{2}\mathrm{e}^{-\phi},$
(108g) $\displaystyle\left\langle
B_{+}(s)a_{k}^{\dagger}B_{+}(0)a_{q}\right\rangle$
$\displaystyle=\left[N_{\mathbf{k}}\delta_{kq}-\left(z_{k}^{*}\tilde{N}_{\mathbf{k}}+2\delta_{k}^{*}\right)N_{\mathbf{q}}z_{q}\right]\kappa^{2}\mathrm{e}^{-\phi},$
(108h) $\displaystyle\left\langle
B_{+}(s)a_{k}B_{+}(0)a_{q}^{\dagger}\right\rangle$
$\displaystyle=\left[\tilde{N}_{\mathbf{k}}\delta_{kq}+\left(-z_{k}N_{\mathbf{k}}+2\delta_{k}\right)\tilde{N}_{q}z_{q}^{*}\right]\kappa^{2}\mathrm{e}^{-\phi},$
(108i)

where $z_{k}\equiv z_{k}(s)$ is given in Eq. (105). Substituting Eqs. (108)
into the two-time correlation function and taking the continuum limit one
obtains

$\displaystyle\left\langle\mathcal{C}(s)\mathcal{C}(0)\right\rangle=\kappa^{2}\mathrm{e}^{-\phi(s)}\left[\Omega_{\mu\mu}\int_{0}^{\infty}\text{d}\nu\
J(\nu)\left(N(\nu)\mathrm{e}^{i\nu s}+\tilde{N}(\nu)\mathrm{e}^{-i\nu
s}\right)-\sum_{n,m\in\\{-1,1\\}}\tilde{\chi}^{\mu\mu}_{nm}(s)\right]+\kappa^{2}|V|^{2}\left(\mathrm{e}^{-\phi(s)}-1\right),$
(109)

where $\tilde{\chi}_{nm}^{ab}(s)$ is defined in Eqs. (95) and (96). Eq. (109)
is only different from Eq. (106) in that $\mathbf{d}_{\bar{\mu}}$ has been
replaced with $\mathbf{d}_{\mu}$.

Using the functions defined in Eqs. (98) we can write the Fourier transform of
Eq. (109) as

$\Gamma^{(\cdot,\cdot)}(\omega)\equiv\int_{0}^{\infty}\text{d}\nu\
\mathrm{e}^{i\omega
s}\left\langle\mathcal{C}(s)\mathcal{C}(0)\right\rangle=\Gamma_{1}^{\mu\mu,-}(\omega)-\Gamma_{2}^{\mu\mu,-}(\omega)+\Gamma_{V,0}^{-}(\omega).$
(110)

###  D.2 Two-time correlation functions $\left\langle
g_{\alpha}^{\dagger}(s)g_{\beta}(0)\right\rangle$

With the expression for
$\langle\mathcal{C}^{\dagger}(s)\mathcal{C}(0)\rangle$,
$\langle\mathcal{C}(s)\mathcal{C}^{\dagger}(0)\rangle$,
$\langle\mathcal{C}^{\dagger}(s)\mathcal{C}^{\dagger}(0)\rangle$ and
$\langle\mathcal{C}(s)\mathcal{C}(0)\rangle$ given in Eqs. (D.1.1), (D.1.2),
(106) and (109) we can now write down the two-time correlation functions
$\langle g_{\alpha}^{\dagger}(s)g_{\beta}(0)\rangle$ for
$\alpha,\beta\in\\{z,+,-\\}$. We will write these expressions in the continuum
limit using the generic operators

$g_{\alpha}=a_{\alpha}\mathcal{C}+b_{\alpha}\mathcal{C}^{\dagger},$ (111)

and the equivalent for $\alpha\to\beta$. One can recover the desired two-time
correlation functions of $g_{-}$, $g_{+}$ and $g_{z}$ by using the
coefficients written in Table 2. One finds that the Fourier transform of the
two-time correlation functions are

$\displaystyle\Gamma_{\alpha\beta}(\omega)\equiv\int_{0}^{\infty}\text{d}s\
\mathrm{e}^{i\omega s}\left\langle
g_{\alpha}^{\dagger}(s)g_{\beta}(0)\right\rangle=a_{\alpha}a_{\beta}\Gamma^{(\dagger,\cdot)}(\omega)+b_{\alpha}b_{\beta}\Gamma^{(\cdot,\dagger)}(\omega)+a_{\alpha}b_{\beta}\Gamma^{(\dagger,\dagger)}(\omega)+b_{\alpha}a_{\beta}\Gamma^{(\cdot,\cdot)}(\omega),$
(112)

where $\Gamma^{(\dagger,\cdot)}(\omega)$, $\Gamma^{(\cdot,\dagger)}(\omega)$,
$\Gamma^{(\dagger,\dagger)}(\omega)$, and $\Gamma^{(\cdot,\cdot)}(\omega)$ are
given in Eqs. (97), (103), (107), and (110), respectively.

$\alpha$ | $-$ | $+$ | $z$
---|---|---|---
$a_{\alpha}$ | $-\sin^{2}\left(\frac{\varphi}{2}\right)$ | $\cos^{2}\left(\frac{\varphi}{2}\right)$ | $\sin^{2}\left(\frac{\varphi}{2}\right)\cos^{2}\left(\frac{\varphi}{2}\right)$
$b_{\alpha}$ | $\cos^{2}\left(\frac{\varphi}{2}\right)$ | $-\sin^{2}\left(\frac{\varphi}{2}\right)$ | $\sin^{2}\left(\frac{\varphi}{2}\right)\cos^{2}\left(\frac{\varphi}{2}\right)$
Table 2: The coefficients used in $g_{\alpha}$ in Eq. (111) to obtain the
polaron frame coupling operators $g_{-}$, $g_{+}$ and $g_{z}$. Recall that
$\cos^{2}(\varphi/2)=(1+\epsilon/\eta)/2$ and
$\sin^{2}(\varphi/2)=(1-\epsilon/\eta)/2$.

##  Appendix E Correlation functions in the displaced frame master equation

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

$H=\frac{\eta^{d}}{2}\tau_{z}^{d}+\sum_{k}\nu_{\mathbf{k}}a_{k}^{\dagger}a_{k}+\sum_{\mu\in\\{z,+,-\\}}g^{d}_{\mu}\tau_{\mu}^{d},$
(113)

where $\eta^{d}=\sqrt{\epsilon^{2}+4|V|^{2}}$ [there is no $\kappa$
renormalisation of $V$] and

$\displaystyle g_{z}^{d}=$ $\displaystyle\
\left[\cos^{2}\left(\frac{\varphi^{d}}{2}\right)-\sin^{2}\left(\frac{\varphi^{d}}{2}\right)\right]\pi_{\Delta\Delta}+\cos\left(\frac{\varphi^{d}}{2}\right)\sin\left(\frac{\varphi^{d}}{2}\right)\left[\pi_{\mu\bar{\mu}}+\pi_{\bar{\mu}\mu}\right]$
(114a) $\displaystyle g_{+}^{d}=$
$\displaystyle-2\cos\left(\frac{\varphi^{d}}{2}\right)\sin\left(\frac{\varphi^{d}}{2}\right)\pi_{\Delta\Delta}+\cos^{2}\left(\frac{\varphi^{d}}{2}\right)\pi_{\mu\bar{\mu}}-\sin^{2}\left(\frac{\varphi^{d}}{2}\right)\pi_{\bar{\mu}\mu},$
(114b)

and $g_{-}^{d}=g_{+}^{d\dagger}$, where

$\pi_{pq}=\sum_{k}\left(p_{k}a_{k}^{\dagger}+q_{k}^{*}a_{k}\right).$ (115)

The Pauli matrices in the eigenbasis are
$\tau^{d}_{\pm}=\ket{\pm^{d}}\bra{\mp^{d}}$ and
$\tau^{d}_{z}=\ket{+^{d}}\bra{+^{d}}-\ket{-^{d}}\bra{-^{d}}$ where

$\begin{pmatrix}\ket{e}\\\
\ket{g}\end{pmatrix}=\begin{pmatrix}\cos\left(\frac{\varphi^{d}}{2}\right)&-\sin\left(\frac{\varphi^{d}}{2}\right)\\\
\sin\left(\frac{\varphi^{d}}{2}\right)&\cos\left(\frac{\varphi^{d}}{2}\right)\end{pmatrix}\begin{pmatrix}\ket{+^{d}}\\\
\ket{-^{d}}\end{pmatrix},$ (116)

with $\cos(\varphi^{d})=\epsilon/\eta^{d}$ and
$\sin(\varphi^{d})=2|V|/\eta^{d}$.

Since Eq. (113) has the same structure as Eq. (12), 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. (78) 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

$\Gamma^{d}_{pq,rs}(\omega)=\int_{0}^{\infty}\text{d}s\ \mathrm{e}^{i\omega
s}\left\langle\pi_{pq}^{\dagger}(s)\pi_{rs}(0)\right\rangle\equiv\frac{1}{2}\gamma_{pq,rs}^{d}(\omega)+iS_{pq,rs}^{d}(\omega).$
(117)

Substituting Eq. (115) into the ECF and using that $\langle
a_{k}^{\dagger}a_{k^{\prime}}\rangle=\delta_{kk^{\prime}}N(\nu_{\mathbf{k}})$,
$\langle
a_{k}a_{k^{\prime}}^{\dagger}\rangle=\delta_{kk^{\prime}}\tilde{N}(\nu_{\mathbf{k}})$
where $\tilde{N}(\nu)=1+N(\nu)$ and that other combinations equal zero, leads
to

$\displaystyle\gamma^{d}_{pq,rs}(\omega)$
$\displaystyle=2\pi\left[\cos\left(\theta_{pr}\right)\Omega_{\bar{p}r}J(\omega)\tilde{N}(\omega)+\cos\left(\theta_{qs}\right)\Omega_{q\bar{s}}J(-\omega)N(-\omega)\right],$
(118a) $\displaystyle S^{d}_{pq,rs}(\omega)$
$\displaystyle=\mathcal{P}\int_{0}^{\infty}\text{d}\nu\
J(\nu)\left[\cos\left(\theta_{pr}\right)\Omega_{\bar{p}r}\frac{\tilde{N}(\nu)}{\omega-\nu}+\cos\left(\theta_{qs}\right)\Omega_{q\bar{s}}\frac{N(\nu)}{\omega+\nu}\right].$
(118b)

To exemplify how these functions relate to the rates in Eqs. (78) we will
derive $\gamma_{--}^{d}(\eta^{d})$ explicitly. This function is given by

$\gamma_{--}^{d}(\eta^{d})=2\Re\int_{0}^{\infty}\text{d}s\
\mathrm{e}^{i\eta^{d}s}\left\langle
g_{-}^{d\dagger}(s)g_{-}^{d}(0)\right\rangle.$ (119)

Substituting in $g^{d}_{-}$ and using Eqs. (118) we can read off that
[ignoring arguments $(\eta^{d})$ on the right-hand-side and temporarily
denoting $\cos(\varphi^{d}/2)=c$ and $\sin(\varphi^{d}/2)=s$],

$\displaystyle\gamma_{--}^{d}(\eta^{d})=$ $\displaystyle\
4c^{2}s^{2}\gamma^{d}_{\Delta\Delta,\Delta\Delta}+2cs^{3}\mathrm{e}^{i\theta}\gamma^{d}_{\Delta\Delta,\mu\bar{\mu}}-2c^{3}s\mathrm{e}^{-i\theta}\gamma^{d}_{\Delta\Delta,\bar{\mu}\mu}$
$\displaystyle+2cs^{3}\mathrm{e}^{-i\theta}\gamma^{d}_{\mu\bar{\mu},\Delta\Delta}+s^{4}\gamma^{d}_{\mu\bar{\mu},\mu\bar{\mu}}-s^{2}c^{2}\mathrm{e}^{-2i\theta}\gamma^{d}_{\mu\bar{\mu},\bar{\mu}\mu}$
$\displaystyle-2c^{3}s\mathrm{e}^{i\theta}\gamma^{d}_{\bar{\mu}\mu,\Delta\Delta}-c^{2}s^{2}\mathrm{e}^{2i\theta}\gamma_{\bar{\mu}\mu,\mu\bar{\mu}}^{d}+c^{4}\gamma^{d}_{\bar{\mu}\mu,\bar{\mu}\mu}.$
(120)

Using the symmetry that $\Omega_{ab}=\Omega_{ba}$ and that
$\mathbf{d}_{\Delta}\in\Re$, we can read off from Eq. (118a) that

$\displaystyle\gamma^{d}_{\Delta\Delta,\Delta\Delta}(\eta^{d})$
$\displaystyle=\Omega_{\Delta\Delta}\gamma_{0}(\eta^{d}),$ (121)
$\displaystyle\gamma^{d}_{\bar{\mu}\mu,\mu\bar{\mu}}(\eta^{d})$
$\displaystyle=\Omega_{\mu\mu}\gamma_{0}(\eta^{d}),$ (122)
$\displaystyle\gamma^{d}_{\mu\bar{\mu},\bar{\mu}\mu}(\eta^{d})$
$\displaystyle=\Omega_{\bar{\mu}\bar{\mu}}\gamma_{0}(\eta^{d}),$ (123)
$\displaystyle\gamma^{d}_{\bar{\mu}\mu,\Delta\Delta}(\eta^{d})$
$\displaystyle=\gamma^{d}_{\Delta\Delta,\mu\bar{\mu}}(\eta^{d})=\cos\left(\theta_{\mu\Delta}\right)\Omega_{\mu\Delta}\gamma_{0}(\eta^{d}),$
(124) $\displaystyle\gamma^{d}_{\mu\bar{\mu},\Delta\Delta}(\eta^{d})$
$\displaystyle=\gamma^{d}_{\Delta\Delta,\bar{\mu}\mu}(\eta^{d})=\cos\left(\theta_{\mu\Delta}\right)\Omega_{\bar{\mu}\Delta}\gamma_{0}(\eta^{d}),$
(125) $\displaystyle\gamma^{d}_{\mu\bar{\mu},\mu\bar{\mu}}(\eta^{d})$
$\displaystyle=\gamma^{d}_{\bar{\mu}\mu,\bar{\mu}\mu}(\eta^{d})=\cos\left(\theta_{\mu\Delta}\right)\Omega_{\mu\bar{\mu}}\gamma_{0}(\eta^{d}),$
(126)

where
$\gamma_{0}(\omega)=2\pi[J(\omega)\tilde{N}(\omega)+J(-\omega)N(-\omega)$].
Using these we find that

$\gamma_{--}^{d}(\eta^{d})=k^{0}_{--}\gamma_{0}(\eta^{d}),$ (127)

where the rate coefficient is

$k^{0}_{--}=\left[c^{4}+s^{4}\right]\Omega_{\mu\bar{\mu}}+4c^{2}s^{2}\Omega_{\Delta\Delta}-s^{2}c^{2}\left[\Omega_{\bar{\mu}\bar{\mu}}+\Omega_{\mu\mu}\right]-2sc\left[c^{2}-s^{2}\right]\left[\Omega_{\bar{\mu}\Delta}+\Omega_{\mu\Delta}\right]\cos\left(\theta_{\mu\Delta}\right).$
(128)

This is rather complicated but only because we have complex transition dipole
moments. If we assume real transition dipoles
($\mathbf{d}_{\bar{\mu}}\to\mathbf{d}_{\mu}$) we find
$k^{0}_{--}=\Omega_{--}=(8\pi/3)(\mathbf{d}_{-}\cdot\mathbf{d}_{-})$ where we
have defined a new dipole vector

$\mathbf{d}_{-}=\left[c^{2}-s^{2}\right]\Re[\mathbf{d}_{\mu}]+2cs\mathbf{d}_{\Delta}.$
(129)

Likewise, if we choose purely imaginary transition dipoles
($\mathbf{d}_{\mu}\to i\Im[\mathbf{d}_{\mu}]$) we find that
$k^{0}_{--}=\Omega_{--}$ but with the dipole vector now defined as

$\mathbf{d}_{-}=\left[c^{2}-s^{2}\right]\Im[\mathbf{d}_{\mu}]+2cs\mathbf{d}_{\Delta}.$
(130)

A final limit worth checking is when the eigenstates fully localise,
$\ket{+}\to\ket{e}$ and $\ket{-}\to\ket{g}$, i.e. $c\to 1$ and $s\to 0$. In
this case, $k_{--}^{0}\to\Omega_{\mu\bar{\mu}}$, and so the decay rate is
$\gamma_{--}^{d}(\eta^{d})=\gamma_{--}^{d}(\epsilon)=\Omega_{\mu\bar{\mu}}\gamma_{0}(\epsilon)$,
which is the decay rate in the standard optical master equation.

##  Appendix F Derivation of the emission spectrum

In this appendix, we derive Eq. (52), 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, 59, 60, 61, 63].
Due to the permanent dipoles, there are additional light-matter interaction
terms in the Hamiltonian in Eq. (II), 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(\omega)=\lim_{t\to\infty}\Re\int_{0}^{\infty}\text{d}\tau\
\mathrm{e}^{i\omega\tau}\left\langle\mathbf{E}_{-}\left(\mathbf{R},t\right)\cdot\mathbf{E}_{+}\left(\mathbf{R},t+\tau\right)\right\rangle,$
(131)

where $\mathbf{R}$ is the position of the detector and the positive frequency
component of the electric field is

$\mathbf{E}_{+}\left(\mathbf{R},t\right)=i\sum_{k}\mathbf{e}_{k}f_{\mathbf{k}}a_{k}\mathrm{e}^{-i\mathbf{k}\cdot\mathbf{R}},$
(132)

and $\mathbf{E}_{-}(\mathbf{R},t)=\mathbf{E}_{+}(\mathbf{R},t)^{\dagger}$. Our
aim is to express the expectation value $\langle a_{k}^{\dagger}a_{q}\rangle$
in terms of dipole operators $\sigma_{\pm}$, which is achieved through the
Heisenberg equation of motion

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

where $H$ is the lab frame Hamiltonian in Eq. (1) and equivalently in Eq.
(II). In doing so we will have arrived at the form of the emission spectrum in
Eq. (52).

Using $[a_{k},a_{q}^{\dagger}]=\delta_{kq}$, one can show that

$\frac{\partial}{\partial
t}a_{k}(t)=-i\nu_{\mathbf{k}}a_{k}(t)-\Delta_{k}\sigma_{z}(t)-\mu_{k}\sigma_{+}(t)-\mu_{k}^{*}\sigma_{-}(t),$
(134)

where $\sigma_{\alpha}(t)=U(t)^{\dagger}\sigma_{\alpha}U(t)$ and
$U(t)=\exp(-iHt)$, which can be solved to yield

$a_{k}(t)=a_{k}(0)\mathrm{e}^{-i\nu_{\mathbf{k}}t}-\int_{0}^{t}\text{d}\tau\
\mathrm{e}^{i\nu_{\mathbf{k}}\left(\tau-t\right)}\left[\Delta_{k}\sigma_{z}(\tau)+\mu_{k}\sigma_{+}(\tau)-\mu_{k}^{*}\sigma_{-}(\tau)\right].$
(135)

The first term in Eq. (135) 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. (135) into Eq. (132) and taking the continuum limit we
obtain

$\mathbf{E}_{+}\left(\mathbf{R},t\right)=\int_{0}^{\infty}\text{d}\nu\
\sqrt{J(\nu)}\int_{0}^{t}\text{d}\tau\
\mathrm{e}^{i\nu_{\mathbf{k}}\left(\tau-t\right)}\left[\mathbf{O}_{z}\sigma_{z}(\tau)+\mathbf{O}_{\mu}\sigma_{+}(\tau)+\mathbf{O}_{\bar{\mu}}\sigma_{-}(\tau)\right],$
(136)

where

$\mathbf{O}_{p}=\int_{d\Omega_{\mathbf{k}}}\
\text{d}\Omega_{\mathbf{k}}\sum_{\lambda}\mathbf{e}_{k}\left(\mathbf{d}_{p}\cdot\mathbf{e}_{k}\right)\mathrm{e}^{i\mathbf{k}\cdot\left(\mathbf{r}-\mathbf{R}\right)}.$
(137)

To progress analytically, we need to know how the emitter operators
$\sigma_{\alpha}(t)$ evolve in time. However, this is very complicated.
Instead, we make the so-called harmonic decomposition (see Chapter 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

$\sigma_{+}(\tau)\approx\sigma_{+}(t)\mathrm{e}^{i\epsilon\left(\tau-t\right)},\quad\
\sigma_{-}(\tau)\approx\sigma_{-}(t)\mathrm{e}^{-i\epsilon\left(\tau-t\right)},\quad\text{
and}\quad\ \sigma_{z}(\tau)\approx\sigma_{z}(t).$ (138)

Making the harmonic decomposition in Eq. (136) yields

$\mathbf{E}_{+}(\mathbf{R},t)=\int_{0}^{\infty}\text{d}\nu\
\sqrt{J(\nu)}\left[\mathbf{O}_{z}\sigma_{z}(t)j_{0}(\nu,t)+\mathbf{O}_{\mu}\sigma_{+}(t)j_{+}(\nu,t)+\mathbf{O}_{\bar{\mu}}\sigma_{-}(t)j_{-}(\nu,t)\right],$
(139)

where

$j_{\alpha}(\nu,t)=\int_{0}^{t}\text{d}\tau\
\mathrm{e}^{i\left(\tau-t\right)\left(\nu+\alpha\epsilon\right)}=-i\frac{1-\mathrm{e}^{-it\left(\nu+\alpha\epsilon\right)}}{\nu+\alpha\epsilon}.$
(140)

The function $j_{\alpha}(\nu,t)$ is dominated by the contribution near to
$\nu=-\alpha\epsilon$ and so we approximate it as a delta function,
$j_{\alpha}(\nu,t)\approx 2\pi\delta(\nu+\alpha\epsilon)$ [2]. After
performing this approximation, one finds that the delta function corresponding
to $j_{+}(\nu,t)$ lies outwith the integration domain $\nu\in[0,\infty]$ and
so the term going as $\sigma_{+}(t)$ in Eq. (139) vanishes. Moreover, the term
proportional to $j_{0}(\nu,t)$ provides a delta function at zero frequency,
leading to the term going as $\sigma_{z}(t)$ in Eq. (139) 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. (52). Thus, the only surviving term in Eq. (139) is
proportional to $\sigma_{-}(t)$, and so we obtain

$\mathbf{E}_{+}\left(\mathbf{R},t\right)\approx
2\pi\sqrt{J(\epsilon)}\mathbf{O}_{\bar{\mu}}\left(\mathbf{r},\mathbf{R}\right)\sigma_{-}(t),$
(141)

where we have made the dependence of $\mathbf{O}_{\bar{\mu}}$ on $\mathbf{r}$
and $\mathbf{R}$ explicit. After substituting Eq. (141) and its Hermitian
conjugate into Eq. (131), one obtains Eq. (52), and an explicit expression for
the Green’s function $\alpha_{\text{prop}}(\mathbf{r},\mathbf{R},\omega)$.

##  Appendix G 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 associated temperature $T=1/\beta$. In the polaron
framework, the environment is not in such a convenient form. As such, we
derive here an effective Hamiltonian that encapsulates the polaron thermalised
state but has the effective environment in a free Gibbs state.

We start with the lab frame Hamiltonian,

$H=\frac{\epsilon}{2}\sigma_{z}+V\sigma_{+}+V^{*}\sigma_{-}+E_{\text{dip}}+\sum_{k}\nu_{\mathbf{k}}a_{k}^{\dagger}a_{k}+\pi_{DD}\mathcal{I}+\pi_{\Delta\Delta}\sigma_{z}+\pi_{\mu\bar{\mu}}\sigma_{+}+\pi_{\bar{\mu}\mu}\sigma_{-},$
(142)

where

$\displaystyle\pi_{pq}$
$\displaystyle=\sum_{k}\left(p_{k}a_{k}^{\dagger}+q_{k}^{*}a_{k}\right),$
(143) $\displaystyle p_{k}$
$\displaystyle=i\left(\mathbf{d}_{p}\cdot\mathbf{e}_{\mathbf{k}}\right),$
(144)

with $p,q\in\\{\mu,\bar{\mu},D,\Delta\\}$. In our calculations, the polaron
frame initial state is

$\rho_{p}^{0}=\left|g\right>\left<g\right|\otimes\rho_{E},$ (145)

where
$\rho_{E}=\exp[-\beta\sum_{k}\nu_{\mathbf{k}}a_{k}^{\dagger}a_{k}]/\mathcal{Z}_{E}$
is a thermal state. Performing the inversion of the unitary transformations to
go from the displaced polaron frame to the displaced frame, the initial state
is $\rho_{d}^{0}=U^{\dagger}\rho_{p}^{0}U$ where the polaron transformation is
$U=B(\delta)\left|e\right>\left<e\right|+B(-\delta)\left|g\right>\left<g\right|$
and $\delta_{k}=\Delta_{k}/\nu_{\mathbf{k}}$ and
$B(\alpha)=\exp[\sum_{k}(\alpha_{k}a_{k}^{\dagger}-\alpha_{k}^{*}a_{k})]$.
Therefore,

$\rho_{d}^{0}=\left|g\right>\left<g\right|\otimes\eta(\delta),$ (146)

where we have defined a displaced thermal state as
$\eta(\alpha)=B(\alpha)\rho_{E}B(-\alpha)$. We can then obtain the lab frame
initial state via $\rho^{0}_{l}=B(-d)\rho_{d}^{0}B(d)$ where
$d_{k}=D_{k}/\nu_{\mathbf{k}}$, leading to

$\rho^{0}_{l}=\left|g\right>\left<g\right|\otimes\eta(\delta-d).$ (147)

The expression for $\eta(\alpha)$ can be rewritten by making use of the
identity,

$\exp\left(\mathrm{e}^{S}X\mathrm{e}^{-S}\right)=\mathrm{e}^{S}\mathrm{e}^{X}\mathrm{e}^{-S},$
(148)

which holds if $\mathrm{e}^{S}\mathrm{e}^{-S}=\mathcal{I}$. Using this
identity yields,

$\eta(\alpha)=\exp\left[-\beta\sum_{k}\nu_{\mathbf{k}}\left(a_{k}^{\dagger}-\alpha_{k}^{*}\right)\left(a_{k}-\alpha_{k}\right)\right].$
(149)

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, $\eta(\delta-d)$. The initial environment state in the effective lab
frame will be
$\bar{\rho}_{E}=\exp[-\beta\sum_{k}\nu_{\mathbf{k}}b_{k}^{\dagger}b_{k}]/\mathcal{Z}_{E}$
and so by comparison with Eq. (149), we know that

$b_{k}=a_{k}-(\delta_{k}-d_{k}).$ (150)

Substituting this into our actual lab frame Hamiltonian in Eq. (142) 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},$ (151)

where

$\displaystyle\bar{H}_{0}$
$\displaystyle=\left(\frac{\epsilon}{2}+G_{\Delta\Delta}\right)\sigma_{z}+V\sigma_{+}+V^{*}\sigma_{-}+G_{\mu\bar{\mu}}\sigma_{+}+G_{\bar{\mu}\mu}\sigma_{-}+\sum_{k}\nu_{\mathbf{k}}b_{k}^{\dagger}b_{k},$
(152) $\displaystyle\bar{H}_{I}$
$\displaystyle=\pi^{\prime}_{\Delta\Delta}\left(\mathcal{I}+\sigma^{z}\right)+\pi^{\prime}_{\mu\bar{\mu}}\sigma_{+}+\pi^{\prime}_{\bar{\mu}\mu}\sigma_{-},$
(153)

and

$\displaystyle\pi^{\prime}_{pq}$
$\displaystyle=\sum_{k}\left(p_{k}b_{k}^{\dagger}+q_{k}^{*}b_{k}\right),$
(154) $\displaystyle G_{pq}$
$\displaystyle=\sum_{k}\left(p_{k}\delta_{k}^{*}+q_{k}^{*}\delta_{k}\right).$
(155)

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. (151) for TEMPO calculations.