Conditional Wave Function Theory: A Unified Treatment of Molecular Structure and Nonadiabatic Dynamics.

We demonstrate that a conditional wave function 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 wave function formally recasts the full interacting wave function of a closed system as a set of lower-dimensional (conditional) coupled "slices". We formulate a variational wave function ansatz based on a set of conditional wave function 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 wave functions and address paradigmatic nonequilibrium processes including strong-field molecular ionization, laser-driven proton transfer, and nuclear quantum effects induced by a conical intersection. This work paves the road for the application of conditional wave function theory in equilibrium and out-of-equilibrium ab initio molecular simulations of finite and extended systems.

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 not only be designed to monitor and/or control the structure of molecules[1−7] and solids[8−12] but also form light–matter hybrid states with new physical properties.[13−21] In view of these exciting developments, accurate first-principles theoretical techniques are also needed 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 unfavorable 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 so-called 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−27] and nuclear quantum effects,[28−31] as well as of systems driven out of equilibrium.[32−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 wave function 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−40] Of these approaches, perhaps the most popular are the 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, using a variety of wave function ansatz,[43−48] semiclassical techniques,[49,50] the quantum-classical Liouville equation,[51−53] path-integral methods,[54,55] or methods based on the exact factorization,[56−58] allows for further accuracy with increased computational effort.

While advances in the 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,43] as the memory required to store the information contained in the BOPESs grows rapidly with the number of correlated degrees of 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 favorable 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 wave function (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 wave function 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 wave function approach (ICWF),[73,74] a recently introduced method for performing quantum dynamics simulations that is multiconfigurational 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 time-independent CWF framework is formulated in real space, and it is an exact decomposition of the time-independent wave function of a closed quantum system that yields a set 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 structures of correlated electron–nuclear systems and eventually extended systems. Importantly, the combination of the static version of the ICWF method using a time-dependent 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 Section , we define the mathematical structure of the time-independent version of the CWF framework. Based on these results, we put forth an imaginary-time version of the ICWF technique in Section 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 (1D). In Section , 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 H2 model system. In Section , we revisit the exact time-dependent CWF framework, and in Section , we present the dynamical ICWF (dyn-ICWF) approach to the time-dependent 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 nuclear quantum effects in the dynamics through a conical intersection (CI). A summary of the main results of this work and an outlook on future directions are offered in Section .

Conditional Eigenstates

We begin by considering a closed system with n electrons and N nuclei, collectively denoted by = (, ). We use the position representation for both subsets; lowercase symbols will be used for the electronic subsystem, e.g., = {1s1, ..., s}, and uppercase symbols = {1σ1, ..., σ} for the nuclear subsystem. Hereafter, electronic and nuclear spin indices, respectively, s and σ, will be made implicit for notational simplicity, and, unless otherwise stated, all expressions will be given in atomic units.

The time-independent CWF can be constructed starting from the nonrelativistic time-independent Schrödinger equation in the position representationwhere Ψγ() is an eigenstate of the molecular Hamiltonian Ĥ with label γ and the corresponding energy eigenvalue Eγ. The molecular Hamiltonian operator Ĥ in eq can be written aswhere the kinetic energy operators are , m and z being the characteristic mass and charge of particle j, respectively. The full electron–nuclear potential energy of the system is W() (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 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 groups. 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 two-particle 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.a This is, in general, not a trivial task, and hereafter, we will assume that eq 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 Ψγ() in terms of single-particle conditional eigenstates of either of the two subsystems, which are defined as followsHere, the index α denotes the particular conditional slice and = (1, ..., , , ..., ) are the coordinates of all degrees of freedom in the system except . Similarly, α = (1α, ..., α, α, ..., α) are some particular positions of all system degrees of freedom except . As shown schematically in Figure , the conditional eigenstates in eq represent one-body slices of the full many-body eigenstates Ψγ() taken along the coordinate of the ith degree of freedom. The particle placement α 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 important sampling techniques.

Figure 1

Schematic representation of the CWF approach to the time-independent Schrödinger equation for one electron and one nucleus in one dimension, i.e., = (r, R). (a) The full ground-state Ψ0(r, R) is plotted together with a pair of conditional ground states ϕα,0(r) for the electronic degree of freedom (in red) and χα,0(R) for the nuclear degree of freedom (in blue) for a given position of the full configuration space {rα, Rα}. Contour plots of the molecular wave function are also shown for clarity. (b) The exact solution of the time-independent Schrödinger equation in eq can be reconstructed provided a sufficiently large ensemble of sampling points xα = {rα, Rα}. This can be done by applying the reassembling transformation or (whose definition can be found in Appendix A) to the ensemble of electronic ϕα,0(r) or nuclear χα,0(R) conditional eigenstates, respectively.

Evaluating eq at α by applying the integral operator in eq yields conditional eigenstates that are the solutions of the following eigenvalue problemwhere we introduced Wα() = W(,α), with W() being the full electron–nuclear interaction potential appearing in the Hamiltonian of eq . In addition, ηα,γ() are the kinetic correlation potentials given by

Provided a large enough collection of CWFs satisfying eq , an exact solution of eq can be reconstructed by undoing the conditional decomposition of eq (see Figure b).[69] That is, given a set of conditional slices that sufficiently span the support of Ψγ, then the corresponding conditional eigenstates can be used to reassemble the full electron–nuclear wave functionusing the transformations , which are discussed in more detail in Appendix A. This expression, eq , can be used to evaluate the kinetic correlation potentials in eq . In this way, the generalized one-body eigenvalue problem in eq can be understood as an exact decomposition and recasting of the eigensolution of the full electron–nuclear system, which yields a set of coupled, non-Hermitian, eigenvalue problems.

Time-Independent Hermitian Approximation

An approximate solution to eq can be formulated by expanding the kinetic correlation potentials around the sampling coordinates α 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 Ψγ() = ∏ψ() into eq . By making this approximation, the eigenvalue problems in eq are restored to a Hermitian formThe 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.

Static Properties with Conditional Eigenstates

In general, the higher-order terms in the Taylor expansion of the kinetic correlation potentials are non-negligible. However, one can still take advantage of the simple Hermitian form of the conditional eigenvectors (hereafter referred to as conditional wave functions (CWFs)) in eq to design an efficient many-body eigensolver by utilizing them as bases for electronic and nuclear degrees of freedom in a variational wave function ansatz.

While there is a diverse literature spanning decades on different forms for variational electron–nuclear wave function ansatz, for illustrative (and practical) purposes, we employ a sum-of-product form, which in the language of tensor decompositions is referred to as the canonical format.[77] For each degree of freedom , we utilize a given electronic or nuclear CWF, respectively, coming from solutions to eq , to approximate the γth full system exact excited state as followswhere in the second line, we have rearranged the sum over particle position λ ∈ {1, ..., N} and excited CWF ν ∈ {1, ..., M} into a single index α = λ + N(ν – 1), such that α ∈ {1, ..., NM}. The particle placement α defining the conditional potentials Wα 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 ) as being in canonical format because we do not mix all possible CWFs ψλ,ν for all possible degrees of freedom , as one does with a single-particle function bases across the different system degrees of freedom in the Tucker format employed in the multiconfigurational time-dependent Hartree (Fock)—MCTDH(F)[43] and multiconfigurational electron–nuclear dynamics ansatz.[78] 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 hierarchical Tucker formats, which when employed in the multilayer extension[79,80] of MCTDH allow for an increase in efficiency for certain problems. However, since the time dependence of the ansatz in eq 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 antisymmetry 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 the static-basis ICWF (or sta-ICWF) ansatz. With this ansatz in hand, we then consider a solution of eq based on the imaginary-time propagation technique,[81] i.e.whereand P̂ζ = ΨζΨζ† are projectors used to remove the wave functions Ψζ 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.

By introducing the ICWF ansatz of eq into eq , we find an equation of motion for the coefficients Cγ = {C1γ, ..., Cγ}where , and the matrix elements of and arewhere 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, may be nearly singular, but its inverse can be approximated by the Moore–Penrose pseudo-inverse.

Based on solving the system of equations in eq 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 Ô(). Given an approximate solution to the eigenfunction Ψγ(), the expectation value of readswith the matrix elements of being given by an analogous expression to eq .

Example I: Ground and Excited BOPESs of H2

As an illustrative example, we now calculate the BOPESs of a model for the H2 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.[82,83] In this model, the relevant coordinates are the internuclear separation, R, and the electronic coordinates, r1 and r2. The Hamiltonian, written in terms of these coordinates, iswhere for M being the proton mass, μe = M/(2M + 1) is the reduced electronic mass and μn = M/2 is the reduced nuclear mass. In eq , the electron–electron repulsion and the electron–nuclear interaction are represented by soft-Coulomb potentialsi.e., the Coulomb singularities are removed by introducing smoothing parameters ϵ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.[84−86] Moreover, it has provided valuable information in the investigation of electron correlation effects.[87−89]

For this model, the BOPESs are defined by the following electronic eigenvalue problemwhere , and {Φγ(r1, r2; 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 –14 along with a simplified version of the ansatz in eq that is specialized to this particular case of parametric nuclear dependence. A thorough description of the numerical procedure, as well as the convergence behavior of the sta-ICWF method for this model can be found in Appendix B.1.

In Figure , we show the first five BOPESs calculated via the sta-ICWF approach using (N, 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 the excited BOPES landscape over this energy range. As a point of comparison, in the bottom panel of Figure , 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.

Figure 2

Exact first five BOPESs of the one-dimensional H2 model system (solid black lines). sta-ICWF results for (N, M) = (32, 5) are shown in the top panel (solid gold lines). Hartree–Fock and CIS results for the ground-state and excited-state BOPESs, respectively, are shown in the bottom panel (dashed lines) alongside exact results (solid lines) and color-coordinated via calculated excited states.

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 equationwhere Ψ(, t) is the electron–nuclear time-dependent wave function, 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 wave function is the correlated electron–nuclear ground state, i.e., Ψ(, 0) = Ψγ=0(), and some nonequilibrium dynamics is triggered by the action of an external driving field (hereafter, we omit the superscript γ for clarity). We can then decompose the time-dependent many-body wave function as in eq by restricting it to the case of γ = 0. We choose to restrict, for the moment, the time dependence of our ansatz to the expansion coefficients Cα. Although in this formulation the basis remains static, by choosing sufficient excited CWF states, γ > 0 in eq , for α covering some anticipated range of motion for the dynamics, we can expect to capture the support of Ψ(t). The equations of motion for Cα can be obtained either by inserting eq directly into eq or by utilizing the Dirac–Frenkel variational procedurebIn eq , the matrix elements of and are identical to those defined in eqs and 14, with the Hamiltonian’s time dependence coming from any external fields and the wave function decomposed into single-particle CWFs for the nuclear and both electronic degrees of freedom. 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 . In this way, the combination of the imaginary-time and real-time sta-ICWF methods yields a “closed-loop” algorithm for the structure and dynamics of molecular systems that does not require explicit BO state information as an input to the method. For the interested reader, a detailed flowchart of the resulting sta-ICWF method can be found in Appendix D.

Example II: Optical Absorption Spectrum of H2

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 Section . We utilize the “δ-kick” method of Yabana and Bertsch,[90] where an instantaneous electric field E(t) = κδ(t) with perturbative strength κ ≪ 1 au–1 couples to the dipole moment operator μ = r1 + r2 and thereby produces an instantaneous excitation of the electronic system to all transition dipole allowed states. The resulting (linear) absorption spectra can then be 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 that smoothly vanishes at the final simulation time Tf.

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 of the H2 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 behavior of the sta-ICWF method for this model, can be found in Appendix B.3. The reader can also find a detailed flowchart of the (real and imaginary) sta-ICWF method in Appendix D.

For the H2 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 Figure , we show vibronic spectra calculated both with sta-ICWF for the absorption from S0 to S2 in comparison with the numerically exact results also calculated via the δ-kick approach. For sta-ICWF, we found that N = 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 furthermore, it not only captures 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 and M, and can, in principle, be eliminated at large enough values of these parameters (see Appendix B.3).

Figure 3

S2 ← S0 spectra of ICWF-kick (gold) and multitrajectory Ehrenfest δ-kick (MTEF-kick) (blue) compared to the exact peak placement overlaid as a black line, showing that while mean-field theory is unable to capture qualitatively the correct vibronic line shape spacing and intensity, the sta-ICWF approach accurately captures the exact spectrum.

For comparison, we also show mean-field, semiclassical results for the vibronic spectra. Specifically, we calculated the absorption spectrum with the multitrajectory Ehrenfest δ-kick (MTEF-kick) method,[61] overlaid as dashed blue lines. The electronic subsystem was solved exactly as a two-particle wave function over the real-space grid for each independent nuclear trajectory. 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.

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, ultrafast laser pulse for this same H2 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) = E0Ω(t) sin(ωt), with E0 = 0.005 au 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 wave function. A thorough description of the numerical procedure, as well as the convergence behavior of the sta-ICWF method for this model, can be found in Appendix B.4, as well as in Appendix D.

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 Figure . 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 Figure . In the MTEF mean-field 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, M) = (4096, 3) was sufficient to obtain highly accurate results for both the expectation value of the electronic dipole moment (top panel of Figure ) and the expectation value of the internuclear separation (bottom panel of Figure ). Further details can be found in Appendix B.4.

Figure 4

Top panel: evolution of the expectation value of the dipole operator ⟨μe⟩ for the 1D H2 model system for N = 4096 (from bottom-up) and M = 3. Bottom panel: evolution of the expectation value of the nuclear interseparation ⟨R⟩ for the 1D H2 model system for N = 4096 and M = 3.

Time-Dependent Conditional Wave Functions

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 wave function, 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 wave functions in the ICWF ansatz. Hence, we take advantage of the time-dependent version of the CWF framework introduced in ref (69), which relies on decomposing the exact many-body wave function, Ψ(, 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 at α(t), one can show that the CWFs in eq obey the following equations of motionwhere Wα(, t) = W(,α(t), t), and we remind that W() is the full electron–nuclear interaction potentials that appear in the Hamiltonian of eq . In eq , ηα(, t) are time-dependent complex potentials containing kinetic correlations and advective terms, i.e.

As in the time-independent CWF framework, the conditional wave functions in eq represent slices of the full wave function taken along single-particle degrees of freedom of the two disjoint subsets. Each individual CWF constitutes an open quantum system, whose time evolution is nonunitary, due to the complex potentials ηα(, t), which now include advective terms due to the inherent motion of the trajectories α(t), which evolve according to Bohmian (conditional) velocity fields[69]An exact solution to eq can be then constructed provided we use a sufficiently large number of slices {α(t)} that explore the full support of |Ψγ(, t)|2 (in analogy with Figure b), i.e.where the transformations can be found in Appendix A. The one-body equations of motion in eq can be then understood each as a coupled set of nonunitary 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 comoving frame in which all coordinates except the ith move attached to the electronic and nuclear flows and only the ith coordinate is kept in the original inertial frame. Within the new coordinates, the convective motion of all degrees of freedom except for the ith coordinate is described by a set of trajectories of infinitesimal fluid elements (Lagrangian trajectories), while the motion of the ith degree of freedom is determined by the evolution of the CWFs in a 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 localized where the full molecular wave function has a significant amplitude.

Time-Dependent Hermitian Approximation

In general, the effective potentials in eq exhibit discontinuous steps, which could introduce instabilities in a trajectory-based solution of the many-body dynamics based on eq . 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 α(t), such thatIn 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 is approximated aswhile the trajectories α(t) are constructed according to eq .

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 single-particle equations of motion. Despite the crudeness of the approximation in eq , the set of equations of motion in eq has found numerous applications, e.g., in the description of adiabatic and nonadiabatic quantum molecular dynamics[69,71] and quantum electron transport.[91−95] In ref (69), for example, results using eq for an exactly solvable model system showed a great degree of accuracy of the time-dependent Hermitian approximation in capturing nonadiabatic dynamics. Alternatively, in ref (71), the set of equations in eq was used to describe the adiabatic double proton transfer for an exactly solvable model porphine, showing great promise in capturing quantum nuclear effects. Regarding the comparison of the time-dependent Hermitian approach in eq with conventional mean-field methods, in ref (92), it was shown that quantum electron transport simulations using eq represent an improvement with respect to time-dependent (Hartree-type) mean-field simulations. Similar conclusions were reported in ref (96), where a simplified semiclassical method based on eq was compared with classical mean-field results.

Methods based on eq , however, are known to fail to describe important nonadiabatic processes such as the splitting of the time-dependent reduced nuclear density with influences from different BOPESs.[69] This type of dynamics has been commonly associated with decoherence effects that neither the Hermitian approximation in eq nor other mean-field methods such as Ehrenfest or Tully’s surface hopping dynamics are able to capture.

Simulating Far-from-Equilibrium Dynamics with Conditional Wave Functions

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 to design an efficient many-body wave function propagator. For that, we expand the full electron–nuclear wave function using the ansatzwhere the coefficients Cα(t) and the CWFs ψα(, t) are initialized using the sta-ICWF method and propagated afterward using the approximated equations of motion in eq along with trajectories obeying eq .

The time evolution of the coefficients C(t) can be then obtained by inserting the ansatz of eq into eq where the matrix elements of , are defined as in eqs and 14, with the time dependence coming from external fields in the Hamiltonian and the time-dependent CWFs, while arewhere hα(t) are the Hermitian Hamiltonians in eq and Ĥ(t) is the full time-dependent Hamiltonian in eq .

Obtaining these matrix elements is straightforward, involving a sum across single-body operators in eqs and 32 and all sums of two-body interactions across each degree of freedom in eq . Note that any operator involving only a single species, e.g., the kinetic energy, is canceled out, and thus the evolution of C is governed exclusively by matrix elements of operators, which either fully (through ) or conditionally (through ) correlate the degrees of freedom.

Equations , 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 ⟨Ô()⟩ as given in eqs with dyn-ICWF by simply taking into account that ψα(t) are now time-dependent CWFs.

The above dyn-ICWF method was first put forth in ref (73). At the time of publishing the work in ref (73), however, there was no theory sustaining the construction of the initial conditional wave function basis ψα(,t) without relying on an exact solution of the time-independent Schrödinger equation. That has been the main limitation of the method thus far. Here, instead, we have shown that the imaginary-time sta-ICWF method (derived in Section ) not only allows us to solve accurately the time-independent Schrödinger equation but also serves as a method to define an optimal set of conditional wave function basis ψα(,0). Therefore, the dyn-ICWF in combination with imaginary-time sta-ICWF provides a self-consistent approach to describe observables that are relevant to equilibrium, as well as far-from-equilibrium processes. An example combining these two methods will be shown in the example of Section , where an initial ground state is prepared using imaginary-time sta-ICWF and a later dynamics, triggered by a laser pulse, is described using dyn-ICWF. The interested reader can find a complete flowchart of the combined method in Appendix D.

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.[97,98] We here study a model system of electron–hydrogen scattering that can be exactly solved numerically.[99] In atomic units, the Hamiltonian of this one-dimensional two-electron model system readswhereare, respectively, the soft-Coulomb interaction and the external potential that models the H atom located at r = 10 au. The initial interacting wave function is taken to be a spin singlet, with a spatial partwhere ϕH(r) is the ground-state hydrogen wave function and ϕWP(r) is an incident Gaussian wavepacketwith α = 0.1 representing an electron at r = −10 au, approaching the target atom with a momentum p.

The time-resolved picture presents scattering as a fully nonequilibrium problem, where the system starts already in a nonsteady state, and so, the imaginary-time sta-ICWF cannot be applied here to prepare the initial wave function. Instead, we stochastically sample the initial probability density |Ψ0(r1, r2)|2 with N trajectories {r1α(0), r2α(0)} that are used to construct CWFs ϕ1α(r1, 0) and ϕ2α(r2, 0), as defined in eq . A thorough description of the numerical procedure, as well as the convergence behavior of the dyn-ICWF method for this model can be found in Appendix C.1. See also Appendix D for a description of the corresponding workflow.

We study the dynamics of the electron–hydrogen scattering by evaluating the time-dependent one-body density, ρe(r1, t) = 2∫|Ψ(r1, r2, t)|2  dr2, for two different initial momenta, viz., p = 0.3 and 1.5 au. For p = 0.3 au, the energy is lower than the lowest excitation of the target (which is about ω = 0.4 au) 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 = 128 results both capture the correct dynamics accurately (see Figure ). In approaching the target atom with the larger momentum p = 1.5 au, 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 = 256, although a lower number of CWFs already reproduces qualitatively the dynamics (see Appendix C.1 and Figure ).

Figure 5

Top panel: reduced electron density at t = 1.8 fs for p = 0.3 au and N = 128. Bottom panel: reduced electron density at t = 0.85 fs for p = 1.5 au and N = 256 and Nin = 10.

Example V: Laser-Driven Proton-Coupled Electron Transfer

We now show dyn-ICWF results for a prototypical photoinduced proton-coupled electron transfer reaction, using the Shin–Metiu model.[100] The system comprises donor and acceptor ions, which are fixed at a distance L = 19.0a0, 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 iswhere 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 readswhere erf() is the error function. The parameter regime studied for this model (R = 5a0, R = 4a0, and R = 3.1a0) is chosen such that the ground-state BOPES, ϵBO1, is strongly coupled to the first excited adiabatic state, ϵBO2, around the mean nuclear equilibrium position Req = −2a0. The coupling to the rest of the BOPESs is negligible.

We set the system to be initially in the full electron–nuclear ground state obtained from the imaginary-time propagation method described above, i.e., Ψ(r, R, 0) = Ψ0(r, R) (the interested reader can find a general workflow of the simulation in Appendix D). We then apply an external strong electric field, E(t) = E0Ω(t) sin(ω t), with E0 = 0.006 au, Ω(t) = sin(πt/20)2, and ω = ϵBO1(Req) – ϵBO0(Req). 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 Figure , where we monitor the BO electronic state populations P(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 Figure , where, as a measure of decoherence, we use the indicator D(t) (whose definition can be found in Appendix C.2). As nonadiabatic transitions occur, the system builds up a degree of coherence that subsequently decays as the system evolves away from the coupling region.

Figure 6

Top panel: population dynamics of the first two adiabatic electronic states P0,1(t). Solid black lines correspond to exact numerical results. Solid blue and red lines correspond to dyn-ICWF results with (N, M) = (256, 1) for the ground and first excited adiabatic populations, respectively. Dashed blue and red lines correspond to mean-field MTEF results. Bottom panel: decoherence dynamics between the ground state and first excited adiabatic electronic states, i.e., D01. Solid black lines correspond to exact results. The solid blue line corresponds to dyn-ICWF results with (N, M) = (256, 1). The dashed blue line corresponds to mean-field MTEF results.

As shown in Figure , the dyn-ICWF method reaches quantitative accuracy for (N, M) = (256, 1) and vastly outperforms the multitrajectory 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 of these aspects of this problem are qualitatively well described by the dyn-ICWF method using only (N, M) = (16, 1) (these results can be found in Appendix C.2).

Example VI: Interference Effects Near a 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[101] and extended further by Schaupp and Engel.[102] The model consists of a quantized electron and proton that can move in two Cartesian directions, along with two fixed “classical” protons, 1, 2. A CI occurs in this model when (treating the quantized proton as a BO parameter) the protons are in a D3 geometry. The potential energy isand we use the parameter values a = 0.5, b = 10, R0 = 1.5, 1 = (−0.4√3, 1.2), and 2 = (0.4√3, 1.2).

We initialize the total system wave function as a direct product of the first excited electronic BO state and a nuclear Gaussian state centered at c = (0, 0.4) with standard deviation σ2 = 5. For this placement of 1, 2, the CI occurs at the origin and, in the BO picture, the initial nuclear wavepacket “falls toward” the CI (see Figure in Appendix C.3). In this picture, the nuclear motion occurs on a single BOPES and the two portions of the nuclear wavepacket around the CI (i.e., the clockwise and anticlockwise components) cause an interference pattern to develop when they do recollide (see Figure ).

Figure 14

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 is initialized as a Gaussian centered at = (0, 0.4) on the lower surface.

Figure 7

Exact and dyn-ICWF reduced nuclear density showing the interference pattern after having traversed the conical intersection at the origin.

While the interference pattern described in Figure can be understood as the adiabatic circular motion around the position of a conical intersection, it is important to emphasize that the concept of CI makes sense only when the adiabatic picture, i.e., the Born–Huang basis expansion, is used to represent the molecular wave function. However, any observable effect that can be explained on the adiabatic basis must arise also in any other picture such as the diabatic picture or the full real-space grid picture used by the dyn-ICWF method. Therefore, while not depending on the BO picture (beyond defining the initial state), the dyn-ICWF method is able to capture the correct CI curvature effects, as well as any interference pattern that forms in the fully reduced nuclear densitySee Appendices C.3 and D for further details on the dyn-ICWF calculation.

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 wave function 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 nonperturbatively driven response properties of an exactly solvable one-dimensional H2 model system that standard mean-field theories fail to describe.

We have also considered a broader class of conditional wave functions that was formally introduced through time-dependent conditional wave function theory, yielding a set of coupled single-particle equations of motion. An approximated set of these time-dependent conditional wave functions are utilized as time-dependent basis of a stochastic wave function 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.

The sta- and dyn-ICWF methods are wave function-based approaches. Therefore, while the simple sum-of-product forms that we have employed for our ansatz in eqs and 30 can be made more efficient, by introducing a tensor network representation for the expansion coefficients such as matrix product states or hierarchical Tucker formats, for example, an exponential scaling with respect to the number of correlated degrees of freedom is expected unless approximations are introduced. That being said, we want to emphasize that the ICWF method is fundamentally different from wave function methods that rely on the Born–Huang expansion of the molecular wave function. Alternatively, the ICWF method describes electronic and nuclear degrees of freedom on the same mathematical footing, viz., the real-space grid picture. It is this particular trait that makes the ICWF an original starting point for developing novel, unexplored, approximations that could eventually yield a significant computational advantage compared to methods that rely on the Born–Huang expansion.

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 wave functions 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.[103] 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 thermalization,[104−109] correlated lattice dynamics,[110−112] and phase transitions.[113−115]