A hybrid method is proposed to propagate system-bath quantum dynamics that use both basis functions and coupled quantum trajectories. In it, the bath is represented with an ensemble of Bohmian trajectories while the system degrees of freedom are accounted through reduced density matrices. By retaining the Hilbert space structure for the system, the method is able to capture interference processes that are challenging to describe in Bohmian dynamics due to singularities that these processes introduce in the quantum potential. By adopting quantum trajectories to represent the bath, the method beats the exponential scaling of the computational cost with the bath size. This combination makes the method suitable for large-scale ground and excited state fully quantum molecular dynamics simulations. Equations of motion for the quantum trajectories and reduced density matrices are derived from the SchrΓΆdinger equation and a computational algorithm to solve these equations is proposed. Through computations in two-dimensional model systems, the method is shown to offer an accurate description of subsystem observables and of quantum decoherence, which is difficult to obtain when the quantum nature of the bath is ignored. The scaling of the method is demonstrated using a model with 21 degrees of freedom. The limit of independent trajectories is recovered when the mass of bath degrees of freedom is much larger than the one of the system, in agreement with mixed quantum-classical descriptions.

## I. INTRODUCTION

There is growing interest in dynamical processes that occurs in complex environments, such as non-adiabatic chemical dynamics, proton coupled electron transfer, and atomic collisions with surfaces.^{1β5} Often these processes can be understood using a system plus bath (environment) perspective, in which the system refers to the particular degrees of freedom (DOFs) that one is interested in, while the remaining degrees of freedom are relegated to the bath. This perspective is also central to understanding decoherence and dissipation in quantum mechanics, which arises due to entanglement and energy exchange between the system and bath.

Explicitly propagating the system-bath quantum dynamics is very challenging as the computational resources needed for exact quantum treatment scale exponentially with the system size. For systems coupled to a bath of independent harmonic oscillators, important analytical results have been obtained.^{6,7} For more complex baths, simulations are often made possible by taking a classical limit of the bath.^{8,9} These quantum-classical methods are appealing due to their computational efficiency and because they can be employed to describe baths with chemical detail. However, one of the problems of these approximate methods is quantum decoherence.^{10} Quantum decoherence cannot be accurately described by a classical bath because it arises due to entanglement. Further, there is evidence that quantum decoherence plays an important role for a variety of systems of chemical relevance.^{11β15} At the short time limit, the frozen Gaussian approximation for the bath wavefunction, where a classical trajectory is replaced by a Gaussian wavefunction parameterized by its position and momentum, can be used to alleviate the problem.^{11} In the full quantum regime, the multi-configuration time-dependent Hartree (MCTDH) method,^{16β18} even when it does not solve the problem of exponential scaling with the system size, has been applied to study system-bath dynamics in systems with dozens of degrees of freedom.

As a possible strategy to avoid the exponential scaling with the dimensionality of the system, there is growing interest in the hydrodynamic or quantum trajectory formulation of the time-dependent SchrΓΆdinger equation (TDSE).^{19β28} Quantum trajectories (QTs) offer the most efficient representation of the time-dependent wavefunction because their distribution follows the time-dependent probability. The main numerical difficulty stems from the approximation to the quantum potential that arises in this representation. In the course of wavepacket dynamics, if there are no nodes generated in the wavefunction, the quantum potential can often be appropriately approximated using a least-squares fit to a polynomial function. When nodes or even quasi-nodes are formed, the required computational scheme becomes much more complicated^{29} and a direct application of the quantum trajectory method (QTM) is not convenient. As an example, this problem arises when the subsystem potential has a more complicated structure than the one of the bath. This occurs, for instance, in a typical proton transfer reaction where the reaction coordinate is of double-well shape while the environment, the collective vibrations of the molecule, can often be described by a set of coupled (an)harmonic oscillators.^{30} Another case in which the QTM cannot be directly used occurs when the subsystem is represented by a finite-dimensional Hilbert space, e.g., a spin coupled to a bath of harmonic oscillators.

In this paper, we propose a method that partially circumvents such difficulties by combining the conventional basis set method with quantum trajectories. Essentially, the method treats two coupled quantum systems using different frameworks of quantum mechanics (wave mechanics and Bohmian mechanics). The strategy is analogous to other hybrid strategies for system-bath dynamics such as mixed quantum-classical Liouville dynamics.^{31} The method takes advantage of the compact representation of quantum dynamics achieved by the quantum trajectories to allow dealing, fully quantum mechanically, with a bath of hundreds of degrees of freedom (DOFs). By incorporating a basis, in turn, it allows treating not only a subsystem represented by a finite-dimensional Hilbert space but also cases where interference is important for subsystem dynamics. Interference is straightforward to capture using a set of basis functions but difficult for trajectories. Below, through numerical simulations in low-dimensional models, we show that such a strategy is computationally tractable and accurately describes decoherence dynamics.

The paper is organized as follows: Section II introduces the theoretical strategy and equations of motion for the quantum trajectories and associated reduced density matrices. Section III discusses strategies to solve the equations of motion and contrasts the behavior of the method against exact quantum results in two-dimensional models with a bilinear and nonlinear system-bath coupling. Our main findings are summarized in Section IV.

## II. FORMALISM

The Hamiltonian for a system-bath problem can be written as

where $Hs=\beta \x88\x91k=1Nspk22mk+Vs(\pi \x9d\x92\x92)$ is the Hamiltonian for the system, $Hb=\beta \x88\x91\Xi \u038c=1NbP\Xi \u038c22M\Xi \u038c+Vb(\pi \x9d\x91\u0388)$ is the one for the bath, and $Vsb\beta \x89\u2018Vsb(\pi \x9d\x92\x92,\pi \x9d\x91\u0388)$ represents the interaction between them. Here, *N*_{s} is the number of degrees of freedom for the system and *N*_{b} for the bath. The quantities $\pi \x9d\x92\x92,\pi \x9d\x91\u0388$ denote system and bath coordinates with conjugate momenta $\pi \x9d\x92\x91=(p1,p2,\beta \x80\xa6)$ and $\pi \x9d\x91\xb7=(P1,P2,\beta \x80\xa6)$β , respectively.

As a starting point for the partial hydrodynamic representation, the wavefunction of the composite system is written as

where $\Omicron \x88(\pi \x9d\x92\x92,\pi \x9d\x91\u0388,t)$ is the subsystem conditional wavefunction and $\Omicron \x95(\pi \x9d\x91\u0388,t)$ is the bath wavefunction. Such exact decomposition always exists^{32} but is not unique. To uniquely determine the evolution of $\Omicron \x88(\pi \x9d\x92\x92,\pi \x9d\x91\u0388,t)$β , one has to specify the initial conditions and the time evolution for $\Omicron \x95(\pi \x9d\x91\u0388,t)$β . Substituting this form for the wavefunction into the TDSE $i\beta \x88\x82\beta \x88\x82t|\Xi \xa8\beta \x9f\xa9=H|\Xi \xa8\beta \x9f\xa9$ gives

where $\beta \x88\x87\Xi \u038c\beta \x89\u2018\beta \x88\x82\beta \x88\x82Q\Xi \u038c$ is the derivative with respect to a bath coordinate $\Xi \u038c$β . For simplicity, we have dropped the explicit arguments of the wavefunction in Eq. (3). Atomic units are used throughout, $\beta \x84\x8f=1$β .

As an ansatz, let $\Omicron \x95(\pi \x9d\x91\u0388,t)$ satisfy the effective bath SchrΓΆdinger equation

with the effective Hamiltonian

Here $\Xi \x9b(\pi \x9d\x91\u0388,t)$ is an unspecified function that couples $\Omicron \x88(\pi \x9d\x92\x92,\pi \x9d\x91\u0388,t)$ and $\Omicron \x88(\pi \x9d\x91\u0388,t)$β , and that can be chosen based on numerical convenience. Substituting Eq. (4) into Eq. (3) and dividing both sides by $\Omicron \x95(\pi \x9d\x91\u0388,t)$ gives

As initial conditions for Eqs. (4) and (6), $\Omicron \x95(\pi \x9d\x91\u0388,t)$ is chosen such that the probability $|\Omicron \x95(\pi \x9d\x91\u0388,t)|2$ coincides with the true bath probability function $\beta \x88\xab|\Xi \xa8(\pi \x9d\x92\x92,\pi \x9d\x91\u0388,t)|2d\pi \x9d\x92\x92$β . This together with Eq. (2) uniquely defines $\Omicron \x88(\pi \x9d\x92\x92,\pi \x9d\x91\u0388,t)$ and $\Omicron \x95(\pi \x9d\x91\u0388,t)$ up to a phase factor. Note that because we do not impose the partial normalization condition $\beta \x88\xab|\Omicron \x88(\pi \x9d\x92\x92,\pi \x9d\x91\u0388,t)|2d\pi \x9d\x92\x92=1$β , $|\Omicron \x88(\pi \x9d\x92\x92,\pi \x9d\x91\u0388,t)|2$ does not generally coincides with the conditional probability in Ref. 33 unless this additional constraint is imposed.

In exact treatments, the choice of Ξ is immaterial because if Eqs. (4) and (6) are solved exactly they are equivalent to solving the TDSE [Eq. (3)] which is independent of Ξ. In approximate treatments, Ξ can be chosen to enhance the convergence properties of the method as discussed in Section II D. In the proposed method, $\Omicron \x88(\pi \x9d\x92\x92,\pi \x9d\x91\u0388,t)$ is represented in a basis while the effective bath SchrΓΆdinger equation [Eq. (4)] is solved using the QTM as described below.

### A. Bath dynamics using quantum trajectories

The bath dynamics is solved within the Bohmian formulation of quantum mechanics. In this quantum trajectory framework, the bath wavefunction $\Omicron \x95(\pi \x9d\x91\u0388,t)$ is represented by an ensemble of the so-called quantum trajectories. The equation of motion for the quantum trajectories is similar to the Newtonian mechanics except for an additional quantum potential.

Specifically, the Bohmian formulation of quantum mechanics or hydrodynamics formulation is based on the polar form of the wavefunction $\Omicron \x95(\pi \x9d\x91\u0388,t)=A(\pi \x9d\x91\u0388,t)eiS(\pi \x9d\x91\u0388,t)$β , where $A(\pi \x9d\x91\u0388,t)$ and $S(\pi \x9d\x91\u0388,t)$ are real functions. Substituting this form into the time-dependent SchrΓΆdinger equation [Eq. (4)] yields

Here $\Xi \xb7(\pi \x9d\x91\u0388,t)=A2(\pi \x9d\x91\u0388,t)$ is the probability density, and

is the so-called quantum potential that enters into the equations of motion together with the external time-dependent classical potential $V(\pi \x9d\x91\u0388,t)=Vb(\pi \x9d\x91\u0388)+\Xi \x9b(\pi \x9d\x91\u0388,t)$β . We have made *β* explicit in Eq. (9) to highlight its full quantum nature. Once the gradient of the wavefunction phase is associated with the momentum of the quantum trajectory,

then Eq. (8) becomes an equation of continuity of the wavefunction density $\Xi \xb7(\pi \x9d\x91\u0388,t)$β . Equation (7) is the quantum Hamilton-Jacobi equation, which differs from the classical Hamilton-Jacobi equation by the quantum potential term. The quantum potential is responsible for all quantum-mechanical effects. It is nonlocal in ** Q** and suggests a formal transition to classical mechanics when $U\beta \x86\x920$β . In this limit, the potential determining the dynamics becomes local and the coupled quantum trajectories become independent.

Differentiation of Eq. (7) with respect to the spatial coordinates defines the time-evolution of ** P**,

where $\pi \x9d\x92\x97=(v1,v2,\beta \x80\xa6,vNb)$ and ${v\Xi \u038c=P\Xi \u038c/M\Xi \u038c,\Xi \u038c=1,2,\beta \x80\xa6,Nb}$β . Transferring Eq. (11) into the Lagrangian frame-of-reference where

gives Newtonβs equation of motion for the quantum trajectory,

### B. System dynamics in a basis

Assume there is a convenient time-independent basis set ${\Omicron \x86n(\pi \x9d\x92\x92)}$β , $H0|\Omicron \x86n\beta \x9f\xa9=En|\Omicron \x86n\beta \x9f\xa9$β , where *H*_{0} is a Hermitian operator in the system subspace, to describe the system. In this basis,

where *n*_{b} is the number of basis functions used to represent the wavefunction. When the system and bath are entangled, the expansion coefficients will depend on the bath coordinates no matter what $\Omicron \x95(\pi \x9d\x91\u0388,t)$ and the basis are chosen. The basis can be chosen to be static or adapted dynamically. A static basis is convenient when the Hilbert space of the subsystem is of finite dimension such as the spin-boson model. In other cases such as in the studies of non-adiabatic molecular dynamics, the adiabatic basis, $(Hs+Vsb(\pi \x9d\x92\x92,\pi \x9d\x91\u0388))\Omicron \x86n(\pi \x9d\x92\x92;\pi \x9d\x91\u0388)=En(\pi \x9d\x91\u0388)\Omicron \x95n(\pi \x9d\x92\x92;\pi \x9d\x91\u0388)$β , is often more convenient.

Substituting Eq. (14) into Eq. (6), left multiplying by $\Omicron \x86m*(\pi \x9d\x92\x92)$ and integrating over ** q** gives

where $V\Beta \u2015sb,mn(\pi \x9d\x91\u0388)=\beta \x88\xab\Omicron \x86m*(\pi \x9d\x92\x92)Vsb(\pi \x9d\x92\x92,\pi \x9d\x91\u0388)\Omicron \x86n(\pi \x9d\x92\x92)d\pi \x9d\x92\x92$β . For a system-bath interaction potential of the form $Vsb=\beta \x88\x91\Xi \xb1S^\Xi \xb1\beta \x8a\x97B^\Xi \xb1$β , where $S^\Xi \xb1$ and $B^\Xi \xb1$ are operators in the subsystem and bath space, respectively, $V\Beta \u2015sb,mn=\beta \x88\x91\Xi \xb1S\Xi \xb1,mnB^\Xi \xb1$β . As the bath wavefunction is represented by an ensemble of trajectories, the above equation can be discretized in terms of these trajectories ${\pi \x9d\x91\u0388(\Xi \xb1)(t),\Xi \xb1=1,\beta \x80\xa6,Ntraj}$β , where $Ntraj$ is the total number of trajectories. Realizing that

where $R\Xi \u038c$ is the so-called nonclassical momentum,^{34β36} and writing Eq. (15) in the moving frame of reference [Eq. (12)], the equations of motion for ** c** is

where the classical momentum $P\Xi \u038c$ has been absorbed into the total time derivative. The first two terms in the right-hand side of Eq. (17) represent terms that can be determined by a single trajectory while the last term involves the spatial derivative of the expansion coefficients. The latter represents the coupling between trajectories because calculating gradients involves determining how the coefficients change among trajectories.

In this way, the time-dependent SchrΓΆdinger equation has been transformed into a discretized version in both subsystem space by using a set of basis functions and also bath space by using an ensemble of trajectories (Eqs. (13) and (17)). Representing the bath configuration with an ensemble of trajectories allows Monte Carlo sampling of bath configurations, thus opening the possibility of treating larger number of DOFs for the bath. The exact quantum dynamics is obtained in the limit of a large number of quantum trajectories and basis functions.

### C. Density matrix and dynamical observables

Below we describe how usual quantum mechanical quantities are calculated in the partial hydrodynamic representation. Consider first the normalization of the total wavefunction in the hybrid picture,

Here we have defined the weight for each trajectory as $w\Xi \xb1\beta \x89\u2018|\Omicron \x95(\pi \x9d\x91\u0388(\Xi \xb1))|2\Xi \u0384\pi \x9d\x91\u0388(\Xi \xb1)$β , where $\Xi \u0384\pi \x9d\x91\u0388(\Xi \xb1)$ is the volume element around $\pi \x9d\x91\u0388(\Xi \xb1)$β . If the dynamics is performed in the Lagrangian frame of reference defined by Eq. (12), the weight is a constant for each trajectory.^{23} If importance sampling from the initial probability density $|\Omicron \x95(\pi \x9d\x91\u0388,t=0)|2$ is used, the weight is the same for each trajectory, $w\Xi \xb1=1/Ntraj$β .

In this hybrid representation, the reduced density matrix of the subsystem is given by

where

is the reduced density matrix associated with each quantum trajectory $\pi \x9d\x91\u0388(\Xi \xb1)(t)$ and $TrB[\beta \x8b\u2015]$ (β $TrS[\beta \x8b\u2015]$β ) represents the partial trace over bath (system) DOFs. Note that $\Omicron \x81(\Xi \xb1)$ does not necessarily satisfy $TrS[\Omicron \x81(\Xi \xb1)]=1$β , as probability is allowed to exchange between neighboring QTs. The expectation value of subsystem observables $S^$ is given by

where $Smn=\beta \x9f\xa8\Omicron \x86m|S^|\Omicron \x86n\beta \x9f\xa9\pi \x9d\x92\x92$β . Here $\beta \x9f\xa8\beta \x8b\u2015\beta \x9f\xa9\pi \x9d\x92\x92$ denotes integration over system coordinates only.

### D. Choosing Ξ

In the limit of large number of QTs, the method is independent of the choice of Ξ. When a finite number of trajectories are employed, a careful choice of $\Xi \x9b(\pi \x9d\x91\u0388,t)$ can enhance the convergence of the method with respect to the number of trajectories. In essence, this term couples the system state $\Omicron \x88(\pi \x9d\x92\x92,\pi \x9d\x91\u0388,t)$ to $\Omicron \x95(\pi \x9d\x91\u0388,t)$β , and it can be used to ensure that the wavefunction $\Omicron \x95(\pi \x9d\x91\u0388,t)$ always has spatial support close to the true bath distribution function $\beta \x88\xab|\Xi \xa8(\pi \x9d\x92\x92,\pi \x9d\x91\u0388,t)|2d\pi \x9d\x92\x92$β . Guaranteeing such spatial support enhances the convergence properties of the method with respect to the number of trajectories. In an ideal situation, one would like to choose Ξ such that $|\Omicron \x95(\pi \x9d\x91\u0388,t)|2\beta \x89\x88\beta \x88\xab|\Xi \xa8(\pi \x9d\x92\x92,\pi \x9d\x91\u0388,t)|2d\pi \x9d\x92\x92$ for all times as it optimizes the sampling. When the condition is satisfied exactly, the method reduces to the exact factorization method of the nuclear-electronic wavefunction in Refs. 33 and 37 for cases in which the effect of the vector potential in that formalism can be avoided or neglected. In this limit, $\beta \x88\xab|\Omicron \x95(\pi \x9d\x92\x92,\pi \x9d\x91\u0388,t)|2d\pi \x9d\x92\x92=1$β , which is the partial normalization condition used in Ref. 33 and Ξ is akin to the exact time-dependent potential energy surface (TDPES). Since this TDPES is challenging to obtain (see Ref. 38 on recent developments on how to do so), and the exact nuclear quantum dynamics may still be numerically challenging on the TDPES, a practical alternative is required. Here we choose Ξ to be *V*_{sb} averaged over subsystem coordinates,

The factor $Z\Xi \xb1$ is used to account for the fact that the reduced density matrix of a single trajectory does not always have a unit trace. If an adiabatic basis is used, this will be the mean-field potential used in Ehrenfest dynamics.^{39} This choice is useful because the Ehrenfest potential is the leading term in the TDPES.^{40}Β

### E. Independent trajectories with Ehrenfest potential

From a computational perspective, it is always advantageous to have independent trajectories. However, the coupled nature of quantum trajectories due to the quantum potential requires propagating an ensemble of them simultaneously, as dynamic properties, such as the quantum potential, can only be extracted from ensemble averages. The equations derived so far are exact. If independent trajectories are desired, an approximation is required to decouple them.

To obtain independent trajectories, one needs to assume that (i) the quantum potential is zero and (ii) the system and bath are not entangled such that the total wavefunction can be written as a product state $|\Xi \xa8(\pi \x9d\x92\x92,\pi \x9d\x91\u0388,t)\beta \x9f\xa9=|\Omicron \x88(\pi \x9d\x92\x92)\beta \x9f\xa9\beta \x8a\x97|\Omicron \x95(\pi \x9d\x91\u0388)\beta \x9f\xa9$β . The latter condition is required such that the terms $\beta \x88\x87\Xi \u038c2cn,\beta \x88\x87\Xi \u038ccn$ coupling the trajectory-dependent reduced density matrices vanish. In this limit, the coupled trajectory picture is reduced to an ensemble of independent trajectories, each carrying a density matrix for the subsystem. The equation of motion for ** c** in this independent trajectory approximation (ITA) is

Here the dynamics for the reduced density matrix associated with a single trajectory is unitary. The non-unitary dynamics or decoherence for the subsystem now comes from the averaging of bath configurations. Physically, this limit is obtained when the mass of the bath is much larger than the system one (see Ref. 41 for a detailed study of this limit).

Note that when the mass of the bath is finite, the ITA does not preserve the invariance of the method with the choice of Ξ. Such invariance is strictly recovered in the limit of $M\Xi \u038c\beta \x86\x92\beta \x88\x9e$β . As shown in Sec. III, the ITA usually does not give an accurate description of decoherence and even observables of the subsystem.

## III. IMPLEMENTATION AND MODEL COMPUTATIONS

### A. Approximate quantum potential and $cn(\pi \x9d\x91\u0388)$

To solve the equations of motion for the quantum trajectories (Eq. (13)), the quantum potential has to be determined using approximations.^{9,20,34,42β44} This is because computing it requires reconstructing the full wavefunction from the quantum trajectories. Here we use the linear quantum force (LQF) method^{34,45} because it conserves energy, is computationally practical, and admits a well defined variational optimization procedure. The details of the method can be found in Ref. 34. Its essence is summarized below.

The quantum potential [Eq. (9)] can be written in terms of the nonclassical momentum $R\Xi \u038c$ [Eq. (16)],

In the LQF, $R\Xi \u038c$ is approximated as an expansion in a linear basis of bath coordinates ** Q**, $\pi \x9d\x92\x87=(1,Q1,Q2,\beta \x80\xa6,QNb)$β ,

Minimization of $\beta \x9f\xa8(R\Xi \u038c\beta \x88\x92R\beta \x88\u038c\Xi \u038c)2\beta \x9f\xa9$ with respect to the expansion coefficients **A** is achieved by solving the linear matrix equation,

where the elements of the overlap matrix **M** and of the right hand-side matrix **B** are

The quantum potential and quantum force can be readily obtained using Eq. (24) after obtaining the analytical expression for $R\Xi \u038c$β . The LQF method has been shown to be able to capture basic quantum mechanical effects, such as wavepacket bifurcation, moderate tunneling, and zero-point energy.^{23,30}

Any existing method to approximate the quantum potential can also be applied to approximate $cn(\pi \x9d\x91\u0388)$β . This is because both cases require the computation of derivatives of a function on unstructured grids determined by the trajectories ${\pi \x9d\x91\u0388(\Xi \xb1)}$β . Here a global least-squares fitting^{23} is used to obtain ${\beta \x88\x87\Xi \u038c2cn,\beta \x88\x87\Xi \u038ccn}$β , which uses a polynomial basis function up to a certain order. For example, in a quadratic fitting, the expansion coefficient *c*_{n} is approximated by $c\beta \x88\u038cn=Cn(0)+\beta \x88\x91\Xi \u038cCn\Xi \u038c(1)Q\Xi \u038c+\beta \x88\x91\Xi \u038cCn\Xi \u038c(2)Q\Xi \u038c2$β . All the fitting coefficients are obtained by minimizing the function $I(n)=\beta \x88\x91\Xi \xb1w\Xi \xb1|cn(\pi \x9d\x91\u0388(\Xi \xb1))\beta \x88\x92c\beta \x88\u038cn(\pi \x9d\x91\u0388(\Xi \xb1))|2$β . A quartic polynomial fitting was found to be adequate for the simulations below.

### B. Model computations

As an exemplifying model, we consider a quartic oscillator coupled with a bath of harmonic oscillators with Hamiltonian

To benchmark the method with exact quantum results, we focus first on the two-dimensional model case where the bath is represented by a single oscillator, i.e., $H(x,Q)=p22m+12m\Omicron \x8902x2+\Xi \xbb4x4+P22M+12M\Omicron \x89Q2+V(x,Q)$β . We consider both the bilinear *V*(*x*,*Q*)=*gxQ* and nonlinear coupling cases *V*(*x*, *Q*)=*gx ^{2}Q^{2}*/4. Then we apply the method to a bath with 20 harmonic oscillators bilinearly coupled to the system oscillator. The nonlinearly coupled system is considerably more complicated than the bilinear one because it does not admit simplifications that apply for harmonic systems. In fact, for quadratic systems, it is possible to eliminate the bath variables and account for the influence of the environment to the subsystem by influence functionals.

^{7}The initial wavefunction is set as a product of Gaussian wavepackets

where $\Xi \xb1=1.0\beta \x80\x89a0\beta \x88\x922$ and $\Xi \xb2\Xi \u038c=16.0\beta \x80\x89a0\beta \x88\x922$β . The initial momentum $P\Xi \u038c$ of the trajectories is 0 because the initial wavefunction is chosen to be real. The remaining parameters for the model and initial wavefunction are shown in Table I. The basis functions used in the computation are the eigenstates of the harmonic part of the system Hamiltonian. The system and bath oscillators are chosen to have comparable vibrational frequency such that the dynamics is not in the adiabatic limit. To characterize the dynamics, we follow the expectation value of the subsystem position and the purity $P=Tr[\Omicron \x81S2]$β . The purity is a sensitive quantity as it depends on all elements of the reduced density matrix of the subsystem and quantifies the degree of entanglement between the system and bath. Exact quantum results are obtained with the split-operator method.^{46} The number of trajectories for all simulations is 1024.

.Β | V(x, Q)
.Β | n_{b}
.Β | m [a.u.] .Β | M [a.u.] .Β | x_{0} [a_{0}]
.Β | Q^{0} [a_{0}]
.Β | g .Β | Ξ» [$Eha0\beta \x88\x924$]
.Β | $\Omicron \x890$ .Β | $\Omicron \x891$ .Β |
---|---|---|---|---|---|---|---|---|---|---|

Model 1Β | gxQΒ | 8Β | 1.0Β | 8.0Β | 0.4Β | 0.3Β | 0.5 [$Eha0\beta \x88\x922$]Β | 0Β | 1.0Β | $5/2$Β |

Model 2Β | $g4x2Q2$Β | 8Β | 1.0Β | 8.0Β | 0.4Β | 0.3Β | 1.0 [$Eha0\beta \x88\x924$]Β | 0.4Β | 1.0Β | $5/2$Β |

.Β | V(x, Q)
.Β | n_{b}
.Β | m [a.u.] .Β | M [a.u.] .Β | x_{0} [a_{0}]
.Β | Q^{0} [a_{0}]
.Β | g .Β | Ξ» [$Eha0\beta \x88\x924$]
.Β | $\Omicron \x890$ .Β | $\Omicron \x891$ .Β |
---|---|---|---|---|---|---|---|---|---|---|

Model 1Β | gxQΒ | 8Β | 1.0Β | 8.0Β | 0.4Β | 0.3Β | 0.5 [$Eha0\beta \x88\x922$]Β | 0Β | 1.0Β | $5/2$Β |

Model 2Β | $g4x2Q2$Β | 8Β | 1.0Β | 8.0Β | 0.4Β | 0.3Β | 1.0 [$Eha0\beta \x88\x924$]Β | 0.4Β | 1.0Β | $5/2$Β |

Consider first the case of two bilinearly coupled harmonic oscillators, which is an analytically solvable problem. Figure 1 shows the expectation value of the position of the subsystem using the proposed method, the proposed method in the independent trajectory approximation, and also the exact quantum results. For completeness, in Figs. S1 and S2 in the supplementary material, we have included a plot of the quantum potential and expansion coefficients *c*_{n}(*Q*) at selected times of the dynamics for the coupled trajectory case. It is clear that the coupled trajectory method is indistinguishable with the exact quantum results. In turn, the results obtained with the independent trajectory approximation are only accurate at short times and then exhibit a phase shift of the periodic motion and overestimate the amplitude of motion, which implies an incorrect energy flow from the bath oscillator to the subsystem. Figure S3 in the supplementary material shows the error in $\beta \x9f\xa8x\beta \x9f\xa9$ when a reduced number of trajectories is employed. As shown, the method offers accurate results even with a quarter of the number of trajectories (256).

Figure 2 shows the purity of the subsystem with time for the harmonic model (model 1 in Table I). The decoherence when the ITA is used is coming from the average inherent to the initial sampling of the bath wavefunction. As can be seen from this figure, the purity computed with the ITA exaggerates the decoherence, whereas the results obtained with coupled trajectories capture all the fine features albeit with a small shift in the amplitude.

The same quantities are also computed for the model with an anharmonic oscillator nonlinearly coupled with a harmonic oscillator (model 2 in Table I), and the results are shown in Figs. 3 and 4. For model 2, the decoherence rate is slower than model 1, but it shares the main features of model 1. The expectation value of the subsystem position in the ITA agrees well with the exact one at short times but deviates at longer times. However, unlike model 1, for model 2 the independent trajectories underestimate the decoherence rate. Thus the effects of using independent trajectories on the decoherence rate depend on the particular system. The results with coupled trajectories follow closely the exact position and exact purity. However, it shows larger oscillations in purity at later times, which may come from inaccuracies in the fitting of the expansion coefficients.

To demonstrate the favorable scaling of the method with respect to the bath size, Fig. 5 shows a simulation for a harmonic oscillator bilinearly coupled with 20 bath oscillators. The oscillator frequencies for the bath are chosen to be uniformly distributed in the range [0.2, 2] and the coupling $g\Xi \u038c=0.2\Omicron \x89\Xi \u038c$β . The system observes decoherence and damping of the oscillations. The simulations take about 2 h using a single central processing unit (CPU) if a global linear basis for the fitting of the coefficients is used, in comparison to around 6 min for the two-dimensional models. Further, parallelization can speed up the computation as illustrated in Ref. 23. Note, however, that the accuracy of the results in Fig. 5 has not been established as exact quantum results are unreachable by the split-operator method.

## IV. FINAL REMARKS

In this paper, we have developed a partial hydrodynamic method to follow quantum dynamics of a system interacting with an environment. The method treats both the subsystem and bath quantum mechanically. The bath wavefunction is represented by an ensemble of coupled quantum trajectories, while the subsystem retains a Hilbert space structure. In this manner, the method takes advantage of the ability of Hilbert space dynamics to capture interference and the effects of strong anharmonicity in the potential. Further, the hydrodynamic strategy for the bath incorporates favorable scaling with the bath size that partially circumvents the limitation of other full quantum treatments such as a split-operator method, path integral, and MCTDH. In this way, it extends the regime of problems that can be dealt with using full quantum dynamics.

The numerical difficulty associated with the method is in computing the terms that couple quantum trajectories. They appear (i) in the approximation of the quantum potential and (ii) in the equation of motion for the expansion coefficients of the subsystem wavefunction in a basis (Eq. (17)). The latter involves derivatives of the coefficients with respect to the bath coordinates. Nevertheless, there is a hierarchy of approximation schemes that can be used to systematically improve the results. Specifically, the method can be systematically improved by using a higher-order polynomial basis or local least-squares fitting at the expense of higher computational cost.^{19,34} Further, the implementation can be readily parallelized and has modest memory requirements because it does not require storing the full wavefunction of the composite system. In addition, the formalism offers a rigorous platform to think about decoherence in terms of trajectories.

The method was applied to two-dimensional model systems which show that the coupled trajectories can recover all the essential features of the exact quantum results. By contrast, the ITA gives an adequate description of subsystem dynamics at short times but deviates from exact results at long times. Also the ITA does not offer an accurate description of the decoherence dynamics, which implies that the coupling terms among trajectories are important for entanglement between the system and bath. The ITA naturally arises when the mass of the bath is large. In this limit, with the choice of the adiabatic basis and Ehrenfest potential for Ξ, the method reduces to an ensemble of Ehrenfest trajectories. Notice that although the ITA used here is not accurate enough at long times, the independence of trajectories is still an attractive feature from a numerical perspective. An improved version of the ITA will be highly useful for high-dimensional systems.

Future prospects include understanding decoherence in non-adiabatic molecular dynamics using this hybrid strategy and exploring the potential of the method to follow the dynamics of large quantum systems.

## SUPPLEMENTARY MATERIAL

See supplementary material for quantum potential and expansion coefficients at selected times, and convergence with respect to the number of trajectories for model 1.

## ACKNOWLEDGMENTS

This material is based on the work supported by the National Science Foundation under No. CHE-1553939.