In the inquiry of quantum open system theory, the pursuit of complete and exact system dynamics, particularly those that incorporate the backflow of information from the environment, remains a paramount and universal research topic [1, 2]. Memory effects in structured environments can be significant in the general dynamics of open quantum systems, rendering the Born and/or Markov approximations ineffectual [3-7]. The systemic evolution in this case corresponds to non-Markovian dynamics whose untapped potential in quantum information science has been realized gradually in recent years [8]. Given the potential of the backflow to counteract environment-induced dissipation and decoherence, a series of representative schemes playing the advantages of non-Markovian dynamics are proposed, for examples, entanglement (coherence) preservations [9-13], precise detections [14, 15] and ground state cooling [16, 17]. Recently, the backflow effect has been recognized as a controllable resource in certain noise-induced quantum control schemes [18-20].

The burgeoning interest in applying non-Markovian dynamics has prompted the pursuit of a universal theoretical framework capable of addressing the complexities of these systems. For practical systems, the evolution is usually unmanageable due to the memory effects, which necessitate the recording of the system's entire history. Few common models can be solved analytically so that numerical simulation of non-Markovian dynamics has also become an important and challenging issue rapidly for years [21, 22]. Until now, non-Markovian dynamic theories have been well formulated in the Schrödinger picture, such as time-convolutionless [23, 24] and Nakajima-Zwanzig master equations [25, 26], Feynman-Vernon influence functional techniques [27, 28], and so on. Corresponding simulation methods have been proposed, employing stochastic simulation techniques with tailored distributions and quantum jump processes [29, 30]. The algorithms have been improved, for faster convergence and higher accuracy, by multitudinous follow-up researches in recent years [31-39].

Although above stochastic simulation approaches are based on the generalized quantum open system theory, it still remains a difficult problem because the computing difficulty will be further upgraded once the system complexity is increased. In quantum information processing (QIP), practical systems often have specific requirements, comprising multiple subsystems, with distinctive structures [5], prohibiting rotating wave approximation (RWA) [17, 40] as well as coupling environment with finite temperature [16, 17], and so on. To solve the state evolution of complex systems, simplifying approximations are often necessary, which inevitably compromise precision [16, 41]. This has redirected research efforts towards exploring quantum Brownian motion within the Heisenberg picture [2, 42]. Only when the observable quantities are considered or the quantum states are confined to a particular framework can exact Heisenberg-Langevin equations be deduced very simply, even for complex quantum systems [14, 17, 43-46].

The operator differential equations have been simulated by translating them into c-number equations with designed stochastic processes [17, 46-49]. In this paper, we simplify this idea so that the mean value and covariance matrix of a continuous variable (CV) Gaussian system can be calculated even if the systemic Hamiltonian and spectral density of the system are arbitrary. Algorithm complexity is almost independent of the structure of the environment, system-environment interaction and the dynamics of the system itself, as the only mathematical deduction required is to obtain the memory kernel function. Correspondingly, other relations and approximations, even the noise correlation functions, are not necessary. Unlike previous works introducing stochastic terms and corresponding integral terms into differential equations [47-50], here all the c-number dynamic equations can be determined with a further simplified algorithm, and the only random elements are the corresponding initial values. In order to demonstrate the validity directly, the accuracy analysis is proved by calculating a linear model and comparing the results with those given by exact master equations [28, 51].

In the Heisenberg picture, the non-Markovian dynamics of an open quantum system can be described by a set of Heisenberg-Langevin equations of the system operators ${\hat{o}}_{i}$. Usually, the standard approach for obtaining system dynamics is to use Green's function technique [8]. Therefore Langevin equations are formally written as The above Langevin equation obtained by solving Heisenberg equations correspond to environmental degrees of freedom formally. ${ \mathcal U }({\hat{o}}_{i}(t))=-{\rm{i}}[\hat{H},{\hat{o}}_{i}(t)]+{ \mathcal G }({\hat{o}}_{i}(t))$ characterize the unitary evolution parts of the system dynamics, and the first term is decided by the systemic Heisenberg equations and the second one is an environmental induced renormalization term. The environmental-system interaction induces the non-Markovian memory effect by taking into account the time nonlocal integral kernel ${ \mathcal K }(t-\tau )$. $\hat{\xi }$ is the quantum input noise operator caused by the environment free energy.

In general, the environment is regarded as a bath or a thermal reservoir so that $\langle \hat{\xi }\rangle =0$ holds [2], which enlightens us to the that some ${\hat{o}}_{i}$ whose expectation values $\langle {\hat{o}}_{i}\rangle $ are not affected by the noise operator can always be found. Here, $\langle$. . . $\rangle$ refers to taking the expectation value with respect to the density matrix of the quantum system. This condition seems harsh, fortunately, the eigen-operators of harmonic oscillator satisfy this requirement in some typical non-Markovian dynamic models. For example, we consider a harmonic oscillator coupled, via a linear beam splitter or Brownian motion interaction, to a noninteracting bosonic environment. Then its eigen-operators, that is, creation and annihilation operators, satisfy It encourages us to expand the relevant quantum operators as the sum of their respective expectation values and corresponding fluctuation operators, that is, ${\hat{o}}_{i}=\langle {\hat{o}}_{i}\rangle +\delta {o}_{i}\equiv {o}_{i}+\delta {o}_{i}$. In a linear system, this processing is rigorous and it is also a commonly used approximation for systems with weak nonlinearity [11, 17, 48, 52, 53]. In nonlinear systems, this processing is inaccurate. The Wigner function exhibits negative values only when the nonlinear system energy is of the order of a few quanta. Fortunately, when the dimension of the Hilbert space is very large, the Wigner function of the system is typically positive. In the case of a Wigner function that is always positive, it can be mapped to a classical distribution, and a c-number Langevin equation can be used to reconstruct this distribution. In the case of nonlinear systems, the radius canonical equation is highly consistent with the full quantum [54]. The quantum state of a CV system is generally considered as a Gaussian state whose all quantum properties can be obtained by expectation values and corresponding covariance matrix of quadrature operators, so the non-Markovian dynamic evolution can be calculated by the statistic result of multiple quantum trajectories [54]. In particular, the initial state ρ₀ of the system is simulated by the Gaussian complex numbers ${o}_{i}^{j}(0)\unicode{x0007E}{ \mathcal N }({o}_{i},\delta {o}_{i})$ with ${o}_{i}=\,\rm{Tr}\,(\hat{{o}_{i}}{\rho }_{0})$ and $\delta {o}_{i}={(\,\rm{Tr}\,({\hat{o}}_{i}^{2}{\rho }_{0})-{(\,\rm{Tr}\,({\hat{o}}_{i}{\rho }_{0}))}^{2})}^{1/2}$. Note that such kind of complex numbers can be achieved by ${o}_{i}^{j}=({{ \mathcal Z }}_{1i}^{j}+{\rm{i}}{{ \mathcal Z }}_{2i}^{j})/\sqrt{2}$ with ${{ \mathcal Z }}_{1i,2i}^{j}\unicode{x0007E}{ \mathcal N }({o}_{i},\delta {o}_{i})$ and ${{ \mathcal Z }}_{1i,2i}^{j}\in { \mathcal R }$. Then every trajectory satisfies a c-number equation where superscript j denotes the jth trajectory. Hence, the expectation values and corresponding covariance matrix in non-Markovian dynamics at arbitrarily time t can be obtained by and Equations (3) and (4) are equivalent to equation (2) if the total Hamiltonian including system and environment is linear. For the system with weak nonlinear interaction in $\hat{H}$, the statistical simulation is thought to be, proved by [47], a more accurate technique compared with the linearization method. Note that the right term in equation (5) loses its non-commutativity because ${o}_{i}^{j}$ and ${o}_{k}^{j}$ are already c-numbers. Therefore, commutation relation $[{\hat{o}}_{i},{\hat{o}}_{k}]$ is also needed as a modification to calculate terms like $\langle {\hat{o}}_{i}(t){\hat{o}}_{k}(t)\rangle $.

A key task for adopting this approach to investigate non-Markovian dynamics is to simulate color noise $\hat{\xi }$ as a c-number ξ^j with some randomness. To be more specific, the beam splitter or Brownian motion interaction with a harmonic oscillator can be generally expressed by the Hamiltonian [27, 28] where the operators $\{\hat{b},{\hat{b}}^{\dagger },\hat{q}\}$ are annihilation, creation, position operators of the system, and the operators $\{{\hat{b}}_{k},{\hat{b}}_{k}^{\dagger },{\hat{q}}_{k},{\hat{p}}_{k}\}$ are annihilation, creation, position, momentum operators of the kth reservoir mode with the eigenfrequency ω_k. ν_k, ${\nu }_{k}^{* }$ and Υ_k are the coupling coefficients between the system and reservoir. Corresponding to the above Hamiltonian, the input color noises can be deduced as [14, 43, 46]

Here a simple replacement ν_k = ω_kΥ_k is applied in equation (7c). Considering the initial state of a bath or a reservoir is usually a thermal equilibrium state with Gaussian distribution, so ${\hat{b}}_{k}$ is also simulated as ${b}_{k}^{j}\sim { \mathcal N }(0,\sqrt{{\bar{n}}_{bk}+0.5})$ since $\langle {\hat{b}}_{k}^{\dagger }{\hat{b}}_{k}+{\hat{b}}_{k}{\hat{b}}_{k}^{\dagger }\rangle =2{\bar{n}}_{bk}+1$ with ${\bar{n}}_{bk}={[\exp (\hslash {\omega }_{k}/{k}_{{\rm{B}}}T)-1]}^{-1}$ [47]. Following the above general procedure, ${\hat{q}}_{k}(0)$ and ${\hat{p}}_{k}(0)$ can also be simulated, so that the only pending problem is coupling coefficients between the system and each reservoir mode. Previous statistical and experimental researches have pointed out that the environment can be specified by its spectral density ${ \mathcal J }(\omega )={\sum }_{k}{\nu }_{k}{\nu }_{k}^{* }\delta (\omega -{\omega }_{k})$ and the following relation exists: for the continuous distributed spectrum and environment.

The dynamic equations of expectation values and second-order mechanical quantities possess U(1) symmetry under the transform ν_k → ν_ke^iθ, meaning that the system dynamics just depend on |ν_k|² and one can always assume ${\nu }_{k}\in { \mathcal R }$, which is strictly established in a linear system, for convenience. Therefore, the interaction coefficients can be obtained by the spectral density: where ∆ω = ω_i+1 - ω_i is the frequency difference when we divide the continuous frequency distribution as discrete distribution. In equation (9), all the information in the spectral density is assigned to the coupling coefficient, implying that the density of the reservoir itself should be a uniform distribution with ∆ω = ω_cut/M, where ω_cut is a large enough cut-off frequency in order to handle the integral ${\int }_{0}^{\infty }$ and M is also a large enough number whose physical meaning is the total number of environmental modes. So far, the color input noise, corresponding to ${\hat{\xi }}_{{\rm{bs}}}$ in equation (7a) as an example, can be expressed in a discrete finite difference form: Hence, a general non-Markovian dynamic can be simulated by connecting equations (3), (4), (5) and (10) with two additional parameters ω_cut and M.

By now, it is evident that our simulation approach disallows the presence of any initial correlations within the system, including among the subsystems or between system and reservior. Fortunately, this restriction is not overly limiting, as the system always starts its evolution from a separable equilibrium state in an actual physical process. Beyond this constraint, our simulation is not limited by other restrictions. In especial, the spectral density and the system-reservoir interaction can be arbitrary here.

To show the accuracy of the non-Markovian simulation and give some more detailed discussions, we need to consider and analyze a general dynamic process as contrast. The condition ’general’ here has two meanings: (i) The system contains both Markovian and non-Markovian environments simultaneously. (ii) The system is influenced by a reservoir with a finite temperature. To satisfying these three requirements simultaneously, a normal mechanical system is selected as our model, as shown in figure 1. The mechanical oscillator is coupled to a general non-Markovian reservoir. Then the Hamiltonian of the system can be written as ${\hat{H}}_{s}\,={\omega }_{m}{\hat{b}}^{\dagger }\hat{b}$ [55, 56].

Here $\hat{b}$ (${\hat{b}}^{\dagger }$) is the annihilation (creation) operator of the mechanical mode with corresponding frequency ω_m. The environment Hamiltonian is selected as ${\hat{H}}_{{xx}}$ in equation (6b), so that the Heisenberg-Langevin equations of motion for the annihilation operators of the system are given by [11, 53] By solving ${\hat{b}}_{k}$ as we have where $f(t)=2{\rm{i}}{\sum }_{k}{\nu }_{k}^{2}\sin ({\omega }_{k}t)=2{\rm{i}}{\int }_{0}^{\infty }{\rm{d}}\omega { \mathcal J }(\omega )\sin (\omega t)$.

If we do the RWA on the system-reservoir coupling terms, the system will degenerate into a simple model with a solitary oscillator coupled to a reservoir with a beam splitter interaction ${\hat{H}}_{{\rm{bs}}}$. Then the corresponding integro-differential Heisenberg equation is where $f(t)={\int }_{0}^{\infty }{\rm{d}}\omega { \mathcal J }(\omega ){{\rm{e}}}^{-{\rm{i}}{\omega }_{k}t}$. Equations (13) and (14) are the basis for the following calculations.

To be more specific, let us first examine the non-Markovian dynamics. By substituting equation (14) into equation (3), the trajectory equation corresponding to this model is and it can be solved numerically by some mature algorithms such as the Runge-Kutta algorithm or linear multistep algorithm. For convenience, the reservoir temperature here is characterized by ${\bar{n}}_{b}({\omega }_{m})\equiv {n}_{b0}={[\exp (\hslash {\omega }_{m}/{k}_{{\rm{B}}}T)-1]}^{-1}$, that is, the mean phonon number at the oscillator frequency. So the phonon distribution used to produce ${{ \mathcal Z }}_{1,2k}^{j}$ can be further obtained as

Auxiliary parameters ω_cut and M are restricted by the spectral structure under the reservoir. Corresponding to the Brownian motion interaction, the common used spectral density ${ \mathcal J }(\omega )$ is of the following form [57]where η is a dimensionless coupling constant and ω_s is the cut-off frequency. The parameter s classifies the environment according to the ω dependence of in the low-frequency region. Specifically, 0 < s < 1, s = 1 and s > 1 correspond to the sub-ohmic, the ohmic and the super-ohmic baths, respectively. In figure 2, we show the high frequency mode of the reservoir has little contribution to the systematic dynamics because the high frequency part of the Ohmic-type spectrum and the phonon number distribution approach to zero. The integral truncation is therefore set as ω_cut≥20ω_s in the following simulations. Moreover, figure 2(b) shows that M ≥ 1000 is large enough for generating an approximate continuous distribution. The insets in figure 2(b) illustrate that the fluctuation of only 100 statistical results is almost consistent with the analytic phonon number distribution, and the curve becomes smoother by adding the statistical results to 3000.

Besides the stochastic statistical method, this system can also be solved exactly by using Zhang's theory, and the exact master equation reads [28]: where $\hat{\tilde{H}}(t)=\tilde{\omega }(t){\hat{b}}^{\dagger }\hat{b}$ is the renormalized Hamiltonian. The renormalized frequency, time-dependent dissipation coefficient and the fluctuation coefficient are given by where the functions $u(t,{t}_{0})=\langle {[\hat{b}(t),{\hat{b}}^{\dagger }({t}_{0})]}_{\mp }\rangle $ and $v(t,t)\,=\langle {\hat{b}}^{\dagger }(t)\hat{b}(t)\rangle $ are related to the nonequilibrium Green's functions of the system in the Schwinger-Keldysh nonequilibrium theory [58], and u(t, t₀), v(t, t) satisfy the integro-differential equations of motion with the boundary conditions u(t₀, t₀) = 1 and v(t₀, t) = 0 under t₀≤Γ≤t. Here the self-energy correction $f(\tau ,\tau ^{\prime} )\,\equiv f(\tau -\tau ^{\prime} )$ has the same form with that in equation (14). Therefore $\tilde{f}(\tau ,\tau ^{\prime} )$ is expressed as: In this expression, $\bar{n}(\omega )$ is the Bose-Einstein distribution mentioned above. The boundary condition v(t₀, t) = 0 ensures that v(Γ, t) can be solved analytically as Recently, the exact analytic solution of u(Γ, t₀) is also reported in [28, 51] expect v(Γ, t). Thus, the system can be calculated exactly with the help of equations (20) and (22).

To check the validity, we first show intuitively some contrast results by calculating main physical observables with two methods. Physical observables $\langle \hat{b}\rangle $ and $\langle \hat{b}\hat{b}\rangle $ can be calculated basing on the arithmetic means of trajectories b^j and ${({b}^{j})}^{2}$. The phonon number modified by commutation relation is also obtained by $\langle {\hat{b}}^{\dagger }\hat{b}\rangle =\sum | {b}^{j}{| }^{2}/N-0.5$. On the other hand, substituting $\langle \hat{o}\rangle (t)=\,\rm{Tr}\,(\hat{o}\rho (t))$ into the master equation (18), the relations between physical observable and Green's function can be found as follows: In figures 3(a)-(c), we present the time evolutions of the physical observables with cryogenic limit n_b0 = 0. The coincidence curves indicate that, for both the first-order and second-order observables, the stochastic simulation results are in good agreement with those obtained by the exact master equation at any time. Note that the second-order observables contain quantum correlations, implying that the stochastic simulation is effective even in the quantum regime. It can also be known from figures 3(a)-(c) that the simulation method will not lose its validity and accuracy even though the coupling strength of system-reservoir interaction increases from the weak interaction regime η = 0.1, across the strong interaction boundary η = 0.3 to the critical value η = ω_m/[ω_sΓ(s)] [59]. Finally, we want to explain that all above-mentioned properties are still true even under higher reservoir temperature. The time evolutions of $\langle {\hat{b}}^{\dagger }\hat{b}\rangle $, since other observables are not affected by temperature, are plotted in figure 3(d) corresponding to n_b0 = 12 and the consistency between the simulation result with that of the exact solution can also be found as we expect.

The accuracy of the simulation is further discussed by considering its N dependence. In figure 4, we show that the result by averaging over N = 100 runs deviates from the exact dynamics in a long-time scale. By increasing the number of trajectories to 1000 and 10000, better results appear, such as the black and red lines shown in figure 4(a). To further validate the accuracy with quantitative discussion, the difference of state evolutions according to the two methods is measured here by the Gaussian fidelity [60-62]: The detailed explanations of its definition and computing method will be provided in the appendix. The inset in figure 4(a) shows the fidelity can be higher than 96%, 99.8% and 99.98% corresponding to N = 100, N = 1000 and N = 10000, respectively. The relation for increasing N to improve simulation accuracy can be further explained intuitively by figure 4(b) in which the time-averaged $\bar{{ \mathcal F }}={\int }_{2}^{10}{ \mathcal F }(t){\rm{d}}t/8$ is plotted. Here $\bar{{ \mathcal F }}$ presents a rapid increase in the interval N ∈ [200, 3000] and its trend becomes stationary, leading $\bar{{ \mathcal F }}\gt 99.9 \% $ to be always satisfied. This is the reason why we set N = 3000 in previous simulations.

Besides the expectation value and the average phonon number, some more complex physical observables can also be calculated with simple statistics. As an example, we calculate the two-time correlation function $\langle {\hat{b}}^{\dagger }(t)\hat{b}(t+\tau )\rangle $, an important component of the non-Markovianity measure in recent work [63], under a semi-classical approximation: and compare the simulation results with an exact two-time correlation It can be seen from figure 5(a) that for the two-time correlation function, the result given by the stochastic simulation has a small error at a short-time scale, but it will be in good agreement with that obtained by the exact expression at a long-time scale. Physically, the neglected commutation relations are weakened by the dissipation effect so that the semi-classical approximation becomes more precise after a long-time evolution. It is also worth noting in figure 5(a) that the error of semi-classical approximation can also be eliminated by a larger initial energy.

One may wonder if the accuracy depends on the spectral density selected. To clarify this point, we finally plot figure 5(b) to show the good consistency between the simulation result and the exact solution. It illustrates that the non-Markovian dynamics precisely correspond to both sub-ohmic and super-ohmic reservoirs. To sum up, we have already proved that our simulation method can be in good agreement with the exact solution.

In summary, we have proposed a general simulation approach to calculate Heisenberg-Langevin equations corresponding to non-Markovian dynamics. This method simulates observable quantities, defined by their mean values and covariance matrix, using definitive c-number equations and stochastic initial conditions. Although this idea has its limitation, that is, the algorithm can only handle the Gaussian state without any initial correlation, it still leads to the following advantages compared to previous works: (a) Determined differential equations allow us to adopt a simple numerical simulation method in order to avoid solving the randomized procedure. (b) In the process of simulation, the only additional derivation is to obtain memory kernel function. (c) The algorithm is of expansibility and the complexity is almost independent of the structure of environment, system-environment interaction and the dynamics of the system itself. These three merits ensure our simulation approach can calculate many of the common non-Markovian systems in recent QIP schemes, even for the special spectral structure in experiment. We believe that this algorithm can be extended to more complex quantum optomechanical systems with nonlinear effect and applied for spectral analysis.

The validity and accuracy analysis are discussed by calculating a linear model and comparing the results with those given by exact master equation. In particular, the time-averaged fidelity satisfies $\bar{{ \mathcal F }}\gt 99.9 \% $ with 3000 trajectories. Given the scalability of our algorithm, we are convinced that our approach can more accurately simulate complex nonlinear effects and Brownian motion. We believe that this unsurpassable advantage holds extensive potential applications in QIP domain, such as the calculations of the cooling and measurement limitation without being bound by the approximate conditions.

In the last part of our paper, we want to give a reasonable prediction, that is, if the system has an initial correlation but not a system-environment correlation, it can still be simulated by adopting our approach with a prior-simulation to prepare the initial states from the uncorrelated ones. This open problem can be further improved by the follow-up researches.