Non-Markovian decoherence of a two-level system in a Lorentzian bosonic reservoir and a stochastic environment with finite correlation time

In this paper we investigate non-Markovian evolution of a two-level system (qubit) in a bosonic bath under influence of an external classical fluctuating environment. The interaction with the bath has the Lorentzian spectral density, and the fluctuating environment (stochastic field) is represented by a set of Ornstein-Uhlenbeck processes. Each of the subenvironments of the composite environment is able to induce non-Markovian dynamics of the two-level system. By means of the numerically exact method of hierarchical equations of motion, we study dependence of the steady states of the two-level system, the reduced density matrix evolution and the equilibrium emission spectrums on frequency cutoffs and coupling strengths of the subenvironments. Additionally we investigate the impact of the rotation-wave approximation (RWA) used for the interaction with the bath on accuracy of the results.


I. INTRODUCTION
Almost every quantum system interacts with its surroundings in a way that makes the system nearly impossible to be fully isolated. The interaction gives rise to processes of decoherence and appears to be harmful in some circumstances, e.g. in quantum information processing [1][2][3][4], and to be a valuable resource in others [5][6][7]. In some cases it is possible to capture evolution of the quantum system by means of a master equation of Lindblad form [8][9][10]. Systems of that type are called Markovian, and their evolution has the form of a quantum Markov process. Markovian systems are common in quantum optics, where a quantum system is often weakly coupled to an environment characterized by negligibly small correlation time.
Markovian systems assume their environments memoryless, however, an environment always changes under influence of a quantum system, which causes the memory effects to appear. Systems characterized by memoryful environments belong to the much wider category of non-Markovian systems [9,11]. Due to increased experimental and computational capabilities, non-Markovian systems are of great interest today. Among them there are such well-known systems as quantum dots [12][13][14][15], micromechanical resonators [16,17], superconducting qubits [18][19][20], and many others. Non-Markovian effects are ubiquitous in physics, chemistry, and biology and for systems interacting with either bosonic or fermionic reservoirs, like photosynthetic systems [21][22][23][24][25], molecular aggregates [26], molecular magnets [27], and solar cells [28]. Recently non-Markovian environments started to gain attention in quantum information processing [29][30][31], where attempts are made to utilize the backflow of information from the environment.
Accurate description of non-Markovian systems is a more complicated problem in comparison with description of Markovian systems. There are a plenty of methods known, but none of them are generally applicable, i.e. each one has its strong and weak sides, and the num-ber of systems it can describe efficiently is often limited. Analytical methods are represented mostly by perturbative expansions for some special parameter regimes, e.g. effective weak coupling theories or the projection operator techniques [32,33]. Numerically exact methods include the ones based on enlarging of the system state space, e.g. by extending the system space by the most relevant environmental modes [34][35][36], the ones utilizing tensor network approximations in propagation of influence functionals [37][38][39] and in mappings on effective 1D fermionic and bosonic chains [40][41][42], and etc. One of the most well-established numerically exact methods is the method of hierarchical equations of motion (HEOM) [43][44][45]. HEOMs utilize infinite systems of recurrent differential equations to encode the memory kernel of systemenvironment interaction and are able to handle a great variety of environmental spectral densities.
Switching from the Markovian approximation to a full non-Markovian description reveals many interesting phenomena. The most common one is the emergence of oscillations [46][47][48]. Non-Markovianity is known to be able to affect steady states (equilibrium states) of a system, e.g. it causes non-canonical steady states to appear [49], and also it affects the correlation functions. Actually, if non-Markovianity of a quantum process is sufficiently high, the quantum regression theorem (QRT) stops giving reliable correlation functions [50][51][52], which leads, for example, to significant differences between the predicted and the actual emission spectrums. It is a common practice to use the rotating wave approximation (RWA) for interaction with Markovian environments. Often the RWA is still used when the evolution become non-Markovian and can cause problems if non-Markovianity of the evolution is sufficiently high [11,53,54]. For example, wrongly used RWA may lead to incorrect shifts of the system frequencies [55] or cutoff all non-Markovianity [54].
When an environment is composite, i.e. consists of several subenvironments, it is possible to utilize one of the subenvironments to control decoherence of a quantum system. The case of a stochastic subenvironment as a control tool is rather popular in literature [56][57][58][59][60][61] and has similarities with dynamical decoupling schemes that alter the environment spectral density via filtering functions realized in sequences of laser impulses [62].
In the paper we investigate non-Markovian evolution of a two-level system (TSL) in a bosonic bath under influence of an additional external fluctuating environment. The bath spectral density function is chosen to be Lorentzian [8,63,64]. The Lorentzian spectral density is suitable, for example, for interaction between a Jaynes-Cummings cell and a zero-temperature bosonic bath [65]. The stochastic environment is chosen to be represented by a set of Ornstein-Uhlenbeck random processes. Following [66], we derive a HEOM capable of handling both RWA and non-RWA couplings with the bath equally accurate. Assuming the bath and the fluctuating environment to be independent, we analyze steady states of the TLS, evolution of the reduced density matrix, and equilibrium emission spectrums. We investigate the dependence on frequency cutoffs and coupling strengths of the subenvironments spectral densities in both RWA and non-RWA cases and provide a comparison with the ones obtained in the Markovian approximation [67].
The paper is organized as follows. In Sec. II we introduce the model, in Sec. III we present the hierarchical equations of motion. Next, we study the TLS evoution numerically. In Sec. IV we study steady states of the TLS, in Sec. V we investigate evolution of the reduced density matrix, and in Sec. VI we investigate the emission spectrums. Finally, we draw conclusions in Sec. VII.

II. MODEL
The full Hamiltonian for the system can be written aŝ whereĤ A = ω 0σ+σ− is the Hamiltonian for the free TLS, ω 0 is the TLS frequency,σ + andσ − are the rising and the lowering operators of the TLS, respectively;b + k andb k form a set of creation and annihilation operators describing the modes of the bosonic bath;Ĥ IF is the Hamiltonian for the interaction of the TLS with the fluctuating environment (stochastic field) andĤ IB describes the interaction between the TLS and the bath. The Hamiltonian for interaction with the stochastic field is defined by the next expression where Ω(t), ξ(t) are random functions,ξ(t) is the complex conjugation of ξ(t). The interaction gives rise to two decoherence channels, a dephasing channel and a relaxation channel, originating from the first and the second terms in (2), respectively. The random function Ω(t) is a real random process and ξ(t) is a complex random process. We assume that Ω(t) and ξ(t) are Markov processes of Ornstein-Uhlenbeck (OU) type [68], and consider the real and imaginary parts of ξ(t) as two independent real OU processes ξ 1 (t) and ξ 2 (t), respectively. Correlation functions of the random processes have the same form where ν ∈ {Ω, ξ 1 , ξ 2 }, ∆ ν is the standard deviation of the OU process and defines the coupling strength with the stochastic field, 1/γ ν is the correlation time of the OU process and γ ν has the physical meaning of the cutoff frequency of the environment. The HamiltonianĤ IB is used in two forms, the form corresponding to the full electric-dipole interaction (non-RWA) and the form representing the interaction in the rotating wave approximation (RWA). Introducing a new auxiliary TLS operatorâ, we can combine both forms of H IB in one expression where g k are the TLS-bath coupling constants,â + is the adjoint ofâ. Ifâ =σ + +σ − , (4) describes the full interaction, forâ =σ − it corresponds to the RWA interaction. For path integral methods it is more naturally to define an interaction with environment in the continuous form via the spectral density function, instead of utilizing the coupling constants g k directly. In the paper we consider the Lorentzian spectral density [8,[63][64][65] where ∆ B defines the coupling strength, γ B is the bath spectral width, also called the environment cutoff frequency. For the noise induced by the bath, 1/γ B represents the correlation time. The parameters have close relation to the corespondent parameters of the stochastic field ∆ ν and γ ν (3) and generally have the same physical meaning in terms of impact on the TLS dynamics.

III. HIERARCHICAL EQUATIONS OF MOTION
Let us suppose that before the initial moment of time the TLS does not interact with the bath and the stochastic field and both of the subenvironments are at equilibrium. Then the total density matrix at t = t 0 has the factorized form where P eq (Ω, ξ 1 , ξ 2 ) is the factorizable Gaussian equilibrium distribution function of the stochastic field,ρ (A) (t 0 ) denotes the initial density matrix of the TLS, and ρ (B) eq (t 0 ) is the equilibrium bath density matrix at zero temperature.
For the factorized initial conditions (6) and the Lorentzian spectral density (5), we can obtain the HEOM by the steps presented in [66], where we have to replace the high-temperature Drude spectral density with the Lorentzian spectral density and to take the limit β → ∞ for the bath subenvironment. The HEOM expression has the same form and can be written as where m combines the field and the bath indexes into one composite index. We assume that the first three components of m index the recursion relations for the stochastic field in the order {Ω, ξ 1 , ξ 2 }. The last two components index the bath recursion relations. The special notation m| k+1 is used for the index m with the k-th component increased by 1  (1) F,k belong to the stochastic field part of the recursion relations and can be expressed via parameters of the stochastic field and the field coupling operators where ν k denotes the k-th element in {Ω, ξ 1 , ξ 2 }, e.g.
The remaining coefficients α B,k define the bath part of the recurrence relations and depend on the bath spectral density. If the spectral density is Lorentzian (5), we have the next expressions for the coefficients whereĉ 1 =â andĉ 2 =â + , the superscript "x" denotes the commutator superoperator defined earlier, and the superscript "R" denotes the superoperator for action from the right, i.e.ĉ R 2ρ =ρĉ 2 . From (11) it follows that the bath part of the HEOM cannot be transformed into the one-indexed form, in contrast with the stochastic field part [44] and the case of non-RWA interaction with the high-temperature Drude bath [66], because there are two unequal constants α k2 . While the infinite-temperature Drude bath appears to be quite similar to the stochastic field as an environment for the TLS [43], we expect that the zerotemperature Lorentzian bath and the stochastic field act on the TLS in qualitatively different ways.

IV. STEADY STATES
The steady states of the TLS reachable from the selected initial state (6) can be obtained by propagating the state forward in time until the reduced density matrix stops changing. Typical response of the steady states on changes of the environment parameters is presented in Fig. 1.
For a wide range of frequency cutoffs and coupling strengths tested, the stochastic field acting alone brings the TLS to the steady state where both the excited and the ground states are equally possible. Qualitatively different picture is observed for the TLS in the Lorentzian bath. Here, if one of the approximations is used, either the RWA or the Markovian, the TLS equilibrates in its ground state. Otherwise, in case of the non-RWA interaction with the bath, the probability of the excited state is above zero and gradually rises with the frequency cutoff, starting near zero and then approaching 1/2 from below ( Fig. 1(a)). Also in the non-RWA case the stationary state exhibits a weak dependence on the coupling strength, the excited state probability rises from zero to 1/2 from below, but much slower, then it is in the frequency cut-off case. Thus, the RWA, same as the Markovian approximation that intrinsically utilizes the RWA, alters the stationary states behavior drastically.
When the TLS interacts with both subenvironments ( Fig. 1(b)-1(e)), the difference between the stationary states for each of the subenvironments becomes visible. Because the difference between the steady states may be big, e.g. for low bath frequency cutoffs in the non-RWA case or for any parameters in the RWA case, the impact of the stochastic field can be significant.
If we fix parameters of the stochastic field and start increasing the bath frequency cutoff γ B (Fig. 1(b)) or the coupling strength with the bath ∆ B (Fig. 1(c)) starting from zero, the bath contribution in a steady state grows and the excited state population in the steady state decreases, because the steady states for decoherence in the bath always lie lower. Because in the RWA case the steady states are all completely unexcited ( Fig. 1(a)), the RWA curve is always below the non-RWA one. The distance between them constantly increases as the non-RWA steady state for decoherence in the bath shifts up with both the bath frequency cutoff γ B and the bath coupling strength ∆ B . The RWA curves exhibit saturation in both figures, but the non-RWA curves reach minimums, go up and approach the Markovian curves from below (more pronounced in Fig. 1(b)). In the Markovian approximation the steady states seem insensitive to any changes of the bath spectral density parameters. Now let us fix the bath parameters and change the stochastic field instead. If we begin to gradually increase the stochastic field frequency cutoff from zero changing all the random processes simultaneously γ ν = γ F (Fig. 1(e)), the field contribution in the resulting steady state starts growing and the steady excited state population of the TLS starts growing either, because the steady states for decoherence in the stochastic field always lie higher. At some value of γ F we reach the maximum, and after it the excited state population starts decreasing. The RWA applied to the interaction with the bath shifts the steady state down with respect to the one of the non-RWA curve, while the Markovian approximation gives much higher steady excited states populations for small frequency cutoffs and significantly overestimates the rate at which they decrease. The situation is similar if we change the field coupling strength in the same way ∆ ν = ∆ F (Fig. 1(d)), but there is no maximum, and the curves continues to rise, approaching the Markovian curve from below.
The observed behavior can be explained via the magni-tude of the environment spectral density (for the stochastic field it is the spectrum of corresponding random processes) in the vicinity of the TLS resonant frequency. The OU random processes and the Lorentzian bath have spectral densities with one peak. The OU process peak is located at ω = 0 and becomes wider and lower when the correspondent frequency cutoff increases. In Fig. 1(d) we see how the stochastic field contribution in the steady states grows at first, causing the rise of the steady states curves, because the spectral density peak widens and the magnitude of the spectral density near the resonant frequency increases, then the contribution falls, because the process of the spectral density peak declining becomes dominating, and the curves go down too. The switch from the growth to the decline of the field spectral density near the resonance frequency determines the maximum of the RWA and the non-RWA curves.
In contrast, the peak of the bath spectral density is always located at resonance and widens when the bath frequency cutoff increases. In Fig. 1(b) we see how the widening increases the magnitude of the bath spectral density in the vicinity of ω 0 and the contribution of the bath in the steady states start growing, causing the curves to go down, because the TLS steady states for interaction only with the bath are located lower. Then the magnitude reaches its maximum and the contribution saturates. The subsequent rise of the non-RWA curves can be explained by the rise of the non-RWA curve for decoherence in the bath only ( Fig. 1(a)).
Coupling strengths of both subenvironments affect only heights of the corresponding spectral density peaks, when a coupling strength grows, the peak grows either. It results in gradual increase of the magnitude of a spectral density near the TLS resonant frequency and can be seen in Figs. 1(c) and 1(e). In Fig. 1(c) the bath contribution increases, shifting the curves down, while in Fig. 1(e) the field contribution increases, shifting the curves up.
The situation becomes more complex if we stop using the restrictions γ ν = γ F with ∆ ν = ∆ F and allow arbitrary changes for each of the underlying random processes of the stochastic field. The detailed study of impact of each of the processes on the steady states will be presented elsewhere.

V. DENSITY MATRIX EVOLUTION
An initially excited TLS placed in an equilibrium non-Markovian environment loses coherence due to the interaction with the environment, but the process is not monotone. At some point during the evolution the information backflow from the environment begins restoring the coherence, then the backflow weakens and the TLS starts losing coherence again. As a result, the oscillatory behavior emerges.
The evolution of the reduced density matrix from the factorized initial state (6) for different parameter regimes is presented in Figs. 2 and 3. The curves corresponding to the non-Markovian evolution, both for the RWA and the non-RWA types of coupling with the bath, exhibit rapidly vanishing oscillations, which are more evident for the RWA curves, including all the non-Markovian curves for decoherence in the stochastic field (Fig. 2). The amplitude of the oscillations has a clear relation to shapes of the subenvironments spectral densities in the vicinity of the TLS resonance frequency ω 0 . If the environment spectral density is flat in the vicinity of ω 0 , which is the case of large frequency cutoffs, there are no oscillations. For example, the Markovian approximation curves in Figs. 2(a) and 3(a) exhibit no oscillations, because the Markovian approximation assumes that environment correlation times are small, which corresponds to large frequency cutoffs. If the environment spectral density is not flat in the vicinity of the TLS resonance, the oscillations appear. The dependence of the oscillations amplitude on the frequency cutoff value in Figs. 2(a) and 3(a) is more evident in case of decoherence in the bath, because the peak of its spectral density is located at the TLS resonance frequency. The coupling strength of the environment impacts the amplitude of the oscillations in the opposite way (Figs. 2(b) and 3(b)).
The evolution becomes faster if the coupling strength increases, i.e. the minimums are located closer and the steady states are reached earlier (Fig. 2(b) and 3). The frequency cutoff impacts the speed of the evolution in a more complex way, there is a cutoff frequency for which the evolution speed is maximum (not shown).
The main difference between the RWA and the non-RWA curves in Fig. 3 resides in values of the minimums. The RWA curves have its minimums placed at the horizontal axis where also the stationary values are located, while the minimums of all the non-RWA curves are placed strictly above the corresponding stationary values. The stochastic field curves for sufficiently large couplings have their first minimum located below the stationary value (Fig. 2), which resembles the behavior of the non-RWA curves for decoherence in the bath. In overall, the non-RWA evolution of the reduced density matrix reminds the smoothed version of the RWA evolution.
When the TLS interacts with both subenvironments simultaneously (Fig. 4), the evolution becomes faster due to the presence of an additional decoherence channel, the first minimum is moved to the left and the stationary value is reached earlier. Because the steady state is shifted up by the stochastic field, the curves lie above the corresponding curves for decoherence in the bath only. Also the stochastic field significantly damps the oscillations, which is evident even for rather weak coupling strengths with the field in comparison with the coupling strength with the bath. The Markovian approximation shows the fastest decoherence among all the curves.

VI. EMISSION SPECTRUM
We obtain the equilibrium emission spectrums of the TLS by applying the Fourier transform to the two-time correlation function σ + (t 2 )σ − (t 1 ) , where t 2 > t 1 . The time t 1 is selected sufficiently big for the reduced density matrix evolution to reach its stationary phase, thereby the correlation function may be considered stationary. We calculate the stationary correlation function in the following way. First we propagate the initial state to the steady state, then apply the operatorσ − to all density matricesρ (A) m (t 1 ), lying in the TLS subspace, next the result is propagated to t 2 , whereσ + is applied.
The equilibrium emission spectrums for interaction with the Lorentzian bath can be obtained only in the case of the non-RWA coupling, because the steady states in the Markovian and the RWA approximations are com- pletely unexcited and do not emit (Fig. 5). If the frequency cutoff is large, the non-RWA emission spectrum has one peak, which shifts to the right when the cutoff becomes smaller (not shown). At some cutoff value the top of the peak becomes a plateau and splits in two practically non-distinguishable peaks of non-equal intensity which are placed symmetrically with respect to ω = ω 0 (not shown). Then the peaks move in the opposite directions, slightly declining, but stop at some value of the cutoff and start moving backwards while becoming more and more distinct (Fig. 5(a)). At the same time the overall intensity of the spectrum gradually decreases. As a result the two-peaked spectrums are fairly weak in comparison with the one-peaked spectrums. In Fig. 5 we use the normalization by the maximum value, so the actual spectrum intensity is not shown. Potentially the appearance of the two peaks is explained by the presence of the two complex-conjugated coefficients α The dependence on the coupling strength with the bath is slightly more simple. If we increase it, the main peak widens and moves to the right, then splits in two peaks of unequal intensity (Fig. 5(b)). If we increase the coupling strength further, the right peak, which is the main peak, decreases, and the left (the side peak) grows, while moving in the opposite directions, the peaks resolution becomes better. For large coupling strengths the side peak becomes the dominant and approaches ω = 0.
For comparison, we show the equilibrium emission spectrums for decoherence in the stochastic eld in Fig. 6. The spectrum differs qualitatively from the case of decoherence in the Lorentzian bath. The stochastic environment spectral density has a peak located at ω = 0 for any value of the frequency cutoff and the coupling strength.
The peak becomes more distinct if the frequency cutoff lowers or the coupling strength rises. The resonance is clearly visible in Fig. 6, where a side peak located at ω = 0 appears if the cutoff frequency is sufficiently small or the coupling strength is sufficiently large. At the same time, the main peak shifts to the right and widens. The spectrum energy redistributes from the main peak to the side peak, and the side peak grows while the main peak decreases.
The Markovian approximation effectively considers the environmental spectral density flat (large frequency cutoffs), or, equivalently, it considers only a small region of the spectral density in the vicinity of the TLS resonance. If the spectral density resonance is located sufficiently far from the TLS resonance, the Markovian approximation loses the essential information about the peak existence. In Fig. 6 it results in wide one-peaked spectrums, which peaks are located in the vicinity of the TLS resonance frequency.
In Fig. 7 we show the impact of the stochastic field on the equilibrium emission spectrums for decoherence in the bath. The stochastic field causes the emergence of the zero-frequency peak, like it does in Fig. 6, so the spectrum obtains the three-peaked form. Also it lowers the left peak of the doublet (the side peak in Fig. 5), widens it and shifts the main peak (the rightmost peak) to the right. The field smoothes the negative frequencies spectrum, making the zero point disappear, and increases intensity of the RWA curve so that it can be observed. The Markovian approximation is clearly inaccurate in the parameter regions selected.

VII. CONCLUSION
We have studied non-Markovian evolution of a TLS in a composite environment consisting of two subenvironments, a zero-temperature bosonic bath characterized by the Lorentzian spectral density and a stochastic field of the Ornstein-Uhlenbeck type, and analyzed the impact of the rotating-wave approximation used for the interaction with the bath.
It was shown that the steady states for decoherence in the bath depend on the coupling type used and the full interaction leads to different steady states in comparison with the cases when either the RWA or the Markovian approximation is used. We investigated the joint influence of the subenvironments on the steady states and found connections with the shape of the environment spectral density in the vicinity of the TLS resonant frequency.
We demonstrated the dependence of the reduced density matrix evolution on the frequency cutoff and the coupling strength of the environment. In all cases, except the cases involving the Markovian approximation, the reduced density matrix exhibits the oscillatory behavior, the amplitude of which can be explained via shapes of the subenvironments spectral densities in the vicinity of the TLS resonance frequency. We showed that increasing the frequency cutoff smooths the oscillations and shifts the first minimum location to the right while increasing the coupling strength acts in the opposite way. Comparing the cases of the full and the RWA couplings to the bath, we found that the minimums of the oscillation in the RWA can be located below the stationary value in contrast to the case of the full interaction where this is not observed.
We investigated the dependece of the TLS equilibrium emission spectrums on the frequency cutoffs and the coupling strengths of the subenvironments and found that, if the TLS interacts only with the bath, the spectrums can have the doublet form in the positive frequencies domain and the doublet form for the negative frequencies domain at the same time, but the intensity of the spectrum is low. For interaction with both subenvironments, the spectrum can have three distinct peaks, one of which is located at the zero frequency, and the other two are located at the opposite sides of the TLS resonance.