{'fraction_non_alphanumeric': 0.04067568703252374, 'fraction_numerical': 0.017144297840154636, 'mean_word_length': 4.0919982884039365, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Non-Markovian Quantum Jump with Generalized Lindblad Master Equation

The Monte Carlo wave function method or the quantum trajectory/jump approach is a powerful tool to study dissipative dynamics governed by the Markovian master equation, in particular for high-dimensional systems and when it is difficult to simulate directly. In this paper, we extend this method to the non-Markovian case described by the generalized Lindblad master equation. Two examples to illustrate the method are presented and discussed. The results show that the method can correctly reproduce the dissipative dynamics for the system. The difference between this method and the traditional Markovian jump approach and the computational efficiency of this method are also discussed.


I. INTRODUCTION
Since the pioneering work of Albert Einstein, who explained the phenomenon of dissipation and Brownian motion in his annus mirabilis of 1905 by use of statistical methods, a rich variety of methods to tackle quantum fluctuations and quantum dissipation in open systems has been proposed[1, 2]. Among them the quantum master equation (QME) approach and the quantum Langevin description (QLE) [3] are two of powerful functional integral techniques for the study of time evolution of open quantum systems. The quantum master equation can be divided into two categories: Markovian and non-Markovian. The Markovian master equation [4] (especially in the Lindbald form) can be derived with the weak coupling limit(or the Born approximation) and the Markovian approximation. It can be solved analytically [5] for some special cases, but for most cases we have to solve and simulate it numerically by the Monte Carlo wave function method or quantum trajectory/jump approach [6,7,8,9,10,11]. This method is very effective for qubit systems even with large number of qubits, say n = 24 [9].
However, the dynamics of an open system is not always Markovian. Strong systemenvironment couplings, correlation and entanglement in the initial state and structured reservoirs may lead the dynamics far from Markovian. Many methods have been proposed to describe the non-Markovian process, including the Lindblad equation with time dependent decay rates [12], generalized Lindblad equation [13] obtained from the correlated projection superoperator techniques [14,15], phenomenological memory kernel master equation [16,17] and the post-Markovian master equation [18,19,20]. The first two methods are local in time while the last two involve an integral of time. For the first method, the only difference from the Markovian master equation is that the decay rates in the equation are time-dependent.
These decay rates may take not only positive values but also negative ones. When decay rates are positive, the Markovian Monte Carlo wave function method can directly be used.
However, the mehtod is not available when the decay rates are negative. This problem was solved in Ref. [21] by introducing reversed jumps.
The generalized Lindblad master equation can well describe the dynamics of an open system beyond the Markovian limit, especially it is very effective for an environment composed of spins [22,23,24] and structured reservoirs [25]. However the extension of the Monte Carlo simulation to this equation remains untouched. In this paper, we will explore the unravel-ing and quantum trajectory approach for the generalized Lindblad equation. The structure of this paper is organized as follows. In Sec.II we briefly review the generalized Lindblad

II. GENERALIZED LINDBLAD MASTER EQUATION
The equation that governs the dynamics of an open quantum system can be derived by means of the projection superoperator technique [12,14]. The form (Markovian or non-Markovian) of the master equation crucially depends on the approximation used in the derivation, reflecting in the projection superoperator chosen. When we project the total system state into a tensor product, we can obtain the Markovian master equation, whereas a non-Markovian master equation can be obtained when we use a correlated projection.
The following is the master equation derived by this method and it is called the generalized Lindblad master equation [13], where H m are Hermitian operators and R λ mn are arbitrary system operators depending on the form of system-environment interactions. If we have only a single component ρ S = ρ 1 , this equation obviously reduces to the ordinary Markovian master equation. In this paper we will focus on the case where we have at least two components. The state of the reduced system in this case is ρ S = m ρ m , we remind that Trρ m < 1.

III. QUANTUM JUMP
For clarity, we define the jump operators W λ mn = R λ mn and non-jump operators W 0 mm = I − iH m dt, where the non-Hermitian effective Hamiltonian is given by There are two subscripts and one superscript for the operator W λ mn . The first subscript m denotes the index of component where the system in, while the second subscript n denotes the index of component for the operation acting on; the superscript λ represents the jump mode. Initially we assume that each operator ρ m (t 0 ) can be written is a non-normalized wave function. After an infinitesimal time dt, it evolves into the following state where the new states are defined by and with probabilities respectively. In Eqs. (3) and (4), is the weight for the component ρ m that satisfies Note that the jumps for ρ m depend on the other components ρ n , (n = m) of the reduced density matrix ρ. This makes our method different from the traditional quantum jump method.
(1) can be simulated numerically by the so-called Monte Carlo wave function approach according to the unraveling given above. We start the time evolution from the state where ρ m (m = 1, 2, 3, ...) are the components for ρ. At time t 0 + dt, where dt is much smaller than the timescale relevant for the evolution of the density matrix, a random number ǫ which is randomly distributed in the unit interval [0, 1] is used to determine the jump. Note that all the components are controlled by this random number. For each component |ψ m , if 0 ≤ ǫ ≤ dp 1 m1 , it jumps to |ψ 1 m1 , if dp 1 m1 < ǫ ≤ dp 1 m1 + dp 2 m1 , it jumps to |ψ 2 m1 , and so on. These jumps are all operated on the component ρ 1 ; if λ dp λ m1 < ǫ ≤ λ dp λ m1 + dp 1 m2 , it jumps to the component 2, namely |ψ 1 m2 . Jumps to the other components can be established in a similar way. If ǫ > nλ dp λ mn , a non-jump takes place and the state ends up in |ψ 0 mm . This operation is acted on the component ρ m itself. We define a generalized jump superoperator W i , which denotes all jumps for all the components controlled by this random number. We repeat this process as many times as n = ∆t/dt for all the components, where ∆t is the total evolution time. We call this single evolution a generalized quantum trajectory. This trajectory contains all the components of the density matrix. Given an operator A, we can write its mean value A (t) = Tr(Aρ(t)) as an average over N trajectories as

IV. APPLICATION
In this section, we use the model and the generalized master equation given in Refs. [25] and [23] as two examples to illustrate our method. First consider a two-state system coupled to an environment. The environment consists of a large number of energy levels which are arranged into two energy bands with the same energy spacing(see Fig.1). The lower energy band contains N 1 levels while the upper one N 2 levels. This model can be understood as a "many level" environment or "container", of which only the relevant parts of the spectrum enter the model. For details of this model, we refer the reader to [26,27]. The total Hamiltonian for a qubit coupled to such an environment in Schrödinger picture is with(we seth = 1) δǫ N 1 n 1 |n 1 n 1 | + n 2 (ω + δǫ N 2 n 2 )|n 2 n 2 |, V = λ n 1 n 2 c(n 1 , n 2 )σ + |n 1 n 2 | + H.c., where the index n 1 denotes the levels of lower energy band and n 2 denotes the levels of upper band, σ z and σ ± are Pauli operators. λ is the overall strength of the interaction, c(n 1 , n 2 ) are coupling constants, they are independent of each other and are identically distributed, satisfying c(n 1 , n 2 ) = 0, c(n 1 , n 2 )c(n ′ 1 , n ′ 2 ) = 0, c(n 1 , n 2 )c * (n ′ 1 , n ′ 2 ) = δ n 1 ,n ′ 1 δ n 2 ,n ′ 2 .
According to H 0 , one can transform the problem into the interaction picture and, with the help of projection superoperator technique, obtain the non-Markovian evolution equation as d dt ρ (1) where γ i = 2πλ 2 N i δǫ (i = 1, 2).
With definitions of Π 1 = n 1 |n 1 n 1 | and Π 2 = n 2 |n 2 n 2 |, Π 1 + Π 2 = I E , the two non-normalized density matrixes can be obtained by ρ (i) S = Tr E (Π i ρ T ), i = 1, 2, where ρ T is the total density matrix for the system and environment. The reduced density matrix for the system is then given by ρ = ρ (1) S . We note that in Eq.(10), there are no environment operators other than the two (c-number) parameters γ 1 , γ 2 . The initial state of the environment is taken into account by means of the distribution of initial ρ In this model, there is only one jump operator for each component, i.e. W 1 12 = √ γ 1 σ + and W 1 21 = √ γ 2 σ − , and non-jump operators W 0 mm = I − iH m dt with H 1 = − 1 2 γ 2 σ + σ − and H 2 = − 1 2 γ 1 σ − σ + . We consider two types of initial condition in the following simulation. First, only the lower band of the environment is populated, i.e. ρ (2) S = 0. Under this condition, the reduced system can be solved analytically. Another case is, the two bands of the environment are all populated. With this initial condition, we solve the master equation numerically. In both cases, we choose an initial states |φ(0) = |e and |φ(0) = 1 2 (|e + |g ) for the system, where |e and |g denote the excited state and ground state, respectively. We compare the analytic solution and the numerical simulation(solve the equation by Runge-Kutta method) to the results obtained from the quantum jump/trajectory approach in Figs.2 and 3. The trajectory number in this quantum jump approach is N = 400. We can see from the figures that the quantum trajectory approach correctly reproduces the system evolution. The errors are sufficiently small, although we choose a small number of trajectories, showing that this method is efficient.
Another example is a qubit coupled to a spin bath [23]. The full system consists of a central spin interacting with a bath of N spins. Such a system can be described by where σ denotes the Pauli matrix for the central spin that is the system we are interested in, and σ k stands for the k-th spin in the bath. After defining an unperturbed part H 0 = ω 2 σ z + 2σ z K z , where K z = 1 2 N k=1 α k σ k z , the Hamiltonian can be transformed into the interaction picture. Assuming the parameters are real and time independent, the master equation reads where ρ m = Tr B (ρ T Π m ), ρ T is the density matrix for the total system(the central spin plus the bath), Π m is a projection superoperator that projects the z-component of the bath angular momentum into an eigenvector with eigenvalue m. We take N = 2 as an example, then the density matrix of the central spin has three components, denoted by ρ 1 , ρ 0 , ρ −1 , respectively. Each component has two jump operators which act on the other two components, and a non-jump operator, which acts on itself. The comparison between directly numerical simulations (by Runge-Kutta) and quantum trajectory method is shown in Fig.4. Here the trajectory number is chosen to be N = 4000. We can find that as the number of jump operators and components increases, the number of quantum trajectory, with that we can obtain a correct result, increases.

V. CONCLUSION AND DISCUSSION
In this paper, we have developed an efficient unraveling for the generalized Lindblad Here J is the number of jump operators for each component(including the non-jump operator