Cavity Born-Oppenheimer Approximation for Correlated Electron-Nuclear-Photon Systems

In this work, we illustrate the recently introduced concept of the cavity Born-Oppenheimer approximation for correlated electron-nuclear-photon problems in detail. We demonstrate how an expansion in terms of conditional electronic and photon-nuclear wave functions accurately describes eigenstates of strongly correlated light-matter systems. For a GaAs quantum ring model in resonance with a photon mode we highlight how the ground-state electronic potential-energy surface changes the usual harmonic potential of the free photon mode to a dressed mode with a double-well structure. This change is accompanied by a splitting of the electronic ground-state density. For a model where the photon mode is in resonance with a vibrational transition, we observe in the excited-state electronic potential-energy surface a splitting from a single minimum to a double minimum. Furthermore, for a time-dependent setup, we show how the dynamics in correlated light-matter systems can be understood in terms of population transfer between potential energy surfaces. This work at the interface of quantum chemistry and quantum optics paves the way for the full ab-initio description of matter-photon systems.


INTRODUCTION
Recent experimental progress has made it possible to study light-matter interactions in the regime of strong and ultrastrong light-matter coupling. Experiments from exciton-polariton condensates, 2,3 near-field spectroscopy, 4,5 plasmon-mediated single-molecule strong coupling, 6 superconducting qubit circuits, 7 quantum information, 8 direct measurements of vacuum fluctuations, 9 and chemistry in optical cavities 10−12 open now the path to shape the emerging correlated lightmatter interactions with the goal toward a new control of material properties. In this new field that has been driven in particular by experiment, traditional theoretical methods from either quantum chemistry or quantum optics can lose their applicability. On the one hand, traditional quantum chemistry concepts such as the Born−Oppenheimer (BO) approximation 13,14 or electronic structure methods such as Hartree−Fock theory, 15 coupled-cluster theory, 16 or density-functional theory (DFT) 17 have been originally designed to treat approximately correlated electron−nuclear problems but are not capable to correctly account for the quantum nature of light. On the other hand, concepts from quantum optics typically describe the quantum nature of the light field in great detail but fail in describing more complex dynamics of matter due to the often employed simplification to a few levels. 18,19 To fill this gap, in this work, we generalize a well-established concept from quantum chemistry, namely the Born−Oppenheimer approx-imation, to the realm of correlated light-matter interactions for systems in optical high-Q cavities.
First theoretical studies in similar direction, e.g. the modification of the molecular structure under strong lightmatter coupling, 20 the nonadiabatic dynamics of molecules in optical cavities, 21,22 or the cavity-controlled chemistry, 23 have already been conducted.
Since the complexity of an exact ab initio description of such correlated many-body systems that contain electronic, nuclear, and photonic (Fermionic and bosonic) degrees of freedom scales exponentially with system size, approximate descriptions have to be employed for any realistic system. Recently, the concept of DFT has been generalized to electron-photon problems and was termed quantum-electrodynamical densityfunctional theory. 24−27 This theory maps the complicated many-body problem into a set of nonlinear equations for the electronic and photonic degrees of the densities/currents that facilitates the treatment of such complex systems, similarly as standard DFT has done over the years to deal with correlated electronic systems. Still for this theory to be applicable, accurate functionals for combined light-matter systems have to be developed to calculate approximate effective potentials and observables. In this work, we use an alternative approach, the cavity Born−Oppenheimer (CBO) 1 approximation, that allows to construct approximate wave functions to the exact eigenstates for such problems. The cavity Born−Oppenheimer approximation has recently been introduced in ref 1, and in this paper we derive the theory in a complete manner and give explicit examples to highlight its applicability for general electron−nuclear-photon systems. This work is structured into three sections: (i) First, the theoretical framework is introduced where we demonstrate how the concept of the Born− Oppenheimer approximation can be generalized to matterphoton coupled systems. (ii) We apply this theoretical framework to study a prototypical electron-photon system, where the photon couples resonantly to an electronic transition. (iii) The last section is devoted to a model system of an electron, a nuclei, and photons, where a photon mode couples to a vibrational excitation. In what follows and without loss of generality, we describe the electron−nuclear-photon problem in Coulomb gauge, dipole approximation, and the Power-Zienau-Woolley frame. 28,29 Our system of interest contains n e electrons, n n nuclei, and n p quantized photon modes, e.g. the matter is located in an optical high-Q cavity. Strong light-matter coupling is obtained, once the light-matter coupling is stronger than the dissipation of the system due to e.g. cavity losses. For simplicity, we neglect dissipative channels in the following. [Since in this work dissipation is neglected, we find modifications of the eigenstates of the matter-photon system with respect to the bare matter eigenstates for all nonvanishing matter-photon coupling strengths.] The original derivation of the Born− Oppenheimer approximation is outlined e.g. in ref 14 where Tn and Ŵn are the nuclear kinetic energy and nuclear interaction, respectively. The electron−nuclear interaction Hamiltonian Ĥe n is given by  (4) and the cavity photon Hamiltonian Ĥp with n p quantized photon modes of frequency ω α takes the form The displacement field operators ω =̂+̂ℏ )/ 2 / consist of the usual photon creation and annihilation operators and [q̂α, p̂α ′ ] = iℏδ α,α′ . Furthermore, the q̂α are directly proportional to the electric displacement field operator of the α-th photon mode 30,31 at the charge-center of the system by the connection D̂α = ϵ 0 ω α λ α qα and the pα are proportional to the magnetic field. In eq 5, the sum runs from 1 to 2n p , to correctly account for the two possible polarization directions of the electromagnetic field. The last three terms in eq 1 describe the light-matter interaction Hamiltonian. The first term is the explicit electron-photon interaction in the dipole approximation with the total electronic dipole moment X e = −∑ i=1 n e er i and the matter-photon coupling strength λ α . 25, 31 The second term gives the explicit nuclear-photon interaction, again in the dipole approximation with the total nuclear dipole moment X n = ∑ i=1 n n Z i eR i , and the last term describes the quadratic dipole-self-interaction term where X now describes the total dipole moment of the system, i.e. X = X e + X n . We then introduce the following abbreviations  (9) In general, we are interested in calculating eigenstates Ψ i (r, R, q) and eigenvalues E i of the particular problem. These states then give us access to any observable of interest. To calculate these quantities, we have to solve the full Schrodinger equation of the correlated electron−nuclear-photon problem that is given bŷ i ii (10) where the Hamiltonian Ĥis given by eq 1. Obtaining general solutions to the Schrodinger equation of eq 10 is an ungrateful task.
[We note that in free space eq 10 has no square-integrable eigenstates in the charge neutral case due to its translational invariance. Hence one either has to go into a comoving frame, e.g., a center-of-mass frame, and consider the corresponding reduced Hamiltonian, or one has to use a confining potential to

Journal of Chemical Theory and Computation
Article localize the molecule.] In practice, the Schrodinger equation is barely solved exactly but only approximately. One of such approximate methods is the cavity Born−Oppenheimer approximation 1 that is capable to partially decouple the electronic degrees of freedom from the nuclear and photonic degrees of freedom. In electron−nuclear problems, such an adiabatic decoupling procedure is commonly assumed 14 and well justified for low lying states, e.g. the ground state. However, severe limitations are known that require going beyond the adiabatic treatment by including nonadiabatic electron−nuclear terms, e.g. at conical intersections. 32 In this work we decouple the electronic degrees of freedom from the nuclear and photon degrees of freedom. This allows us, on the one hand, to simplify the problem much more than if we decoupled the nuclear from the electronic and photonic degrees of freedom, as has been done in refs 20 and 33, and the additional photonic degree of freedom becomes formally equivalent to a nuclear (phononic) degree of freedom. The latter can be understood as follows: in eq 5, the photon degree of freedom is written in terms of a quantum harmonic oscillator that contains a kinetic energy term Tp and a potential term Ŵp that are both connected via the virial theorem. 34 In this sense, we can regard the description of the photon modes as formally equivalent to the description of the nuclei. As a consequence, we can apply traditional methods to solve the electron−nuclear problem to the generalized electron−nuclear-photon problem. One of these methods is the Born−Oppenheimer approach that in the cavity accumulates an additional photonic degree of freedom, reminiscent of the nuclear degree of freedom. Thus, the same arguments for the validity of the usual Born− Oppenheimer approximation that apply in the case of nuclear motion also apply for any extended system, as they do not depend on the details of the interactions that produce the potential-energy surfaces. In practice, the main problem for the standard Born−Oppenheimer approximation is to solve the resulting electronic equation, while simple approximations to the nuclear equation, such as harmonic approximations, are often sufficient. On the other hand, a decoupling of the electronic degrees of freedom provides most flexibility for the applications that we consider, e.g. a single electron coupled to one mode. From a physical perspective, however, this decoupling scheme seems counterintuitive at a first glance. The usual simplified argument for the decoupling of the nuclear from the electronic degrees of freedom is that the nuclei move "slowly" compared to the electrons, i.e., the kinetic-energy contribution is negligible, and hence a classical approximation seems reasonable. In the case of quantized photons, the term Tp in eq 5 is related to the square of the magnetic field operator, thus pα is proportional to the magnetic field. Therefore, the magnetic field can be interpreted as an analogue to the nuclear velocity in real-space, although the conjugate momentum = ∂ α ℏ α p i q is defined in the q α -space of the harmonic oscillator. The coordinate q α describes the displacement of the harmonic oscillator of the photon mode with specific energy ω α . In this sense, while the usual Born−Oppenheimer approximation is justified by "slow" nuclei, we can justify the cavity Born− Oppenheimer approximation, if the magnetic field in the photon mode is "small". This is in particular the case for all eigenstates, due to ∂ t q α = p α . Along these lines, we conclude that the cavity Born−Oppenheimer approximation is applicable, if p α remains small, thus the magnetic field remains "small". If this is the case, the time-derivative of q α remains "small", thus the electric displacement field changes only "slowly" over time, and the electrons can adapt "quasiinstantaneously" to these "slow" changes of the electric displacement field. That this approach can indeed give highly accurate results will be demonstrated in the following.
2.2. Cavity Born−Oppenheimer Approximation. In this section, we derive the approximate cavity Born−Oppenheimer states to eq 10. This goal is achieved in three successive steps. First, we solve the electronic part of the eq 10, where we consider explicitly all terms containing an explicit electronic contribution. This electronic Schrodinger equation has only a parametric (conditional) dependence on the nuclear and field degrees of freedom, or alternatively nuclear and field coordinates enter the electronic equation as c-numbers. In principle, the electronic Schrodinger equation has to be solved for every possible combined nuclear and photon-field configuration, and the eigenvalues of the electronic Schrodinger equation then enter the nuclear and photon-field Schrodinger equation through the emerging potential-energy surfaces. Having solved both equations, we can then construct the approximate cavity Born−Oppenheimer states in a factorized manner. To obtain the approximate cavity Born−Oppenheimer states, as a first step, we solve the electronic Schrodinger equation e e n p e p e n j j j (11) for each fixed set of nuclear coordinates R and photon displacement coordinates q. For each fixed set of (R, q), the electronic eigenfunctions of eq 11 {ψ j (r, R, q)} form a complete basis in the electron many-particle Hilbert space. In the electronic Schrodinger equation of eq 11, (R, q) enter the electronic cavity Born−Oppenheimer Hamiltonian as (classical) parameters, thus the eigenvalues ϵ j also parametrically depend on R, q. For each fixed set of (R, q), we can then expand (also known as the Born-Huang expansion 35 ) the exact many-body wave function Ψ i (r, R, q) that is a solution to the full Schrodinger equation of eq 10 as Here, the exact wave function is decomposed into sums of product states consisting of an electronic wave function ψ j (r, R, q) and a nuclear-photon wave function χ ij (R, q). The latter is obtained by solving the following equation where T̂n(R) and T̂p(q) are given by eqs 3 and 5, respectively.
The eigenvalues E i of eq 13 are the exact correlated eigenvalues of eq 10. The term in the second line of eq 13 describes the nonadiabatic coupling between cavity Born−Oppenheimer potential energy surfaces (PES). The cavity Born−Oppenheimer approximation now neglects the offdiagonal elements in the nonadiabatic coupling terms of eq 13. Then eq 13 can be rewritten in a much simpler form (14) where the newly generalized cavity PES V j (R, q) are given explicitly by

Journal of Chemical Theory and Computation
The first two terms are the nuclear and the photon potentials of eqs 3 and 5, and all anharmonicity in the PES can be attributed to the electron-photon, electron−nuclear, nuclear−nuclear, and nuclear-photon interaction contained in eq 1. Furthermore, the eigenvalues E i of eq 14 are an approximation to the exact correlated eigenvalues and provide by the variational principle an upper bound. With this reformulation, we have the advantage that we can solve the electronic Schrodinger equation of eq 11 and the nuclear-photon Schrodinger eq 14 separately. The ground-state Ψ 0 in the cavity Born− Oppenheimer approximation then becomes and accordingly for the excited states. In Born−Oppenheimer calculations for systems that only contain electrons and nuclei often the harmonic Born−Oppenheimer approximation is carried out 14 that can be realized by expanding V j (R, q) around its minimum value and in this way even simplifies the problem further. In the harmonic approximation, we have to solve eq 11 not for all possible values of (R, q), but only at the minimum of ϵ j (R, q). However, in this work, we do not apply the harmonic approximation to correctly demonstrate the full capacity of the cavity Born−Oppenheimer concept. Before we introduce our examples, let us comment on the expectable accuracy of the cavity Born−Oppenheimer states when decoupling electronic from photonic and nuclear degrees of freedom. Our simplified physical arguments for the decoupling scheme so far have been that the nuclei are "slow" and the magnetic-field contribution small, such that we can neglect the corresponding kinetic terms in the equation for the electronic subsystem. However, the decisive quantities that indicate the quality of this approach are the nonadiabatic coupling elements of eq 13 and the distance between the potential-energy surfaces. If these elements are small and the potential-energy surfaces are far apart, we can expect a good quality of the approximate cavity Born−Oppenheimer states. This argument is similar to standard Born−Oppenheimer treatment that loses its validity at crossing of eigenvalues, i.e. conical intersections.

DISCUSSION AND RESULTS
In the following, we now want to illustrate the concept of the cavity Born−Oppenheimer approximation for two specific setups. We numerically analyze first a model system consisting of a single electron coupled resonantly to a photon mode. In this example, the nuclei can be understood as frozen, leading to an external potential acting on the electronic degrees of freedom. This model will allow us to study the decoupling mechanism introduced for the correlated electron-photon interaction in detail. In the second example, we then analyze a model system that contains electron−nuclear-field degrees of freedom. Here, potential-energy surfaces emerge that have nuclear-photon nature.

Light-Matter Coupling via Electronic Excitation.
In this section, we illustrate the concept of the cavity Born− Oppenheimer approximation for a simple coupled electronphoton model system. The system of interest is a model system for a GaAs quantum ring 36 that is located in an optical cavity and thus coupled to a single photon mode. 31 The model features a single electron confined in two-dimensions in realspace (r = r x e x + r y e y ) interacting with the single photon mode with frequency ℏω α = 1.41 meV and polarization direction e α = (1,1). The polarization direction enters via the electron-photon coupling strength, i.e. λ α = λ α e α and depends on the specific experimental setup. The photon mode frequency is chosen to be in resonance with the first electronic transition. We depict the model schematically in Figure 1 (a). The bare electron ground-state n λ=0 (r) has a ringlike structure shown in Figure 1 (b) due to the Mexican-hat-like external potential that is given by with parameters ℏω 0 = 10 meV, V 0 = 200 meV, d = 10 nm, and m 0 = 0.067m e 36 and shown in Figure 1 (c). For the single electron, we employ a two-dimensional grid of N = 127 grid points in each direction with Δx = 0.7052 nm. In contrast, we include the photons for the exact calculation in the photon number eigenbasis, where we include up to 41 photons in the photon mode.
For the cavity Born−Oppenheimer calculations, we calculate the photons also on an uniform real-space grid (qrepresentation) with N = 41 with Δq = 6.77 aJ fs 2 and construct the projector from the uniform real-space grid to the photon number states basis explicitly. This projector can be calculated by employing the eigenstates of the quantum harmonic oscillator in real-space. For a more detailed discussion of the model system, we refer the reader to refs 31 and 36. Since this model can be solved by exact diagonalization in full Fock space, 37 all exact results shown in the following have been calculated employing the full correlated electron-photon Hamiltonian. 1,25,30,31 For this model, the potential-energy surfaces from eq 15 can be calculated explicitly as

Article
In Figure 2 (a), we show the PES surfaces V j ({q α }) for the weak-coupling regime of λ α = 0.0034 meV 1/2 /nm. We find that all PES have a strong harmonic nature, due to the dominant q̂α 2 term in eq 18. The eigenvalues ϵ j and the integral in the last line of eq 18 are the corrections to the harmonic potential. In this case, both are rather small for all excited-state surfaces in the weak-coupling regime, i.e. for the ground-state surface adiabatic term in the last line of eq 18 is around 2 orders of magnitude smaller than ϵ 0 . In general, a harmonic correction that can be obtained by calculating the second derivative at the minimum value will shift the frequency of the photon mode. We define as harmonic approximation to eq 18 where q j,0 is the minimum value of the j-th PES eq 18. In the weak-coupling regime, we find ω̃α ≈ ω j,α . All corrections beyond the second derivative of these terms are then called the anharmonic corrections. We find the lowest cavity PES that is the ground-state PES shown in black, well separated from the first and second excited cavity PES that are shown in solid red and dotted blue. The first and second excited cavity PES are close to being degenerate. This 2-fold degeneracy has its origin in the two-dimensional external potential, similar to the s/p degeneracy in the hydrogen atom. In Figure 2 (b), we show the cavity PES surfaces in the strong-coupling regime with λ α = 0.134 meV 1/2 /nm. [In this work, the weak and strong coupling regimes are defined by the relation of the Rabi splitting to the photon frequency ω α , thus accordingly to the definitions in the Rabi model, see e.g. ref 38 and references therein.] While the second PES shown in blue and the fourth potential energy surface shown in yellow keep the harmonic shape, in the lowest cavity PES shown in black and the third cavity PES shown in solid red, two new minima with a double-well structure appear. [Note that if we would like to express this electron-dressed photon system in terms of the original creation and annihilation operators, we will need new combinations of these operators, i.e., photon-interaction terms. Physically these interaction terms describe the coupling between photons mediated via the electron.] The minima of the cavity PES are strongly shifted away from the equilibrium position at the origin. This electron-dressed potential for the photon modes induces a new vacuum state with two maxima. Since the cavity PES is symmetric, the vacuum state still has a displacement observable of ⟨q α ⟩=0, i.e., we have a stable vacuum with zero field. However, with respect to the bare vacuum the other observables, e.g., the vacuum fluctuations, will clearly change. Furthermore, we find for the harmonic approximation in the ground-state cavity PES, ω̃0 ,α ≈ 0.8 ω α , hence an effective softening of the photon mode in the groundstate cavity PES with the strong displacement of q 0,0 = 18.85 aJ fs 2 . A similar behavior has been observed before in the context of polaron physics in the Holstein Hamiltonian. 39,40 We further analyze this transition in Figure 3. In Figure 3 (a), we show how the ground-state PES depends on the electronphoton coupling strength λ α . We find that for absent and weak coupling, the ground-state surface can be well described by a single harmonic potential that has the minimum at q α = 0. If we increase the electron-photon coupling to strong coupling, we find around λ α = 0.044 meV 1/2 /nm the splitting of the singlewell structure to a double-well structure. For strong coupling, e.g. λ α = 0.1342 meV 1/2 /nm this double-well structure becomes strongly pronounced. In Figure 3 (b) and (c), we plot the corresponding electron density n λ (r) = ∫ dq α Ψ 0,λ * (r,q α ) ×Ψ 0,λ (r,q α ) of the exact correlated ground state Ψ 0,λ (r,q α ) for different values of λ. In the weak-coupling regime, shown in Figure 3 (b), we find that the electron is only slightly distorted in comparison to the ringlike structure of the bare electron ground state 31 shown in Figure 1 (b). In contrast, in the strong coupling regime, shown in Figure 3 (c), the electron density becomes spatially separated and localized in direction of the polarization direction of the quantized photon mode.
The consequences of the ground-state transition identified in Figure 3 become also apparent if we study the difference of the correlated and bare electron density. Let us define the bare electron density. Here, we refer to the electron density that is the ground-state of the external potential without coupling to the photon mode or alternatively λ α = 0, thus n λ=0 (r). This density is shown in Figure 1 (b). Then we define Δn λ (r) = n λ (r) − n λ=0 (r). In Figure 4, we plot Δn λ (r) as a function of the electron-photon coupling strength λ α . In the weak-coupling limit, shown in Figure 4 (a) for λ α = 0.0034 meV 1/2 /nm, we find that the electron density is slightly distorted such that in the correlated density more density is accumulated perpendicular to the polarization direction of the photon mode compared to the bare electron density. However, once the strong-coupling regime is approached, we also identify a transition in Δn λ (r). In the strong coupling regime, that is entered in Figure 4

Journal of Chemical Theory and Computation
Article the ground-state electron density is reoriented until ultimately in Figure 4 (e) the electron density is arranged in direction of the polarization direction of the photon mode, up to higher strong-coupling regions shown in Figure 4 (f).
The additional insights from the ground-state transition can be obtained by evaluating the exact correlated electron-photon eigenvalues. In Figure 5, we plot the exact eigenvalues from the weak-to the strong-coupling regime. The ground-state energies are plotted by the black line and are increasing for stronger coupling. 1 For the first excited state in the case of λ α = 0 coupling, we find a 3-fold degeneracy that is split once the electron-photon coupling is introduced. For strong coupling the first-excited state (shown in blue) and the ground-state become close leading to the splitting of the electron-density shown in Figure 4. [We emphasize that this behavior is similar to what is in molecular systems known as static correlation for e.g. stretched molecules. 41 ] Higher-lying states show energy crossings that are typical for electron-photon problems and have been previously observed e.g. in the Rabi model. 38,42,43 We find not only allowed level crossings at λ α ≈ 0.031, 0.067, 0.113 meV 1/2 /nm but also an avoided level crossing at λ α ≈ 0.055 meV 1/2 /nm between the fifth and sixth eigenvalue surface. In the Rabi model, level crossings are used to define transition from the weak, strong, ultrastrong, 44 and deep-strong coupling regime. 45 Similarly to the Rabi model, 42 we find in the strong coupling regime a pairing of states in terms of the energy. Two states each with different parity become close to degeneracy. Since in the strong-coupling regime the interaction terms in the Hamiltonian become dominant and we apply the interaction in dipole coupling, the eigenstates of the full Hamiltonian become close to the eigenstates of the dipole operator that are the parity eigenstates. We can expect a different behavior beyond the dipole coupling, e.g. if electric quadrupole and magnetic dipole coupling, or higher multipolar coupling terms are also considered. In contrast to the Rabi model, 42 we find an overall increase of the ground-state energy for increasing coupling strength. This behavior is due to the inclusion of the quadratic dipole self-interaction term of eq 8. In Figure 5, we indicate by the dashed line, the ground-state transition discussed before. In the coupling region indicated by (I), we find a single minimum in the PES and Δn is located perpendicular to the polarization direction, while in the coupling regime (II), we find two minima and a double-well structure in the PES and Δn are located along the direction of the polarization of the photon mode. The quality of the cavity Born−Oppenheimer approximation is shown in Table 1 in terms of overlaps ⟨Ψ j |Ψ j,CBO ⟩ 2 between approximate and exact states. If the eigenenergies shown in Figure 5 are well separated as in the strong coupling regime for λ α = 0.1342 meV 1/2 /nm, then the cavity Born−Oppenheimer approximation is well justified. For states that are close to degeneracy, as e.g. the states #2 and #4 in the weak-coupling for λ α = 0.0034 meV 1/2 / nm, we find a lower quality. However, this low quality could be improved by symmetry considerations. Overall, we find a very high and sufficient quality of the approximate energies and states in comparison to its corresponding exact values.
The remaining part of this section is concerned with the time-dependent case. Here, we employ the full correlated Figure 4. Difference of the correlated ground-state electron density to the bare electron density (Δn λ = n λ − n λ=0 ) from the weak-to the strong-coupling limit. The red color refers to surplus density regions, while the blue color refers to regions with reduced density.

Journal of Chemical Theory and Computation
Article electron-photon Hamiltonian and choose as initial state a factorized initial state that consists of the bare electronic ground state and a bare photon field in a coherent state with ⟨â †a⟩ = 4 where λ α = 0.0034 meV 1/2 /nm. This example is also the first time-dependent example studied in ref 31. To numerically propagate the system, we use a Lanczos scheme and propagate the initial state in 160000 time steps with Δt = 0.146 fs. In Figure 6, we briefly analyze this setup by evaluating the dipole moment ⟨x̂+ y⟩ in Figure 6 (a), the purity γ = Tr(ρ ph 2 ) that contains the reduced photon density matrix ρ ph and the Mandel Q parameter 46 that is defined as  (20) in Figure 6 (b) and the photon occupation ⟨â †a⟩ in Figure 6 (c). In the case of the dipole moment of this example shown in Figure 6 (a), we find first regular Rabi oscillations up to the maximum at t = 5 ps and around t = 10 ps, and we find the necklike feature 47 typical for Rabi oscillations. In Figure 6 (b), we show the purity γ in dashed black lines. The purity γ is a measure for the separability of the many-body wave function into a product of an electronic and a photon wave function. We find that γ is close to 1 up to t = 5 ps, which means that the many-body wave function is close to a factorizable state. After t = 5 ps, γ deviates strongly from 1, and the system is not factorizable anymore. This dynamical buildup of correlation has also an effect on the nonclassicality of the light-field visible in the Mandel Q-parameter shown in Figure 6 (b) in solid black lines. While initially Q ≈ 0 that indicates the coherent statistics of the photon mode, after t = 5 ps also this observable deviates from 0 and nonclassicality shows up. From Figure 6 (c), where we plot the photon number, we see that until t = 5 ps a photon is absorbed that is later re-emitted, and, after t = 15 ps, we again observe photon absorption processes. In the following, we analyze this dynamics of the correlated electron-photon problem in terms of population in the cavity Born− Oppenheimer surfaces calculated in Figure 2 (a). In Figure 7, we show the occupation of the photon number states in the first cavity PES in (a) and the third cavity PES in (b). The values (P1,P3) give the population of the first cavity PES and the third cavity PES, respectively. All other cavity PES have populations which are an order of magnitude smaller, since P1+P3 is close to 1 for all times. In Figure 7 (a), we find that at the initial time t = 0 ps, the first cavity PES is populated with a photon state, which has a coherent distribution with ⟨aα † aα⟩ = 4, which is in agreement with our initial condition. During the time propagation, we observe a transfer of population from the first cavity PES to the third cavity PES. In the first cavity PES, we see until t = 9.3 ps a depletion of population, while in the third cavity PES (Figure 7 (b)), we observe an increase of the population. After this time, the population is again transferred back from the third cavity PES to the first cavity PES (Rabi oscillation). However, not only the amplitude of the population is changing but also the center of the wave packets. In principle, if the same photon state would be populated in the two different cavity PES, the system could still be factorizable. For small times, up to t = 5 ps the center of the wave packet in the first cavity PES remains close to its initial value. Later it changes to smaller photon numbers, which indicates photon absorption. We can conclude that the dynamics of the many-body system is dominated by the population transfer from the first cavity PES to the third cavity PES and vice versa. While for this example, a good approximate description may be a two-surface approximation reminiscent of the Rabi model, 42 we expect a different behavior for more complex cavity Born−Oppenheimer surfaces e.g. in many-electron problems, multiphoton modes, or strongcoupling situations.

Light-Matter
Coupling via Vibrational Excitation. The second system that we analyze is the Shin-Metiu model 48,49 coupled to cavity photons. Without coupling to photon modes, this system exhibits a conical intersection between Born−Oppenheimer surfaces and has been analyzed heavily in the context of correlated electron−nuclear dynamics, 50 exact forces in nonadiabatic charge transfer, 51 or nonadiabatic effects in quantum reactive scattering, 52 to mention a few. In our case, we place the system, consisting of three nuclei and a single electron into a optical cavity, where it is coupled to a single mode that is in resonance with the first vibrational excitation. The outer two nuclei are fixed, and the free electron and the nuclei are restricted to one-dimension. The model is schematically depicted in Figure 8. The Hamiltonian of such a system is given by 48 where V n (R) is the Coulomb interaction of the free nuclei with the two fixed nuclei, r is the electronic coordinate, and R is the   Figure 6.

Journal of Chemical Theory and Computation
Article nuclear coordinate. V e is the sum of the electron interaction with the three nuclei, i.e. three terms each of which is of the following form 49 where x is the electron−nuclear distance, and erf describes the error-function. We fix the nuclear mass M to the mass of a hydrogen atom, choose Z = 1, and set the length L = 10 Å. Furthermore, we use the dipole operators X e = −er and X n = eR.
In the nuclear dipole moment operator X n , the two outer nuclei cancel each other due to their fixed positions at ±L/2. Further R c can be used to tune the energy difference Δ between the ground-state and the first-excited state potential energy surface. For the cavity Shin-Metiu model, we represent the electron on a grid of dimension N r = 140 with Δr = 0.4233 Å, and the nuclear coordinate on a grid of dimension N R = 280 with ΔR = 0.0265 Å, while the photon wave function is expanded in the photon number eigenbasis, where the mode can host up to 81 photons in the photon mode. To get first insights on how the light-matter coupling is capable of changing the chemical landscape of the system, in Figure 9, we calculate the ordinary PES surfaces of eq 15 for the case of q α = 0. The solid red line shows the ground-state energy surface, while the blue line shows the excited state energy surface for R c = 1.5 Å with ℏω α = 72.5 meV and R c = 1.75 Å with ℏω α = 69.3 meV. In both examples, the photon frequencies ω α correspond to the first vibrational transition of the exact bare Hamiltonian. Next, we tune the matter-photon coupling strength λ α from the weakcoupling regime to the strong-coupling regime. The corresponding cavity PES are shown in gray in Figure 9. The inset in the figures shows the energy gap Δ depending on the matterphoton coupling strength λ α . In the left figure, we choose the value R c = 1.5 Å, and, in the case of λ α = 0, we find well separated cavity Born−Oppenheimer surfaces. The matterphoton coupling (chosen here from λ α = 0 to λ α = 82.55 eV 1/2 / nm with a Rabi splitting Ω R = (E 5 −E 3 )/ℏω α = 43.81%) opens the gap significantly, as shown in the inset. Additionally, for R c = 1.5 Å, we find that the double-well structure visible in the first-excited state becomes more pronounced for stronger lightmatter coupling. The right figure shows the results for R c = 1.75 Å, where in the field-free case a much narrower gap Δ is found. Introducing the matter-photon coupling in the system from λ α = 0 to λ α = 84.48 eV 1/2 /nm with Ω R = 64.04% also opens the gap significantly, and we find a similar qualitative behavior as in the previous example with the notable difference, that we observe in the present example a similar single-well to double-well transition but now in the first-excited state. However, since we restricted ourselves to a specific cut in the full two-dimensional cavity Born−Oppenheimer surface by choosing q α = 0, Figure 9 does not show the full picture. Therefore, in Figure 10, we show the full two-dimensional cavity PES for R c = 1.75 Å. In the figure, the x-axis shows the nuclear degree of freedom (R), while the y-axis shows the photonic degree of freedom q α . In the case of λ α = 0, that is the upper panel in the figure, we find that the photonic degree of freedom introduces harmonicity into the surface. We also indicate the minima in the surfaces by white crosses. In agreement with Figure 9, we find a double minimum for the ground-state cavity PES and a single minimum for the excited state cavity PES. In the case of strong-coupling that is shown in the lower panel of the figure, we observe new emerging normal modes. These new normal modes are caused by the entanglement of the matter and photon degrees of freedom and are manifest in the displacement of the minima out of the equilibrium positions. In the first-excited state surface in strong coupling, we also observe a single-well to double-well transition, as observed in the coupling to the electronic excitation and discussed in the first

Journal of Chemical Theory and Computation
Article part of this work. Here, we find that now two minima appear in the first-excited state surface. If we adopt an adiabatic picture we can conclude that now two new reaction pathways are possible from the first excited state surface to the ground-state surface.
To conclude, we have seen how the photonic degrees of freedom alter considerably chemical properties in a model system containing electronic, nuclear, and photonic degrees of freedom. We have identified the change of traditional Born− Oppenheimer surfaces, gap opening, and transitions from single-well structures to double-well structures in the firstexcited state surface from first principles. The gap opening can be connected to recent experiments, 53 where a reduction in chemical activity has been observed for vibrational strong coupling.

SUMMARY AND OUTLOOK
In this paper, we introduced the concept of the cavity Born− Oppenheimer approximation for electron−nuclear-photon systems. We used the cavity Born−Oppenheimer approximation to analyze the ground-state transition in the system that emerges in the strong-coupling limit. During this transition the ground-state electron density is split, and the ground-state cavity PES obtains a double-well structure featuring finite displacements of the photon coordinate. Furthermore, we illustrated for a time-dependent situation with a factorizable initial state how the complex correlated electron-photon dynamics can be interpreted by an underlying back-and-forth photon population transfer from the ground-state cavity PES to an excited-state cavity PES. In the last section, we have demonstrated how this transition can also appear in case of strong-coupling and vibrational resonance. Here, we find that the first-excited state surface can obtain a double-well structure leading to new reaction pathways in an adiabatic picture. In future studies toward a full ab initio description for cavity lightmatter systems, where solving the electronic Schrodinger equation of eq 11 by exact diagonalization is not feasible, the density-functional theory for electron-photon systems can be used. 25, 26 The discussed methods can be still improved, e.g. along the lines of a more accurate factorization method such as the exact factorization 54−56 known for electron−nuclear problems, or trajectory based methods 50,57 can be applied to simulate such systems dynamically. This work has direct implications on more complex correlated matter-photon problems that can be approximately solved employing the cavity Born−Oppenheimer approximation to better understand complex correlated light-matter coupled systems.