Stochastic Schrödinger equation derivation of non-Markovian two-time correlation functions

We derive the evolution equations for two-time correlation functions of a generalized non-Markovian open quantum system based on a modified stochastic Schrödinger equation approach. We find that the two-time reduced propagator, an object that used to be characterized by two independent stochastic processes in the Hilbert space of the system, can be simplified and obtained by taking ensemble average over one single noise. This discovery can save the cost of computation, and accelerate the converging process when taking the average over noisy trajectories. As a result, our method can be widely applied to many open quantum models, especially large-scale systems and extend the quantum regression theory to the non-Markovian case. In the short-time simulations, it is observed a significant difference between Markovian and non-Markovian cases, which can be applied to realize the environmental spectrum detection and enhance the measurement sensitivity in varying open quantum systems.

www.nature.com/scientificreports/ results of Lindblad master equations. Therefore, it is natural to extend the SSE approach to calculate the two-time correlation functions (TTCF) or multi-time correlation functions (MTCF).
In the system-environment framework, the interaction Hamiltonian can be generally described by the Lindblad coupling operator L, and the total Hamiltonian is where H S(R) is the Hamiltonian of the system (reservoir). b k and b † k are the annihilation and creation operators of the bosonic environment, is a suitable perturbation parameter that can take one. g k (g * k ) and ω k are coupling constants and the eigen frequency of the kth mode of the environment, respectively.
In this paper, we derive a general formula for two-time correlation functions (TTCF) in the non-Markovian regime, where the TTCF is defined as where A and B are two arbitrary operators in the Hilbert space of the system, and t > t ′ is always satisfied. Moreover, we assume that system and environment are decoupled at t = 0 , and the initial density matrix can be rewritten in the form of a pure state ρ tot (0) = |ψ 0 ��ψ 0 | ⊗ |0��0| , where |ψ 0 � and |0� are the initial states of the system and the environment respectively. In an OQS, there are two major sources of noise, quantum fluctuations and thermal fluctuations. Recent experiments have successfully cooled the system down to the few milli-Kelvin regime and suppressed the thermal noises 43 . For simplicity, we study a zero-temperature environment and focus on the quantum fluctuations only.
The aim of the paper is twofold. Firstly, we derive a non-Markovian stochastic Schrödinger equation (SSE) to characterize the dynamics of TTCFs in accordance to the generalized system-environment Hamiltonian. The derived evolution functions for quantum mean values can be analytically governed by a single Gaussian noise. The subsequent process of numerical simulation for the ensemble average of all possible trajectories can be significantly simplified. Secondly, we investigate several typical applications of TTCF, e.g. two-level systems, and second-order correlation function g (2) . In particular, we study how a TTCF is influenced by various non-Markovian environments, by adjusting the eigen frequencies of system and coupling strength. We discuss the physical mechanism in the model and the experimental feasibility of our proposed model.

Background
In this section, we briefly review the stochastic Schrödinger equation (SSE) approach and compare it with the quantum regression theorem (QRT).
Stochastic Schrödinger equation approach. The evolution of a general quantum open system (QOS), characterized by the Hamiltonian (1), is simply governed by the Schrödinger equation Since the reservoir consists of a large number of bosonic oscillators, H R = k ω k b † k b k , so it is natural to choose the Bargmann coherent state basis to represent the environment, |z� = ⊗ k |z k � , a tensor product of all the environmental oscillators, where |z k � is defined as b k |z� = z k |z k � for the kth mode in the environment. Therefore, the state vector of the system in the Bargmann basis, |ψ t (z * )� = �z|�(t)� , can be expressed as a stochastic trajectory. In addition, its evolution, in the interaction picture, is governed by the following equation (setting = 1) where the complex Gaussian noise z * t = −i k g k z * k e iω k t , satisfying the auto-correlation function and dµ(z)|z��z| = I . Moreover, the evolution operator of system-environment U(t), defined as |�(t)� = U(t)|�(0)� , can be explicitly shown in the interaction picture U(t) = e iH R t e −iH tot t , which obeys the evolution equation For the cases that the system and environment are entangled initially, we can decompose the density matrix into a product state basis. Therefore, for simplicity, we assume the system and the environment are separate at t = 0 , and the environment is in the vacuum state, |�(0)� = |ψ 0 � ⊗ |0� . It is essential to calculate the term, Here, we compare two equivalent methods, the Heisenberg equation approach and the quantum-state diffusion approach, to compute it. Firstly, let us consider the observable b k (t) = U † (t)b k U(t) in the Heisenberg picture and its evolution equation acquires (6) �z|b k |�(t)� = �z|b k U(t)|0�|ψ 0 �. where the correlation function is α(t, s) = k |g k | 2 e −iω k (t−s) for the zero-temperature environment. Secondly, we apply the chain rule of the functional derivative of the Gaussian stochastic process z * t and the Eq. By comparing the results of Eqs. (10) and (12), the initial condition of O operator must satisfy SSE approach can be used to solve many OQS models, including discrete qubit systems, continuous variable systems, and multi-layer mixes systems in the presence of non-Markovian environments 44,45 . It can also be extended to deal with both zero-temperature and finite temperature environments.
Quantum regression theorem. Quantum regression theorem indicates that one-time expectation values of system operators and TTCFs (MTCFs) have the same form of evolution equations [25][26][27] . To compute the expectation values, we introduce the Liouville operator with respect to time, L(t) and get the reduced density matrix as It is important to mention that this is assuming that there are no non-linear dynamics involved in the environment. From this equation, we can find the mean value of a system observable A is determined by For the two-time correlation functions, we keep the first operator as A(t) and the second operator is denoted as B(t ′ ) . Assuming the time t ′ is fixed and t > t ′ , the TTCF is expressed as Under the Born-Markov approximation, ρ tot (t) is always given by ρ S (t) ⊗ ρ R (0) , and consequently, ρ tot,B (t ′ , τ ) = Bρ S (t ′ + τ ) . It indicates that, for a given t ′ , the new reduced density matrix ρ tot,B evolves in the similar manner as ρ S , Therefore, the TTCF is downgraded to a single operator mean value, that www.nature.com/scientificreports/ The SSE approach can derive the same result. Formally, the reduced density matrix of the system can be obtained by taking ensemble average over all possible trajectories, either analytically or numerically,

Consequently, the formal master equation can be expressed as
If the environment is Markovian and the correlation function, α(t, s) = Ŵ 2 δ(t, s) , the formal master equation can be expressed in the Lindblad form (the temperature is zero, for simplicity): Therefore, the evolution equation of one-time expectation values of system operators is For TTCFs, assuming the time t ′ is fixed and t > t ′ , the evolution equation is expressed Following the similar method, the evolution equation of the term P S (τ ) is In the Markoivan regime, it is easy to show that the term P S (t) follows the same operator equation form as the master equation, Consequently, we can conclude that the TTCF evolution equation can be expressed as the same as the single-time expectation value evolution equation However, from our derivations, it is clear to show that the QRT is not valid in general conditions . Particularly, when the environment is structured where the Markovian assumption is no longer applicable, the evolution equation of TTCFs needs to compute in the language of non-Markovian OQS carefully. In some previous work, the SSE approach had been applied to compute the stochastic propagator where two independent stochastic processes are employed. From our experience, the convergence speed usually depends on the type of systemenvironment coupling and the size of the system's Hilbert space. The operation that computes the ensemble average over multiple stochastic processes will increase the computational burden significantly. To solve this issue, our approach can generally calculate TTCFs with only applying one stochastic noise.

Method
We consider a general OQS coupled to a bosonic environment with a general Hamiltonian (1) in the interaction picture, and the corresponding SSE of the system's state vector, |ψ t (z * )� = �z|�(t)� , is The evolution functions of the TTCFs are formally expressed as Here, we need to point out that U(t)U † (t ′ ) � = U(t − t ′ ) in the interaction picture. And we define the two terms www.nature.com/scientificreports/ which can be simulated independently with the same set of stochastic noise z * t . To contact with the previous work, particularly comparing with the method in the ref. 18 , we insert the identity of a second independent noise w * t into |η(t 1 , t 2 , z * )� and turn it into where dµ(w) = d 2 we −|w| 2 /π and the term �z|U(t)U † (t ′ )|w� is a stochastic propagator. Once being expanded into two independent noise-basis, every element of the propagator can be independently characterized as a stochastic function of noises z * t and w * t . However, the computation complexity is also doubled. And the consequent TTCF (29) evolves For those simple OQSs, e.g. consisting of few qubits, this approach is helpful in calculating TTCFs and MTCFs. Our approach will only employ one stochastic noise z * t to simulate the TTCF, which can be easily extended to a complicate-structured quantum model. Applying the similar derivation in the Eq. (12), we arrive at the timeconvolutionless evolution equation of η(t, t ′ , z * )�, where the functional deferential chain rule in Eq. (11) is applied again. Compute the last term of the Eq. (32), In the first term of the above expression, it should be noted that b k (t) = U(t)b k U † (t) cannot apply the analytical solution of b k (t) in Eq. (8) directly. In order to apply it, we need to slightly adjust the form of b k (t) to where the term e iω k (t+s) is zero, due to the rotating wave approximation. Consequently, the above expression can be simplified as After repeating the above approach, we can obtain the complete form of the evolution equation of the term |η(t, t ′ , z * )�, Short-time limit. In the short time-scale, a regime where non-Markovian effects are dominant, we consider the t ′ = 0 condition and the subsequent TTCF is simplified and expressed as In the following derivation, set the parameter = 1 for simplicity. In accordance with the short-time assumption, the evolution equation of |η� in the Eq. (36) can be further simplified as We notice that the evolution equation of |η(t, z * )� has the same operator form as |ψ t (z * )� , in the Eq. (10). In addition, the evolution equation can be alternatively written in the form of SSE, with the ansatz O operator. Comparing the two equivalent approaches, we find out that   12), with the same initial value O(t, s = t) = L and the same operator form evolution equation (13). But |η(t, z * )� and |ψ t (z * )� still characterize two different sets of trajectories, due to the different initial values of |η(t = 0, z * )� = B|ψ 0 � and |ψ 0 � respectively. Now, we can conclude that the TTCF in the short-time limit can be computed by taking the ensemble average over two trajectories where R(t) = |η(t, z * )��ψ t (z)| is the equivalent stochastic density matrix. Similar to the master equation approach, M(R(t)) can be determined by an evolution equation Here, the subsequent evolution equation of the short-time limit TTCF is Generally, the O operator can be resolved in a sum of operator basis with the ascending order of noises. From the solution, we can prove that QRT is not valid in a general OQS. However, in the Markovian limit, the O operator equals to the Lindblad operator L, and the correlation function α(t, s) = Ŵ 2 δ(t, s) . Consequently, the operator M(R(t)Ō † ) can be explicitly expressed as Ŵ General discussion. In order to study the genreal TTCF of OQSs in the non-Markovian regime, we derive the general formula of TTCFs from Eq. (32) by expanding the |η(t, t ′ , z * )� in a linear sum of ascending order of noise terms By substituting the expansion into Eq. (32), we can equal the terms in two sides of the equation with the same order of noise, and obtain a series of deferential equations for each noise-independent |η k �, In addition, we notice the terms containing z * t in the two sides of the Eq. (32) are always equal to each other, so that the initial conditions of |η k � can be obtained Up to now, we apply the SSE approach and employ only one noise to calculate the non-Markovian TTCFs by generating evolution equations. This method starts from the microscopic OQS Hamiltonian, and is tolerant to many factors, including the central system Hamiltonian, the system-environment coupling operator, the correlation functions of the non-Markovian environment, etc. So it can be used to systematically solve TTCFs of OQS.

Numerical results
To examine the derived TTCF equations, we consider a dissipative two-level system H s = ω 2 σ z , in which the dissipative type of coupling is L = σ − . In refs. 13,14,46 , the O operator has been discussed in detail, that O(t, ) = f (t, s)σ − , where the coefficient factor f(t, s) is governed by the initial condition f (t, s = t) = 1 and the evolution function where F(t) = t 0 dτ α(t, τ )f (t, τ ) and α(t, s) is the correlation function. The subsequent SSE and master equation are explicitly expressed www.nature.com/scientificreports/ We now calculate the non-Markovian evolution equations of TTCFs, by applying the commutation relations between the Pauli matrices. First, we consider the evolution of the mean averages �σ i �, i = (x, y, z, +, −) from the derived master Eq. (48), which results in a group of linear differential equations: Then, we calculate the TTCFs by applying the results of short-time limit (42) and general derivation respectively and compare the difference.
Case I: short-time limit. In the short-time limit, set the time index t ′ = 0 , and TTCFs can be obtained by the following set of differential equations: and similarly, In order to study the impacts on TTCFs induced by the non-Markovian environment, we use the Ornstein-Uhlenbeck style correlation function, where Ŵ is the global coupling strength and the factor γ indicates the spectral width of the coupling. In addition, 1/γ scales the memory time of the environment. Therefore, if γ takes infinity, the environment does not have memory effect and the system-environment interaction is close to the Markov limit. Note that our method works for arbitrary spectral functions and the corresponding correlation functions of the environment, including Ohmic, super-and sub-Ohmic, etc. Without loss of generality, the initial state of the system is chosen to be |ψ 0 � = (|0� + |1�)/ √ 2 . In Fig. 1, we select γ = 0.5 and 2 to simulate the non-Markovian dynamics and compare with the Markovian limit when γ = 10 . Two TTCFs, σ x (t)σ y (0) and �σ + (t)σ − (0)� , are investigated as the example. We observe considerable difference between the non-Markovian and Markovian dynamics that the non-Markovian dynamics oscillate deeper and last longer time to achieve the steady value, while the Markov dynamics quickly decay to its long-time limit value. In addition, from the microscopic perspective, the Markovian environment indicates the equal interaction strength to every mode, a flat spectral function. This can be observed in Fig. 1b and (d), that the green dotted line is very close to a flat line when the coupling is Markovian. Since the selected Ornstein-Uhlenbeck correlation function can gradually show the transition from Markov to non-Markovian regime, we witness the peaks of the Fourier spectrum are higher and sharper at the resonant location ±ω . It indicates the strong non-Markovian effects occur in the short-time regime. Moreover, some sub-peaks of the spectrum can only be observed in the non-Markovian regime.
Case II: general TTCFs. Here, we will continue to discuss the TTCF in a more general manner, where the t ′ takes a non-zero value. By expanding the noisy term |η� into a polynomial consisting of ascending orders of noise z * s k , as shown in Eq. (45), our method can be applied to systematically solve many models. For this two-level dissipative system, we have proved that the order of noise in the trajectory will not exceed 1 in our previous study. With the exact |η 0 � and |η 1 � , governed by www.nature.com/scientificreports/ and the boundary condition the noisy trajectory |η(t, t ′ , z * )� can be simulated as When t = t ′ , the initial condition is |η� = B|ψ t (z * )� . Then we can solve the TTCF �A(t)B(t ′ )� numerically stating from arbitrary t ′ . Here, we take the TTCF �σ x (t)σ z (t ′ )� as an example. In Fig. 2a, we display the evolution of �σ x (t)σ z (t ′ )� as the function of the time difference τ = t − t ′ , and choose t ′ = 5 . In Fig. 2c, the t ′ is changed to 10. The Fourier spectral functions of the two TTCFs are shown in Fig. 2b,d respectively. We also compare the dynamics in the Markovian (blue dashed line) and non-Markovian (red solid line) regimes. In the non-Markovian regime, there is a significant difference between the two TTCFs: �σ x (τ + 5)σ z (5)� (Fig. 2a) and �σ x (τ + 10)σ z (10)� (Fig. 2c), which indicates that the QRT and the time translation invariance are not valid. Particularly, in Fig. 2a, we observe that the dynamics does not change gradually with the change of γ , but a sudden change. In contrast, the two Markovian dynamics in (a) and (c) are similar to each other, which indicates that the QRT can be valid with the Born-Markov approximation. Moreover, the absolute magnitude of �σ x (τ + 10)σ z (10)� in (c) is much smaller, comparing to the value in (a). It is because the dissipative system loses energy into the environment gradually and the correlation for the long time difference vanishes quickly. Therefore, the short-time approximation is a reasonable way to simplify the problem. However, for some other complicated models, such as the spin-boson model, the TTCF will show significant differences even in the long-time limit.

Conclusion
In this paper, we have derived the exact evolution equation of the non-Markovian two-time correlation function for general open quantum system models, using the stochastic Schödinger equation approach. This method allows us to turn the hard open quantum system problems into a numerical ensemble average problem over all possible trajectories. Starting from a fully microscopic Hamiltonian, we can obtain a time-local SSE for the trajectory. In addition, our method employs one Gaussian noise to simulate two different trajectories and reduces ∂ t η 0 = −iH s η 0 − L †  www.nature.com/scientificreports/ the computational complexity significantly, which makes it possible to solve complicated and large-size models. Moreover, our method can be extended to study the multiple time correlation functions. Our numerical results disclose a fact that the non-Markovian TTCF is far different from the Markov limit, especially when the time is short, a strong non-Markovian regime. Our work can help better understand the dynamics of system in the presence of a non-equilibrium environment in a full quantum framework, and extend the research of non-Markovian quantum effects to more complicated chemical and biological systems.