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

## I. INTRODUCTION

In general, the interaction of atomic systems with light through transition dipoles is well understood, and the optical master equation describing exciton creation and annihilation through photon emission and absorption, respectively, is derived in many introductory texts dedicated to open quantum systems.^{1–4} Atomic systems have highly symmetric electron orbitals and so possess negligible permanent dipoles. However, many physical systems do not share this property and can possess permanent dipoles stronger than their transition dipole moments. Such systems include molecules with parity mixing of the molecular state,^{5–12} quantum dots with asymmetric confining potentials,^{13–23} nanorods with non-centrosymmetric crystallographic lattices,^{24–26} and superconducting circuits.^{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} second-harmonic generation,^{37,38} and bathochromic shifts.^{39,40} Understanding the role of strong permanent dipoles and their influence on control schemes, such as external driving fields, is highly relevant for the design of new quantum technologies, as decoherence induced by the local electromagnetic field limits the capabilities of quantum computation implementations. Studying the impacts of permanent dipoles may provide additional channels for control over such quantum systems, as well as advancements in quantum chemistry that may uncover novel biochemical processes.

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

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

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

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

## II. THE MODEL AND POLARON FRAME

^{56–58}$ \Pi = i \u2211 k e k f k ( a k \u2020 \u2212 a k )$ is the electric displacement field, where $ k = { k , \lambda }$ is a four-vector representing both the wavevector

**k**and polarization state

*λ*of the mode with polarization vector $ e k$, energy $ \nu k$, and coupling strength $ f k = \nu k / ( 2 V )$, where $V$ is the field volume.

*a*and $ a k \u2020$ are the field mode annihilation and creation operators, and the dipolar self-energy term is

_{k}*H*require that $ d p \u2208 \mathbb{R} 3$ for $ p \u2208 { \Delta , D}$ and $ d \mu \u2208 \u2102 3$.

At this point, it is often assumed that either $ | d \Delta | \u2248 | d D | \u2248 0$, which leads to the standard optical master equation (SOME),^{1–4} or that $ | d D | \u2248 | d \mu | \u2248 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.

*d*denotes the displaced frame, and displacement operators are given by

*exactly*with the terms that result from displacing the light–matter interactions, such that Eq. (10) is independent of

*D*. This is a remarkable result, as it tells us that the standard perturbative master equation would be considerably inaccurate for most values of $ d D$ as this term would have a non-negligible contribution to the dynamics.

_{k}Throughout this paper, we make the standard assumption that the transition dipole is small enough to permit an accurate second-order perturbative expansion in its magnitude, $ | d \mu |$. Even for perturbative transition dipole moments, the master equation derived using *H _{d}* in Eq. (10), referred to as the

*displaced frame master equation*(DFME), will become inaccurate if $ | d \Delta |$ is large.

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

*H*, which ensures that the perturbation theory yields the Redfield equation.

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

^{2}We have defined dimensionless dipole vectors $ d \u0303 p = d p / d ref$ measured against an arbitrary reference value which for convenience we take to be $ d ref = 1 \u2009 D \u2248 9 \xd7 10 \u2212 6 eV \u2212 1$. The purpose of introducing $ d ref$ is to ensure that the spectral density has dimensions of energy. For free space and unpolarized light,

## III. REALISTIC PARAMETERS

^{12}Conversely, the multipolar-gauge free-space spectral density

*ν*is large. For simplicity, in this paper, we use a free space model with permanent dipole strengths inspired by realistic molecules embedded in solvents and cavities.

_{c}As we derive in Appendix C, $ S = \u222b 0 \u221e d \nu J ( \nu ) / \nu 2 = s 0 ( \nu c d ref ) 2 / [ 2 ( 2 \pi ) 3 ]$ is a dimensionless Huang–Rhys parameter where *s*_{0} is a Purcell enhancement factor. $ s 0 = 1$ corresponds to free space, and $ s 0 > 1$ to restricted environmental geometries and solvents. Multiplying the spectral density by *s*_{0} is accurate within the weak field regime.^{59} The cutoff frequency in Eq. (24) is required to enforce the electric dipole approximation implicit in Eq. (1), such that $ \nu c \u223c 1 / r d$, where $ r d \u2273 50 2 a 0$ is the size of the dipole, giving $ \nu c \u223c [ 1 , 10 ] \u2009 eV$.^{57} Unless otherwise stated, we use $ \nu c = 1 \u2009 \u2009 eV$.

In all our numerical results—except when comparing to TEMPO which we discuss in that section—we take parameter values using the gamma-globulin protein as an example: $ | d \Delta | \u223c 100 \u2009 D , \u2009 | d \mu | = 0.1 \u2009 D$, and emitter energy $ \u03f5 = 3.16 \u2009 \u2009 eV$.^{32,60} Taking the lifetime of an individual dipole—i.e., when $ | d \Delta | = 0 = V$—to be $ 1 \u2009 ns$, requires that $ s 0 = 5.705 \xd7 10 6$ for the values introduced so far.^{61} These values fix $ S = 9.341 \xd7 10 \u2212 7$ and so the field remains weak, thus justifying our introduction of the Purcell enhancement factor. Moreover, these values predict that polaron formation ( $ \kappa \u226a 1$) will occur at $ | d \Delta | \u223c 100 \u2009 D$. Finally, we choose the inverse temperature as $ \beta = 2 \u2009 \u2009 eV \u2212 1$ corresponding to $ T \u2248 5800 \u2009 K$.

To summarize, the important parameters in this model are the dipole strengths $ | d \mu |$ and $ | d \Delta |$, their relative angle $ \theta \mu \Delta $, the eigenenergy *η*, the driving strength $ | V |$, the relative complex phase $ \u03d1 \mu V$, and the cutoff frequency *ν _{c}*.

## IV. EFFECTS OF STRONG PERMANENT DIPOLES

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

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

*γ*, and oscillations in the coherence at frequency $ \eta \xaf$. The transition and decoherence rates are

_{d}^{1}

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

*V*are induced by the driving. The other three Fourier transforms of the $ C ( s )$ two-time correlation functions also depend on these four contributions, except that $ \u27e8 C \u2020 ( s ) C \u2020 ( 0 ) \u27e9$ and $ \u27e8 C ( s ) C ( 0 ) \u27e9$ 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 amplitude $ | V |$, the relative phase $ \u03d1 \mu V$, and dipole magnitudes Ω

_{pq}, as

For additional context, when $ | V | \u226a \u03f5$—the parameter regime in which polaron theory is valid—the eigenbasis angle is $ \phi \u2248 0$ and so $ g z \u2248 0$ and $ g + \u2248 C$. In this limit, the secular approximation used in the master equation in Eqs. (25) becomes exact, and $ \gamma \u2191 \u2248 2 \u211c [ \Gamma ( \u2212 \eta ) ]$ with $ \Gamma ( \omega )$ in Eq. (30) and $ \gamma \u2193 \u2248 \u222b 0 \u221e d s \u2009 e i \eta s \u27e8 C ( s ) C \u2020 ( 0 ) \u27e9$ which, as we show in Appendix E, evaluates to $ \gamma \u2193 \u2248 2 \u211c [ \Gamma \xaf ( \eta ) ]$ where $ \Gamma \xaf ( \omega )$ is equal to Eq. (30) but with $ \Gamma V , 1 ( \omega ) \u2192 \u2212 \Gamma V , 1 ( \omega )$.

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

### A. One photon processes, **Γ**_{1}(*ω*)

_{1}(

*ω*)

^{1}—our model in the absence of permanent dipoles and driving—which is

We now return to Eq. (32). In Ref. 62, an analytic solution to Eq. (33) was found by exploiting the fact that $ K ( \epsilon )$ only depends on $ J ( \nu )$ through its moments $ \mu m = 4 \Omega \Delta \Delta \u222b 0 \u221e d \nu [ J ( \nu ) / \nu 2 ] \nu m$ for $ m = 1 , 2 , \u2026 , \u221e$, and that moments of lower order contribute relatively more. We can then evaluate $ K ( \epsilon )$ by replacing $ J ( \nu )$ with a truncated spectral density $ J \u2032 ( \nu ) = \u2211 k = 1 N * | f \u2032 k | 2 \delta ( \nu \u2212 \nu k \u2032 )$, as long as we choose the coupling strengths $ { f k \u2032}$ and energies $ { \nu k \u2032}$ of the modes such that $ J \u2032 ( \nu ) / \nu 2$ has the same lowest moments as *μ _{m}* for $ m = 1 , 2 , \u2026 , 2 N *$. The solution converges rapidly for increasing $ N *$, and a single-mode truncation is often very accurate, with truncation mode parameters: $ \nu \u2032 1 \u2261 \nu s = \mu 2 / \mu 1$ and $ f \u2032 1 \u2261 f s = \mu 2 1 / 2 \u221d \Omega \Delta \Delta 1 / 2$. A single-mode truncation captures all of the essential physics described by $ K ( \omega )$. Thus, for the analytical analysis in this work, we use a single-mode truncation to evaluate Eq. (33); however, in all of the simulations in this paper, we use converged values.

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

### B. Two-photon processes, **Γ**_{2}(*ω*)

_{2}(

*ω*)

*σ*and $ \sigma \xb1$ interactions do not commute, and so are unique to permanent dipole interactions. The rate function is

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

Before we move onto the processes induced by driving, in Fig. 3, we plot the decay rate $ \gamma \u2193$ as a function of *ν _{c}* and $ | d \Delta |$ for

*V*= 0. Since

*V*= 0, the eigenbasis coincides with the $ { | e \u27e9 , | g \u27e9}$ basis and the decay rate is $ \gamma \u2193 = 2 \u211c [ \Gamma ( \u03f5 ) ]$, where $ \Gamma ( \omega )$ is in Eq. (30).

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

*ν*and $ | d \Delta | \u223c 100 \u2009 D$.

_{c}### C. Driving-induced one-photon processes, **Γ**_{V,1}(*ω*)

_{V,1}(

*ω*)

Importantly, the prefactor in Eq. (45) is proportional to $ cos ( \u03d1 \mu V ) \u2009 cos ( \theta \mu \Delta )$, which allows one to control whether this function suppresses or enhances the rates and Lamb shifts. As we discuss later, this control could be easy to achieve in practice by tuning the driving phase.

### D. Driving-induced zero-photon processes, **Γ**_{V,0}(*ω*)

_{V,0}(

*ω*)

### E. Numerical analysis of the rates

In Fig. 4, panel (a) shows the total decay rate for different driving strengths *V* as a function of $ | d \Delta |$, and panel (b) shows the four contributions to the rate with $ V = 10 \u2212 4 \u2009 eV$. Due to detailed balance, the excitation rate can be obtained from Fig. 4 by $ \gamma \u2191 = e \u2212 \beta \eta \gamma \u2193$ and, provided pure dephasing is negligible (which is true for the parameters in Fig. 4 because $ 2 \kappa | V | \u2aa1 \u03f5$), the decoherence rate can be obtained from Fig. 4 by $ \gamma d = ( 1 + e \u2212 \beta \eta ) \gamma \u2193 / 2$. The key points resulting from the analysis of Fig. 4 are as follows: (1) one can both increase and decrease rates through the physical mechanisms enabled by the presence of permanent dipoles and driving, particularly by controlling $ \Gamma V , 1 ( \omega )$ processes, and (2) all physical processes found in the polaron ECFs contribute non-negligibly to the total rates.

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

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

## V. COMPARISON TO EXACT DYNAMICS

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

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

## VI. EMISSION SPECTRUM

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

^{64–68}$ \alpha prop ( r , R , \omega )$ accounts for propagation and filtering of the light from the dipole to the detector at positions

**r**and

**R**, respectively.

^{64}The polarization spectrum is

*ρ*. Using $ \rho l = B ( \u2212 D / \nu ) U \u2020 \rho p U B ( D / \nu )$, where

_{l}*U*is the polaron transformation in Eq. (11), $ B ( \xb1 D / \nu )$ are displacement operators, and

*ρ*is the polaron frame density operator, we have

_{p}^{1,69}The QRT utilizes the cyclicity of the trace to rewrite the two-time expectation value as

Due to the polaron transformation, Eq. (54) captures the photonic sideband. However, the QRT contains an implicit Born approximation^{69} as well as any approximations used in the derivation of the master equation. Moreover, the QRT produces spectral lines at the Lambshifted eigenenergies of *H _{S}* because the ground-excited coherence

*ρ*oscillates at these frequencies [see Eq. (25b)]. This means that the QRT is unable to capture spectral lines associated with multiphoton transitions. However, multiphoton transitions scale with the square of the spectral density, and so for weak light–matter coupling, we expect the associated spectral lines to be much smaller than the photonic sideband.

_{eg}*κ*, which in turn provides an effective measurement of $ \Omega \Delta \Delta $ and, therefore, of the strength of the $ | d \Delta |$ dipole moment. In measurements, $ I p , \xd7 ( \omega )$ cannot be separated from $ I p , 0 ( \omega )$, but as we show in the following, $ \u222b \u2212 \u221e \u221e d \omega I p , \xd7 ( \omega )$ can be accurately approximated as the integral of the measured spectrum over the Lorentzian parts of the dominant peak(s).

^{70}

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

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 in Table I. We show that it is possible to evaluate $ \rho e e ( \u221e )$ and *κ* from $ I p , 0 ( \omega )$ using Eqs. (57) and (58) and, provided that the temperature of the experiment is known, these values can be used to obtain *ϵ* and $ | V |$. If the spectral density is also known, it is also possible to obtain $ | d \Delta |$ from *κ*.

Param. . | Analytic . | Meas. . |
---|---|---|

$ \kappa 2$ | 0.632 | 0.621 |

$ \rho e e ( \u221e )$ | 0.450 | 0.450 |

$ \u03f5 / eV$ | 0.100 | 0.100 |

$ | V | \u2009 / eV$ | 0.010 | 0.010 |

$ | d \Delta | / \u2009 D$ | 100 | 101 |

Param. . | Analytic . | Meas. . |
---|---|---|

$ \kappa 2$ | 0.632 | 0.621 |

$ \rho e e ( \u221e )$ | 0.450 | 0.450 |

$ \u03f5 / eV$ | 0.100 | 0.100 |

$ | V | \u2009 / eV$ | 0.010 | 0.010 |

$ | d \Delta | / \u2009 D$ | 100 | 101 |

*κ*can be approximately found using Eq. (58) with $ \u222b \u2212 \u221e \u221e d \omega I p , \xd7 ( \omega ) \u2248 \u222b R d \omega \u2009 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. In the example in Fig. 7, this region is between the vertical dashed lines in the inset. One could increase the accuracy of the $ \kappa 2$ estimation in Fig. 7 by integrating over the other two peaks, but these contributions are small. If the spectral density of the photon bath is known, the relation $ \kappa = exp \u2009 [ \u2212 \varphi ( 0 ) / 2 ]$ can be inverted to obtain $ | d \Delta |$. The remaining parameters of the system,

*ϵ*and $ | V |$, can be estimated by assuming that the system thermalizes with respect to $ H S = \eta \tau z / 2$, yielding the steady state

*ϵ*. Then, using $ \eta \xaf \u2248 \eta = ( \u03f5 2 + 4 \kappa 2 | V | 2 ) 1 / 2$, one can estimate $ | V |$.

## VII. IMPORTANCE OF INITIAL CONDITIONS

*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 ) \u2297 \rho E$ in the polaron frame models an initial state in the lab frame that is more similar to $ \rho \beta $. For example, in our calculations in Fig. 6, $ \rho S ( 0 ) = | g \u27e9 \u27e8 g |$ and so $ \rho S ( 0 ) \u2297 \rho E$ in the polaron frame becomes $ \rho l ( 0 ) = | g \u27e9 \u27e8 g | \u2297 B ( [ D + \Delta ] / \nu ) \rho E B ( \u2212 [ 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 thermalized state, and so is a good initial state in which to model excitation.

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

*π*is given in Eq. (6), $ \u03f5 \u0303 = \u03f5 + 2 G \Delta \Delta , \u2009 V \u0303 = V + G \mu \mu \xaf$, where $ G p q = \u2211 k ( p k \Delta k * + q k * \Delta k ) / \nu k$. We note that because the displacement direction of the polaron transformation is state dependent, the necessary effective Hamiltonian depends on $ \rho S ( 0 )$.

_{pq}By comparing the Hamiltonians in Eqs. (61) and (5), we see stark differences if both are assumed to have the initial state $ \rho S ( 0 ) \u2297 \rho E$. This includes a renormalization of the transition energy $ \u03f5 \u2192 \u03f5 \u0303$ and the driving term $ V \u2192 V \u0303$. 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 D D I \u2032$ in Eq. (5), have a similar effect. As we showed in Eq. (10), this type of interaction can be removed through a displacement transformation and so, ultimately, the role of an identity interaction is to change the initial environment state from a thermal Gibbs state to a displaced one.

## VIII. CONCLUSION

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

We have shown three key results. (1) In Sec. IV, we showed that transition and decoherence rates can be engineered for practical application by exploiting permanent dipoles, for example, to eliminate optical decoherence and transitions. By using the novel physical processes explicit in the polaron rate equations we derive, one can design systems to exploit these effects. The novel physical processes arising from the permanent dipoles are as follows: a photon sideband with a relative contribution to the emission spectrum scaling as $ \kappa 2$, two-photon processes scaling as $ \Omega \mu \Delta 2$, and a term linear in the driving amplitude scaling as $ | V | \Omega \mu \Delta \u2009 cos ( \u03d1 \mu \Omega )$. (2) In Sec. V, we proved that the optical polaron description provides a much more accurate master equation by comparing it to the DFME and TEMPO. (3) In Sec. VI, we indicated distinguishable features of permanent dipoles in emission spectra and described possible measurements to obtain $ \rho e e ( \u221e )$, *κ*, $ \eta \xaf$ and by extension the bare energy splitting *ϵ*, the driving strength $ | V |$, and the permanent dipole magnitude $ | d \Delta |$. The distinguishable features are the photonic sideband and the altered widths and positions 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 $ d \mu $ and $ d \Delta $ are perpendicular and are unique to pure dephasing interactions that are induced by permanent dipoles, as opposed to vibrational environments. The linear driving term, $ \Gamma V , 1 ( \omega )$, may be particularly useful for quantum technologies. As we have shown, one could engineer devices in which the optical decoherence and transition rates can be varied from essentially zero to faster than those in the same system but without permanent dipoles, by simply changing the driving phase by *π*. This has potential applications in quantum computing and energy storage.

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

## ACKNOWLEDGMENTS

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

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

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

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

## DATA AVAILABILITY

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

### APPENDIX A: DISPLACEMENT OPERATORS

### APPENDIX B: CALCULATION OF $ \u27e8 C \u27e9$

*a*annihilates with $ | 0 \u27e9$. One finds that

*C*given in Eq. (B1) yield

### APPENDIX C: THE HUANG–RHYS PARAMETER

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

*s*

_{0,}which is justified for weak field interactions.

### APPENDIX D: NON-SECULAR MASTER EQUATIONS

Notice that if the coherences are initially zero ( $ \rho + \u2212 ( 0 ) = \rho \u2212 + ( 0 ) = 0$), then if $ k \u2212 = k + = 0$ as well, the coherences will be zero at all times *t*. This occurs in the polaron frame if $ V \u2192 0$ since $ g z \u221d \u2009 sin \u2009 ( \phi ) \u2192 0$, i.e., the eigenstates fully localize, but this does not happen in the displaced frame in the same limit because $ g z d$ does not become zero.

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

##### 1. Preliminary calculations

All $ \u27e8 g \alpha \u2020 ( s ) g \beta ( 0 ) \u27e9$ depend on four two-time correlation functions: $ \u27e8 C \u2020 ( s ) C ( 0 ) \u27e9 , \u2009 \u27e8 C ( s ) C \u2020 ( 0 ) \u27e9 , \u2009 \u27e8 C ( s ) C ( 0 ) \u27e9$, and $ \u27e8 C \u2020 ( s ) C \u2020 ( 0 ) \u27e9$. Each of these in turn depends on either seven or nine unique two-time correlation functions involving *a _{k}*, $ a k \u2020$, and $ B \xb1$. In this appendix, we list the results of all necessary two-time correlation functions, which are each derived using the same mathematics as in the explicit examples in Appendix B.

###### a. $ \u27e8 C \u2020 ( s ) C ( 0 ) \u27e9$

###### b. $ \u27e8 C ( s ) C \u2020 ( 0 ) \u27e9$

###### c. $ \u27e8 C \u2020 ( s ) C \u2020 ( 0 ) \u27e9$

###### d. $ \u27e8 C ( s ) C ( 0 ) \u27e9$

To calculate this two-time correlation function, we require

##### 2. Two-time correlation functions $ \u27e8 g \alpha \u2020 ( s ) g \beta ( 0 ) \u27e9$

*g*by using the coefficients written in Table II. The Fourier transforms of the two-time correlation functions are

_{z}α
. | − . | + . | z
. |
---|---|---|---|

$ a \alpha $ | $ \u2212 sin 2 ( \phi 2 )$ | $ \u2009 cos 2 ( \phi 2 )$ | $ \u2009 sin 2 ( \phi 2 ) \u2009 cos 2 ( \phi 2 )$ |

$ b \alpha $ | $ \u2009 cos 2 ( \phi 2 )$ | $ \u2212 sin 2 ( \phi 2 )$ | $ \u2009 sin 2 ( \phi 2 ) \u2009 cos 2 ( \phi 2 )$ |

α
. | − . | + . | z
. |
---|---|---|---|

$ a \alpha $ | $ \u2212 sin 2 ( \phi 2 )$ | $ \u2009 cos 2 ( \phi 2 )$ | $ \u2009 sin 2 ( \phi 2 ) \u2009 cos 2 ( \phi 2 )$ |

$ b \alpha $ | $ \u2009 cos 2 ( \phi 2 )$ | $ \u2212 sin 2 ( \phi 2 )$ | $ \u2009 sin 2 ( \phi 2 ) \u2009 cos 2 ( \phi 2 )$ |

### APPENDIX F: CORRELATION FUNCTIONS IN THE DISPLACED FRAME MASTER EQUATION

*κ*renormalization of

*V*] and

A final limit worth checking is when the eigenstates fully localize, $ | + d \u27e9 \u2192 | e \u27e9$ and $ | \u2212 d \u27e9 \u2192 | g \u27e9$, i.e., $ c \u2192 1$ and $ s \u2192 0$. In this case, $ k \u2212 \u2212 0 \u2192 \Omega \mu \mu $, and so the decay rate is $ \gamma \u2212 \u2212 d ( \eta d ) = \gamma \u2212 \u2212 d ( \u03f5 ) = \Omega \mu \mu \gamma 0 ( \u03f5 )$, which is the decay rate in the standard optical master equation.

### APPENDIX G: DERIVATION OF THE EMISSION SPECTRUM

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

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

*H*is the lab frame Hamiltonian in Eq. (1). By obtaining $ \u27e8 a k \u2020 a q \u27e9$ in terms of $ \sigma \xb1$, we will have arrived at the form of the emission spectrum in Eq. (51).

^{2}After performing this approximation, one finds that the delta function corresponding to $ j + ( \nu , t )$ lies out with the integration domain $ \nu \u2208 [ 0 , \u221e ]$, and so, the term going as $ \sigma + ( t )$ in Eq. (G9) 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. (G9) to be proportional to $ 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. (51). Thus, the only surviving term in Eq. (G9) is proportional to $ \sigma \u2212 ( t )$, and so, we obtain

**r**and

**R**explicit. After substituting Eq. (G11) and its Hermitian conjugate into Eq. (G1), one obtains Eq. (52), and an explicit expression for the Green's function $ \alpha prop ( r , R , \omega )$.

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

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

### APPENDIX I: EFFECTIVE HAMILTONIAN

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

*b*, which we will relate to the

_{k}*a*to ensure we model our desired initial state, $ \eta ( \delta \u2212 d )$. The initial environment state in the effective lab frame will be $ \rho \xaf E = exp \u2009 [ \u2212 \beta \u2211 k \nu k b k \u2020 b k ] / Z E$, and so by comparison with Eq. (I8), we know that

_{k}## REFERENCES

*The Theory of Open Quantum Systems*

*Quantum Interference and Coherence: Theory and Experiments*

*Quantum Optics*