Light-matter interactions via the exact factorization approach

The exact factorization approach, originally developed for electron-nuclear dynamics, is extended to light-matter interactions within the dipole approximation. This allows for a Schrödinger equation for the photonic wavefunction, in which the potential contains exactly the effects on the photon field of its coupling to matter. We illustrate the formalism and potential for a two-level system representing the matter, coupled to an infinite number of photon modes in the Wigner-Weisskopf approximation, as well as to a single mode with various coupling strengths. Significant differences are found with the potential used in conventional approaches, especially for strong couplings. We discuss how our exact factorization approach for light-matter interactions can be used as a guideline to develop semiclassical trajectory methods for efficient simulations of light-matter dynamics.


Introduction
The interaction of light with matter involves the correlated dynamics of photons, electrons, and nuclei. Even at a non-relativistic level the solution of Schrödinger's equation for the coupled subsystems is a daunting computation. In a given situation however, one is often measuring properties of only one of these subsystems. For example, one might be wanting to know how the electrical conductivity of a molecule is affected by the photons, as in the recent experiment showing the increased conductivity of organic semiconductors due to hybridization with the vacuum field [1]. On the other hand, one might want to understand how molecular dissociation after electronic excitation is affected in the presence of light, as in the recent study of light-induced versus intrinsic nonadiabatic dynamics in diatomics [2]. Or, one might want to measure the superradiance from a collection of atoms [3]. In each of these three cases, the observable of interest involves one of the subsystems alone, electronic, nuclear, and photonic, respectively, yet to capture the dynamics of the relevant subsystem, clearly the effects of all subsystems are needed. The question then arises: can we write a Schrödinger equation for one of the subsystems alone, such that the solution yields the wavefunction of that subsystem? The potential appearing in the equation would have to incorporate the couplings to the other subsystems as well as to any externally applied fields. Hardy Gross, with co-workers, in fact already answered exactly these questions [4][5][6] for the case of coupled electronic and nuclear subsystems in the presence of a classical light field neglecting the magnetic field contribution. That is, for systems of electrons and nuclei, interacting with each other via a scalar potential (usually taken as Coulomb), and in the presence of an externally applied scalar potential, such as the electric field of light, it was shown that one can exactly factorize the complete molecular wavefunction into a wavefunction describing the nuclear system, and a wavefunction describing the electronic system that is conditionally dependent on the nuclear subsystem [4][5][6][7]: Ψ (r, R, t) = χ(r, t)Φ R (r, t), where r = r 1 , . . . , r Ne and R = R 1 , . . . , R Nn represent all electronic and nuclear coordinates respectively. The equation for the nuclear subsystem has a Schrödinger form, with scalar and vector potentials that completely account for the coupling to the electronic system. One can reverse the roles of the electronic and nuclear subsystems, to instead get a Schrödinger equation for the electronic system, which is particularly useful when one is most interested in the electronic properties [8], e.g. in field-induced molecular ionization.
Recently rapid experimental and theoretical advances have however drawn attention to fascinating phenomena that depend on the quantization of the light field in its interaction with matter. This includes few-photon coherent nonlinear optics with single molecules [9], direct experimental sampling of electric-field vacuum fluctuations [10,11], multiple Rabi splittings under ultrastrong vibrational coupling [12], exciton-polariton condensates [13,14], polaritonically enhanced superconductivity in cavities [15], or frustrated polaritons [16] among others. Optical cavities can be used to tune the effective strength of the light-matter interaction, and, in the strong-coupling regime in particular, one finds for example non-radiative energy transfer well beyond the Förster limit between spatially separated donors and acceptors [17], strong coupling between chlorosomes of photosynthetic bacteria and confined optical cavity modes [18], photochemical reactions can be suppressed with cavity modes [19], the position of conical intersections can be shifted or they can be removed [2,20], or state-selective chemistry at room temperature can be achieved by strong vacuum-matter coupling [21]. Strong vacuum-coupling can change chemical reactions, such as photoisomerization or a prototypical deprotection reaction of alkynylsilane [21,22]. This has given rise to the burgeoning field now sometimes called "polaritonic chemistry" [20,[23][24][25][26][27]. In addition, novel spectroscopies have been proposed which explicitly exploit correlated states of the photon field. For example the use of entangled photon pairs enables one to go beyond the classical Fourier limit [28,29], or correlated photons can be used to imprint correlation onto matter [20,27,30,31].
In this paper, we extend the exact factorization approach to non-relativistic coupled photon-matter systems within the dipole approximation. We focus particularly on finding the potential driving the photonic system in the present study. One motivation is towards developing mixed quantum-classical methods for the light-matter system. The observation that in a matter-free system, the photonic Hamiltonian is a sum over harmonic Hamiltonians for each mode of the radiation field suggests that a classical treatment of the photonic system would be accurate: if the system begins in a Gaussian wavepacket, classical Wigner dynamics exactly describes the motion [32]. Coupling to matter within the dipole approximation where the coupling operator is linear in the photonic variable preserves the quadratic nature of the Hamiltonian, and one might then think that again a classical Wigner treatment would be exact. However, although accurate, it is not exact. This implies that the true potential driving the photonic motion is in fact not quadratic. The exact factorization approach defines exactly what this potential should be. In this paper we explain the formalism and give some examples of this potential, that clearly show deviations from harmonic behaviour throughout the dynamics.
The theory is described in Section 2, presenting the Hamiltonian that we will consider, and the formalism of the factorization approach. Section 3 demonstrates the approach on two examples, that we choose as the simplest cases for this initial study. The matter system is described by a two-level system while the photonic system is chosen to either be an infinite number of modes treated within the Wigner-Weisskopf approximation, or a single cavity mode chosen to be resonant with the spacing of the two levels, explored over a range of coupling strengths. We find and interpret the potential driving the photonic system, which depends significantly on whether the initial state of the system is chosen correlated or fully factorized. Finally in Section 4 we summarize and discuss the relevance of this approach for future investigations of light-matter dynamics.

QED-Hamiltonian
In this work, we consider the non-relativistic limit of a system of N e electrons, N n nuclei, and N p quantized photon modes, treated within the dipole approximation in Coulomb gauge [27,33,34]. For now, we do not consider any classical external fields, and neglect spin-coupling. The Hamiltonian of this coupled system is then defined by [20,[35][36][37][38] H(q, r, R) =Ĥ p +Ĥ e +Ĥ n +Ĥ ep +Ĥ np +Ĥ en +Ĥ pen , (1) which operates in the space of: r = {r 1 ..r i ..r Ne } representing all electronic spatial coordinates, R = {R 1 ..R I ..R Nn } representing all nuclear coordinates, and q = {q 1 ..q α ..q Np } representing all photonic displacement coordinates. The first term characterizes the cavity-photon Hamiltonian Hereq α = α 2ωα (â + α +â α ) defines the photonic displacement coordinate for the αth mode, with creation(a + ) and annihilation(a) operators [35,36], and the commutation relation [q α ,p α ] = ı δ α,α . The photonic displacement coordinate is directly proportional to the modeprojected electric displacement operator,D α = 0 ω α λ αqα , whilep α is proportional to the magnetic field [36,37]. The αth mode has frequency ω α = k α c = απc/V , with k α the wavevector and V the quantization volume. The electron-photon coupling strength is given by where S α denotes the mode function, e.g. a sine-function for the case of a cubic cavity [36,39], k α the wave vector, and X the total dipole of the system. In particular, we emphasize at this point that the mode functions introduce a dependence of the coupling constants on the quantization volume of the electromagnetic field. By confining this volume, for example with an optical cavity, one can tune the interaction strength. Finally, we note that the sum in equation (2) goes from 1 to 2N p , to take the two polarization possibilities of the electromagnetic field into account. The second term of equation (1)

denotes the electronic Hamiltonian
where m e defines the electronic mass,p i the electronic momentum operator conjugate tor i . The third term in equation (1) denotes the nuclear Hamiltonian with analogous identifications to the electronic Hamiltonian and eZ I here being the nuclear charge. The remaining terms in equation (1) denote the couplings between the subsystems. The electron-nuclear coupling appears as the usual Coulombic interaction: the electron-photon coupling, in dipole approximation, (where e is the magnitude of the electronic charge) bilinearly couples the total electric dipole moment with the electric field operator for each mode of the photonic field. Similarly, the nuclear-photon coupling iŝ Finally, H pen represents the dipole self-energy of the matter in the radiation field: This self-energy term is essential for a mathematically well defined light-matter interaction. Without this term the Hamiltonian is not bound from below, and loses in addition translational invariance (in case of a vanishing external potential) [40].
The dynamics of such a coupled system is given by the solution of the time-dependent Schrödinger equation where Ψ (r, R, q, t) is the full matter-photon wavefunction, that contains the complete information of the coupled system. However it is difficult to obtain an intuitive understanding and interpretation of such a coupled system from the high-dimensional Ψ (r, R, q, t), and moreover, we may not be interested in all the information as we might be interested in one of the subsystems. If one of these subsystems varies on a much slower time-scale than the others (in particular the nuclei), what is often done in coupled electron-nuclear systems is a Born-Oppenheimer (BO) adiabatic approximation where the faster time-scale subsystem (in particular the electrons) are assumed to instantaneously adjust to the positions of the nuclei, and hence if they begin in an eigenstate, they remain in an eigenstate parameterized by the nuclear coordinate. The eigenenergy maps out a BO potential energy surface (PES) which provides the potential for the nuclear dynamics. These potential-energy surfaces are clearly an approximation within the adiabatic ansatz, but in fact an exact PES can be defined quite generally without the need for any adiabatic approximation, which brings us to the main point of this paper. For the electronnuclear problem, these arise from the exact factorization approach mentioned earlier in the introduction. In the next section we will extend the idea of the exact factorization for electron-nuclei systems to coupled photon-matter systems. Before moving to this, we note that equation (1) is the most general form of Hamiltonian that we will consider in the present work. In later sections, in particular in the explicit examples, we will simplify to just a two-level electronic system interacting with the photonic field in a cavity. In that case, many of the terms in equation (1) are zero, and we simplify the remaining terms even further to a model Hamiltonian Here σ i are the Pauli matrices. The first term is the two-level system that replaces the electronic Hamiltonian (including the dipole self-energy, which simplifies to a constant energy shift for a two-level system), where the energy-level difference is ω 0 , and d eg , appearing in the third term, is the dipole moment of the transition. The second term describes the free photon field, as in equation (2), while equation (8) reduces to the third term with λ α as the coupling strength evaluated at the position of the atom in the cavity. The TDSE also simplifies, to where we use the notation − → Ψ (q, t) being a 2-vector defined at every q and t. A cartoon of the problem is given in Figure 1.

Exact factorization approach
The exact factorization (EF) may be viewed as a reformulation of the quantum mechanics of interacting coupled systems where the wavefunction is factored into a marginal amplitude and a conditional amplitude [4][5][6][7]41]. With non-relativistic electron-nuclear systems in mind, the equations for these amplitudes were derived for Hamiltonians of the form whereV is a scalar potential that includes coupling between the electrons and nuclei (usually Coulombic) and any externally applied fields. HereT e,n are kinetic energy operators of the electronic and nuclear systems, just as in equations (4) and (6), that have the form of − i(I) ∇ 2 i(I) /2m i(I) (that is, no vector potential). The EF then proves that the exact full molecular wavefunction can be factored as The equation for the nuclear amplitude χ has a TDSE form [5,6,42,43], equipped with a time-dependent scalar potential (R, t) and a time-dependent vector potential A I (R, t) that include entirely the effects of coupling to the electronic system as well as external fields. The equation for the conditional electronic amplitude Φ R has a more complicated form, involving a coupling operatorÛ en , that acts on the parametric dependence of Φ R . The factorization is unique, up to a gauge-like transformation, provided Φ R satisfies the "partial normalization condition" (PNC), dr|Φ R (r, t)| 2 = 1; under such a transformation, and A transform as scalar and vector potentials do in electrodynamics. The nuclear N n -body probability density and current-density can be obtained in the usual way from the nuclear amplitude χ(R, t), so in this sense, χ can be identified as the nuclear wavefunction of the system. The form of equation (15) is similar to the BO approximation, however with the important difference that equation (15) is an exact representation of the wavefunction, not an approximation, and further that it is valid for time-dependent systems, with time-dependent external fields, as well. The BO approximation assumes that the electronic system remains always in the instantaneous ground (or eigen)-state associated with the nuclear configuration R, and therefore misses all the physics associated with non-adiabatic effects, including wavepacket branching and decoherence. These effects are contained exactly in the coupling terms in the EF equations: the scalar and vector potentials and the coupling operator of the electronic equation. It is important to note that there is no assumption of different timescales in the EF approach, in contrast to the BO approximation.
As the scalar potential plays a role analogous to the BO PES, but now for the exact system, it is denoted the time-dependent potential energy surface (TDPES), while the vector potential (TDVP) is an exact time-dependent Berry connection. The gauge-freedom is a crucial part of the EF approach: in particular, whether a gauge exists in which the vector potential can be transformed into part of the TDPES has been explored in some works [44][45][46][47], especially since the common understanding is that Berry phases appear only out of an adiabatic separation of timescales, while the EF is exact and does not assume any such separation. Further, equally valid is the reverse factorization [8], Ψ (r, R, t) = χ(r, t)Φ r (R, t), which is particularly useful when one is interested in the electronic system, since in this factorization, the electronic system follows a TDSE in which the potentials can be analysed and interpreted.

Exact factorization approach for QED
Here, we extend the exact factorization to systems of coupled photons, electrons, and nuclei. Since all the kinetic operators in the Hamiltonian within the dipole approximation, equation (1), are of similar form to those that were considered in the original EF approach, equation (14), the mathematical structure of the equations and coupling terms will be similar when we make a factorization into two parts.
There are three possibilities for such a factorization, and we expect each to be useful in different contexts. One possibility, which is perhaps the most natural extension of the factorization of references [4][5][6][7]41], is to take the nuclear system as the marginal one, with the PNC for every nuclear configuration R at each time t. This would yield a TDSE for the nuclear system, much like in the original EF approach, but now the TDPES and TDVP includes not only the effects on the nuclei of coupling to the electrons, but also to the photons. This would be a particularly useful factorization for studying light-induced non-adiabatic chemical dynamics phenomena, when the quantum nature of light is expected to play a role. In fact, an approximation based on the normal BO approximation for the electron-ion dynamics has been used to study the cavity-induced changes in the potential energy surfaces in the strong coupling regime [48]. This would be a particularly useful factorization for studying light-induced non-adiabatic phenomena, when the quantum nature of light is expected to play a role. A second possibility is the natural extension of the reverse factorization [8], where the electronic system is the marginal amplitude with the PNC dqdR|Φ r (q, R; t)| 2 = 1, for all t and every electronic configuration r, which would yield a TDSE for electrons, with the e-TDPES and e-TDVP now incorporating the full effects on the electrons of coupling to the nuclei as well as the photons. This could be particularly useful for studying, for example, the impact of the vacuum field on electrical conductivity in a molecule or semiconductor. This leaves the third possibility, where the photonic system is chosen as the marginal: with the PNC for each field-coordinate q and all times t. This is the factorization we will focus on in the present paper: it gives a TDSE for the photonic system, within which the scalar potential, which we call the q-TDPES, and vector potential, the q-TDVP, contain the feedback of the matter-system on the radiation field. In free space, the potential acting on the photons is quadratic as is evident from equation (2), however, in the presence of matter, the potential determining the photonic state deviates from its harmonic form due to interactions with matter. The cavity-BO approach introduced in reference [34] has demonstrated these deviations within the BO approximation. The EF approach now renders this concept exact, beyond any adiabatic assumptions. The equations for each of these three factorizations follow from a straightforward generalization of the original EF equations, as the non-multiplicative operators (the kinetic operators) have the same form; hence the derivation proceeds quite analogously to that given in references [5,6,42]. In particular, for the factorization equation (19), we obtain where the matter HamiltonianĤ m is given bŷ defined in an analogous way to the BO Hamiltonian, but now for the photonic system. The electron-photon coupling potentialÛ ep is given bŷ the q-TDPES by and the q-TDVP by The factorization (19) is unique up to a gauge-like transformation, provided the PNC, equation (20) is satisfied. The gauge-like transformation has the structure of the usual one in electromagnetism, except here the scalar and vector potentials arise due to coupling, rather than due to external fields, and they are potentials on the photonic system, not on the matter system. The equations are form-invariant under the following transformation: Further, one can show that the displacement-field density represented by χ reproduces that of the full wavefunction, i.e. |χ(q; t)| 2 = drdR|Ψ (q, r, R; t)| 2 , and that the phase of χ together with the q-TDVP provide the displacement-field probability current in the natural way: where χ(q; t) = |χ(q; t)| exp(iS(q, t)). This means that observables associated with multiplication by q can be obtained directly from χ(q, t), for example, the electric field while the magnetic field is

Exact factorization for simplified model hamiltonian: two-level system in radiation field
For our exploration of the QED factorization in this paper, we will turn to the simplified model Hamiltonian of equation (12), where the matter system's Hamiltonian is a 2×2 matrix. First, it is useful to write equation (12) HereĤ qBO is analogous to the BO Hamiltonian in the usual electron-nuclear case. We can define q-BO states as normalized eigenstates: and these can be used as a basis to expand the fully coupled wavefunction, i.e.
which would be analogous to the Born-Huang expansion but now for the cavity-matter system. Now in the EF approach, the fully coupled wavefunction is instead factorized as a single product: where the PNC becomes and holds for every q and each time t.
We note that there are two useful bases for this problem. One is obtained from diagonalizing the field-free two-level system, i.e. that defined by eigenvectors of the Pauli-σ z matrix. The other basis is the q-BO basis, defined by the eigenvectors ofĤ qBO , as in equation (35).
The EF equations follow directly from equations (21)- (27) but with the much simplifiedĤ qBO above, and all integrals over r and R replaced by 2×2 matrixmultiplication.

Photonic time-dependent potential energy surface
Unlike the original electron-nuclear factorization, the q-TDVP can always be chosen to be zero due to the onedimensional nature of each photon-displacement mode. This means that one can always transform to a gauge in which the q-TDPES contains the entire effect of the coupling of the matter system on the radiation field, i.e. it is the only potential that is driving the photonic dynamics. For the matter system, both the q-TDPES and the photon-matter coupling operator incorporate the effect of the photonic system on the matter. In the original electron-nuclear factorization, the exact TDPES proved to be a powerful tool to analyze and interpret timeresolved dynamics of the system in cases ranging from dynamics of molecules in strong fields [5,6,8,[49][50][51][52], to non-adiabatic proton-coupled electron-transfer [53][54][55], to nuclear-velocity perturbation theory [56,57] and dynamics through a conical intersection [58,59]. It provides an exact generalization of the adiabatic BO-PES.
In the present work, we will study the q-TDPES (q, t) of equation (26) for the case of the radiation field coupled to a two level-system, using the model Hamiltonian (12). Given a solution − → Ψ (q, t) for the coupled system, found from equation (13), we will extract the exact q-TDPES by inversion.
To do this, we first ensure that we work in the gauge where A α = 0. Similarly to previous work [6,54], this gauge can be fixed by choosing the phase S(q, t) of the photonic wavefunction, χ(q, t) = |χ(q, t)| exp(iS(q, t)), to satisfy So, from the given solution − → Ψ (q, t), we compute and insert into the q-TDPES where we have identified here the first term in equation (41) as wBO , a "weighted q-BO" surface, weighted by the probabilities of being in the q-BO eigenstates: using the expansion equation (36), The second term arises from kinetic effects from the parametric dependence of the conditional matter wavefunction, hence we denote it as kin , and it originates from the electron-photon coupling operator: in this gauge the only term that contributes to the electron-photon coupling operator expectation value is Both wBO and kin are invariant under different gauge choices, while the last term in equation (41) is gaugedependent, hence its name GD .

Results and discussion
We will consider two extremes within the simplified model Hamiltonian equation (12). The first is the Wigner-Weisskopf limit where the two-level system is coupled to an infinite number of cavity modes. This is the classic model for spontaneous emission: an atom initially in an excited state in a vacuum decays to the ground-state by spontaneously emitting a photon. The second system we study is the two-level system coupled to a single resonant mode. The Hamiltonian is then the same as the Jaynes-Cummings one, but as we will begin with a photonic vacuum and excited atom, we will not see the famous collapses and revivals, but we will see Rabi oscillation type behavior for weak coupling. In both cases, our central question is what are the structures and features of the q-TDPES potential that drives the photonic system away from its vacuum state? As initial condition, we take the photon modes in the vacuum state, and the two-level system in the excited state. For the single resonant mode case, we will compare the effect of starting in a fully factorized matter-photonic state with that of starting in a q-BO state. The fully factorized initial state would be the physical one when an excited atom is instantaneously brought into a closed cavity and just then its dynamics is studied, while the excited q-BO state results when there is initially an external dissipative coupling together with an applied resonant field to maintain the atom in that excited state before the dynamics is examined.
We notice that the dipole matrix element and coupling parameter appear only together as a product in this model equation (12), d eg λ. Physically, these are fixed by the problem at hand, specifically the volume of the cavity and the dipole coupling between the two levels in the atom, apart from fundamental constants. But here, in this model we choose them arbitrarily, and compare dynamics for different d eg λ that range from relatively weak coupling to strong coupling.

Wigner-Weisskopf limit
We first consider the Wigner-Weisskopf limit, in which our two-level system is coupled to an infinite number of modes. In this limit, the accepted well-known approximate solution for the coupled system is known analytically, which makes the q-TDPES particularly straightforward to find. We begin by briefly reviewing this solution.
The solution for − → Ψ of the coupled problem can be found in the standard literature [60]. The initial state is taken to be a purely factorized state of the electron in the excited state and all photon modes in their ground states, i.e. where which follows from the harmonic nature of the free photon field. The coupling in the off-diagonal elements of equation (12) then cause Ψ to evolve in time, as under the reasonable assumption that the coefficients of the two-photon and higher states are negligible. Here the one-photon states of the photonic system are The coefficients a(t) and b α (t) can be found by substituting equation (46) into the TDSE equation (13). After making the Wigner-Weisskopf approximations (taking the continuum limit so V → ∞, taking a(t) to change with a rate much slower than the resonant frequency ω 0 and performing a Markov rotating-wave approximation, and neglecting a divergent Lamb shift), we arrive at where g α = πωα 2 λ α d eg and the decay (spontaneous emission rate), Γ = (d eg λ) 2 ω 2 0 V c 3 . The Wigner-Weisskopf solution is accurate for weak coupling, so that in this limit the solution also generates accurate q-TDPES.
With this Wigner-Weisskopf solution, we can then find the corresponding "exact" q-TDPES, equation (41), using equations (40) and (39). However, this yields an infinite dimensional surface, since q = (q 1 ..q α ...q ∞ ), which is challenging to visualize. Instead, we plot some onedimensional cross-sections of the q-TDPES, along the ith mode, setting q α =i = 0. In the following, we useq to denote all modes not equal to q i . We will abbreviate quantities such as (q i ,q = 0, t) by (q i , t), understood to be looking at the cross-section where the displacementcoordinate of all other modes is zero. We will choose two different modes to look along: one resonant with ω 0 , and the other slightly off-resonant. With this choice of cross-sections through the origin of all modes but one, it can be shown that the phase of the nuclear wavefunction that satisfies the zero-q-TDVP condition, equation (27), S(q i , t) ≡ 0. This leads to some simplification in the components of (q i , t).
Before we discuss the q-TDPES, in Figure 2 we plot the autocorrelation function as this gives an indication of what time-scales to expect in the behavior of the q-TDPES (q i , t) for different coupling strengths d eg λ = {0.01, 0.1, 0.4}. In the upper panel, we have chosen to plot the q-TDPES along the mode of the radiation field that is resonant with the two-level system. In this case, the decay of the autocorrelation depends primarily on (d eg λ) 2 , through Γ , i.e. A Φ (t) ∝ e −Γ t , although there are some small polynomial corrections.
In the slightly off-resonant case, we have chosen ω i = 0.41 while ω 0 = 0.4. In fact, we observe partial revivals in the autocorrelation function for very long times in the case of the weakest coupling shown (d eg λ = 0.01), as shown in the inset, with the amplitude decreasing with each revival. However the initial decay follows a similar d eg λ-scaling pattern to that of the on-resonant section. In either case, the dynamics of the decay is essentially the same for all coupling strengths, provided the time is scaled appropriately, and their q-TDPES's also map on to each other at the corresponding times. In the following then, we will focus on the case d eg λ = 0.01, for both the crosssection taken along the on-resonant mode and the slightly off-resonant mode.
In Figures 3 and 4 we show the different components of the q-TDPES (q i , t) for the different time snapshots indicated by the colored dots in the decay-plot in the top left panel, for the on-resonant and off-resonant sections respectively. The displacement-field density, |χ(q i , t)| 2 = |χ(q i ,q = 0, t)| 2 at these time-snapshots is shown in the top middle panel, and we observe the gradual evolution from the vacuum state towards the state with one photon during the decay. This is also seen in the conditional probability amplitudes shown in the top right panel, which we obtain from These are the coefficients of expansion of Φ qi (t) in the q-BO basis and are related to the coefficients in equation (36) and C (2) (q i , t) begin close to 0 and 1, respectively, as expected, and as the coupling kicks in and the atom decays, one might expect them to evolve to 1 and 0, respectively. This is in fact correct for almost all q i , however non-uniformly in q i . As expected from the nature of the bilinear coupling Hamiltonian equation (12), the conditional electronic amplitude associated with larger photonic displacements q i couple more strongly than those associated with smaller ones, so the conditional amplitude on the upper surface falls away from 1 starting on the outer edges and then moving in. In fact, the conditional amplitude at q = 0 remains forever stubbornly at the upper surface, unaffected by the coupling to the field. This non-uniformity is reflected in the q-TDPES (q i , t), plotted in the middle panel, and leads to a strong deviation from the harmonic form it has in the absence of matter. The potential, driving the photonic motion, loses its harmonic form in the initial time steps as the decay begins, peeling away starting from the outer q i . The potential nearer q i = 0 remains harmonic for the initial stages, but as time goes on, more of the surface peels away from the upper surface, while a peak structure develops near q i = 0 that gets increasingly localized and increasingly sharp as the atom decay process completes and the photon is fully emitted. It is this peak structure in the potential driving the photonic system that excites the system from the zero-photon state towards the one-photon state.
We turn now to the components of this exact surface. In the weighted BO surface, wBO (q i , t) that is plotted in the middle right panel, we see the same peeling away from the outer edges, but sticking resolutely to the original upper surface at q i = 0. As the decay occurs, wBO (q i , t) gradually melts to the lower surface everywhere except for a shrinking region near the origin that sticks to the upper surface. The peak seen in the full q-TDPES on the other hand comes from kin (q i , t), plotted in the lower left panel, which gets sharper and sharper as the photon is emitted. Mathematically, this structure follows from the change in the conditional-dependence of Φ q near q i = 0, as the electronic state associated with q i = 0 remains on the upper q-BO surface while away from q = 0, in a shrinking region, the electronic state is associated with the lower surface. This gets sharper as χ(q i = 0, t) gets smaller and smaller there. One can show from the analytic solution, that, in the long-time limit, the surface at q i = 0 grows exponentially with t at a rate determined by Γ , while for q = 0, kin (q i = 0, t → ∞) → 0.
These features of wBO and kin are very similar for both the cross-section that cuts along the resonant mode (Fig. 3) and the section that cuts along the slightly off-resonant (Fig. 4). The remaining component of the q-TDPES, GD is much smaller than the other components, and has a different structure in the two cases. In fact, it is straightforward to show from the analytic solution that GD (q i = 0, t) is independent of t, and that uniformly shifting GD (q i , t) so that GD  there is a symmetric step-like feature in GD , of the size of the difference in the mode frequency of interest and the resonant mode, and as t gets larger, this feature sharpens. Thus, we can see that in the Wigner-Weisskopf limit, the potential driving the photonic modes deviates significantly from its initial harmonic form during the decay, although once again becoming harmonic almost everywhere (except at q = 0) in the long-time limit. The atomphoton correlation is required to capture these effects, and if one wanted to model this exact q-TDPES, the conditional dependence of the electronic amplitude is crucial to include.

Two-level system coupled to a single resonant cavity photon mode
We now turn to the other limit, tuning the cavity so that there is just one mode that couples appreciably to the two-level atom, with a mode frequency that is resonant with the atomic energy difference. The q-BO surfaces can be easily found by diagonalizing H qBO of equation (34), keeping only one mode with ω α = ω 0 in the field: For couplings λd eg 1/2, the q-BO surfaces are approximately parallel and harmonic except at large q (see also Ref. [34]). So in this case if the initial photonic state is a vacuum, then the ensuing dynamics is driven by a largely harmonic potential, without much perturbation from the atom, except at larger q. Deviations from parallel harmonic surfaces, and hence non-q-BO behavior, occurs at larger q and as the coupling increases. We will investigate the q-TDPES driving the photonic dynamics for three different coupling strengths, (d eg λ = {0.01, 0.1, 0.4}) and will include a plot of the two q-BO surfaces with our results for comparison with the exact q-TDPES.
In Figures 5-7, we plot the exact q-TDPES for coupling strengths d eg λ = 0.01, 0.1 and 0.4, respectively, beginning with the atom in the excited q-BO level, multiplied by the photonic ground-state. On the upper panel (left) we plot the autocorrelation function to indicate the approximate periodicity of the system dynamics. Comparing A Ψ (t) in these three figures, we find a decrease of the approximate period with the increase of coupling strength until the periodicity breaks down for the strong coupling d eg λ = 0.4. The weakest coupling strength we have chosen is on the borderline of being in the Rabi regime [61], while the strongest is far from it. The photonic distribution (middle) and conditional electronic coefficients (right) are shown in the top panel of Figures 5-7. The initial coefficients are C (1) (q, 0) = 1 and C (2) (q, 0) = 0. After some time we see a transition of the electron from the excited state to the ground-state as indicated by these coefficients. We observe that the transfer begins earlier for higher values of q and then is followed by lower q-values, and again the conditional amplitude at q = 0 sticks to the upper surface at all times in all cases as there is no coupling for q = 0. The q-dependence of these coefficients has a significant role in shaping the structure of the q-TDPES that we will shortly discuss. At the same time, the probability of photon emission increases, as indicated by the morphing of the initial gaussian in χ(q, t) towards its first-excited profile. For the weakest coupling strength (Fig. 5), after a half period the system begins to move back approximately to its initial state, as the photon is reabsorbed and atom becomes excited again. For strong coupling d eg λ = 0.4 (Fig. 7), the periodic character is lost and we find more wells and structure appearing in the displacement-field density profile. With such strong coupling the q-BO surfaces are quite distorted from a pure harmonic, as evident in the plot (blue lines in the middle panel), and the anharmonicity brings more frequencies into play. A one-photon state that is associated with the lower q-BO surface has a wider profile with density maxima further out than a one-photon state associated with the upper surface would have, for example. In fact the character of the coupled cavity-matter system becomes quite mixed, as is evident from the conditional electronic coefficients shown on the right, and as one goes along the photonic coordinate q one associates with different superpositions of the electronic states. This leads to interesting structure in the exact q-TDPES, that, when decomposed in terms of the q-BO surfaces, has components that vary a lot with q (i.e. not just approximately piecewise-in-q).
The q-TDPES for initial state prepared in the upper q-BO state begins with the weighted q-BO component, wBO (middle right panel) on top of the upper q-BO surface as expected. For the weakest coupling, d eg λ = 0.01, wBO (q, t) then melts down to the lower surface over half a cycle, peeling away from the outer higher q-values first, in a similar way to what was seen in the Wigner-Weisskopf limit. This potential approaches the lower surface before returning back to the upper q-BO surface, but the region near q = 0 remains bound to the upper surface. The time-dependent double-well structure in the potential is again important in driving the photon emission. A similar trend is seen for the stronger coupling 0.1 in Figure 6, but for the strongest coupling d eg λ = 0.4 of Figure 7, wBO (q, t) shows a more complicated correlation in q, with structures mirroring those in the displacement-field density discussed above. As for the kinetic component, for the weaker couplings, a peak structure in kin (q, t) (lower left panel) develops that grows and narrows during the photon emission stage, similar to what was seen in Wigner-Weisskopf, but this then reverses during the reabsorption here. Again for the stronger coupling, the structure is more complicated, mirroring the more complicated dynamics. The gauge-dependent part, GD (lower right panel) is generally a smaller contribution to the total q-TDPES compared to the other components, but again we see step-like features for the weaker couplings, and more complicated dynamics for the strongest coupling.
The dynamics depends significantly on whether the initial state is the correlated q-BO state of Figures 5-7, or a fully factorized one, and now we turn to the surfaces, conditional probabilities, and displacement-field densities for the latter case, plotted in Figures 8-10. The initial coefficients were C (1) (q, 0) = 1 and C (2) (q, 0) = 0 when beginning in the q-BO states, but when beginning in the fully factorized state, these coefficients deviate from these uniform values, especially for larger q, with deviation increasing with the coupling strength. Although the photonic field still begins in the vacuum state, the electronic state is not purely in the upper q-BO surface; the electronic state associated with larger q already has some component in the ground-state. So at these larger values of q, the initial wBO (q, 0) surface dominates the q-TDPES and is anharmonic from the very start, lying intermediate between the upper and lower q-BO surface. In the weak show that the q-TDPES has a tamer structure for the fully-factorized initial state than for the correlated q-BO initial state, especially at larger q; this is likely because less energy is available at these larger q for the system to exchange between the atomic and photonic systems because the atomic state correlated with large q is not completely in its excited state initially.
To summarize: at time zero the exact q-TDPES starts on the upper q-BO-surface, which, depending on coupling strength and choice of initial state, ranges from lying directly on top of the upper q-BO surface (weaker coupling and with q-BO initial state), to in between the two q-BO surfaces with deviations from the upper being larger for larger q (stronger coupling, or fully factorized initial state). After some time the potential starts to melt down onto the lower BO-surface, first starting at higher q-values and then followed by lower q-values, with peak structures developing in the interior region. Around q = 0 the kinetic-component dominates, which leads to an increasing and after half a period decreasing peak. For stronger coupling we observe several peak features in the potential and significant deviations from the curvature of the q-BO surfaces throughout q small contribution below the lower BO-surface; the deviations at larger q arise from the gauge-dependent component.

Summary and outlook
We have introduced an extension of the exactfactorization approach, originally derived for coupled electron-nuclear systems, to light-matter systems in the non-relativistic limit within the dipole approximation. We have presented different possible choices for the factorization but in this work have focussed on the one where the marginal is chosen as the photonic system and the matter system is then conditionally-dependent on this. This choice is particularly relevant when one is primarily interested in the state of the radiation field since the exact factorization yields a time-dependent Schrödinger equation for the marginal, while the conditional is described by an equation with an unusual matter-photon coupling operator. The equation for the marginal is, in a sense, simpler than that in the electron-nuclear case, since the vector potential, q-TDVP, appearing in the equation can always be chosen to be zero, so only a scalar potential remains, the q-TDPES. We have studied the potential appearing in this equation in a gauge where the q-TDVP is zero, for a two-level system coupled to an infinite number of modes in the Wigner-Weisskopf approximation, and for a two-level system coupled to a single photonic field mode with a range of coupling strengths. In all cases we find a very interesting structure of the potential that drives the photonic dynamics, and in particular, large deviations from the harmonic form of the free-photon field. These deviations completely incorporate the effect of the matter system on the photonic dynamics. We also studied the effect of beginning in an initially purely factorized light-matter state, compared to a q-BO initial state, finding significant differences for larger coupling strengths in the ensuing dynamics, implying that in modelling these problems a careful consideration of the initial state is needed.
To use the exact factorization for realistic light-matter systems, approximations will be needed, since solving the exact factorization equations is at least as computationally expensive as solving the Schrödinger equation for the fully coupled system. The success of such an approximation depends on how well the q-TDPES is modelled. The components of the exact q-TDPES beyond the weighted BO depend significantly on the q-dependence of the conditional probability amplitude; approximations that neglect this dependence (Ehrenfest-like) will likely lead to errors in the dynamics. It has been shown recently that mixed quantum-classical trajectory methods that are derived from the exact factorization approach can correctly capture decoherence effects [62][63][64]. Since photons are intrinsically non-interacting and therefore even simpler to treat than nuclei, we expect in analogy to the electronnuclear case that semiclassical trajectory methods derived from systematic and controlled approximations to the full exact factorization of the light-matter wavefunction will be able to capture decoherence effects beyond the Ehrenfest limit for light-matter coupling. This will be subject of future investigations.