Capturing Vacuum Fluctuations and Photon Correlations in Cavity Quantum Electrodynamics with Multi-Trajectory Ehrenfest Dynamics

TL;DR

This paper demonstrates that multi-trajectory Ehrenfest dynamics can effectively simulate vacuum fluctuations and photon correlations in cavity QED systems, providing a promising semi-classical approach for quantum light-matter interactions.

Contribution

It extends mixed quantum-classical Ehrenfest methods to model quantum vacuum effects and photon correlations in cavity QED, showing accurate results compared to exact quantum benchmarks.

Findings

Accurately captures spontaneous emission dynamics.

Reproduces quantum correlations and vacuum fluctuations.

Offers a promising semi-classical simulation approach.

Abstract

We describe vacuum fluctuations and photon-field correlations in interacting quantum mechanical light-matter systems, by generalizing the application of mixed quantum-classical dynamics techniques. We employ the multi-trajectory implementation of Ehrenfest mean field theory, traditionally developed for electron-nuclear problems, to simulate the spontaneous emission of radiation in a model quantum electrodynamical cavity-bound atomic system. We investigate the performance of this approach in capturing the dynamics of spontaneous emission from the perspective of both the atomic system and the cavity photon field, through a detailed comparison with exact benchmark quantum mechanical observables and correlation functions. By properly accounting for the quantum statistics of the vacuum field, while using mixed quantum-classical (mean field) trajectories to describe the evolution, we identify a surprisingly accurate and promising route towards describing quantum effects in realistic correlated light-matter systems.

Equations

$$\hat{H}={\hat{H}}_A+{\hat{H}}_F+{\hat{H}}_{AF}.$$

$${\hat{H}}_A=\sum _k{\epsilon }_k|k⟩⟨k|.$$

$${\hat{H}}_F=\frac{1}{2}\sum _{\alpha =1}^{2N}\left({\hat{P}}_{\alpha }^2+{\omega }_{\alpha }^2{\hat{Q}}_{\alpha }^2\right).$$

$${\hat{Q}}_{\alpha }$$

$${\hat{P}}_{\alpha }$$

$$\hat{H}_{AF}=\sum_{\alpha=1}^{2N}\Big{(}\omega_{\alpha}\hat{Q}_{\alpha}(\lambda_{\alpha}\cdot\hat{\mu})+\frac{1}{2}\Big{(}\lambda_{\alpha}\cdot\hat{\mu}\Big{)}^{2}\Big{)},$$

$$\frac{\partial }{\partial t}{\hat{\rho }}_W(X,t)=-i\mathcal{L}{\hat{\rho }}_W(X,t).$$

$${\hat{\rho }}_W(Q,P)=\frac{1}{(2\pi \hslash {)}^{2N}}\int dZ{e}^{iP\cdot Z}⟨Q-\frac{Z}{2}|\hat{\rho }|Q+\frac{Z}{2}⟩.$$

$$i\mathcal{L}\cdot =\frac{i}{\hslash }[{\hat{H}}_W,\cdot ]-\frac{1}{2}(\{{\hat{H}}_W,\cdot \}-\{\cdot ,{\hat{H}}_W\}),$$

$${\hat{\rho }}_W(X,t)={\hat{\rho }}_A(t){\rho }_{F,W}(X,t),$$

$$\hat{\rho}_{A}(t)=\text{Tr}_{F}\Big{(}\hat{\rho}(t)\Big{)}=\int dX\hat{\rho}_{W}(X,t),$$

$$\frac{d}{dt}\hat{\rho}_{A}(t)=-i\Big{[}\hat{H}_{A}+\hat{H}_{AF,W}(X(t)),\hat{\rho}_{A}(t)\Big{]},$$

$${\rho }_{F,W}(X,t)=\sum _j^{N_{traj}}{w}^j\delta (X-X^j(t))$$

$$\frac{\partial Q_{\alpha }}{\partial t}=\frac{\partial H_{F,W}^{Eff}}{\partial P_{\alpha }},\frac{\partial P_{\alpha }}{\partial t}=-\frac{\partial H_{F,W}^{Eff}}{\partial Q_{\alpha }}.$$

$$H^{Eff}_{F,W}=\frac{1}{2}\sum_{\alpha}\Big{(}P_{\alpha}^{2}+\omega_{\alpha}^{2}Q_{\alpha}^{2}-2\omega_{\alpha}\lambda_{\alpha}Q_{\alpha}\mu(t)\Big{)},$$

$$⟨O(t)⟩={\text{Tr}}_A\int dX{\hat{O}}_W(X,t){\hat{\rho }}_W(X,t=0).$$

$${\rho }_{F,W}(X,0)=\prod _{\alpha }\frac{1}{\pi }\exp \left[-\frac{P_{\alpha }^2}{{\omega }_{\alpha }}-{\omega }_{\alpha }{Q}_{\alpha }^2\right].$$

$$:\hat{N}_{pt}:=\frac{1}{2}\sum_{\alpha}\Big{(}\frac{\hat{P}^{2}_{\alpha}}{\omega_{\alpha}}+\omega_{\alpha}\hat{Q}^{2}_{\alpha}-1\Big{)},$$

$$\hat{E}(r,t)=\sum _{\alpha }\sqrt{2{\omega }_{\alpha }}{\zeta }_{\alpha }(r){\hat{Q}}_{\alpha }(t)$$

$$\zeta_{\alpha}(r)=\sqrt{\frac{\hbar\omega_{\alpha}}{\epsilon_{0}L}}\sin\Big{(}\frac{\alpha\pi}{L}r\Big{)}.$$

$$:{\hat{E}}^2(r,t):=:\hat{I}(r,t):=2\sum _{\alpha }{\omega }_{\alpha }{\zeta }_{\alpha }^2(r){\hat{Q}}_{\alpha }^2(t)-\sum _{\alpha }{\zeta }_{\alpha }^2(r).$$

$$:g^2(r_1,r_2,t):=\frac{⟨:{\hat{E}}^{+}(r_1,t){\hat{E}}^{+}(r_2,t)\hat{E}(r_2,t)\hat{E}(r_1,t):⟩}{⟨:\hat{I}(r_1,t):⟩⟨:\hat{I}(r_2,t):⟩}.$$

$$:$$

$$\hat{H}=$$

$$+$$

Full text

I Introduction

Profound changes in the physical and chemical properties of material systems can be produced in situations where the quantum nature of light plays an important role in the interaction with the system Feist2015 ; Schachenmayer2015 ; Cirio2016 . A few notable recent examples of such effects are few-photon coherent nonlinear optics with single molecules maser2016 , direct experimental sampling of electric-field vacuum fluctuations riek2015 ; moskalenko2015 , multiple Rabi splittings under ultra-strong vibrational coupling george2016 , exciton-polariton condensates byrnes2014 ; kasprzak2006 , and frustrated polaritons Schmidt2016 . These exciting developments have been strongly driven by experimental efforts, thus exposing the immediate need for the development and improvement of theoretical approaches that can bridge the gap between quantum optics and quantum chemistry ruggenthaler2018 .

Due to the similarity of the electron-photon and the electron-nuclear problems, simulation methods that have traditionally been of use in the quantum chemistry community, such as semiclassical and mixed quantum-classical methods, offer a potentially interesting avenue to bridge this gap. In particular, the family of trajectory-based quantum-classical methods has the advantage of providing a very intuitive, qualitative understanding of nonadiabatic molecular dynamics. Further, these techniques typically do not exhibit the pernicious exponential scaling of computational effort inherent in grid-based quantum calculationsThoss2004 . Available techniques in this family of exact and approximate approaches are Ehrenfest mean field dynamics, fully linearized and partially linearized path integral methods, forward-backward trajectory methods Hsieh2012 ; he12 ; SKR18 , and trajectory surface-hopping algorithms KM13 . All these techniques have some ability to describe essential quantum mechanical effects such as tunnelling, interference, zero-point energy conservation.

Recently, Subotnik and co-workers have performed investigations of light-matter interactions where an adjusted Ehrenfest theory based method is used to simulate spontaneous emission of classical light CLSNS18 ; CLSNS218 ; LNSMCS18 . Here, in contrast with these works, we focus on the description of quantized light fields. We then generalize the well established multi-trajectory Ehrenfest method to treat quantum mechanical light-matter interactions. We highlight the possibilities and theoretical challenges of this method in comparison to the exact treatment of the quantum system, by applying this approach to investigate spontaneous emission for a model atom in an optical cavity.

The remainder of this work is divided into three sections: in Sec.II we briefly review general interacting light-matter systems, and the multi-trajectory Ehrenfest dynamics method. In this framework, we then introduce a one-dimensional model system comprising a single (two or three level) atomic system coupled to a multi-mode quantum electrodynamical (QED) cavity. In Sec.III we investigate the performance of multi-trajectory Ehrenfest (MTEF) dynamics in describing the process of spontaneous emission. We conclude our results in Sec.IV and discuss some prospects for future work.

II Theory

II.1 Quantum Mechanical Light-Matter Interactions

To begin, we describe a general coupled field-matter system using Coulomb gauge and the dipole approximation Faisal1987 ; Flick2017a . The total Hamiltonian for the system is Tokatly2013 ; Flick2017 ; Pellegrini2015 ; Flick2015 ; Craig1998

[TABLE]

The first term, $\hat{H}_{A}$, is the atomic Hamiltonian, which may be generally expressed in the spectral representation,

[TABLE]

Here $\{\varepsilon_{k},|k\rangle\}$ are the atomic energies and stationary states of the atomic system in absence of coupling to the cavity. The second term is the Hamiltonian of the uncoupled cavity field $\hat{H}_{F}$,

[TABLE]

The photon-field operators, $\hat{Q}_{\alpha}$ and $\hat{P}_{\alpha}$, obey the canonical commutation relation, $[\hat{Q}_{\alpha},\hat{P}_{\alpha^{\prime}}]=\imath\hbar\delta_{\alpha,\alpha^{\prime}}$, and can be expressed using creation and annihilation operators for each mode of the cavity field,

[TABLE]

where $\hat{a}^{\dagger}_{\alpha}$ and $\hat{a}_{\alpha}$ denote the usual photon creation and annihilation operators for photon mode $\alpha$. The coordinate-like operators, $\hat{Q}_{\alpha}$, are directly proportional to the electric displacement operator, while the conjugate momenta-like operators, $\hat{P}_{\alpha}$, are related to the magnetic field Pellegrini2015 ; Flick2015 . The upper limit of the sum in Eq.(3) is $2N$, as there are (in principle) two independent polarization degrees of freedom for each photon mode, however in the 1D cavity models presented here only a single polarization will be considered.

The final term in Eq.(1) represents the coupling between the atom and the cavity field,

[TABLE]

where we denote ${\hat{\mu}}$ as the electronic dipole moment vector of the atomic system, and ${\lambda}_{\alpha}$ as the electron-photon coupling vector Tokatly2013 ; Flick2015 . In the case of a two-level electronic system the quadratic term in the atom-field coupling Hamiltonian simply results in a constant energy shift and hence has no effect on observables SRR18 , and we neglect this term in the case of the three level model system. Furthermore, we note that this Hamiltonian can easily be extended to include nuclear degrees of freedom, however this has been omitted in the present work.

II.2 Multi-Trajectory Ehrenfest Dynamics

In this section we apply the well-known multi-trajectory Ehrenfest method, traditionally introduced to study electron-nuclear systems Ehrenfest1927 ; Makri1987 ; C.Tully1998 , to coupled light-matter systemsMcLachlan1964 ; C.Tully1998 ; Flick2017a .

A particularly simple and instructive route to derive the MTEF mean field theory is via the quantum-classical Liouville (QCL) equationqcle . This equation of motion for the density matrix is formally exact for an arbitrary quantum mechanical system that is bilinearly coupled to a harmonic environment, as is the case in the atom-field Hamiltonian studied here. The QCL equation can be written in a compact form as

[TABLE]

It describes the time evolution of $\hat{\rho}_{W}(X,t)$, which is the partial Wigner transform of the density operator taken over the photon field coordinates, which are thus represented by continuous phase space variables $X=(Q,P)=(Q_{1},Q_{2},...,Q_{2N},P_{1},P_{2},...,P_{2N})$. The partial Wigner transform of the density operator, $\hat{\rho}$, is defined as

[TABLE]

The QCL operator is defined as

[TABLE]

where $\hat{H}_{W}$ denotes the Wigner transform of $\hat{H}$, $[\cdot,\cdot]$ is the commutator, and $\{\cdot,\cdot\}$ is the Poisson bracket in the phase space of the environmental variables.

At this point, one may arrive at MTEF equations by assuming that the total density of the system can be written as an uncorrelated product of the atomic and photonic reduced densities at all times,

[TABLE]

where the reduced density matrix of the atomic system is

[TABLE]

and the Wigner function of the cavity field is $\rho_{F,W}(X,t)=Tr_{A}(\hat{\rho}_{W}(X,t))$. If one seeks solutions to the QCL equation of this form, the Ehrenfest mean-field equations of motion for the atomic system are obtained:

[TABLE]

where $\hat{H}_{AF,W}$ denotes the Wigner transform of $\hat{H}_{AF}$. The evolution of the Wigner function of the photon field can be represented as a statistical ensemble of independent trajectories, with weights $w^{j}$,

[TABLE]

that evolve according to Hamilton’s equations of motion,

[TABLE]

The effective photon field Hamiltonian is,

[TABLE]

where $\mu(t)=\text{Tr}_{A}(\hat{\rho}_{A}(0)\hat{\mu}(t))$.

The exact expression for the average value of any observable, $\langle O(t)\rangle$, can be written as

[TABLE]

We note here that for this class of systems the Ehrenfest equations of motion for the photon field coordinates correspond to a mode resolved form of Maxwell’s equations. In applying the MTEF dynamics method numerically, we use the above expressions in the following manner:

We first perform Monte Carlo sampling from the Wigner transform of the initial density operator of the photon field $\hat{\rho}_{F,W}(X,0)$ to generate an ensemble of initial conditions, for the trajectory ensemble $(Q_{\alpha}^{j}(0),P_{\alpha}^{j}(0))$. In this work we used uniform weights $w^{j}=\frac{1}{N_{traj}}$, however other importance sampling schemes could be employed as the only requirement is that the sum of the weights is normalized, $\sum_{j}w^{j}=1$. 2. 2.

We evolve each initial condition independently according to the Ehrenfest equations of motion, producing a trajectory. In the following we refer to such a solution as an ensemble of independent trajectories. 3. 3.

Average values are constructed by summing over the entire trajectory ensemble, and normalizing the result with respect to $N_{traj}$, the total number of trajectories in the ensemble,

\langle O(t)\rangle=\sum_{j}^{N_{traj}}\text{Tr}_{A}\Big{(}\hat{O}_{W}(Q^{j},P^{j},t)\hat{\rho}_{A}(0)\Big{)}/N_{traj}.

Here $\rho_{F,W}(X,0)$ is Wigner transform of the zero temperature vacuum state

[TABLE]

II.3 Observables and Normal Ordering

Before we proceed with a discussion of our simulation results, we must note that the Wick normal ordered form for operators (denoted $:\hat{O}:$ for some operator $\hat{O}$) is used when calculating average values in this study. The reason for using the normal ordered form, in practice, is to remove the effect of vacuum fluctuations from the results, which ensures that both $\langle E\rangle=0$ and $\langle I\rangle=0$, irrespective of the number of photon modes in the cavity field, when the field is in the vacuum state. The effect of this operator ordering is particularly evident for the photon number operator,

[TABLE]

where normal ordering produces a constant shift due to the zero-point energy term.

The quantized electric field operator is defined as

[TABLE]

with

[TABLE]

The corresponding normal-ordered electric field intensity operator is given by

[TABLE]

The effect of normal ordering on this quantity is shown in Fig. 2, where the intensity of the electric field is plotted in both its canonical and normal ordered forms. In addition to a constant shift with respect to the normal ordered quantity, which is identically zero, the canonical average field intensity also displays additional oscillations near the boundaries and the atomic position, corresponding to the vacuum fluctuations for this system.

We also consider the second order correlation function for the photon field Glauber ,

[TABLE]

This function is frequently used in quantum optics to discriminate between classical light and non-classical states of the photon field that exhibit photon bunching $(g^{2}>1)$ or photon anti-bunching $(g^{2}<1)$. The normal-ordered form of the numerator in $g^{2}$, also referred to $G^{2}(r_{1},r_{2},t)$, is

[TABLE]

The partial Wigner transforms of the polynomial functions of the bath coordinate operators are simply polynomial functions of the continuous bath coordinates, \Big{(}\hat{Q}_{\alpha}^{n}(t)\Big{)}_{W}=(Q_{\alpha}(t))^{n} Hillery . The same is also true for the corresponding momenta, and thus the average values of the preceding operators can be easily calculated using mean-field trajectories.

II.4 Model System

Following previous workBuzek ; Flick2017 , we investigate a model atomic system in a one-dimensional electromagnetic cavity, as depicted in Fig 1.

[TABLE]

where the upper limit of the first and last summation $m$ denotes the number of atomic energy levels. In the case of a two-level atomic system, this corresponds to a special case of the spin-boson model. With the position of the atom fixed at $r_{A}=\frac{L}{2}$ in this study, half of the $2N$ cavity modes decouple from the atomic system by symmetry. We adopt the same parameters as in Ref. Flick2017 ; Su1991 , which are based on a 1D Hydrogen atom with a soft Coulomb potential (in atomic units): $\{\varepsilon_{1},\varepsilon_{2}\}=\{-0.6738,-0.2798\}$, $\lambda_{\alpha}(\frac{L}{2})=0.0103\cdot(-1)^{\alpha}$, $L=2.362\cdot 10^{5}$ and $\mu_{12}=1.034$. For the three-level atom, we adopt all the same parameters for the field and the atom-field coupling as for the two-level case. The atomic energies for the three level model are $\{\varepsilon_{1},\varepsilon_{2},\varepsilon_{3}\}=\{-0.6738,-0.2798,-0.1547\}$ and as before the numerical parameters are based on the 1D soft-Coulomb Hydrogen atom. The dipole moment operator only couples adjacent states, such that, the only nonzero matrix elements are $\{\mu_{12},\mu_{23}\}=\{1.034,-2.536\}$ and their conjugates.

III Results and Discussion

We now investigate the performance of the MTEF method in the context of cavity-bound spontaneous emission. In all calculations shown below we use 400 photon modes to represent the cavity field. We choose the atom to be initially in the excited state, and the cavity field is in the vacuum state at zero temperature. In all simulations reported here we use an ensemble of $N_{traj}=10^{4}$ independent trajectories, sampled from the Wigner transform of the initial field density operator given in the previous section. This level of sampling is sufficient to converge the atomic observables to graphical accuracy, however observables and correlations functions of the photon-field would require a slightly larger trajectory ensemble for graphical convergence. All observables shown below correspond to their normal ordered forms. For our benchmark numerical treatment we solved the time-dependent Schrödinger equation by using a truncated Configuration Interaction (CI) expansion. More precisely, the photon field state-space is truncated at two excitations per photon mode, whereas for the atomic system a two and three state discrete variable representation is used in each case. Numerical convergence is checked to ensure that the CI basis that we employ is complete for the models and parameter regimes studied in this work.

III.1 Two Level Atom: One Photon Emission Process

In Fig. 3 we show the intensity of the cavity field along the axis of the cavity, at four different times. As the spontaneous emission process proceeds, a photon wavepacket with a sharp front is emitted from the atom and travels toward the boundaries where it is reflected, and then travels back to the atom (e.g. panel (c) of Fig. 3). The emitted photon is then absorbed and re-emitted by the atom, which results in the emergence of interference phenomena in the electric field. This produces a photonic wavepacket with a more complex shape (e.g. panel (d) of Fig. 3). We observe that the MTEF simulations capture the qualitative character of the spontaneous emission process extremely accurately, as well as the wavepacket propagation through the cavity. However, MTEF dynamics fails to reproduce the interference phenomena in the field due to re-emission. We do note, however, that the MTEF simulations are capable of describing the remaining field intensity at the atomic position (e.g panel (e) of Fig. 3). This feature corresponds to a bound electron-photon state, or polariton, which is an emergent hybrid state of the correlated light-matter system.

We also plot the excited state population of the atomic system, and the average value of the photon number for the field, in Fig. 4. Again, MTEF is able to capture the qualitative behaviour of both of these quantities very nicely. However, it fails to quantitatively reproduce the correct values for the emitted photon number and atomic population transfer, as these quantities are underestimated. Furthermore, as a result of this loss in accuracy, only a part of the subsequent re-excitation and re-emission processes are captured.

In Fig. 5 we investigate the normalized second order correlation function, $g^{2}(r_{1},r_{2},t)$ for the cavity photon field. The unperturbed vacuum state, which is coherent, corresponds to $g^{2}(r_{1},r_{2},t)=1$, given by the black background seen in Fig. 5. The vacuum state is disturbed by the emitted wavepacket, corresponding to anti-bunched light with $g^{2}(r_{1},r_{2},t)<1$. The simplicity of the one-dimensional/one-photon process is quite clear in Fig. 6, where we show the associated one-dimensional cuts of $g^{2}$, along with projections of $g^{2}(r_{1},r_{2},t)$ along the positive and negative diagonals, $r_{\pm}=(r_{1}\pm r_{2})/\sqrt{2}$. Here we find similar to the intensity a nice qualitative agreement between MTEF and the exact result for the first three time snapshots. However for the last time-snap-shot the exact solution shows a broader correlation than MTEF, which corresponds to the fact that MTEF is not able to accurately capture re-emission. Furthermore, as we only consider a one-photon process in this case, the correlation is symmetric in $r_{+}$ and $r_{-}$.

III.2 Three Level Atom: Two Photon Emission Process

We now investigate the three-level system, for the same observables as the previous section. The initial state for the atomic system is now the second excited state. The photonic initial state remains the zero temperature vacuum state.

In Fig. 7 we show the intensity of the cavity field during the two-photon emission process. Similar dynamics are observed compared with the two-level case. However, due to the additional intermediate atomic state, we now observe a double-peak feature in the emitted photonic wavepacket. This feature corresponds to the emission of two photons, as the excited atom initially drops to the first excited state emitting one photon, and then further relaxes to the ground state, emitting a second photon. The polariton peak (the central feature in the field intensity profile) is overestimated in the MTEF simulations. This overestimation is due to the incomplete relaxation of the second excited state within the Ehrenfest description.

In Fig. 8, we show the time evolution of the atomic state populations and total photon number. Again, the emitted photonic wavepacket moves through the cavity, is reflected at the mirrors, and returns to the atom. The first and second excited state are then repopulated due to stimulated absorption. A second spontaneous emission process ensues, and the emitted field again takes on a more complex profile due to interference. For the intensity, as well as the atomic population and photon number, we observe that MTEF displays qualitatively correct short-time dynamics. However it fails to describe the correct spatial structure of the (re)emitted two photon wavepacket, as well as the correct amplitude for the observables, in accordance with what was observed previously in the two-level case.

In Fig. 9 we show $g^{2}(r_{1},r_{2},t)$ for the two-photon emission process. The energy level spacing in the three-level truncation of the 1D soft-Coulomb hydrogen atom is uneven, such that the two emitted photons are of different frequencies. Hence, in contrast to the one-photon process, we expect to observe asymmetric features in the second order correlation function. In the exact result we observe that the vacuum state is locally disturbed by a structured, anti-bunched photon wavepacket. The fine, multi-lobed, spatial structure of the photon wavepacket is blurred into a single, rather narrow, feature in the MTEF result. However, MTEF dynamics indeed show the correct spatial asymmetry that is expected in $g^{2}(r_{1},r_{2},t)$. In the corresponding one-dimensional cuts of $g^{2}(r_{1},r_{2},t)$, shown in Fig. 10, we show in further detail the comparison of MTEF dynamics and the exact results in this more complex two-photon case.

IV Summary and Outlook

In this work we have adapted the multi-trajectory Ehrenfest method (MTEF) to simulate correlated quantum mechanical light-matter systems. We applied this mixed quantum-classical dynamics method, which is traditonally applied to electron-nuclear dynamics problems, to two and three level model QED cavity bound atomic systems, and in order to simulate observables and correlation functions for both the atomic system, and the photon field. We find that MTEF dynamics is able to qualitatively characterize the correct dynamics for one and two photon spontaneous emission processes in a QED cavity. However, MTEF dynamics does suffer from some quantitative drawbacks. Furthermore, we also observed that MTEF dynamics simulations can, in fact, capture some quantum mechanical features such as bound polariton states and second order photon correlations. Moreover, as experimental advances drive the need for realistic descriptions of light-matter coupled systems, trajectory-based quantum-classical algorithms emerge as promising route towards treating more complex and realistic systems. In particular, combining the ab initio light-matter coupling methodology recently presented in Jestädt et. al. JRORA18 with our multi-trajectory approach, provides a computationally feasible way to simulate photon-field fluctuations and correlations in realistic three-dimensional systems. Work along these lines is already in progress. Furthermore, an alternative to the independent trajectory-based approach employed here is the conditional wavefunction approach, which allows one address nonadiabatic dynamics problems in complex systems with higher accuracy than MTEF dynamicsGuille18 , and opens up an interesting potential route for mixed quantum-classical methods in correlated light-matter systems.

Acknowledgements

We would like to thank Niko Säkkinen and Johannes Flick for insightful discussions and acknowledge financial support from the European Research Council (ERC-2015-AdG-694097). AK acknowledges support from the National Sciences and Engineering Research Council (NSERC) of Canada.