Piecewise ensemble averaging stochastic Liouville equations for simulating non-Markovian quantum dynamics

Here we present a novel stochastic Liouville equation with piecewisely correlated noises, in which the inter-piece correlation is rigorously incorporated by a convolution integral involving functional derivatives. Due to the feature of piecewise correlation, we can perform piecewise ensemble average and serve the average of the preceding interval as the initial condition of the subsequent propagation. This strategy avoids the long-time stochastic average and the statistical errors are saturated at long times. By doing so, we circumvent the intrinsic difficulty of the stochastic simulations caused by the fast increase in the variance of the quantum Brownian motion. Therefore, as demonstrated by the numerical examples, the proposed method enables us to simulate the long-time quantum dissipative dynamics with long memories in the non-perturbative regime.


Introduction
The dynamics of many quantum systems are significantly altered by their environments. To avoid the curse of dimensionality due to the direct consideration of the environments, the quantum Brownian motion paradigm is used in diversified fields of physics, chemistry and biology [1][2][3], ranging from the fundamental issues of quantum measurements [4], quantum statistical mechanics [5,6], quantum phase transition [7], and quantum information [8], to the applied aspects of quantum biology [9,10] and chemical reactions [11].
Driven by the environmental quantum fluctuations, the evolution of the system of interest is neither unitary nor deterministic, and is in principle non-Markovian [12,13]. These features, sometimes interplayed with many-body effects, make the simulation of quantum dissipative dynamics highly challenging. The situation is further complicated by the fact that a deterministic description of an open quantum system requires the complete knowledge of the microscopic structure of the whole world, which is impossible [14]. A natural solution then is a stochastic description by representing the fluctuating environment with stochastic fields, which originated from the Einstein relation [15] and thrived due to the stochastic re-interpretation of the Feynman-Vernon influence functional [16][17][18]. At present, a series of stochastic methods have been developed for simulating the quantum dissipative dynamics, including the quantum Langevin equation [19,20], quantum jump [21], quantum Monte Carlo [18,22], quantum state diffusion and its non-Markovian extension (NMQSD) [23][24][25], and the stochastic Liouville equations (SLE) [26][27][28][29]. While most of the stochastic methods were designed for bosonic environments, they were also extended to fermionic bathes [30,31]. The stochastic scheme not only offers practical simulation tools, but also serves as a methodology development platform, e.g. for deriving quantum master equations [32], hierarchical equation of motion (HEOM) [33][34][35], and various hybrid stochastic-hierarchical equations [36][37][38][39].
Despite the impressive progress, direct stochastic simulations with available methods only work well for short-time dynamics or weak to moderate dissipation. The main obstacle lies in the intrinsic increase of the ensemble variance at long times. Since the workhorse of stochastic methods is the linear stochastic differential equation, the difficulty may be clearly illustrated by the simple colored-noise extension of the geometric Brownian motionẋ with the initial condition x(0) = 1. Here, ζ(t) is a zero-mean colored noise assuming the covariance ζ(t)ζ(s) = γ(t − s) = γ 0 exp(−|t − s|/τ ), with τ being the memory time, κ the fiction coefficient, and γ 0 the dissipation strength. The ensemble average 1} are finite at any finite time. However, in the long-time region (t τ ) the variance will grow exponentially as exp[2(2γ 0 τ − κ)t] in the cases of long memory (τ > κ/2γ 0 ). Nonlinear systems like the spin-boson model, though not analytically solvable, have been numerically demonstrated to assume a similar exponential increase in variance [40]. Consequently, to satisfy the preset error tolerance, the number of stochastic trajectories, i.e., the computational cost, grows exponentially with time.
The basic idea to circumvent the fast error increase is to avoid long time stochastic average. Attempts toward this direction took advantage of the finite memory time. For example, Cerrillo and Cao have suggested a transfer tensor method by extracting non-Markovian dynamics maps from the short-time evolution to predict long-time propagation [41]. Stockburger has split the noises into long-range and short-range ones. While the short-range noise is accounted for in the finite-memory propagation, the long-range noise is sampled directly. The resulting finite-memory stochastic propagation [40] and the finite dephasing-time variant [42] can enhance the computational efficiency by orders of magnitude. However, with a long memory the sampling of the long-range noise drastically degrades the efficiency of the finite-memory propagation and makes it (and other available stochastic schemes) impractical to simulate low-temperature dissipative dynamics.
In this work we propose a piecewise ensemble averaging stochastic equation (PEASLE) for simulating non-Markovian quantum dynamics. To this end we partition the entire time domain into equidistant pieces [T n , T n+1 ] (T n = nT) with T being a time step of an arbitrary length. The covariances of the involving noises have a piecewise structure so that the noises are uncorrelated between different pieces. The stochastic average can therefore be done piecewisely and independently, and the average of the preceding piece may serve as the initial condition of the subsequent piece. At first glance this idea seems impossible because the time correlation in the non-Markovian dynamics will be lost in the independent piecewise average. Even if the inter-piece correlation can be retrieved by modifying the stochastic Liouville equation, one may still think that a time-convolution, like that employed in the Nakajima-Zwanzig projection method [43], has to be used for propagating each stochastic trajectory from the initial time, which would ruin the piecewise average and is practically prohibited for stochastic simulations. Here we will demonstrate the possibility and present the rigorous derivation of PEASLE based on the noise decomposition. Starting from a conventional stochastic scheme, the noises may be decomposed into two parts, a principal part assuming piecewise correlations and an auxiliary part accounting for the inter-piece correlations. The partial ensemble average over the auxiliary noise converts the system-noise correlation to an integral involving functional derivatives, and yields an integro-differential stochastic equation that rigorously takes into account the inter-piece correlation.

One-dimensional Brownian motion
We now illustrate the form of PEASLE and its basic features with the one-dimensional Brownian motion in equation (1). We first split the noise ζ as ζ(t) = ξ(t) + η(t), where the noise ξ is the principal part assuming the covariance with θ(x) being the Heaviside step function (unity for x > 0 and zero for x 0). Within such a scheme, the different pieces of the noise ξ are uncorrelated, and the missing inter-piece correlation has to be incorporated by the auxiliary noise η. To this end, the noise η is of zero mean, uncorrelates with ξ and assumes the non-vanishing covariance To proceed, we take a partial stochastic average over the noises η(t) and define a new stochastic quantityx(t) = x(t) η . The average involves the noise-system correlation η(t)x(t) η , which, according to the Furutsu-Novikov theorem [44], may be converted to an integral over functional derivatives, yielding a new stochastic differential equation, where K = t/T with · denoting the floor function. It includes explicitly the principal noise ξ only, while the influence of the auxiliary part η is incorporated through the integral over the function derivative The corresponding procedure is illustrated in figure 1. As plotted in figure 1(b), the principal covariance, stationary in each square and zero outside, is of a clear piecewise pattern. The auxiliary noise η, as depicted in figure 1(c), is complementary to the principal noise ξ and reproduces the full covariance of the original noise ζ.
Thanks to the elimination of inter-piece correlation, the principal noises within difference pieces can be generated independently and identically. Consequently, the ensemble average may be carried out separately, and the average at t = T n may serve as the initial condition for the propagation of each trajectory in the time interval T n < t < T n+1 .
Note that the choice of T affects the computational efficiency of numerical simulations. For such a one-dimensional model and with a long memory, the variance grows quadratically at short time but exponentially at long time, and there exists a crossover τ c from a power law rise to an exponential growth. Generally the constraint with T < τ c is necessary for an efficient stochastic simulation because a larger value of T in the exponential growth region leads to larger statistic errors and degrades the computational efficiency. Obviously, the extreme choice T = ∞ turns PEASLE into the conventional SLE scheme.
We emphasize that the proposed PEASLE scheme is exact. The seemingly lost correlation between t T n and t > T n in the piecewise average is fully compensated by the convolution integral over the functional derivative. As clearly indicated by the comparison between the PEASLE trajectories and the conventional ones in figure 1(d), the long-time ensemble average is avoided in PEASLE, and the variance of stochastic trajectories is greatly suppressed. Further, because the functional derivative δx(t)/δξ(τ ) is solved exactly, the computational effort of individual trajectories in PEASLE is about the same as that in the conventional scheme.

Quantum dissipative systems
For general many-body systems there are still no reliable methods to simulate their quantum dynamics. But within the framework of quantum Brownian motion, a generic model, proposed by Caldeira and Leggett [16,17], enables practical simulations by simplifying the environment with infinite 'effective' harmonic oscillators. In this paradigm the total Hamiltonian can be written asĤ =Ĥ s +Ĥ b +Ĥ int , with the three terms on the right-hand side describing the modified system including the renormalization energy, the bath, and the system-bath coupling, respectively. To be specific, the bath is harmonic,Ĥ b = j ω jâ † jâ j , and the coupling to the system via a system operatorL is linear,Ĥ int =L j c j (â † j +â j ). Such a linear bath can be exactly solved and its dissipative influence on the system is encapsulated by its correlation function α(t) = ∞ 0 dωJ(ω) coth βω/2 cos(ωt) − i sin(ωt) which is in turn determined by the spectral density function J(ω) = j c 2 j δ(ω − ω j ). Here, β is the reciprocal temperature rescaled by the Boltzmann constant and is the reduced Planck constant.
The scheme illustrated in section 2 can be straightforwardly applied for simulating quantum dissipative dynamics. To illustrate the procedure, we start with the SLE derived from the stochastic decoupling via the dissipative interaction [28,32] or via the influence functional [45], The initial condition for equation (4), ρ s (0), will be chosen according to the process under investigation. In writing equation (4), it is assumed that the whole system evolves from a factorized initial density matrix ρ tot (0) = ρ s (0)ρ b (0) and further the bath starts from a thermal equilibrium state of noninteracting harmonic oscillators ρ b (0) = exp(−βĤ b )/tr b [exp(−βĤ b )]. The system's reduced density matrix is produced with the ensemble averageρ s (t) = ρ s (t) . In equation (4), ζ and ν are zero-mean noises with non-vanishing covariances, where α r/i is the real/imaginary part of the bath correlation function. There exist many ways to decompose the noises so as to realize the covariances specified in equation (5). Such freedom has enabled the development of noise reduction methods, which achieve a substantial reduction in the number of required samples by several orders of magnitude [46,47]. In this work, we adopt a simple way for noise decomposition, i.e., we only decompose ζ into a principal part plus an auxiliary part. The partial average over the auxiliary noise yields the PEASLE where K = t/T ,ρ s is the new stochastic density matrix whose average directly produces the reduced density matrix, i.e.,ρ s (t) = ρ s (t) ξ,ν , and ξ is the principal noise with a zero mean and non-vanishing correlations, Equation (6) is the central result of this work. Compared to the conventional SLE of equation (4), equation (6) is an integro-differential stochastic equation involving functional derivatives. In this sense PEASLE generalizes NMQSD from a wavefunction representation of states to a density matrix representation. The key difference between the two methods lies in the statistical features of the involving noises. Specifically, the noises in PEASLE are piecewisely correlated, which allows for simulations based on piecewise ensemble averages. Consequently, the original difficulty concerning the long-time ensemble average is eliminated, while the remaining challenge is the evaluation of the functional derivatives with respect to stochastic noises. The above derivation thus reveals that the calculation of the functional derivatives in NMQSD actually allows a piecewise ensemble average to speed up the stochastic simulations.
Admittedly, the difficulty of taking the long-time average in the conventional SLE is converted into the calculation of the functional derivatives in PEASLE. The central task now is to determine the functional derivatives, which is equivalent to the derivation of the quantum master equation, and analytical results are only available for a few cases [25,29]. For general systems the function derivatives can be approximated by perturbative expansion [48]. To go beyond the perturbative regime the hierarchical approach is applicable when the correlation function can be well approximated with few exponentials [34]. Note that because the correlation function is a decaying function, it can always be expanded in a series of exponentials via direct curve fitting [49], or with the Padé decomposition as well as pole expansion of its spectrum [50]. For simplicity of formulation and presentation, we assume that the real and imaginary parts of the bath correlation function are of single exponentials, i.e., α(t) = κ r e −γ r t + iκ i e −γ i t . The extension to multi-exponential cases is straightforward and will not be present here. To proceed, define auxiliary density matrices where · ξ,ν;s represents the partial average with respect to the noises ξ(τ ) and ν(τ ) for τ s. Thus, the functional integral in equation (6) is the sum ρ 0,1 (t; T K ) + ρ 1,0 (t; T K ). Functional calculus leads to the stochastic hierarchical equations, The initial values of ρ m,n (t; s) at t = s are determined by solving the following hierarchical equations: Therefore, with t = s equation (8) is nothing but a reexpression of the HEOM auxiliary density matrices in terms of functional derivatives [34]. The initial condition is ρ 0,0 (0; 0) = ρ(0) and ρ m,n (0; 0) = 0 for (m, n) = (0, 0).
Equations (9) and (10) dictate the evolution of auxiliary density matrices ρ m,n (t; s). In practical implementation we first integrate equation (10) to obtain ρ m,n (T k ; T k ) (k = t/T ), and then use them as the initial condition of equation (9) for propagating the stochastic matrices ρ m,n (t; T k ) from t = T k to t = T k+1 .
The above procedure is a similar to the hierarchy of pure state (HOPS) approach [37], and the computational cost of HOPS and PEASLE is thus alike concerning the propagation of individual trajectories, which is about two to four orders of magnitude larger than that of the conventional SLE. Further, although equation (9) is formally analogous to a hybrid stochastic-hierarchical formalism [36], it is still a pure stochastic method because environmental fluctuations are directly represented by noises, while the deterministic hierarchical equations are used only for computing the functional derivatives.

Numerical simulations
We now demonstrate the outstanding efficiency of PEASLE by simulating the zero-temperature dissipative dynamics of the spin-boson model, i.e.,Ĥ s = Δ 2 σ x andL = σ z , where σ x and σ z are the Pauli matrices and Δ is the tunneling amplitude between the two states σ z = ±1. The spin-boson model has diversified applications in physics and chemistry, such as quantum optics, electron transfer reactions, the Kondo problem, and quantum computing [2]. Despite its simplicity and importance, the spin-boson model cannot be solved exactly and serves as an ideal testing platform for methodology development [51]. The initial condition is ρ s (0) = (1 + σ z )/2. The bath spectral density function is J(ω) = π 2 αωω 4 c (ω 2 c + ω 2 ) −2 where α is the Kondo parameter characterizing the dissipation strength and ω c the high frequency cutoff of the bath. With such a spectral density function the imaginary part of the correlation function is α i (t) = π 2 8 ω 3 c t e −ω c t and the real part can be approximated with exponentials α r (t) ≈ j γ j e −κ j t [49]. The imaginary part of the correlation function is included exactly, while the real part is approximated with five exponentials, which are extracted from reference [49] and specified in table 1. In this simulation we will investigate the crossover point from coherent to incoherent motion at α = 0.5 together with a relatively large frequency cutoff ω c = 10Δ. The hierarchy is truncated to the seventh tier and stochastic equations in PEASLE and HOPS are integrated with a third-order strong scheme [52] together with a step size of Δt = 0.005/ω c and T = 5Δt. For efficiency verification the trace distances between the HEOM [ρ heom (t)] and the stochastic [ρ ss (t)] results are adopted as the errors, i.e., The direct simulations for such strong dissipation go beyond the ability of conventional stochastic methods. This point can be analyzed with a limiting case Δ/αω c → 0 in which the variance of the quantum expectation σ z can be approximated as  Because the noise ν cannot be purely imaginary, the covariance of its real part α νν (t 1 − t 2 ) = Rν(t 1 )Rν(t 2 ) is semi-positive definite. If α νν (t) decays faster than 1/t at t → ∞, the double integral in equation (12) approaches to Γt where Γ = ∞ 0 dtα νν (t). Then Var( σ z ) asymptotically increases as e 4Γt . Simulations, shown in figure 2, have verified the above analysis. With 3 × 10 7 trajectories the conventional SLE yields reliable results only for Δ · t < 1.2. The HOPS method, compared to SLE, can only produce 10 5 trajectories with the same computation time, but reduces the errors by more than ten orders of magnitude. Despite the drastic improvement we still witness large errors and their steady long-time increase, which may be unfavorable for simulating long time dynamics. By contrast PEASLE can effectively simulate the dynamics with such a dissipation strength and the errors with 10 5 trajectories are always less than 1.2 × 10 −3 . The errors of PEASLE are two to three orders smaller than that of HOPS, which implies a 4-6 orders of efficiency enhancement. More importantly, the errors, as shown in figure 2(b), are saturated after a fast initial increase, which is highly appreciated for simulating long time dynamics.

Summary
In this work we have rigorously derived a stochastic Liouville equation with piecewisely correlated noises and the ensemble average can therefore be performed independently for different time intervals. By doing so we can successfully avoid the long-time stochastic average and the statistical errors are free from long-time increase. The above examples show that PEALSE is very efficient for simulating long-time dynamics with long memory and strong dissipation.
Here we have developed the PEASLE approach by starting from the conventional SLE. However, PEASLE provides a new direction to enhance the simulation efficiency of stochastic methods and can also be started from the hybrid stochastic-hierarchical equations, NMQSD, the finite-memory propagation, or the transfer tensor method. Further combination with these schemes may stimulate more efficient, versatile methods for studying the quantum Brownian motion.