Conditional wavefunction theory: a unified treatment of molecular structure and nonadiabatic dynamics

We demonstrate that a conditional wavefunction theory enables a unified and efficient treatment of the equilibrium structure and nonadiabatic dynamics of correlated electron-ion systems. The conditional decomposition of the many-body wavefunction formally recasts the full interacting wavefunction of a closed system as a set of lower dimensional (conditional) coupled `slices'. We formulate a variational wavefunction ansatz based on a set of conditional wavefunction slices, and demonstrate its accuracy by determining the structural and time-dependent response properties of the hydrogen molecule. We then extend this approach to include time-dependent conditional wavefunctions, and address paradigmatic nonequilibrium processes including strong-field molecular ionization, laser driven proton transfer, and Berry phase effects induced by a conical intersection. This work paves the road for the application of conditional wavefunction theory in equilibrium and out of equilibrium ab-initio molecular simulations of finite and extended systems.


I. INTRODUCTION
Emerging experimental capabilities in the precise manipulation of light and matter are opening up new possibilities to understand and exploit correlations and quantum effects that can be decisive in the functional properties of molecules and materials. Light-driven states can be not only designed to monitor and/or control the structure of molecules [1][2][3][4][5][6][7] and solids [8][9][10][11][12], but also to form light-matter hybrid states with new physical properties [13][14][15][16][17][18][19][20][21]. In view of these exciting developments, accurate first principles theoretical techniques are also needed in order to help interpret observations, to enable the predictions of simplified models to be scrutinized, and ultimately, to help gain predictive control. Our ability to treat the full correlated quantum structure and dynamics of general electron-ion systems unfortunately remains limited by the unfavourable scaling of the many-body problem.
A standard approach to address this problem in molecular and solid-state systems has been to 'divide-and-conquer' in the sense that the electronic structure and the electron-nuclear interactions are treated separately. Introduced almost a century ago by Born and Oppenheimer [22], the adiabatic approximation, i.e., the assumption that electrons adjust instantaneously to the motion of nuclei, is the cornerstone of this socalled standard approach. The Born-Oppenheimer (BO) approximation has been crucial to the development of a vast majority of approaches in quantum chemistry and condensed-matter theory [23,24], and the concept of ground state Born-Oppenheimer potential-energy surface (BOPES) is the foundation for understanding the properties of systems at thermal equilibrium such as chemical reactivity [25][26][27] and nuclear quantum effects [28][29][30][31], as well as of systems driven out of equilibrium [32][33][34][35].
Accurately describing systems driven away from equilibrium and including nonadiabatic electron-nuclear effects, places even more stringent demands on the development of practical first principles tools.
In the standard approach one directly builds upon the BO approximation by expanding the full molecular wavefunction in the Born-Huang basis [36]. Within this framework, nonadiabatic processes can be viewed as nuclear wavepacket dynamics with contributions on several BOPESs, connected through nonadiabatic coupling terms that induce electronic transitions [37]. In this picture, trajectory based quantum dynamics methods offer a trade-off between physical accuracy and computational cost [38][39][40]. Of these approaches, perhaps the most popular are Ehrenfest mean field theory [41], and Tully's surface hopping dynamics [42]. Both of these approaches consist of an ensemble of uncorrelated trajectories. Reintroducing correlation, for example by using a variety of wavefunction ansätze [43][44][45][46][47][48], semiclassical techniques [49,50], the quantum-classical Liouville equation [51][52][53], path-integral methods [54,55], or methods based on the exact factorization [56][57][58], allows for further accuracy with increased computational effort.
While advances in ab initio electronic structure theory in quantum chemistry and condensed matter have made computing the ground state energies both routinely efficient and rather accurate in many cases, obtaining accurate excited state information remains a challenging problem in its own right. Even in cases where the excited state electronic structure is available, performing fully quantum nuclear dynamics calculations using the standard approach quickly becomes infeasible [35,59] as the memory required to store the information contained in the BOPESs grows rapidly with the number of correlated degrees freedom. In this respect, gaining the ability to rigorously treat selected nuclear degrees of freedom quantum mechanically without incurring an overwhelming computational cost is the goal.
An alternative approach for describing quantum effects in coupled electron-ion systems is using a real space representation of all degrees of freedom. This route might sound less intuitive as it avoids routine concepts such as BOPESs and nonadiabatic couplings that are fundamental in the present description and understanding of quantum molecular dynamics. However, this feature might be turned into an attractive playground from the computational point of view, as these quantities are usually demanding to obtain and fit from ab initio electronic structure calculations. In this framework, one of the leading approximate methods to describe the coupled electron-nuclear dynamics for large systems is time-dependent density functional theory coupled to classical nuclear trajectories through the Ehrenfest method [60]. Due to its favourable system-size scaling, the real-space picture Ehrenfest method has been successful for a great many applications, from capturing phenomena associated with vibronic coupling in complex molecular systems [61] and photodissociation dynamics in small molecules [62], to radiation damage in metals [63]; its efficiency allows calculations on large systems for even hundreds of femtoseconds [64]. It has also been recently combined with the nuclear-electronic orbital method as a way to include quantum effects for selected nuclear degrees of freedom, to study proton transfer processes in molecular excited states [65].
It is well-known, however, that the Ehrenfest approach can be inaccurate due to its mean-field nature. One classic example of this breakdown occurs in photochemical reaction dynamics, where mean field theory can often fail to correctly describe the product branching ratios [39,66]. Generally speaking, the mean field description of any transport property can potentially suffer some deficiency; this is sometimes referred to as a violation of detailed balance [67], but it ultimately stems from the lack of time-translational invariance that is inherent to any approximate method that does not rigorously preserve the quantum Boltzmann distribution [68].
The conditional wavefunction (CWF) framework introduced in Ref. [69] offers a route to go beyond the limits of mean field theory while retaining a real-space picture; it is an exact decomposition and recasting of the total wavefunction of a closed quantum system [70]. When applied to the time-dependent Schrödinger equation the conditional decomposition yields a set of coupled, non-Hermitian, equations of motion [69]. One can draw connections between CWF theory and other formally exact frameworks proposed in the literature to develop novel approximate schemes that provide a completely new perspective to deal with the long-standing problems of nonadiabatic dynamics of complex interacting systems [71,72]. An example is the time-dependent interacting conditional wavefunction approach (ICWF) [73,74], a recently introduced method for performing quantum dynamics simulations that is multi-configurational by construction. Using a stochastic wave-function ansatz that is based on a set of interacting single-particle CWFs, the ICWF method is a parallelizable technique which achieves quantitative accuracy for situations in which mean-field theory drastically fails to capture qualitative aspects of the dynamics, such as quantum decoherence, using orders of magnitude fewer trajectories than the converged mean-field results [73].
In this work we introduce an exact time-independent version of the CWF mathematical framework. The timeindependent CWF framework is formulated in real-space, and it is an exact decomposition of the time-independent wavefunction of a closed quantum system that yields a set of of coupled nonlinear eigenvalue problems and associated conditional eigenstates. Based on this framework, we put forth a static-basis version of the ICWF method, which allows us to establish an efficient and accurate algorithm for calculating the ground-and excited-state structure of correlated electron-nuclear systems and eventually extended systems. Importantly, the combination of the static version of the ICWF method using a timedependent conditional eigenstate basis sets the stage for the implementation of a general purpose ab initio molecular simulator that is formulated in the real-space picture and that self-consistently treats stationary states as well as driven dynamics.
This manuscript has the following structure: in Sec. II we define the mathematical structure of the timeindependent version of the CWF framework. Based on these results we put forth an imaginary time version of the ICWF technique in Sec. III for solving the time-independent Schrödinger equation and the performance of the resulting algorithm is assessed through the calculation of the ground-state and the low-lying excited state BOPESs of the hydrogen molecule in one-dimension. In Sec. IV a real-time extension of this multiconfigurational ansatz is presented, along with an algorithm for solving the time-dependent Schrödinger equation using a stochastic static-basis ansatz. The ability of the resulting algorithm in capturing static and dynamic properties is then assessed by evaluating the absorption spectrum and a laser-induced dynamics of the aforementioned H 2 model system. In Sec. V we revisit the exact time-dependent CWF framework and in Sec. VI we present the dynamical ICWF approach to the timedependent Schrödinger equation. The performance of the time-dependent ICWF method in combination with its imaginary time variation for preparing the initial state is demonstrated for three model systems, viz., a laser driven proton-coupled electron transfer model, an electron-atom scattering process, and an example of the effect of the Berry phase in the quantum nuclear dynamics through a conical intersection. A summary of the main results of this work and an outlook on future directions is offered in Sec. VII.

II. CONDITIONAL EIGENSTATES
We begin by considering a closed system with n electrons and N nuclei, collectively denoted by x = (r, R). We use the position representation for both subsets; lowercase symbols will be used for the electronic subsystem, e.g., r = {r 1 s 1 , ..., r n s n }, and uppercase symbols R = {R 1 σ 1 , ..., R N σ N } for the nuclear subsystem. Hereafter, electronic and nuclear spin indices, respectively s j and σ j , will be made implicit for notational simplicity, and, unless otherwise stated, all expressions will be expressed in atomic units.
The time-independent CWF can be constructed starting from the non-relativistic time-independent Schrödinger equation in position representation, where Ψ γ (x) is an eigenstate of the molecular Hamiltonian H with label γ, and corresponding energy eigenvalue E γ . The molecular Hamiltonian operatorĤ in Eq. (1) can be written as: where the kinetic energy operators areT j = 1 2mj (−i ∇ j − z j A(x j )) 2 , being m j and z j the characteristic mass and charge of particle j respectively. The full electron-nuclear potential energy of the system is W (x) (written in the position basis rather than, say, the BO or Born-Huang basis), and A is the vector potential due to an arbitrary static external electromagnetic field.
Note that the total Hamiltonian in Eq. (1) is invariant under translations and rotations of all particles. This means that the eigenstates of the system will be invariant under transformations by the translation and rotation group. Together with the inversion symmetry, this implies that all one-body quantities such as the electron density or any nuclear reduced density are constant and that twoparticle position correlation functions only depend on the distance between their arguments. This is obviously not a convenient starting point to describe the structure of a quantum system. The solution to this problem relies on transforming the Hamiltonian to a fixed coordinate system that reflects the internal properties of the system [75]. This is, in general, not a trivial task and hereafter we will assume that Eq. (1) already reflects such internal properties, either by exploiting a particular symmetry of the system or by simply introducing a parametric dependence on, e.g., a fixed (heavy) nuclear position.
At this point we can decompose the eigenstates Ψ γ (x) in terms of single-particle conditional eigenstates of either of the two subsystems, which are defined as follows: Here the index α denotes the particular conditional slice andx i = (x 1 , ..., x i−1 , x i+1 , ..., x n×N ) are the coordinates of all degrees of freedom in the system except ., x α n×N ) are some particular positions of all system degrees of freedom except x i . As shown schematically in Fig. 1, the conditional eigenstates in Eq. (3) represent one-body slices of the full many-body eigenstates Ψ γ (x) taken along the coordinate of the i-th degree of freedom. The particle placement x α defining the CWFs has not yet been specified, and although in principle it can be chosen arbitrarily, it will be proven convenient in practice to exploit importance sampling techniques.
Evaluating Eq. (1) atx α i by applying the integral operator in Eq. (3) yields conditional eigenstates that are the solutions of the following eigenvalue problem: where we introduced W α i (x i ) = W (x i ,x α i ), with W (x) the full electron-nuclear interaction potential appearing in the Hamiltonian of Eq. (2). In addition, η α,γ i (x i ) are kinetic correlation potentials given by, Provided a large enough collection of CWFs satisfying Eq. (4), an exact solution of Eq. (1) can be reconstructed by undoing the conditional decomposition of Eq. (3) (see Fig. 1.b) [69]. That is, given a set of conditional slices that sufficiently spans the support of Ψ γ , then the corresponding conditional eigenstates can be used to reassemble the full electron-nuclear wavefunction, using the transformations D xi which are discussed in more detail in Appendix A. This expression, Eq. (6), can be used to evaluate the kinetic correlation potentials in Eq. (5). In this way, the generalized one-body eigenvalue problem in Eq. (4) can be understood as an exact decomposition and recasting of the eigensolution of the full electron-nuclear system, that yields a set of coupled, non-Hermitian, eigenvalue problems.

A. Time-independent Hermitian approximation
An approximate solution to Eq. (4) can be formulated by expanding the kinetic correlation potentials around the sampling coordinates x α using Taylor series, and then truncating at zeroth order, i.e.: At this level the kinetic correlation potentials engender only a global phase that can be simply omitted as expectation values are invariant under such global phase transformations. Note that these approximated kinetic correlation potentials can be alternatively obtained by introducing a mean field ansatz Ψ γ (x) = n×N i=1 ψ(x i ) into Eq. (5). By making this approximation the eigenvalue problems in Eq. (4) are restored to a Hermitian form, The Hermitian limit allows the full many-body problem to be approximated as a set of independent single-particle problems. That is, the superscript γ refers exclusively to the conditional eigenstate excitation number.

III. STATIC PROPERTIES WITH CONDITIONAL EIGENSTATES
In general the higher order terms in the Taylor expansion of the kinetic correlation potentials are nonnegligible. However, one can still take advantage of the simple Hermitian form of the conditional eigenvectors (hereafter referred to as conditional wavefunctions CWFs) in Eq. (8) to design an efficient many-body eigensolver by utilizing them as bases for electronic and nuclear degrees of freedom in a variational wavefunction ansatz.
While there is a diverse literature spanning decades on different forms for variational electron-nuclear wavefunction ansätze, for illustrative (and practical) purposes we employ a sum-of-products form, which in the language of tensor decompositions is referred to as the canonical format [76]. For each degree of freedom x i we utilize a given electronic or nuclear CWF, respectively, coming from solutions to Eq. (8), in order to approximate the γ th full system exact excited state as follows: Where in the second line we have rearranged the sum over particle position λ ∈ {1, . . . , N c } and excited CWF ν ∈ {1, . . . , M} into a single index α = λ + N c (ν − 1), such that α ∈ {1, . . . , N c M}. The particle placement x α defining the conditional potentials W α i has not yet been specified, and in principle it can be chosen arbitrarily, however in practice we choose to sample from initial guesses for the reduced densities of the electronic and nuclear subsystems.
We refer to this ansatz (Eq. 9) as being in canonical format because we do not mix all possible CWFs ψ λ,ν i for all possible degrees of freedom x i , as one does with a single particle function bases across the different system degrees of freedom in the Tucker format employed in the Multi-Configurational Time-Dependent Hartree (Fock) -MCTDH(F) [43,59], and Multi-Configurational Electron Nuclear Dynamics ansätze [77]. In principle this choice can be relaxed, and one can utilize various choices of tensor network representation for the expansion coefficients C, such as Matrix Product States or hierarchichal Tucker formats, which when employed in the Multi-Layer extension [78,79] of MCTDH allow for an increase in efficiency for certain problems. However, since the time dependence of the ansatz in Eq. (9) is entirely within the expansion coefficients, one only needs to calculate the matrix elements at time zero, creating a quite efficient time propagation framework. Note that although we use a simple Hartree product over electronic degrees of freedom, the above ansatz can be straightforwardly extended to have fermionic anti-symmetry via treating the CWFs as the spatial component of spin orbitals in Slater determinants.
Hereafter, and for reasons that will be apparent later, we will call Eq. (9) the static-basis ICWF (or sta-ICWF) ansatz. With this ansatz in hand, we then consider a solution of Eq. (1) based on the imaginary time propagation technique [80], i.e.: wherê andP ζ = Ψ ζ Ψ ζ † are projectors used to remove the wavefunctions Ψ ζ from the Hilbert space spanned byĤ. The first excited state, for instance, is thus obtained by removing the ground state from the Hilbert space, which makes the first excited state the ground state of the new Hamiltonian,Ĥ ζ=1 .
By introducing the ICWF ansatz of Eq. (9) into Eq. (10) we find an equation of motion for the coefficients where C ξ = C ξ C ξ, † , and the matrix elements of H and S are: where again, the α, β indices refer to the index over particle placement and excited CWFs. Obtaining these matrix elements involves a sum over all two-body interactions across each degree of freedom, and a sum across one-body operators. In practice S may be nearly singular, but its inverse can be approximated by the Moore-Penrose pseudoinverse.
Based on solving the system of equations in Eq. (12) for C γ , one already has the ingredients to put forth a time-independent ICWF eigensolver algorithm that will ultimately be used to evaluate the expectation value of generic observablesÔ(x). Given an approximate solution to the eigenfunction Ψ γ (x), the expectation value of O reads: with the matrix elements of O being given by an analogous expression to Eq. (14).

A. Example I: Ground and Excited BOPESs of H2
As an illustrative example we now calculate the BOPESs of a model for the H 2 molecule. We adopt a model where the motion of all particles is restricted to one spatial dimension, and the center-of-mass motion of the molecule can be separated off [81]. In this model the relevant coordinates are the internuclear separation, R, and the electronic coordinates, r 1 and r 2 . The Hamiltonian, written in terms of these coordinates, is where µ n = M/2 and µ e = M/(2M + 1) are the reduced nuclear and electronic masses respectively, and M is the proton mass. In Eq. (16) the electron-electron repulsion and the electron-nuclear interaction are represented by soft-Coulomb potentials, i.e., the Coulomb singularities are removed by introducing a smoothing parameter ee = 2 and en = 1. The above model system qualitatively reproduces all important strong-field effects such as multiphoton ionization, above threshold ionization, or high-harmonic generation [82][83][84].
Moreover, it has provided valuable information in the investigation of electron correlation effects [85][86][87]. For this model, the BOPESs are defined by the following electronic eigenvalue problem: (19) where H el =Ĥ −T nuc , and {Φ γ (r 1 , r 2 ; R)} are the (complete, orthonormal) set of BO electronic states. A parametric dependence on the nuclear coordinates is denoted by the semicolon in the argument. The BOPESs, γ (R), can be calculated using the imaginary time sta-ICWF method described in Eqs. (10)- (14) along with a simplified version of the ansatz in Eq. 9 that is specialized to this particular case of parametric nuclear dependence. A thorough description of the numerical procedure as well as the convergence behaviour of the sta-ICWF method for this model can be found in Appendix B 1.
In Fig. 2 we show the first five BOPESs calculated via the sta-ICWF approach using (N c , M) = (32,5). In the top panel, the exact BOPESs are plotted against the sta-ICWF data, overlaid as solid gold lines. The results demonstrate that the sta-ICWF ansatz used in a variational context captures the entire group of excited BOPES landscape over this energy range. As a point of comparison, in the bottom panel of Fig. 2 we also show the result of mean-field type calculations of the BOPESs for this system. Specifically, we show Hartree-Fock and configuration interaction singles (CIS) data for the ground state and excited state BOPESs respectively, which suffer from well-known inaccuracies in capturing the binding energy and excited state properties of the system.

IV. TIME-DEPENDENT PROPERTIES WITH CONDITIONAL EIGENSTATES
The sta-ICWF eigensolver described above can be easily extended to describe dynamical properties. For that, we consider the time-dependent Schrödinger equation: where Ψ(x, t) is the electron-nuclear time-dependent wavefunction and the Hamiltonian of the systemĤ(t) may contain a time-dependent external electromagnetic field.
In practice, we are interested in situations where the initial wavefunction is the correlated electron-nuclear ground state, i.e., Ψ(x, 0) = Ψ γ=0 (x), and some nonequilibrium dynamics is triggered by the action of an external driving field. We can then decompose the time-dependent many-body wavefunction as in Eq. (9) by restricting it to the case of γ = 0 (hereafter we omit the superscript γ for clarity). The expansion coefficients C α then obey an equation of motion that can be obtained either by inserting Eq. (9) directly into Eq. (20) or by utilizing the Dirac-Frenkel variational procedure [88]: In Eq. (21), the matrix elements of S and H are identical to the ones defined in Eqs. (13) and (14), with the hamiltonian's time dependence coming from any external fields. The values of the coefficients at time t = 0, i.e., C(0), may be obtained from the imaginary time sta-ICWF method of Eq. (12). In this way, the combination of the imaginary time and real-time sta-ICWF methods yield a "closed-loop" algorithm for the structure and dynamics of Here we demonstrate an application of the real-time sta-ICWF approach to simulate the optical absorption spectrum for the molecular Hydrogen model introduced in Sec. III A. We utilize the "δ-kick" method of Yabana and Bertsch [89], where an instantaneous electric field E(t) = κδ(t), with perturbative strength κ 1a.u. −1 couples to the dipole moment operator µ = r 1 + r 2 , and thereby produces an instantaneous excitation of the electronic system to all transition dipole allowed states. The resulting (linear) absorption spectra can then calculated via the dipole response, ∆µ(t) = µ(t) − µ(0 − ) In practice, due to the finite time propagation, the integrand is also multiplied by a mask function M(t) that smoothly vanishes at the final simulation time T f . The system is first prepared in the ground state using the imaginary time sta-ICWF. See Appendix B 2 for a thorough description of the imaginary-time sta-ICWF method and its use for preparing the ground state the H 2 model system. The field-driven dynamics is then generated by applying the kick operator to the relevant degree of freedom. A thorough description of the numerical procedure as well as the convergence behaviour of the sta-ICWF method for this model can be found in Appendix B 3.
For the H 2 model the occupation of excited electronic states and subsequent coupled electron-nuclear dynamics produce a characteristic vibronic peak structure usually explained via the Franck-Condon vertical transition theory. In the top panel of Fig. 3 we show vibronic spectra calculated both with sta-ICWF for the absorption from S 0 to S 2 in comparison with the numerically exact results, also calculated via the δ-kick approach. For sta-ICWF, we found that N c = 4096 and M = 3 was sufficient to obtain accurate results. The results demonstrate that the sta-ICWF ansatz used in a variational context achieves an accurate vibronic spacing, and furthremore it captures not only the electron-nuclear correlation inherent to vibronic spectra, but also solves the electron-electron subsystem accurately. The deviation from the exact results does grow with increasing energy, although this is ameliorated with increasing N c and M, and can in principle be eliminated at large enough values of these parameters (see Appendix B 3).
For comparison, we show also mean-field results for the vibronic spectra. Specifically, we calculated the absorption spectrum with the multi-trajectory Ehrenfest δ-kick method (MTEF-kick) [61], overlaid as dashed blue lines. We see that the vibronic spacing calculated with the MTEF-kick approach fails in capturing the correct peak spacing in addition to showing unphysical spectral negativity.

B. Example III: Laser driven dynamics of H2
The present formalism is not restricted to just perturbative fields and can deal with any arbitrary external field. Going beyond the linear response regime, we investigate the effect of strong driving by a few-cycle ultra-fast laser pulse for this same H 2 model system. The system is first prepared in the ground state using the imaginary time sta-ICWF, and then the field-driven dynamics is generated by applying an electric field of the form E(t) = E 0 Ω(t) sin(ωt), with E 0 = 0.005a.u. and an envelope Ω(t) with a duration of 20 optical cycles. The carrier wave frequency ω = 0.403 is tuned to the vertical excitation energy between the ground and second excited BOPESs at the mean nuclear position of the ground state wavefunction. A thorough description of the numerical procedure as well as the convergence behaviour of the sta-ICWF method for this model can be found in Appendix B 4.
The intense laser pulse creates a coherent superposition of the ground and second excited BO states whereby the bond length of the molecule increases, as shown in the bottom panel of Fig. 4. The nuclear wavepacket then eventually returns to the Franck-Condon region, creating the resurgence of the electronic dipole oscillation seen in the top panel of Fig. 4. In the MTEF meanfield description of this process the short-time limit is rather accurately captured, while the subsequent effects of the laser pulse on the nuclear dynamics and the resurgence in the dipole response are not. These results show that the sta-ICWF method is able to capture the electronic correlations inherent to the electronic dipole moment during the initial laser driven dynamics, as well as the electron-nuclear correlations that arise during the subsequent nonequilibrium dynamics. For this particular problem we found that (N c , M) = (4096, 3) was sufficient to obtain highly accurate results for both the expectation value of the electronic dipole moment (top panel of Fig. 4) and the expectation value of the inter-nuclear separation (bottom panel of Fig. 4). Further details can be found in Appendix B 4.

V. TIME-DEPENDENT CONDITIONAL WAVEFUNCTIONS
While the sta-ICWF method shows promising performance in the examples studied thus far, it faces the same limitations as any method that relies on a static basis. Perhaps the most significant aspect can be framed in terms of capturing the full support of the time-dependent wavefunction, which is exacerbated in cases where the time dependent state strays far from the span of the static basis. One strategy to address these scenarios would be to incorporate time-dependent conditional wavefunctions in the ICWF ansatz. Hence, we take advantage of the timedependent version of the CWF framework introduced in Refs. [69], which relies on decomposing the exact manybody wavefunction, Ψ(x, t), in terms of time-dependent single-particle CWFs of either the electronic or nuclear subsystems as: Evaluating the time-dependent Schrödinger equation in Eq. (20) at x α i (t), one can show that the CWFs in Eq. (23) obey the following equations of motion: where t), and we remind that W (x) is the full electron-nuclear interaction potentials that appears in the Hamiltonian of Eq. (2). In Eq. (24), η α i (x i , t) are time-dependent complex potentials containing kinetic correlations and advective terms, i.e.: As in the time-independent CWF framework, the conditional wavefunctions in Eq. (23), represent slices of the full wavefunction taken along single-particle degrees of freedom of the two disjoint subsets. Each individual CWF constitutes an open quantum system, whose timeevolution is non-unitary, due to the complex potentials η α i (x i , t), which now include advective terms due to the inherent motion of the trajectories x α (t), which evolve according to Bohmian (conditional) velocity fields [69]: An exact solution to Eq. (20) can be then constructed provided we use a sufficiently large number of slices {x α (t)} that explore the full support of |Ψ γ (x, t)| 2 (in analogy with Fig. 1.b), i.e.: where the transformations can be found in Appendix A. The one-body equations of motion in Eq. (24) can be then understood each as a coupled set of non-unitary and nonlinear time-dependent problems. The derivation of the exact time-dependent CWF mathematical framework corresponds to the transformation of the many-body time-dependent Schrödinger equation to the partially co-moving frame in which all coordinates except the i-th move attached to the electronic and nuclear flows and only the i-th coordinate is kept in the original inertial frame. Within the new coordinates, the convective motion of all degrees of freedom except for the i-th coordinate is described by a set of trajectories of infinitesimal fluid elements (Lagrangian trajectories), while the motion of the i-th degree of freedom is determined by the evolution of the CWFs in an Eulerian frame [72]. The purpose of this partial time-dependent coordinate transformation is to propagate all trajectories along with the corresponding probability density flow such that they remain localised where the full molecular wavefunction has a significant amplitude.

A. Time-dependent Hermitian approximation
In general the effective potentials in Eq. (25) exhibit discontinuous steps which could introduce instabilities in a trajectory-based solution of the many-body dynamics based on Eq. (24). Therefore, in a similar manner to the time-independent case, an approximate solution can be formulated by expanding the kinetic and advective correlation potentials around the conditional coordinates x α (t), such that In this limit, the kinetic and advective correlation potentials only engender a global phase that can be omitted, as expectation values are invariant under such global phase transformations. The resulting propagation scheme is restored to a Hermitian form. That is, Eq. (24) is approximated as: while the trajectories x α (t) are constructed according to Eq. (26). This approximation to the time-dependent CWF formalism is clearly a major simplification of the full problem, as it recasts the many-body time-dependent Schrödinger equation as a set of independent singleparticle equations of motion. Despite the crudeness of the approximation in Eq. (28), the set of equations of motion in Eq. (29) has found numerous applications, e.g., in the description of adiabatic quantum molecular dynamics [71] and quantum electron transport [90][91][92][93][94].

VI. SIMULATING FAR-FROM-EQUILIBRIUM DYNAMICS WITH CONDITIONAL WAVEFUNCTIONS
In general circumstances where the kinetic and advective correlation potentials are important, we can make use of the simple Hermitian form of the conditional equations of motion in Eq. (29) to design an efficient manybody wavefunction propagator. For that, we expand the full electron-nuclear wavefunction using the ansatz: where the coefficients C α (t) and the CWFs ψ α i (x i , t) are initialized using the sta-ICWF method and propagated afterwards using the approximated equations of motion in Eq. (29) along with trajectories obeying Eq. (26).
The time evolution of the coefficients C(t) can be then obtained by inserting the ansatz of Eq. (30) into Eq. (20), where the matrix elements of S, H, are defined as in Eqs. (13) and (14), with the time dependence coming from external fields in the hamiltonian and the time dependent CWFs, while H i are: where h α i (t) are the Hermitian Hamiltonians in Eq. (29) andĤ(t) is the full time-dependent Hamiltonian in Eq. (20).
Obtaining these matrix elements is straightforward, involving a sum across single body operators in Eq. (13) and Eq. (32), and all sums of two-body interactions across each degree of freedom in Eq. (14). Note that any operator involving only a single species, e.g. the kinetic energy, is cancelled out, and thus the evolution of C is governed exclusively by matrix elements of operators which either fully (through H) or conditionally (through H i and H A ) correlate the degrees of freedom.
Equations (26), (29), and (31) define a set of coupled differential equations that hereafter will be referred to as the dynamical ICWF (dyn-ICWF) method. One can then evaluate the expectation value of a generic observable Ô (x) as given in Eqs. (15) with dyn-ICWF by simply taking into account that ψ α i (t) are now time-dependent CWFs.

A. Example IV: Impact electron ionization
The theoretical description of electron scattering remains challenging, as it is a highly-correlated problem that generally requires treatment beyond perturbation theory [95,96]. We here study a model system of electron-Hydrogen scattering that can be exactly solved numerically [97]. In atomic units, the Hamiltonian of this one-dimensional two-electron model system reads: where W (r 1 − r 2 ) = 1 v ext (r) = 1 (r − 10) 2 + 1 . are respectively the soft-Coulomb interaction and the external potential that models the H atom located at r = 10a.u. The initial interacting wavefunction is taken to be a spin singlet, with a spatial part where φ H (r) is the ground state hydrogen wavefunction. φ W P (r) is an incident Gaussian wavepacket, with α = 0.1, representing an electron at r = −10a.u., approaching the target atom with a momentum p. The time-resolved picture presents scattering as a fully non-equilibrium problem, where the system starts already in a non-steady state, and so, the imaginary time sta-ICWF cannot be applied here to prepare the initial wavefunction. Instead, we stochastically sample the initial probability density |Ψ 0 (r 1 , r 2 )| 2 with N c trajectories {r α 1 (0), r α 2 (0)} that are used to construct CWFs φ α 1 (r 1 , 0) and φ α 2 (r 2 , 0) as defined in Eq. (23). A thorough description of the numerical procedure as well as the convergence behaviour of the dyn-ICWF method for this model can be found in Appendix C 1.
We study the dynamics of the electron-Hydrogen scattering by evaluating the time-dependent one-body density, ρ e (r 1 , t) = 2 |Ψ(r 1 , r 2 , t)| 2 dr 2 , for two different initial momenta, viz., p = 0.3a.u and p = 1.5a.u. For p = 0.3a.u., the energy is lower than the lowest excitation of the target (which is about ω = 0.4a.u.) and hence the scattering process is elastic. In this regime, mean-field results (here represented by extended time-dependent Hartree-Fock calculations) and dyn-ICWF results with N c = 128 results both capture the correct dynamics accurately. In approaching the target atom with the larger momentum p = 1.5a.u., the incident wavepacket collides inelastically with the target electron at around 0.24 fs, after which, a part of the wavepacket is transmitted while some is reflected back leaving the target partially ionized. In this regime, the mean-field method fails to describe the transmission process quantitatively, and the reflection process even qualitatively due to its inability to capture electron-electron correlation effects. This is in contrast with dyn-ICWF results, which quantitatively capture the correlated dynamics for N c = 256, although a lower number of CWFs already reproduces qualitatively the dynamics (see Appendix C 1).

B. Example V: Laser Driven Proton-Coupled Electron Transfer
We now show dyn-ICWF results for a prototypical photo-induced proton-coupled electron transfer reaction, using the Shin-Metiu model [98]. The system comprises donor and acceptor ions which are fixed at a distance L = 19.0a 0 , and a proton and an electron that are free to move in one dimension along the line connecting the donor-acceptor complex. Based on the parameter regime chosen, this model can give rise to a number of challenging situations where electron-nuclear correlations play a crucial role in the dynamics.
The total Hamiltonian for the system is, where m is the electron mass, and M is the proton mass. The coordinates of the electron and the mobile ion are measured from the center of the two fixed ions, and are labeled r and R, respectively. The full electron-nuclear potential reads: where erf() is the error function. The parameter regime studied for this model (R f = 5a 0 , R l = 4a 0 and R r = 3.1a 0 ) and are chosen such that the ground state BOPES, BO , is strongly coupled to the first excited adiabatic state, BO , around the mean nuclear equilibrium position R eq = −2a 0 . The coupling to the rest of the BOPESs is negligible. We set the system to be initially in the full electronnuclear ground state obtained from out imaginary time propagation method described above, i.e., Ψ(r, R, 0) = Ψ 0 (r, R). We then apply an external strong electric field, E(t) = E 0 Ω(t) sin(ωt), with E 0 = 0.006a.u, Ω(t) = sin(πt/20) 2 and ω = 1 BO (R eq ) − 0 BO (R eq ). The external field induces a dynamics that involves a passage through an avoided crossing between the first two BOPESs, with further crossings occurring at later times as the system evolves. When the system passes through the nonadiabatic coupling region, the electron transfers probability between the ground state and the first excited state. This is shown in the top panel of Fig. 6, where we monitor the BO electronic state populations P n (t) (whose definition can be found in Appendix C 2). As a result of the electronic transition, the reduced nuclear density changes shape by splitting into two parts representing influences from both ground and excited state BOPESs. This can be seen in the bottom panel of Fig. 6, where, as a measure of decoherence, we use the indicator D nm (t) (whose definition can be found in Appendix C 2). As nonadiabatic transitions occur the system builds up a degree of coherence which subsequently decays as the system evolves away from the coupling region. As shown in Fig. 6, the dyn-ICWF method reaches quantitative accuracy for (N c , M) = (256, 1), and vastly outperforms the multi-trajectory Ehrenfest mean-field method in describing both the adiabatic populations and the decoherence measure. More specifically, while both the dyn-ICWF method and MTEF dynamics correctly capture the exact adiabatic population dynamics at short times, the latter breaks down at long times as it fails to capture the qualitative structure of the time-evolving indicator of decoherence. Noticeably, all these aspects of this problem are qualitatively well decribed by the dyn-ICWF method using only (N c , M) = (16, 1) (these results can be found in Appendix C 2).

C. Example VI: Berry phase effects and molecular conical intersection
We next study dynamics around conical intersections (CIs) using a minimal generalization of the above Shin-Metiu model first proposed by Gross and co-workers [99], and extended further by Schaupp and Engel [100]. The model consists of a quantized electron and proton that can move in two Cartesian directions, along with two fixed 'classical' protons, R 1 , R 2 . A CI occurs in this model when (treating the quantized proton as a BO parameter) the protons are in a D 3h geometry. The potential energy is, and we use the parameter values a = 0.5, b = 10, R 0 = 1.5, R 1 = −0.4 √ 3, 1.2 , R 2 = 0.4 √ 3, 1.2 . We initialize the total system wavefunction as a direct product of the first excited electronic BO state and a nuclear Gaussian state centered at R 0 = (0, 0.4) with standard deviation σ 2 = 5. For this placement of R 1 , R 2 the CI occurs at the origin, and in the BO picture, the initial nuclear wavepacket "falls towards" the CI (see Fig.  C.3 in Appendix C 3), while the Berry phase associated with the two possible paths around the CI cause an interference pattern to develop. Using dyn-ICWF and propagating entirely in the real space grid picture, this characteristic interference pattern can also be captured. Therefore, while not depending on the BO picture (beyond defining the initial state) the dyn-ICWF method retains the correct Berry curvature effects. See Appendix C 3 for further details on the dyn-ICWF calculation.

VII. CONCLUSIONS
In this work we have introduced an exact mathematical framework that avoids the standard separation between electrons and nuclei and hence enables a unified treatment of molecular structure and nonadiabatic dynamics without relying on the construction and fit of Born-Oppenheimer potential-energy surfaces and the explicit computation of nonadiabatic couplings.
We have introduced a time-independent conditional wavefunction theory, which is an exact decomposition and recasting of the static many-body problem that yields a set of single-particle conditional eigenstates. Based on the imaginary time propagation of a stochastic ansatz made of approximated conditional eigenstates, the resulting method, called sta-ICWF, is able to accurately capture electron-electron correlations intrinsic to molecular structure. A real-time counterpart of the above method has been also derived following the Dirac-Frenkel variational procedure, and its combination with the imaginary time version yields an accurate method for solving out of equilibrium properties of molecular systems where nonadiabatic electron-nuclear correlations are important. This has been shown by reproducing the exact structural, linear response, and non-perturbatively driven response properties of an exactly solvable onedimensional H 2 model system that standard mean-field theories fail to describe.
We have also considered a broader class of conditional wavefunctions that was formally introduced through time-dependent conditional wavefunction theory, yielding a set of coupled single-particle equations of motion. An approximated set of these time-dependent conditional wavefunctions are utilized as time-dependent basis of a stochastic wavefunction ansatz that is meant to describe observables that are relevant to far-from-equilibrium processes. The resulting propagation technique (called dyn-ICWF) in combination with sta-ICWF provides a fully self-consistent approach and, moreover, the method achieves quantitative accuracy for situations in which mean-field theory drastically fails to capture qualitative aspects of the combined electron-nuclear dynamics.
Importantly, the conditional decomposition holds for an arbitrary number of subsets (up to the total number of degrees of freedom in the system), and applies to both fermionic and bosonic many-body interacting systems. Our developments thus provide a general framework to approach the many-body problem in and out of equilibrium for a large variety of contexts. For example, using conditional wavefunctions in a form compatible with time-dependent density functional theory, in connection with alternative tensor network decompositions, or in combination with classical/semiclassical limits for specified degrees of freedom, are particularly appealing routes to follow, and work in this direction is already in progress [101]. Furthermore, the extension to periodic systems is currently under investigation and should allow the ab initio description of driven electron-lattice dynamics such as for example laser driven heating and thermalisation [102][103][104][105][106][107], correlated lattice dynamics [108][109][110] and phase transitions [111][112][113].
Assuming that the conditional sampling points,x α i are distributed according to a normalized distribution N (x α i ), one can approximately reconstruct the full wavefunction based on the interpolation with a Gaussian function G σ (x i ) with a given width σ as: In this way, the full wavefunction is reconstructed as a Gaussian weighted average: in the numerator of Eq. (A2), the contribution from each conditional slice α is weighted with a Gaussian distribution, and it becomes larger if the evaluated point,x, is closer to the sampling pointx α . To compensate the non-uniform sampling distribution contribution, the interpolation weight is divided by the distribution function N (x α i ). In addition, the denominator of Eq. (A2) ensures normalization of the interpolation weight.
By considering a dense sampling (N c → ∞), the reconstructed wavefunction of Eq. (A2) can be rewritten as: and substituting Eq. (A1) into Eq. (A3) one obtains: x n×N ). Therefore, for a dense sampling, Ψ Rec,γ Nc,σ (x) can be understood as the convolution of the full wavefunction Ψ(x) and the Gaussian weight G σ (x i ). Furthermore, in the narrow Gaussian width limit (σ → 0), G σ (x i ) can be treated as a Dirac delta function and hence Eq. (A4) can be written as: In conclusion, one can exactly reconstruct the full electron-nuclear wavefunction in terms of conditional wavefunctions using the reassembling operator D xi defined as: .

(A6)
Appendix B: Convergence of the real and imaginary time versions of the sta-ICWF method In this section we discuss the convergence of the imaginary-and real-time sta-ICWF method for the examples in Secs. III A, IV A, and IV B. For that we first notice that, due to the stochastic nature of the sta-ICWF method, given a set of sampling points N c and their conditional eigenstates M, we may also consider a number N in of different sets of N c sampled points and their associated M conditional eigenstates. This can be accounted for by rewriting the expectation value of Eq. (15) as: The dispersion of Ō (t) with respect to N in is then quantified through its standard deviation, i.e.:

Ground and Excited BOPESs of H2
We discuss here the convergence of the imaginary time version of the sta-ICWF method in capturing the ground state and excited state BOPESs for the H 2 model system introduced in Sec. III A. Finding the BOPESs for this particular model is equivalent to solving Eq. (19) using the imaginary time evolution technique: where {Φ γ (r 1 , r 2 ; R)} are the (complete, orthonormal) set of BO electronic states and we have definedĤ ζ el aŝ The BO electronic states, Φ γ (r 1 , r 2 ; R), are then expanded in terms of CWFs with the following simplified version of the ansatz in Eq. 9 that is specialized to the particular case of parametric nuclear dependence: Slicing points (r α 1 , r α 2 ) are generated by sampling from reduced one-body electronic densities, which in this case are simply chosen to be Gaussian functions ρ e (r i ) = Ae −r 2 i /10 . The conditional eigenstates φ α,ν i (r i ; R), for ν ∈ {1, . . . , M } are then evaluated on each slice using the Hermitian approximation, i.e.: The coefficient vector C γ is randomly initialized and then propagated in imaginary time until the target state is reached according to Eq. (12) of the main text withĤ being substituted withĤ el .
To achieve converged results, a grid (0, 9]a.u. for the internuclear separation with 181 grid points is chosen for the nuclear degrees of freedom. For the electron coordinates, the grid covers the interval [-35,+35]a.u. with 200 grid points. The fourth-order Runge-Kutta integration method [120,121] was used to propagate the imaginary time sta-ICWF equations of motion with a time-step dτ = 0.01a.u, and the Moore-Penrose pseudoinversion method with a tolerance of 10 −8 was used to approximate the numerical inversion of the overlap matrix in Eq. (13). Importantly, the matrices S and H of Eq.

Ground state of H2
We investigate here the ground-state energy for the model H 2 introduced in Sec. III A as well as the convergence behaviour of the imaginary time version of the sta-ICWF method in capturing it. We aim to solve Eq. (10), which for this particular model system reduces to: whereĤ is the Hamiltonian in Eq. (16). For that, we choose the conditional eigenstates basis by sampling N c points (r α 1 , r α 2 , R α ) from guesses to the reduced electronic and nuclear densities ρ e (r i ) = A e e −r 2 i /10 and ρ n (R) = A n e −(R−2) 2 respectively. These positions are then used to construct and diagonalize the Hermitian Hamiltonians in Eq. (8). In this way we obtain 3 × N c × M conditional eigenstates φ α,ζ 1 (r 1 ), φ α,ζ 2 (r 2 ), χ α,ζ (R) . Given a random initialization of the coefficients vector C, we then evolve it in imaginary time according to Eq. the imaginary time sta-ICWF equations of motion with a time-step dτ = 0.01a.u, and the Moore-Penrose pseudoinversion method was used to approximate the numerical inversion of the overlap matrix in Eq. (13). Importantly, the matrices S and H of Eq. (12) need only be constructed at the initial time.
From the exact symmetric ground-state wave function, we found an equilibrium separation of R = 2.2a.u. and the ground-state energy is E 0 = −1.4843a.u. We then define the relative error of the sta-ICWF calculation with respect to the exact calculation as E r = H 0 − E 0 /|E 0 |, where The error E r is presented in Fig

Optical absorption spectrum of H2
We discuss here the convergence of the real-time version of the sta-ICWF method in capturing the optical absorption spectrum of the H 2 model system introduced in Sec. III A. The simulation starts with the preparation of the ground state coefficients C(0) using the imaginary time version of the sta-ICWF method. The relevant degree of freedom of the kick operator is then applied to each CWF, the Hamiltonian and inverse overlap matrices are reconstructed, and C is propagated to the desired time according to Eq. (21). A kick strength of κ = 10 −4 a.u −1 was sufficient to generate the kick spectra within the linear response regime and a total propagation time of T f = 1500a.u. was used to generate the spectra, alongside the mask function To achieve convergence, a grid [−35, +35]a.u. with 200 grid points is chosen for the electronic coordinates. The fourth-order Runge-Kutta algorithm was used to propagate the imaginary time sta-ICWF equations of motion with a time step dt = 0.01a.u, and the Moore-Penrose pseudo-inversion method with a tolerance of 10 −8 was used to approximate the numerical inversion of the overlap matrix in Eq. (13). Again, the matrices S and H of Eq. (12) need only be constructed at the initial time.
In . In all of these cases, we considered a number of different initial sampling points, which have been used to calculate the associated (standard deviation) error bars as in Eq. (B2). As the number of conditional eigenstates basis elements in the ansatz expansion of Eq. (9) increases, the variational nature of the method ensures convergence to the exact linear absorption lineshape. Similarly, the error bars shrink as the number of conditional eigenstates in the basis N c ×M allows to span the relevant part of the Hilbert space.

Laser driven dynamics of H2
We discuss here the convergence of the real-time version of the sta-ICWF method in capturing the laser driven dynamics of the H 2 model system introduced in Sec. III A. As explained in Sec. IV B of the main text, the system is first prepared in the ground state using the imaginary time sta-ICWF as explained in Sec. B 2, and then the fielddriven dynamics is generated by applying an electric field of the form E(t) = E 0 Ω(t) sin(ωt), with E 0 = 0.005a.u. and an envelope Ω(t) with a duration of 20 optical cycles. The carrier wave frequency ω = 0.403 is tuned to the vertical excitation between the ground BO state and second excited electronic surface.
For the dynamics we used a grid (0, 9]a.u. for the internuclear separation with 181 grid points is chosen  In this section we discuss the convergence behaviour of the dyn-ICWF method for the examples of Secs. VI A, VI B, and VI C. As it happened for the sta-ICWF method, the stochastic nature of the dyn-ICWF method allows us to consider a number N in of different initial sampling points for a given set of parameters (N c , M). This is taken into account by writing expectation values as in Eq. (B1) and its standard deviation as in Eq. (B2).

Impact electron ionization
We discuss here the convergence behaviour of the dyn-ICWF method in capturing the laser driven protoncoupled electron transfer described in Sec. VI A.
The time-resolved picture presents scattering as a fully non-equilibrium problem, where the system starts already in a non-steady state, and so, the imaginary time sta-ICWF cannot be applied here to prepare the initial wavefunction. Instead, we stochastically sample the initial probability density |Ψ 0 (r 1 , r 2 )| 2 with N c trajectories {r α 1 (0), r α 2 (0)} that are used to construct CWFs φ α 1 (r 1 , 0) and φ α 2 (r 2 , 0) as defined in Eq. (23). These CWFs are then used to construct the ansatz in Eq. (30) with an initial C vector that is obtained using C(0) = S −1 G, where G is the vector containing the overlap between the initial wavefunction and the CWFs, i.e.: Given C(0), and φ α 1 (r 1 , 0) and φ α 2 (r 2 , 0) for an ensemble of sampling points {r α 1 (0), r α 2 (0)}, these objects are then propagated according to the dyn-ICWF equations of motion.
To achieve converged results, we choose the size of the simulation box to be 150 × 150a.u 2 with an homogeneous grid consisting of 500 grid points in each direction. The fourth-order Runge-Kutta algorithm was used to propagate the dyn-ICWF equations of motion with a time-step dt = 0.01a.u, and the Moore-Penrose pseudoinversion method with a tolerance of 10 −8 was used to approximate the numerical inversion of the overlap matrix in Eq. (31).
In Fig. C.1, we show the one-body electronic density ρ e (r 1 , t), for two different initial momenta and final times, viz., p = 0.3a.u and p = 1.5a.u and t = 1.8fs and t = 0.85fs. For p = 0.3a.u., a very small number of CWFs ((N c , M) = (16, 1)) is already able to capture the correct dynamics quantitatively. In approaching the target atom with the larger momentum p = 1.5a.u., the conventional mean-field method fails to describe the ionization process due to the lack of electronelectron correlation effects. This is in contrast with dyn-ICWF results, which qualitatively captures the correlated dynamics for a small number of CWFs (N c , M) = (64, 1).

Laser Driven Proton-Coupled Electron Transfer
We discuss here the convergence behaviour of the dyn-ICWF method in capturing the laser driven protoncoupled electron transfer described in Sec. VI B. We suppose the system to be initially seating in the full electron-nuclear ground state, i.e., Ψ(r, R, 0) = Ψ 0 (r, R). This state is prepared using the imaginary time version of the sta-ICWF method with ground state CWFs only (i.e., M = 1). The sta-ICWF provides as output the initial expansion coefficients C(0) and the ground state CWFs φ α i (r i , 0) and χ α J (R J , 0). We then apply an external strong electric field, defined in Sec. VI B of the main text, and the coefficients and the CWFs are propagated using the dyn-ICWF equations of motion.
To achieve converged results, a grid [−9, 9]a.u. with 301 grid points is chosen for the nuclear degrees of freedom. For the electron coordinates, the grid covers the interval [-75,+75]a.u. with 250 grid points. The fourth-order Runge-Kutta algorithm was used to propagate the dyn-ICWF equations of motion with a time-step dt = 0.1a.u, and the Moore-Penrose pseudo-inversion method with a tolerance of 10 −8 was used to approximate the numerical inversion of the overlap matrix in Eq. (31). This very small number of CWFs, even if associated to large deviations across different stochastic particle placements, is able to captured nearly quantitatively both the adiabatic populations and the decoherence indicator of Eq. (C2) and Eq. (C3). This results demonstrate that the dyn-ICWF technique achieves quantitative accuracy for situations in which mean-field theory drastically fails to capture qualitative aspects of the dynamics using three orders of magnitude fewer trajectories than a mean-field simulation.

Berry phase effects and molecular conical intersection
We discuss here some of the technical details of the Berry phase intereference calculation demonstrated in Sec. VI C. As in Ref. [100] we took an electronic spatial grid from -12 to 12 a.u. with 81 grid points and a nuclear grid from -1.5 to 1.5 a.u. with 51 grid points alongside a time step of dt = 0.02a.u. The initial wavefunction was constructed on this grid and the exact dynamics were propagated directly using a fourth-order Runge-Kutta integrator.
The time-resolved picture presents this problem as a fully non-equilibrium problem, where the system starts already in a non-steady state, and so, the imaginary time sta-ICWF cannot be applied here to prepare the initial wavefunction. Instead, we stochastically sample the initial probability density |Ψ 0 (r, R)| 2 with N c trajectories {r α (0), R α (0)} that are used to construct CWFs φ α r (r, 0) and φ α R (R, 0) as defined in Eq. (23). In this process, we respected the symmetry of the underlying initial state by symmeterizing the initial particle placement (and thereby complementarily symmetric slice CWFs) around the R y , r y axes, meaning for each particle R α = (R α x , R α y ) FIG. C.3. BOPESs for the first two excited states with electronic quantum numbers ζ = 1 (lower surface) and ζ = 2.
As mentioned in the main text, the initial nuclear state of is intialized as a gaussian centered at R = (0, 0.4) on the lower surface. we set R α+1 = (−R α x , R α y ). In the dyn-ICWF, the pseudoinverse tolerance for S was set to 10 −8 and the evaluation matrix elements of the electron-nuclear interaction potential term of Eq. (40), W αβ = dRdrφ α * (r)χ α * (R)W en φ β (r)χ β (r), (C4) was accelerated by using an SVD decomposition to break up the 4 index potential W en (r x , r y , R x , R y ) into a sum over electronic and nuclear two index vectors: W en (r x , r y , R x , R y ) = Nσ l=1 σ l u l (r x , r y )v l (R x , R y ). (C5) By tossing out σ l < 10 −4 we found that we were able to retain the accuracy of this potential to within a numerically tolerable limit with a speed up in computation time at a factor between 3.6 and 4.3 depending on hardware. A cubic interpolation to a grid twice as fine was used to smooth the images of the nuclear density In Fig. C.3, we show the first and second excited BOPESs associated to the extended Shin-Metiu model introduced in Sec. VI C. dyn-ICWF results for N c = {1024, 1600, 2400} are shown in Fig. C.4. Due to the finesse of the interference pattern and its fragility with respect to the symmetry of the problem, the number of CWFs required to reproduce quantitatively the exact dynamics is relatively high compared to previous examples in Sec. C 1 and C 2. And yet, note that whilst the N c = 1024 result do not reproduce the interference pattern accurately, they do qualitatively capture the nuclear dynamics by avoiding the forbidden region surrounding the conical intersection. This is in contrast to the meanfield result (figure 7 of Ref. [100]), which fails to capture this qualitative feature of the nuclear dynamics.