Simulating Vibronic Spectra without Born–Oppenheimer Surfaces

We show how linear vibronic spectra in molecular systems can be simulated efficiently using first-principles approaches without relying on the explicit use of multiple Born–Oppenheimer potential energy surfaces. We demonstrate and analyze the performance of mean-field and beyond-mean-field dynamics techniques for the H2 molecule in one dimension, in the later case capturing the vibronic structure quite accurately, including quantum Franck–Condon effects. In a practical application of this methodology we simulate the absorption spectrum of benzene in full dimensionality using time-dependent density functional theory at the multitrajectory Ehrenfest level, finding good qualitative agreement with experiment and significant spectral reweighting compared to commonly used single-trajectory Ehrenfest dynamics. These results form the foundation for nonlinear spectral calculations and show promise for future application in capturing phenomena associated with vibronic coupling in more complex molecular and potentially condensed phase systems.

Simulating vibronic effects from first principles calculations is one of the central goals in theoretical spectroscopy that has implications in chemistry, physics and materials science.
The involvement of nuclear vibrational quantum states during electronic transitions plays a decisive role in determining the spectral features associated with these processes. This has been well established by the utility of the Franck-Condon principle, for example, which represents an early paradigm for the role of nuclear quantum effects in electronically nonadiabatic processes. Describing this interplay between the electronic and vibrational degrees of freedom requires a quantum mechanical description that is both accurate and scalable to relatively large system sizes. One popular method to calculate vibronic spectra is to take a sumover-states approach, where matrix elements of the transition operators between the states involved in the various processes that generate the desired spectral signal are constructed.
In this approach the states of interest are represented using the Born-Oppenheimer (BO) basis; one must already have some a priori knowledge of the BO states that are involved, along with the associated potential energy surfaces and nonadiabatic couplings.
An alternative strategy to the sum over states in the BO basis is to take a coordinate space perspective, and construct the response function for the system of interest from direct time-propagation of the system in that picture. 1,2 This invariably requires some level of approximation in the representation dynamics of the electronic and nuclear degrees of freedom, with different consequences for their coupling depending on the method chosen. The mixed quantum-classical Ehrenfest approach is a practical approximation to the fully quantum mechanical dynamics of the system, and despite it's approximate dynamics, a formally exact representation of the quantum equilibrium structure of the correlated electronic and vibrational degrees of freedom can be included in a multi-trajectory Ehrenfest (MTEF) simulation through the use of the Wigner representation. [3][4][5] In this case, the Wigner transform maps the vibrational quantum states onto phase space distributions of continuous position and momentum coordinates which can be sampled by an appropriate Monte Carlo procedure to capture the quantum equilibrium structure of the problem. The limitations of the Ehrenfest approach, and other independent trajectory semiclassical methods, are well known 6-10 and while there have been many attempts to ameliorate these shortcomings, with some exceptions, 11,12 most rely on the BO framework in their implementation. [13][14][15][16][17] In this work we take a different approach to go beyond mean field theory based on the recently introduced interacting conditional wave function (ICWF) formalism, which is able to capture correlated electronic and nuclear dynamics. [18][19][20][21] We apply MTEF and ICWF dynamics to an exactly solvable one dimensional H 2 model, and show that these methods are able to recover electron-nuclear correlations in linear vibronic spectra without the need to calculate multiple BO surfaces. In addition, we show that the MTEF method can be easily extended to ab initio non-adiabatic molecular dynamics simulations by calculating the vibronic spectra for benzene, where we find good agreement with experimental results.
The linear spectrum of a system is given by the Fourier transform of time correlation function (TCF) C AB (t) = [Â(t),B] of the transition dipole operator,μ, C µµ (t) = μ(t)μ(0) , 1,22 (unless otherwise stated all expressions are in atomic units): where the trace occurs over nuclear and electronic degrees of freedom,ρ eq is the equilibrium density matrix for the coupled system, and we evolveμ(t) in the Hilbert representation.
Traditionally vibronic spectra are explained by invoking the Frank-Condon approximation in the BO picture, where the electronic system is instantly excited, thus promoting the unperturbed ground state nuclear system to a different electronic surface. If one has access to the electronic states involved in a particular spectral range then the contributions to the spectrum due to each electronic transition can be identified by resolving the transition dipole operator in the basis of the electronic states of interest, and the vibronic side peaks of that transition can be calculated by propagating the initial state's nuclear subsystem under the effect of the non-equilibrium electronic occupation. When it is feasible to resolve the nuclear wavefunction dynamics, this can be one of the most accurate methods of calculating molecular vibronic spectra. 23,24 Although resolving eq. (1) in the BO framework is a powerful analysis tool, it is computationally impractical for systems with many nuclear degrees of freedom, particularly when one desires spectra over multiple surfaces. One can bypass this computational bottleneck by representing the system in a real space basis and using the "δ-kick" method, 25 which captures electronic transitions to all dipole-transition allowed states (resolved on the grid) within a single calculation by utilising the dipole response to a perturbative, but impulsive external field H f ield = E(t)μ, i.e. with E(t) = κδ(t) and κ << 1. Using first order perturbation theory, the dipole response ∆µ(t) = µ(t) − µ(0) can be written in powers of the field, 1,2 whereμ I (t) is evolved in the interaction representation. Hence, the linear response spectra may also be obtained via the following relation, provided the strength of the perturbing field κ is sufficiently small. This δ-kick approach only requires the initial state of the full system as input, followed by time propagation for a sufficient interval so as to obtain the desired energy resolution. Importantly, this technique can also serve as a foundation for calculating non-linear optical response spectra. 26 While the methods described above are formally equivalent, difference between the calculated spectra can arise when approximations are made. Here we briefly describe two methods for performing coupled electron nuclear dynamics simulations, the quantum-classical meanfield MTEF method, and the ICWF formalism which was designed to go beyond the mean field limit.
A typical approach to Ehrenfest theory is to assume a separable electronic-nuclear wave function ansatz, take the classical limit of the nuclear portion, and initialise the nuclei at the equilibrium position with zero nuclear momentum. 27,28 This single trajectory Ehrenfest (STEF) method is often employed when a mixed quantum-classical method is needed to couple electronic and nuclear dynamics, 29 in some cases providing a stark difference in electronic dynamics compared to fixed nuclei. 30,31 Although attempts at capturing quantized vibrational effects in STEF with the δ-kick method have been made, 32 they can contain unphysical spectral features (see SI).
An alternative route to Ehrenfest is also possible in the density matrix picture, and proceeds via the quantum-classical Liouville equation. 33 The major difference is that this representation results in a multi-trajectory Ehrenfest picture of the dynamics, where the initial quantum statistics of the correlated system can, in principle, be captured exactly. Here, we outline the evolution equations, and offer more details in the supporting information.
The time evolution of the reduced electronic density is where the subscript W refers to the partial Wigner transform over the nuclei, X = (R, P) is a collective variable for the nuclear position R and momentum P, and the effective electronic mean-field Hamiltonian isĤ Ef f e,W (X(t)) =Ĥ e +Ĥ en,W (X(t)), whereĤ e refers to the electronic portion of the hamiltonian, andĤ en to the electron nuclear coupling. The nuclear dynamics is represented as an ensemble of N independent Wigner phase-space trajectories, ρ n,W (X, t) = 1 N N i δ(X i − X i (t)), that evolve according to Hamilton's equations of motion generated from the effective nuclear mean-field Hamiltonian, ∂R i H Ef f n,W = H n,W (R i (t)) + T r e ρ e (t)Ĥ en,W (X i (t)) .
The average value of any observable, O(t) , can then be written as follows, which can be evaluated by sampling initial conditions fromρ W (X, 0), and evolving the expectation value of the observable with according to the above equations of motion. Using this dynamics method in conjunction with the BO basis representation to evaluate eq.s (1) and (4 -6) ultimately leads to the following equations of motion, with sums over BO states denoted by a, (see SI for details) Where a (R) are the BO surfaces and d aa are the non-adiabatic coupling vectors (NACVs) between states a and a .
In contrast to the previous expression, utilising MTEF in the real space δ-kick approach requires initialising the electronic wave function as the BO eigenstate for each initially sampled nuclear geometry. The δ−kick is applied and the electronic wave function is propagated using the time dependent Schrödinger equation equivalent to eq. (4) alongside the nuclei according to eq. (5). Calculating the spectrum via MTEF dynamics in the BO picture is from here on referred to as MTEF-BO, and calculating it via the δ−kick method is referred to as MTEF-kick.
tical ICWF implementation are recently developed methods which have shown to be able to capture non-equilibrium correlated nuclear-nuclear and electron-nuclear phenomena beyond the mean field limit. [18][19][20][21] This approach is based on taking single-particle slices (the CWFs) of the time-dependent wave function of full system, and approximating the equations of motion for these CWFs by the Hermitian components of the sliced Hamiltonian, and finally, in the ICWF extension, utilising these electron-nuclear CWFs as a basis of Hartree products in a wave function ansatz.
Here we describe an implementation of this approach utilising the static and timedependent variational principles for the expansion coefficients in a static CWF basis. The basis is chosen via sampling electronic and nuclear positions (r α , R α ), α ∈ {1, . . . , N c }, where r and R are understood to be collective position variables, from initial guesses to the electronic and nuclear densities. These are used to construct the Hermitian limit of the CWF propagators 18 for a system with N e electrons and N n nuclear degrees of freedom. Taking eigenstates of h α e (r i ) and h α n (R l ), denoted φ α (r i ) and χ α (R l ) respectively, as our CWF basis we write the following wave function ansatz: where we have taken a Hartree product of electronic and nuclear CWFs for each degree of freedom. While the Hartree product over electronic degrees of freedom has been sufficient for accuracy in applications of ICWF so far, this ansatz can in principle be trivially extended to have fermionic anti-symmetry via inclusion of Slater determinants. We then utilise the Dirac-Frenkel variational procedure [34][35][36] to develop equations of motion for C(t), which leads to the following evolution equation for the expansion coefficients, where In practice S may be nearly singular, but its inverse can be approximated by the Moore-Penrose pseudoinverse. 37 The ground state wave function is obtained from this approach using imaginary time evolution, 38,39 and the δ−kick spectra (ICWF-kick) is calculated by applying the perturbative field to the CWFs at time zero and recalculating the S and H matrices, equivalent to propagating in the interaction representation. This "closed-loop" of initial state preparation and time-propagation ensures that our ICWF approach is a fully self-consistent method that increases in accuracy with increasing N c , and requires no BO state information.
To investigate the performance of the MTEF and ICWF approaches to vibronic spectral lineshapes we studied the vibronic transitions in an exactly solvable one dimensional model system for molecular Hydrogen. [40][41][42] The total Hamiltonian can be written in the center of mass frame in atomic units aŝ where µ n = m p /2 and µ e = 2m p /(2m p + 1) are the reduced nuclear and electronic masses, R is the internuclear separation, and r i are the electronic positions. We take the proton mass to be m p = 1836. The electronic and nuclear degrees of freedom were each resolved on grids for the numerically exact solution and ICWF-kick approaches, while the MTEFkick electronic wave functions were time evolved on the (r 1 , r 2 ) grid, and the MTEF-BO information was calculated by solving the electronic subsystem across the nuclear grid; see the computational methods section for more details. A kick strength of κ = 10 −4 a.u. was sufficient to generate the kick spectra within the linear response regime and, unless otherwise stated, a total propagation time of 10, 000 a.u. ≈ 242f s was used to generate the spectra.  In Fig. 1 we show mean field spectra calculated both with (MTEF-BO), and without (MTEF-kick) the use of multiple BO surfaces for the absorption from S 0 to S 2 in comparison with the numerically exact results. We see that in the BO picture the MTEF method recovers the vibronic absorbtion peak placement quite accurately for the first five peaks, with a broadening occurring for the higher energy peaks that leads to a loss of structure.
This broadening of the spectral signal is due to the well-known fact that the MTEF dynamics does not preserve the correct quantum statistics and thus cannot fully capture the electronnuclear correlation in the problem (see the SI for a detailed discussion of this issue). The pre-peak features in Fig. 1b are also unphysical artefacts of MTEF. The MTEF-BO spectra were converged to within graphical accuracy using N = 50, 000 trajectories.
Focusing on the MTEF-kick results in Fig. 1c, we see that this approach recovers vibronic side peak structures again without any BO surface information, albeit with inaccurate spacing, while STEF-kick captures only the vertical electronic transition from the minimum of the S 0 surface. The MTEF-kick spectra converged to within graphical accuracy using N = 30, 000 trajectories. The average peak spacing in the MTEF-kick spectra is approximately 0.32eV; this corresponds remarkably well with the natural frequency of the harmonic approximation to the ground state surface expanded around the equilibrium geometry, which is also 0.32eV in this case. This result is unsurprising as the electronic kick induces a very small population transfer to the upper surface proportional to the square of the kick strength, which results in the mean forces on the nuclei in MTEF-kick essentially corresponding to those of the initial state. The influence of initial state on the MTEF-kick spectra is further demonstrated by analysing the emission spectra in Fig. 2. The initial state here was chosen by hand as the lowest lying nuclear state on the S 2 surface. Once again we see that MTEF-BO recovers the peak placement quite well, while the MTEF-kick data has a less accurate vibronic spacing. Fitting the MTEF-kick peaks, we find an excellent correspondence between mean spacing of the five lowest energy MTEF peaks and the excited surface natural frequency of 0.21eV.
For ICWF-kick, we found that N c = 4096 and mixing the three lowest energy CWF eigenstates was sufficient to obtain quite accurate results. In Fig. 3 we demonstrate that the    ICWF ansatz used in a variational context achieves a much more accurate vibronic spacing than the MTEF-kick approach, without the failing of peak broadening or unphysical spectral negativity apparent in the MTEF-BO results. The accuracy of these results underscores that the ICWF ansatz is a robust framework to capture the electronic and vibronic quantum dynamics, being accurate for not only the electron-nuclear correlation inherent to vibronic spectra, but also the electronic subsystem itself, which in the MTEF results was solved exactly either on a grid or using explicit BO state information. The deviation from the exact results does grow with increasing energy, although this is ameliorated with increasing N c , and can in principle be eliminated at large enough values of N c (see SI).
Finally we demonstrate the application of MTEF-kick to real 3D molecular systems using the ab initio Octopus 43 real-space time dependent density functional theory (TDDFT) 44 package to calculate the linear vibronic MTEF-kick spectra of Benzene. The initial conditions for the nuclear subsystem were obtained by calculating the normal mode frequencies and dynamical matrix of the molecule, and sampling Wigner transforms of the ground state wave functions in the harmonic approximation; see SI for more details. The adiabatic-LDA functional was used, along with norm-conserving Troullier-Martins pseudo-potentials, and the trajectories were evolved for 201 eV ≈ 132f s. A kick strength of κ = 5x10 −3 a.u. was used to generate the kick spectra within the linear response regime in this case, and the graphical convergence of the MTEF results was found to be achieved with N = 500 trajectories.
In Fig. 4, we compare the MTEF-TDDFT-kick results to its STEF-TDDFT-kick counterpart, each scaled to match the peak intensity of an experimental data set for the optical absorption of benzene digitized from Ref. 45 We see that there is remarkably good agreement across the wide energy range available from experiment, before molecular dissociation pathways become available around 13.8eV. Again, the full spectrum is resolved without resorting to calculations of individual transitions as would be required in a BO state calculation. The three STEF peaks in the 7eV region correspond to the energy range of the doubly degenerate, dipole allowed 1 E 1u ← 1 A 1g , π * ← π transition, 31,45,46 with the energy degeneracy  artificially lifted by the discrete grid. The experimental band preceding the central peak, in the range 6eV to 6.5eV is commonly ascribed to the dipole forbidden, but vibronically [45][46][47] and in the MTEF-TDDFT-kick results we see a low energy tail extending through the 5eV-6.25eV range, well away from the STEF results, eventually transitioning into the broad peak centered around 7eV. It's reasonable to expect that the broadening of the MTEF signal relative to the experimental signal is due to the effects discussed above that arise due to the mean field treatment.
We have demonstrated that semi-classical MTEF simulations can capture vibronic structure with the correct spectral sign in the region of the transition. Moreover, we have shown how this can be achieved without using multiple BO surfaces via the δ-kick method, and that the vibronic spacing predicted with the MTEF-kick approach matches the profile of the initial state. We have addressed these shortcomings by combining the ICWF formalism with the δ-kick method, which provides more accurate vibronic spectra in a computationally efficient and systematically improvable fashion. Finally, we demonstrated that MTEF-kick is easily applied to ab initio molecular systems by simulating the vibronic spectra of benzene and finding good agreement to experimental results.
These linear response results establish a solid basis for further investigations into nonlinear response of field driven molecular systems utilising the practical and efficient MTEF and ICWF techniques along with ab initio electronic structure methods. Work in preparation by the present authors also explores the utility of ICWF with electron-electron and electronnuclear correlated systems, and explores the response of these systems under nonperturbative electric fields. Furthermore we expect that MTEF-kick will improve in accuracy for periodic systems, as changes in the electronic configuration are often to likely produce smaller changes in the nuclear forces than in molecular hydrogen. This makes this method interesting to pursue in periodic systems in particular, where there is a dearth of theoretical frameworks for ab initio, nonpertubrative electron-nuclear coupling. 48 Work in this direction is in progress, as is the implementation of the ICWF method within an ab initio framework for molecular and periodic systems.

Computational Methods
In the 1D H 2 model, the electronic coordinates are each resolved on a 65a 0 wide interval with spacing 0.6a 0 , while the nuclear grid extends to R max = 6.3125a 0 with 0.0625a 0 spacing. For the 1D H 2 MTEF-BO results, the potential energy surfaces a (R) and µ aa W (R) were calculated on a nuclear grid with ∆R = 0.02a 0 up to R max = 8a 0 , fit to a cubic spline function, and interpolated every 0.01∆R. 50 The NACV between S 0 and S 2 in this model is numerically zero. These quantities were resolved for the first allowed dipole transition, between the ground state S 0 and the second excited state S 2 , and the results were found to be well converged within about 5 × 10 4 trajectories.

MTEF Equations of Motion
Starting from a density matrix representation of the full system,ρ, we Wigner transform over the nuclear subsystem, producing a unique mapping onto a nuclear position R and momentum P phase space X = (R, P), where R and P are collective variables R = (R 1 , . . . , R Nn ), P = (P 1 , . . . , P Nn ), with R i , P i ∈ R d . The partial wigner transform is defined for any operator aŝ leaving a Hilbert space operator character over the electronic degrees of freedom, dependent on the continuous nuclear phase space parameters. In general, developing equations of motion forρ W (R, P), (or any operator), requires taking the partial Wigner transformation of the Liouville von-Neumann equation of motion for ρ: Where the final line defines the "Moyal product" also known as the "star product". 51 where {A(R, P), B(R, P)} refers to the normal Poisson bracket.
To derive MTEF equations of motion from the QCLE, one takes the mean field approximation by assuming that the full system can be written as a sum of correlated and uncorrelated parts,ρ and then neglecting the contribution of the correlated part in the dynamics. Note that while the ensuing dynamics do not explicitly treat the effect of subsystem correlation, the initial state generally is correlated, and therefore is implicitly included in the dynamics.
In the equations of motion resulting from inserting this approximation into the QCLE, the evolution of the reduced Wigner density of the nuclear subsystem can be exactly represented, via the method of characteristics, by a sufficiently large ensemble of multiple independent trajectories, ρ n,W (X, t) = 1 N N i δ(X i − X(t)). Each trajectory evolves according to Hamil-ton's equations of motion generated from the mean-field effective Hamiltonian, Where H n,W and H en,W refer to the partially Wigner transformed nuclear and electronnuclear coupling operators, respectively. The electronic density associated with each trajectory , ρ i e (t), evolves according to the following commutator: The exact expression for the average value of any observable, O(t) , can be written as The mean field limit of this expression simple corresponds to evaluating the integral by sampling initial conditions for an ensemble of independent trajectories fromρ W (X, 0), and then generating the time evolution for each trajectory by approximatingÔ W (X, t) by it's mean-field counterpart.
Following the sampling of an initial nuclear condition, X i , from the Wigner distribution associated to the nuclear subsystem wave function, the electronic system is initialised as: i.e. implicitly as the BO electronic state at R i . Under this scheme, the electronic subsystem's initial conditions are implicitly correlated with the nuclear subsystem's quantum statistics.
In cases where the nuclear initial state is impractical to calculate exactly one may utilise the normal modes of the molecular system, or phonon coordinates of a periodic system, to treat the full nuclear wavefunction as a Hartree product of N uncoupled harmonic oscillators, where N is the number of non-rotational and non-translational nuclear degrees of freedom: With c Defining the dynamical matrix, such that we obtain the nuclear Hamiltonian in dimensionless normal mode coodinates: Of course, the simple harmonic wave function solutions to the above Hamiltonian have well known analytical expressions and are trivially Wigner transformed, the ground state harmonic oscillator wavefunction's Wigner function for instance is: 4 We can therefore sample these transforms for (Q, S) and then use eq. (23) to back transform to from normal mode coordinates to cartesian coordinates.

MTEF-BO Equations of Motion in the Born Oppenheimer Basis
In We then utilise a position representation in the nuclear dof by expanding in the space of nuclear position states 1 R = dR |R R|, leading tô For a transition between two electronic states g and e, we can expand in the adiabatic basis |φ a (R) , (a = g, e) which are dependent on the nuclear positions R defined by, Taking the partial Wigner transform of eq. (27) leads tô whereĥ e,W (R) is the normal electronic hamiltonian operator, now dependent on R in the Wigner nuclear phase space. Starting with the separability approximation for the density operator, and neglecting correlations, we haveρ W =ρ e ρ n (R, P), with where T r X . . . = . . . dRdP, and P scalar terms are cancelled by the commutator. We are of course interested in evaluating the dipole-dipole correlation function: whereσ = [μ W ,ρ W ], and we resolve the dipole operator aŝ Where Z R refers to the ionic charge of each nuclei. In practice we can neglect the intra-state R term as we are focused entirely on the transition dipole moment.
Taking the initial state as the ground state, (|Ψ = |χ 0 g φ g ) ρ W (R, P, 0) = ρ n g (R, P) leads toσ And therefore the correlation function becomes We can construct an identical quantity from a different initial condition as a superposition state (|Ψ = 1 √ 2 |χ g (|φ g + i |φ e )) giving, ρ W (R, P, 0) = ρ n g (R, P) For this different initial condition we propagatẽ With this different initial condition, we take the MTEF form of the nuclear density arising from the Monte Carlo integration described above, The subsequent equations of motion for the system are for the electronic density, needed for the nuclear trajectories are: Where in the last two equations we have used the identity d i Note that for transitions like the S 0 /S 2 transition 1D H 2 focused on in the body of this paper, the non-adiabatic coupling vector (NACV) d aa = 0, means that the mean field force acting on the nuclei is at all times a 1 2 superposition of the S 0 and S 2 surfaces. These are propagated alongside the dipole matrix element equations of motion, needed for the correlation function:

Application to Displaced Harmonic Oscillator Model
In order to investigate the limitations of MTEF, we can utilise a model which captures the essential physics of the S 0 /S 2 1D H 2 transition which was focused on in the first portion of the main text. Recall that for this transition, the NACV's between the two electronic adiabatic states are zero, that is φ a (R)|∂ R φ a (R) = 0 ∀ a, a in the BO basis, with a, a restricted to S 0 /S 2 This means that matrix elements for the partially Wigner transformed molecular hamiltonian can be written aŝ As described in detail in the first section of this SI, MTEF is rooted in a mean field approximation to the QCLE, which is itself the first order expansion of the partially Wigner transformed Liouville von-Neumann equation. Taking eq. (13) to second order provides, Which in our model Hamiltonian eq. (41) becomes, Such that the error in time propagation resultant from taking only the first order expansion, compared to the second, is proportional to the difference in energy surface curvature.
If we take the analytically solvable Displaced Harmonic Oscillator (DHO) model 22,52 by using surfaces a (R) = 1 2 ω 2 a (R − D a ) 2 + E a , we see that for identical surfaces ω e = ω g that the 2 nd order and higher terms in the Wigner transformed Liouville von-Neumann equation are zero, rendering the QCLE exact for this case.
To demonstrate the effect of varying surface curvature, we took parameters similar to harmonic surface fits to the BO surfaces in 1D H 2 , and for simplicity, took the FC approximation alongside setting µ aa (R) = (1 − δ aa )a.u.. We solve the exact and MTEF-TCF spectra for the DHO with different values of ω e relative to ω g by propagating for T f = 2 · 10 4 a.u..
In Fig. (7) we see iin the left column that for identical upper and lower surfaces, mean field theory is of course exact, and for varying surfaces, MTEF displays a peak broadening and prepeak features. The origin of this broadening is from an effective damping in the time dependent signal, shown in Fig. (8).  The conditional wave function (CWF) approach can be developed starting from the full molecular wave function for electrons and nuclei, Ψ(r, R, t), which can be formally decomposed in terms of the CWFs of each subsystem: From these definitions one can show that the CWFs, ψ α e (t) and ψ α n (t), obey non-Hermitian equations of motion involving complex potentials which are functionals of the full wave function and cause the time-evolution of the individual CWFs to be non-unitary. 18 The recently developed Interacting-CWF (ICWF) method 21 avoids the direct calculation of these nonlocal complex potentials by positing the following multiconfigurational CWF basis ansatz for the full many-body wave function: The basis functions in this sum are chosen to be single particle CWFs that satisfy the mean-field, or Hermitian, limit of the CWF equations in which the complex potentials trivially vanish. The upper limit of the sum, N c , refers to the total number of configurations, which can be stochastically sampled. Including interactions between the trajectories in the ensemble through the coefficients C(t) = {C 1 (t), ..., C Nc (t)} corrects the Hermitian-CWF evolution. The time evolution of these coefficients is obtained by inserting eq. (46) directly into the TDSE.
As described in the text, for the kick spectra adapted ICWF algorithm, the CWFs are instead selected as eigenstates of the Hermitian propagators, and used as a static basis.
The imaginary and real time equations of motion for the expansion coefficient C are then solved using the respective variational principles, [34][35][36]39 allowing for a completely closed-loop algorithm for wave function preparation and propagation.
To generate the kick spectra, after preparing the ground state C(0), the relevant degree of freedom of the kick operator exp(−iκμ) is applied to each CWF, the Hamiltonian and inverse overlap matrices are reconstructed, and C is propagated to the desired time. This procedure is equivalent to propagating in the interaction representation, withV I (t) = κδ(t)μ.
Since these matrices are only constructed at time zero, this algorithm is extremely efficient, With increasing non-redundant variational parameters, one is guaranteed to better capture the initial state and minimize the error of time dependent propagation. 35 As an example of the convergence properties of ICWF-kick, see where Θ is the Heaviside function, and η was set to 0.1Ha/a 0 for both subsystems.
The electronic CAP cut offs, r l and r r , were placed 10a 0 from the walls, while the nuclear CAP start was set at R 0 = 5.6875a 0 .
Physical Review A 2016,