Environmental non-additivity and Franck-Condon physics in non-equilibrium quantum systems

TL;DR

This paper demonstrates that non-additive interactions between vibrational and electromagnetic environments are crucial for accurately modeling non-equilibrium quantum systems, revealing effects like Franck-Condon physics and population inversion.

Contribution

It introduces a non-additive master equation using a collective coordinate approach to better capture vibrational effects in quantum emitters.

Findings

Photon emission rates decrease with vibrational coupling, aligning with Franck-Condon physics.

Additive models fail to predict vibrational influence on emission.

Non-additivity enables population inversion under incoherent excitation.

Abstract

We show that for a quantum system coupled to both vibrational and electromagnetic environments, enforcing additivity of their combined influences results in non-equilibrium dynamics that does not respect the Franck-Condon principle. We overcome this shortcoming by employing a collective coordinate representation of the vibrational environment, which permits the derivation of a non-additive master equation. When applied to a two-level emitter our treatment predicts decreasing photon emission rates with increasing vibrational coupling, consistent with Franck-Condon physics. In contrast, the additive approximation predicts the emission rate to be completely insensitive to vibrations. We find that non-additivity also plays a key role in the stationary non-equilibrium model behaviour, enabling two-level population inversion under incoherent electromagnetic excitation.

Figures

Equations

$${\Gamma }_{e\to g}=\sum _{n,m}{p}_n\mathcal{J}(\Delta {\omega }_{\tilde{m},n}){\left|⟨\tilde{m}|n⟩\right|}^2.$$

$$H_{\mathrm{I}}^{\mathrm{PH}}=|e⟩⟨e|\otimes \sum _k{g}_k(b_k^{†}+b_k)+|e⟩⟨e|\sum _k\frac{g_k^2}{{\nu }_k},$$

$$H_{\mathrm{I}}^{\mathrm{EM}}=\sum _l(f_l{\sigma }^{†}{a}_l+f_l^{\ast }\sigma a_l^{†}),$$

$$H_{\mathrm{S}}^{\prime }$$

$$H_{\mathrm{I}}^R$$

$$H_{\mathrm{B}}^R$$

$$\mathcal{L}[\rho (t)]=-i\left[H_{\mathrm{S}}^{\prime },\rho (t)\right]+{\mathcal{K}}_{\mathrm{R}}[\rho (t)]+{\mathcal{K}}_{\mathrm{EM}}[\rho (t)],$$

$${\mathcal{K}}_{\mathrm{R}}[\rho (t)]=\left[S,\rho (t)\zeta \right]+\left[{\zeta }^{†}\rho (t),S\right],$$

$$\zeta =\frac{\pi }{2}\sum _{jk}{J}_{\mathrm{R}}({\lambda }_{jk})\left[\coth (\frac{{\lambda }_{jk}}{2k_{\mathrm{B}}{T}_{\mathrm{R}}})+1\right]S_{jk}\text{outerproduct}{\psi }_j{\psi }_k,$$

$${\mathcal{K}}_{\mathrm{EM}}[\rho (t)]$$

$${\Gamma }_{e\to g}=\sum _n⟨g,n|\mathcal{L}[{\rho }_X(0)]|g,n⟩,$$

$$H_{\mathrm{I}}(t)=S(t)\otimes B(t)+S(t{)}^2\sum _m\frac{h_m^2}{{\tilde{\nu }}_m}+\sum _{\alpha }{A}_{\alpha }(t)\otimes E_{\alpha }(t),$$

$${\partial }_t\rho (t)=\mathcal{L}[\rho (t)]=-i\left[H_{\mathrm{S}}^{\prime },\rho (t)\right]+{\mathcal{K}}_{\mathrm{R}}[\rho (t)]+{\mathcal{K}}_{\mathrm{EM}}[\rho (t)],$$

$$\mathcal{K}_{\rm R}[\rho(t)]=-i\sum_{m}\frac{h_{m}^{2}}{\tilde{\nu}_{m}}[S^{2},\rho(t)]-\int_{0}^{\infty}d\tau\big{[}S,\big{[}S(-\tau),\rho(t)\big{]}\big{]}C_{+}(\tau)-\int_{0}^{\infty}d\tau\big{[}S,\big{\{}S(-\tau),\rho(t)\big{\}}\big{]}C_{-}(\tau),$$

$$C_{\pm }(\tau )=\frac{1}{2}⟨B(\tau )B\pm B(-\tau )B⟩,$$

$$C_{+}(\tau )={\int }_0^{\infty }d\omega J_{\mathrm{R}}(\omega )\coth \left(\frac{\omega }{2k_{\mathrm{B}}{T}_{\mathrm{R}}}\right)\cos \omega \tau \text{and}{C}_{-}(\tau )={\int }_0^{\infty }d\omega J_{\mathrm{R}}(\omega )\sin \omega \tau .$$

$${\mathcal{K}}_{\mathrm{R}}[\rho (t)]=\left[S,\rho (t)\zeta \right]+\left[{\zeta }^{†}\rho (t),S\right],$$

$$\zeta =\frac{\pi }{2}\sum _{jk}{J}_{\mathrm{R}}({\lambda }_{jk})\left[\coth (\frac{{\lambda }_{jk}}{2k_{\mathrm{B}}{T}_{\mathrm{R}}})+1\right]S_{jk}\text{outerproduct}{\psi }_j{\psi }_k,$$

$$\begin{array}{rl}{\mathcal{K}}_{\mathrm{EM}}[\rho (t)]=& -\underset{0}{\overset{\infty }{\int }}d\tau \left(\left[{\sigma }^{†},\sigma (-\tau )\rho (t)\right]C_{12}(\tau )+\left[\rho (t)\sigma (-\tau ),{\sigma }^{†}(t)\right]C_{21}(-\tau )\right)\\ & -\underset{0}{\overset{\infty }{\int }}d\tau \left(\left[\sigma ,{\sigma }^{†}(-\tau )\rho (t)\right]C_{21}(\tau )+\left[\rho (t){\sigma }^{†}(-\tau ),\sigma (t)\right]C_{12}(-\tau )\right).\end{array}$$

$$C_{12}(\tau )=\underset{0}{\overset{\infty }{\int }}d\tau \mathcal{J}(\omega )(n(\omega )+1)e^{-i\omega \tau }\text{and}{C}_{21}(\tau )=\underset{0}{\overset{\infty }{\int }}d\tau \mathcal{J}(\omega )n(\omega )e^{i\omega \tau },$$

$${\mathcal{K}}_{\mathrm{EM}}[\rho (t)]=-\left[{\sigma }^{†},{\chi }_1\rho (t)\right]-\left[\sigma ,{\chi }_2\rho (t)\right]+h.c.,$$

$${\chi }_1$$

$${\chi }_2$$

$${\chi }_1^{add.}={\Gamma }_{\downarrow }(ϵ)\sigma ,$$

$${\chi }_2^{add.}={\Gamma }_{\uparrow }(ϵ){\sigma }^{†},$$

$${\mathcal{K}}_{\mathrm{EM}}[\rho (t)]=\frac{{\Gamma }_0}{2}(n(ϵ)+1){\mathcal{L}}_{\sigma }[\rho (t)]+\frac{{\Gamma }_0}{2}n(ϵ){\mathcal{L}}_{{\sigma }^{†}}[\rho (t)],$$

$$H^{\mathrm{IBM}}=\left(ϵ+\sum _k\frac{g_k^2}{{\nu }_k}\right)|e⟩⟨e|+|e⟩⟨e|\otimes \sum _k{g}_k(b_k^{†}+b_k)+\sum _k{\omega }_k{b}_k^{†}{b}_k,$$

$${\rho }_{\mathrm{PH}}(0)=\exp (-\sum _k{\omega }_k{b}_k^{†}{b}_k/k_{\mathrm{B}}{T}_R)/\text{tr}[\exp (-\sum _k{\omega }_k{b}_k^{†}{b}_k/k_{\mathrm{B}}{T}_R)],$$

$${\rho }_{eg}(t)={\rho }_{eg}(0)e^{-iϵt}{e}^{-\Gamma (t)}\text{and}{\rho }_{ge}(t)=({\rho }_{eg}(t){)}^{\ast },$$

$$\Gamma (t)=-{\int }_0^{\infty }d\omega \frac{J(\omega )}{{\omega }^2}\coth (\frac{\omega }{2k_{B}T})(1-\cos (\omega t)).$$

Full text

The Franck-Condon (FC) principle May and Kühn (2004); Nitzan (2006) is an invaluable tool in the study of solid-state and molecular emitters. The principle states that electronic transitions of an emitter occur without changes to the motions of its nuclei or those of its environment. As a result, transition rates become dependent on the overlap between vibrational configurations in the initial and final states, which are generally displaced from one another [see Fig. 1(a)]. This picture provides an intuitive starting point for studying the complex interactions between the electronic and vibrational degrees of freedom of an emitter and its environment, for example through rate equations derived from Fermi’s Golden Rule May and Kühn (2004); Nitzan (2006).

Faithfully representing the full non-equilibrium dynamics of such systems requires moving beyond rate equations and instead employing an explicitly time-dependent approach. This should be non-perturbative in the electron-vibrational coupling and thus capable of capturing the dynamical influence of vibrational displacement on the electronic states. Examples include polaron Jang and Cheng (2012); Nazir and McCutcheon (2016) and collective coordinate Iles-Smith et al. (2014, 2016); Wertnik et al. (2018) master equations, hierarchical equations of motion Tanimura and Kubo (1989); Ishizaki and Tanimura (2005); Ishizaki and Fleming (2009), path integrals Makri and Makarov (1995a, b); Nalbach et al. (2011), and tensor network methods Rosenbach et al. (2016); Schröder and Chin (2016); Strathearn et al. (2018). Nevertheless, it is interactions with the electromagnetic environment that ultimately give rise to the observed electronic (e.g. optical) transitions. Our focus is thus on the important question of how to incorporate electromagnetic interactions into the dynamical formalism, such that they respect the non-perturbative nature of the vibrational coupling.

Given that interactions with the electromagnetic field in free space are weak, it is often assumed that the Markovian dynamics they generate can be added to the equations of motion unmodified due to the presence of vibrations Kreisbeck et al. (2011); Dijkstra and Tanimura (2012); Fassioli et al. (2012); Kaer et al. (2012); Ulhaq et al. (2013); McCutcheon and Nazir (2013); Betzholz et al. (2014); Killoran et al. (2015); Chen et al. (2016); Barth et al. (2016); Qin et al. (2017); Stones et al. (2017); Juhász and Csurgay (2018); Chan et al. (2018). Though justifiable in certain circumstances Scala et al. (2007a, b); McCutcheon and Nazir (2013); Giusteri et al. (2017); Kołodyński et al. (2018), additivity is in general a stringent requirement Giusteri et al. (2017); Kołodyński et al. (2018); Mitchison and Plenio (2018); Hofer et al. (2017); Gonz lez et al. (2017) that can break down even if all environments are weakly coupled to the system Mitchison and Plenio (2018). In fact, we shall show below that the dynamics obtained in this manner can exhibit fundamental flaws, such as disregarding the FC principle. In certain cases, both the vibrational and electromagnetic environments may be treated non-perturbatively Olšina et al. ; del Pino et al. (2018), but this comes at an inevitable cost in terms of computational effort and complexity within the formalism.

Here we seek to retain both the simplicity of the Markovian description of the electromagnetic interactions and a non-perturbative treatment of the electron-vibrational coupling, but without the undesirable additivity restriction. This is made possible through a collective coordinate (CC) transformation Garg et al. (1985); Thoss et al. (2001); Iles-Smith et al. (2014); Wertnik et al. (2018), which incorporates non-perturbative effects of the vibrational environment into an enlarged (augmented) system [see Fig. 1(b),(c)]. This in turn enables a Markovian master equation to be derived in the eigenbasis of the augmented system space, rather than that of the original bare emitter, by tracing out the electromagnetic environment and residual vibrational modes Wertnik et al. (2018). On doing so we find that electromagnetic transitions become sensitive to the non-perturbative vibrational dynamics captured by the CC mapping, and our procedure thus retains the non-additive effects crucial to obtaining quantum dynamics that are consistent with the FC principle [Fig. 1(c)]. If instead we enforce additivity [Fig. 1(b)], the resulting electronic decay dynamics becomes independent of the electron-vibrational coupling. We show that capturing non-additivity is also vital for accurately representing the stationary non-equilibrium behaviour within our model. Specifically, under incoherent electromagnetic excitation the non-additive interplay between the electromagnetic field and vibrations directly enables electronic population inversion for situations impossible within the additive approach.

Before examining non-equilibrium dynamics explicitly, we can illustrate the shortcomings of an additive approximation through arguments based on a simple Fermi Golden Rule calculation. We consider a two-level molecular emitter (a monomer) with electronic excited state $\ket{e}$ and ground state $\ket{g}$, separated by an energy $\epsilon$ ($\hbar=1$ throughout). Coupling to the electromagnetic field induces transitions between the electronic states, which are also assumed to couple with strength $\eta$ to a single (harmonic) vibrational mode of frequency $\Omega$, leading to the formation of a displaced manifold associated to the excited electronic configuration. This is the situation depicted qualitatively in Fig. 1(a), though our considerations here and throughout the rest of the paper also apply in the case of continuum phonon environments, where the discrete mode would be identified as the CC post mapping (see below).

We assume for the purpose of calculating the rate that shortly after excitation the system has relaxed to thermal equilibrium in the excited state manifold, $\rho_{\rm eq}=\sum_{m}p_{m}\ket{e,\tilde{m}}\!\bra{e,\tilde{m}}$, where $p_{m}=e^{-m\Omega/k_{\rm B}T}/\sum_{n}e^{-n\Omega/k_{\rm B}T}$ with temperature $T$, and the displaced vibrational basis is denoted $\ket{\tilde{m}}=D(\eta/\Omega)\ket{m}$ for vibrational Fock state $\ket{m}$ and displacement operator $D(\alpha)$. From Fermi’s Golden Rule the electronic excited to ground state decay rate is then May and Kühn (2004)

[TABLE]

There are two principal components to this expression. One is the overlap between vibrational configurations, $|\!\bra{\tilde{m}}n\rangle\!|^{2}$, which is known as the FC factor. The other is the electromagnetic spectral density $\mathcal{J}(\omega)$. This describes the system-field coupling strength weighted by the electromagnetic density of states, and should be sampled at all energy differences between relevant states in the excited and ground manifolds, $\Delta\omega_{\tilde{m},n}$. In the additive approximation, however, the electromagnetic field coupling is treated in isolation from the vibrational interactions, and the electromagnetic spectral density is then incorrectly sampled only at the single frequency $\epsilon$ corresponding to the bare electronic ground and excited state splitting. The expression for the emission rate then reduces to $\Gamma_{e\rightarrow{}g}\approx\mathcal{J}_{0}\sum_{n}\left|\bra{\tilde{0}}n\rangle\right|^{2}=\mathcal{J}_{0}$, where $\mathcal{J}_{0}=\mathcal{J}(\epsilon)$ and we have used $\sum_{n}\ket{n}\!\bra{n}=\mathds{1}$. Thus, in the additive case the FC factor vanishes, and the transition rate loses its dependence on the electron-vibrational coupling. Note that this reasoning can be used to show that the flat spectral density approximation commonly used in quantum optics theory Carmichael (2009) also fails in regimes of strong coupling to vibrational modes.

We now develop a microscopic description in order to establish the extent to which non-additivity can influence the quantum dynamics of electron-vibrational models beyond the heuristic arguments outlined above. Our Hamiltonian is written as $H=H_{\rm S}+H_{\rm I}+H_{\rm B}$, with system Hamiltonian $H_{\rm S}=\epsilon\ket{e}\!\bra{e}$. The electronic configuration of the emitter molecule is directly influenced by both vibrational and electromagnetic environments, such that $H_{\rm I}=H_{\rm I}^{\rm PH}+H_{\rm I}^{\rm EM}$. Within the harmonic approximation the electron-vibrational coupling is written

[TABLE]

where ${b}_{k}$ is the annihilation operator for the $k^{\rm th}$ phonon mode and the second term shifts the excited state due to the reorganisation energy associated to vibrational displacement. The coupling to the phonon environment is characterised by its spectral density, for which we take the common form $J(\nu)=\sum_{k}|{g_{k}}|^{2}\delta(\nu-\nu_{k})=\alpha\nu_{0}^{2}\gamma\nu/\left[(\nu^{2}-\nu_{0}^{2})^{2}+\gamma^{2}\nu^{2}\right]$. Here $\alpha$ and $\nu_{0}$ define the coupling strength and peak position, respectively, and $\gamma$ controls whether $J(\nu)$ is narrow (underdamped) or broad (overdamped) Dijkstra and Tanimura (2015); Iles-Smith et al. (2016). In addition to the phonon environment, we also have an explicit coupling to the electromagnetic field, given by $H_{\rm I}^{\rm EM}=-{\bf d}\cdot{\bf E}$ in the dipole approximation, where ${\bf d}$ is the emitter dipole operator and ${\bf E}$ is the electric field operator Loudon (2000); Breuer and Petruccione (2002); Agarwal (2013). Ignoring polarisation degrees of freedom and working in the rotating wave approximation, this then takes the form

[TABLE]

where ${\sigma}=\ket{g}\!\bra{e}$ and ${a}_{l}$ is the annihilation operator for the $l^{\rm th}$ mode of the electromagnetic field. The spectral density for the light-matter coupling is defined as $\mathcal{J}(\omega)=\sum_{l}|f_{l}|^{2}\delta(\omega-\omega_{l})=(2\pi\epsilon^{3})^{-1}\Gamma_{0}\omega^{3}$ Loudon (2000); Breuer and Petruccione (2002); Agarwal (2013), where $\Gamma_{0}$ is the spontaneous emission rate for the two-level emitter in the absence of phonons. Finally, $H_{\rm B}=H_{\rm B}^{\rm EM}+H_{\rm B}^{\rm PH}=\sum_{l}\omega_{l}{a}_{l}^{\dagger}{a}_{l}+\sum_{k}\nu_{k}{b}_{k}^{\dagger}{b}_{k}$ is the sum of the internal Hamiltonians for the electromagnetic and vibrational environments.

Applying a CC mapping to the phonon bath allows us to incorporate its influence on the electronic dynamics non-perturbatively, and in particular to capture the resulting dynamical generation of electron-vibrational correlations Garg et al. (1985); Thoss et al. (2001); Iles-Smith et al. (2014, 2016). Our Hamiltonian maps as $H=H_{\rm S}+H_{\rm I}+H_{\rm B}\rightarrow H_{\rm S}^{\prime}+H_{\rm I}^{R}+H^{\rm EM}_{\rm I}+H^{R}_{\rm B}+H_{\rm B}^{\rm EM}$, which leaves the light-matter coupling unchanged. Here, we have introduced the transformed Hamiltonians

[TABLE]

where ${b}+b^{\dagger}=\sum_{k}g_{k}({b}_{k}^{\dagger}+{b}_{k})/\eta$ defines creation and annihilation operators for the CC, ${c}_{m}$ is the annihilation operator for the $m^{\rm th}$ mode of the residual environment to which it couples, and we have expressed the reorganisation energy as $\sum_{k}{g_{k}^{2}}/{\nu_{k}}=\int_{0}^{\infty}d\nu J(\nu)/\nu=\pi\alpha/2$. The CC parameters can be written in terms of the quantities defining the vibrational spectral density: $\eta^{2}=\pi\alpha\nu_{0}/2$ and $\Omega=\nu_{0}$ Iles-Smith et al. (2016). Coupling between the augmented emitter-CC system and the residual phonon environment is described by an Ohmic spectral density $J_{\rm R}(\nu)=\sum_{m}|h_{m}|^{2}\delta(\nu-\tilde{\nu}_{m})=\gamma\nu/2\pi\nu_{0}$ Iles-Smith et al. (2016), and ensures that the vibrational environment still acts as a continuum of modes after the mapping. As in the single mode case discussed earlier, the coupling to the CC leads to the formation of two vibronic manifolds associated to the ground and excited electronic configurations. The coupling to the residual environment induces transitions within each vibronic manifold. This leads both to broadening and to dynamical relaxation of the phonon environment, which typically occurs on a sub-picosecond timescale.

From the mapped Hamiltonian we derive a second-order Born-Markov master equation by tracing over the residual environment and the electromagnetic field Breuer and Petruccione (2002), both of which are assumed to remain in thermal equilibrium at temperatures $T_{\rm R}$ and $T_{\rm EM}$, respectively. The resulting master equation can be written $\partial_{t}\rho(t)=\mathcal{L}[\rho(t)]$ with Liouvillian 111For further details see the Supplemental Material:

[TABLE]

where $\rho(t)$ is the reduced state of the augmented emitter-CC system. Here, $\mathcal{K}_{\rm R}$ is a superoperator representing the action of the residual phonon environment Iles-Smith et al. (2016):

[TABLE]

with $S={b}^{\dagger}+{b}$ and

[TABLE]

where the eigenbasis of the augmented system is defined through $H_{\rm S}^{\prime}\ket{\psi_{j}}=\psi_{j}\ket{\psi_{j}}$, giving $\lambda_{jk}=\psi_{j}-\psi_{k}$ and $S_{jk}=\bra{\psi_{j}}S\ket{\psi_{k}}$. We solve for the eigenvalues $\psi_{j}$ and eigenstates $\ket{\psi_{j}}$ numerically, taking the basis $\{\ket{g},\ket{e}\}$ for the TLE and a Fock (number) state basis for the CC.

The effects of the electromagnetic field interaction are contained within $\mathcal{K}_{\rm EM}$. Importantly, the augmented emitter-CC system Hamiltonian, $H_{\rm S}^{\prime}$, is treated (numerically) exactly within the formalism. This is crucial in capturing non-additive effects of the electromagnetic and vibrational environments, as it means that when we move the electromagnetic interaction Hamiltonian [Eq. (3)] into the interaction picture, we do so with respect to the full augmented system Hamiltonian $H_{\rm S}^{\prime}$ [Eq. (4)]. The mapping thus ensures that the electromagnetic environment is sensitive to the underlying eigenstructure of both the electronic and vibrational states. The superoperator representing the dynamical influence of the electromagnetic environment then takes the form:

[TABLE]

where we have introduced the rate operators $\chi_{1}=\sum_{jk}\sigma_{jk}\Gamma_{\downarrow}(\lambda_{jk})\ket{\psi_{j}}\!\bra{\psi_{k}}$ and $\chi_{2}=\sum_{jk}\sigma^{\ast}_{jk}\Gamma_{\uparrow}(\lambda_{jk})\ket{\psi_{k}}\!\bra{\psi_{j}}$, with transition rates given by $\Gamma_{\downarrow}(\lambda)=\pi\mathcal{J}(\lambda)(n(\lambda)+1)$ and $\Gamma_{\uparrow}(\lambda)=\pi\mathcal{J}(\lambda)n(\lambda)$, for field occupation number $n(\lambda)=(\exp{\lambda/k_{\rm B}T_{\rm EM}}-1)^{-1}$ and $\sigma_{jk}=\bra{\psi_{j}}\sigma\ket{\psi_{k}}$. From these expressions, it is evident that the interaction between the system and the electromagnetic field is dependent on the eigenstructure of the augmented system, and thus on the emitter-vibrational coupling through the identification of the CC, and its coupling to the electronic system. We therefore refer to this theory as being non-additive.

In contrast, within the additive approach $\mathcal{K}_{\rm EM}$ is derived without reference to the vibrational coupling. In the present setting, this amounts to neglecting the modification of the system eigenstructure due to vibrational interactions encoded in the mapped system Hamiltonian $H_{\rm S}^{\prime}$, and instead moving the electromagnetic interaction Hamiltonian [Eq. (3)] into the interaction picture with respect to the original system Hamiltonian $H_{\rm S}$. This results in the standard Lindblad dissipator common in quantum optics theory: $\mathcal{K}_{\rm EM}[\rho(t)]=\frac{\Gamma_{0}}{2}(n(\epsilon)+1)\mathcal{L}_{\sigma}[\rho(t)]+\frac{\Gamma_{0}}{2}n(\epsilon)\mathcal{L}_{\sigma^{\dagger}}[\rho(t)]$, where $\mathcal{L}_{O}[\rho]=2O\rho O^{\dagger}-\{O^{\dagger}O,\rho\}$. It is thus clear that within the additive approximation the electromagnetic field superoperator loses its explicit dependence on the vibrational environment. For vanishing electromagnetic interactions ($\Gamma_{0}\to 0$), the additive and non-additive theories become equivalent and the problem reduces to the independent boson model, for which an exact solution can be obtained. We verify that the CC master equation agrees with this exact solution for a range of parameters in the Supplemental Material.

We are now in a position to investigate the impact of non-additive effects on the dynamics of our model system. We begin by considering the decay of an emitter initialised in its excited state with the collective coordinate in a thermal state set by the residual bath temperature $T_{\rm R}$: $\rho(0)=\ket{e}\!\bra{e}\otimes\rho_{th}$, where $\rho_{th}=\exp(-\Omega{b}^{\dagger}{b}/k_{\rm B}T_{\rm R})/\tr[\exp(-\Omega{b}^{\dagger}{b}/k_{\rm B}T_{\rm R})]$. This approximates a canonical thermal state of the original vibrational Hamiltonian in the unmapped representation at the same temperature, and is thus consistent with rapid (vertical) excitation of the system whereby the electronic state changes suddenly but the vibrational states remain unchanged Note (1). The vibrational environment will subsequently relax towards the displaced thermal state associated to the excited state manifold, captured dynamically within our approach.

Fig. 2(a) shows the emitter excited state population dynamics predicted by the additive (dotted) and non-additive (solid) theories for increasing electron-phonon coupling at ambient temperature. Both theories give rise to exponential decay, with the rate in the additive theory remaining constant across all electron-phonon coupling strengths. The non-additive theory, in contrast, displays a monotonic decrease in the decay rate with increasing phonon coupling. This can be seen explicitly in Fig. 2(b), where we extract the decay rates directly from the master equation. Specifically, the excited to ground state transition rate can be written as

[TABLE]

with the Liouvillian taken to be additive or non-additive depending on which case is under investigation. Here we must modify the initial state to account for the aforementioned rapid residual bath induced relaxation of the CC to a displaced thermal state prior to emission: $\rho_{X}(0)=\ket{e}\!\bra{e}\otimes e^{-X}\rho_{th}e^{X}$, where $X=\Omega^{-1}\eta(b^{\dagger}-b)$. As expected, the rate from the additive theory displays no variation with phonon coupling strength, in line with the simple Golden Rule calculation discussed previously but at odds with the FC principle. This again highlights deficiencies with the phenomenological additive treatment of the electromagnetic field. Conversely, the non-additive theory shows a steady reduction of the emission rate as a function of phonon coupling, consistent with FC physics. As the displacement between the ground and excited state manifolds increases linearly with the electron-phonon coupling strength, this reduces the overlap between the vibrational states and thus suppresses electromagnetic transitions.

It is important to stress that discrepancies between the additive and non-additive treatments in our model extend further than spontaneous emission processes. For example, we now consider situations in which the emitter is driven incoherently via thermal occupation of the electromagnetic environment, which constitutes an important building block of widely used models for natural and artificial solar energy conversion. Fig. 3(a) shows the steady state population of the electronic excited state as a function of electron-phonon coupling strength, where the additive treatment once again displays no variation, simply matching the equilibrium distribution expected in the absence of vibrations. The non-additive treatment, on the other hand, shows a monotonic increase in the steady state population of the excited state manifold. Most strikingly, at large coupling strengths there emerges a steady state population inversion. In the absence of phonons, such an inversion would be impossible, with emission and absorption processes balancing each other in equilibrium. This remains true in the presence of phonons when the electromagnetic field is treated additively, as highlighted in Fig. 3(b). Here, the additive theory approaches, but never exceeds, a maximum steady state population $\langle{\sigma}^{\dagger}{\sigma}\rangle=0.5$ in the limit of very large temperatures. In contrast, within the non-additive theory, cooperative effects between the electromagnetic and vibrational environments lead to non-equilibrium stationary states that display substantial levels of population inversion. That such effects should be possible, even for continuum environments, is made clear from the CC mapping. Within the non-additive theory the electromagnetic field has access to the full vibrational structure of the emitter, providing the necessary states to drive a population inversion. This points to a crucial difference between non-additive and additive treatments, where disregarding the eigenstructure of the combined electronic and vibrational system misses key aspects of the non-equilibrium physics.

In summary, we have demonstrated that for models of electronic systems strongly coupled to vibrational environments, including the electromagnetic field in an additive manner can lead to dynamics inconsistent with the FC principle. By developing a dynamical formalism based on collective coordinate mappings, we capture the impact of non-additive effects to recover both transient and stationary non-equilibrium behaviour consistent with FC physics. Furthermore, we find that for common model assumptions on the forms of vibrational and electromagnetic couplings, non-additive phenomena enable steady-state population inversion under incoherent electromagnetic excitation conditions. Our findings open up a number of avenues for further exploration in both natural and artificial light harvesting, for example whether non-additive effects may be harnessed to enhance work extraction in models of quantum heat engines and solar energy conversion devices.

*Acknowledgements.—*The authors wish to thank Adam Stokes and Neill Lambert for useful discussions. H.M. is supported by the EPSRC. A.N. and J.I.S. are supported by the EPSRC, grant no. EP/N008154/1. J.I.S. also acknowledges support from the Royal Commission for the Exhibition of 1851.

Appendix A Supplemental Material

A.1 Further details of the master equation derivation

A.1.1 Non-Additive

In this supplement we outline further details of the master equations used in the main text. We begin by considering the non-additive case, whereby, after the collective-coordinate (CC) mapping, our interaction Hamiltonian is given by the sum of residual phonon bath and electromagnetic field coupling terms [Eqs. (3) and (5) in the main text], $H_{\rm I}=S\otimes B+S^{2}\sum_{m}{h_{m}^{2}}/{\tilde{\nu}_{m}}+\sum_{\alpha}A_{\alpha}\otimes E_{\alpha}$. Here, $S={b}^{\dagger}+{b}$ and $B=\sum_{m}h_{m}({c}_{m}^{\dagger}+{c}_{m})$ describe the CC-residual bath interaction, while $A_{1}={\sigma}^{\dagger}$, $A_{2}={\sigma}$, $E_{1}=\sum_{l}f_{l}{a}_{l}$ and $E_{2}=\sum_{l}f^{\ast}_{l}{a}_{l}^{\dagger}$ define the coupling of the two-level emitter (TLE) to the electromagnetic field. We now move into the interaction picture with respect to the augmented system Hamiltonian describing the coupled TLE and CC, $H_{\rm S}^{\prime}=H_{\rm S}+\eta\ket{e}\!\bra{e}({b}^{\dagger}+{b}+\pi\alpha/2\eta)+\Omega{b}^{\dagger}{b}$, plus the residual and electromagnetic bath Hamiltonians $H_{\rm B}^{R}=\sum_{m}\tilde{\nu}_{m}{c}_{m}^{\dagger}{c}_{m}$ and $H_{\rm B}^{\rm EM}=\sum_{l}\omega_{l}{a}_{l}^{\dagger}{a}_{l}$, respectively. This gives

[TABLE]

where $S(t)=e^{iH_{\rm S}^{\prime}t}Se^{-iH_{\rm S}^{\prime}t}$, $A_{\alpha}(t)=e^{iH_{\rm S}^{\prime}t}A_{\alpha}e^{-iH_{\rm S}^{\prime}t}$, $B(t)=\sum_{m}h_{m}({c}_{m}^{\dagger}e^{i\tilde{\nu}_{m}t}+{c}_{m}e^{-i\tilde{\nu}_{m}t})$, $E_{1}=\sum_{l}f_{l}{a}_{l}e^{-i\omega_{l}t}$ and $E_{2}=\sum_{l}f^{\ast}_{l}{a}_{l}^{\dagger}e^{i\omega_{l}t}$. Within the interaction picture, we then follow the standard procedure to derive a Redfield master equation, tracing out the residual and electromagnetic environments within the Born-Markov approximations Breuer and Petruccione (2002). Moving back into the Schrödinger picture, the resulting master equation may be written in the general form

[TABLE]

where $\rho(t)$ is the reduced density operator of the augmented system, from which either the TLE or CC dynamics may be obtained by tracing out the relevant degrees of freedom. The superoperators $\mathcal{K}_{\rm R}$ and $\mathcal{K}_{\rm EM}$ encode, respectively, the influence of the residual bath and the electromagnetic field interactions on the augmented system dynamics. Note that due to the Born-Markov approximations there are no mixed terms between the residual phonon bath and the electromagnetic field in the master equation above. Nevertheless, our master equation is still non-additive with respect to the original phonon environment and the electromagnetic field due to the CC mapping, which incorporates non-perturbative vibrational effects into the enlarged augmented system Hamiltonian $H_{\rm S}^{\prime}$ used to move into the interaction picture. This results in an electromagnetic superoperator that has explicit dependence on the form and strength of the system-vibrational coupling, as we shall see below.

Starting with the influence of the residual phonon environment, we find

[TABLE]

with residual bath correlation functions defined as

[TABLE]

where the expectation value is taken with respect to a thermal state at temperature $T_{\rm R}$. These can be found given the form of $B(t)$ above to be

[TABLE]

Here, we have defined the residual bath spectral density $J_{\rm R}(\nu)=\sum_{m}|h_{m}|^{2}\delta(\nu-\tilde{\nu}_{m})$ as in the main text. We now decompose the system operators into the eigenbasis of the augmented TLE-CC Hamiltonian, $H_{\rm S}^{\prime}\ket{\psi_{j}}=\psi_{j}\ket{\psi_{j}}$, giving $S(t)=\sum_{jk}S_{jk}e^{i\lambda_{jk}t}\outerproduct{\psi_{j}}{\psi_{k}}$, where $S_{jk}=\bra{\psi_{j}}S\ket{\psi_{k}}$ and $\lambda_{jk}=\psi_{j}-\psi_{k}$. Inserting this decomposition into Eq. (14) and exchanging the order of integration over time and frequency, we can now perform the integral over $\tau$ using the Sokhotski-Plemelj theorem, $\int_{0}^{\infty}d\tau e^{\pm i\epsilon\tau}=\pi\delta(\epsilon)\pm iP(1/\epsilon)$, where $P$ stands for the Cauchy principal value. After some algebra, and neglecting imaginary contributions as justified in Refs. Iles-Smith et al. (2014, 2016), we find that the residual bath superoperator may be written

[TABLE]

with

[TABLE]

as given in the main text.

Next, we consider the superoperator for the electromagnetic field, which can be written

[TABLE]

Here, the bath correlation functions are defined as $C_{\alpha\alpha^{\prime}}(\tau)=\langle E_{\alpha}(\tau)E_{\alpha^{\prime}}\rangle$ with the expectation value taken with respect to a thermal state at temperature $T_{\rm EM}$. This gives

[TABLE]

where we have defined the electromagnetic field spectral density $\mathcal{J}(\omega)=\sum_{l}|f_{l}|^{2}\delta(\omega-\omega_{l})$ and the field occupation number $n(\omega)=(\exp{\omega/k_{\rm B}T_{\rm EM}}-1)^{-1}$, again as in the main text.

For the non-additive treatment of the electromagnetic field interactions we again decompose the system operators into the eigenbasis of the augmented TLE-CC Hamiltonian such that $\sigma(t)=\sum_{jk}\sigma_{jk}e^{i\lambda_{jk}t}\ket{\psi_{j}}\!\bra{\psi_{k}}$, where $\sigma_{jk}=\bra{\psi_{j}}\sigma\ket{\psi_{k}}$. Inserting this decomposition into Eq. (19) and exchanging the order of the time and frequency integrals, we perform the integral over $\tau$ using the Sokhotski-Plemelj theorem as before. Again, after some algebra and neglecting small imaginary terms, we find

[TABLE]

with rate operators

[TABLE]

where $\Gamma_{\downarrow}(\lambda)=\pi\mathcal{J}(\lambda)(n(\lambda)+1)$ and $\Gamma_{\uparrow}(\lambda)=\pi\mathcal{J}(\lambda)n(\lambda)$ as in the main text. Note that these non-additive rate operators account for the full eigenstructure of the electronic and CC degrees of freedom.

A.1.2 Additive

Within the additive master equation the description of the residual phonon bath is unchanged, as $\mathcal{K}_{\rm R}[\rho(t)]$ has no dependence on the electromagnetic field coupling. However, the electromagnetic field superoperator $\mathcal{K}_{\rm EM}[\rho(t)]$ is altered significantly, as it is no longer sensitive to the (non-perturbative) vibrational coupling.

Specifically, to get $\mathcal{K}_{\rm EM}[\rho(t)]$ in the additive case, one should completely ignore the presence of vibrations when moving the relevant system operators $\sigma$ and $\sigma^{\dagger}$ into the interaction picture. This results in $\sigma(t)=e^{iH_{\rm S}t}\sigma e^{-iH_{\rm S}t}=\sigma e^{-i\epsilon t}$ and $\sigma^{\dagger}(t)=e^{iH_{\rm S}t}\sigma^{\dagger}e^{-iH_{\rm S}t}=\sigma^{\dagger}e^{i\epsilon\tau}$, where $H_{\rm S}=\epsilon\ket{e}\bra{e}$. The impact is a simplification of the rate operators to

[TABLE]

and

[TABLE]

which when inserted into Eq. (21) results in the standard Lindblad form common in quantum optics theory:

[TABLE]

with $\mathcal{L}_{O}[\rho]=2O\rho O^{\dagger}-\{O^{\dagger}O,\rho\}$.

A.2 Independent boson model: comparison of the collective coordinate master equation and exact solutions

In the absence of coupling to an optical environment ($\Gamma_{0}=0$), the model outlined in the manuscript reduces to the independent boson model

[TABLE]

which is exactly solvable Breuer and Petruccione (2002), for example via a polaron transformation of the form $H^{\prime}=U^{\dagger}HU$, where $U=\ket{e}\!\bra{e}\prod_{k}D(-g_{k}/\omega_{k})+\ket{g}\!\bra{g}$. This removes the interaction term, permitting a partial trace to be taken over the phonon modes, which we consider here to be initialised in the thermal state

[TABLE]

with total initial density matrix $\chi(0)=\rho_{\rm S}(0)\rho_{\rm PH}(0)$. Since the interaction term in Eq. (27) commutes with the system Hamiltonian, within the independent boson model the phonon bath causes no transitions between electronic eigenstates and thus the system populations are static, $\rho_{gg}(t)=\rho_{gg}(0)$ and $\rho_{ee}(t)=\rho_{ee}(0)$, where $\rho_{ij}(t)=\langle i|\rho_{\rm S}(t)|j\rangle$ for $i,j\in\{g,e\}$ are system density matrix elements. The system coherences do evolve, and are governed by

[TABLE]

where the decoherence function is defined as

[TABLE]

The dynamics of these exact solutions can be compared to the collective coordinate master equation used in the main manuscript by simply setting the electromagnetic coupling $\Gamma_{0}$ to zero in the latter. In Fig. 4 above we compare the two approaches using the same parameters as in Fig. 2 of the manuscript (though with $\Gamma_{0}=0$) and find that the collective coordinate master equation matches the exact solution well into the strong phonon-coupling regime. The initial condition for the two-level system is taken to be $\rho_{\rm S}(0)=\frac{1}{2}(\ket{e}+\ket{g})\!(\bra{e}+\bra{g})$. For the exact solution, the environment is initialised in the thermal state $\rho_{\rm PH}(0)$ of Eq. (28), as stated before, while for the collective coordinate master equation the CC is initialised in the thermal state $\rho_{th}=\exp(-\Omega{b}^{\dagger}{b}/k_{\rm B}T_{\rm R})/\tr[\exp(-\Omega{b}^{\dagger}{b}/k_{\rm B}T_{\rm R})]$ with the residual environment held in thermal equilibrium at temperature $T_{R}$ throughout the dynamics.