Non-Markovian dynamics of the driven spin-boson model

The non-Markovian dynamics of the driven spin-boson model at zero and finite temperature is investigated theoretically. Based on a unitary transformation and convolution theorem, the analytical expression of population difference in Laplace space is obtained. This method provides a tool to help us to understand the interference between driving and dissipation within the non-Markovian frame. We study the transient and the steady-state non-Markovian dynamics of population difference P(t) in the case of different temperature, driving and spin–bath coupling (SBC) strength. It has been shown that it is necessary to use non-Markovian methods when SBC is strong. At zero temperature, the curve of steady-state amplitude P lt versus Rabi driving amplitude Ω is N-like shaped. While at finite temperature, the curve becomes Λ-like shaped. These curves come from the interference between driving and dissipation. Also, from the interference pattern above, we know that large coherent oscillations can be obtained by optimizing the parameters: Rabi driving amplitude, temperature and SBC strength.


Introduction
The dissipation-induced decoherence in a quantum microscopic system can be effectively modeled by a spin-boson model (SBM) [1,2], which describes the interaction between a quantum two-level system and a bosonic bath. The controllable dynamics of a driven SBM (DSBM) is at the core of many vastly different state-of-the-art technologies, especially solid-state implementations of individually addressable two-level-system [3]. The DSBM has been widely studied in experiment and theory, which is related to various physical and chemical processes [1,2,4,5]. Recently, great progress has been made in the experiment of controlled coherent dynamics of qubits in various devices, such as superconducting devices based on Josephson tunnel junctions [6], optically and electrically controlled single spins in quantum dots [7], individual charge in quantum dots [8], and nitrogen vacancy center in diamond [9]. The Hamiltonian of the DSBM reads (set = 1): with H s (t) = − 1 2 Δσ x + Ω cos(ω L t)σ z , and H s (t) is the Hamiltonian of the driven spin and is widely used to describe the controlled qubit with the harmonic driving Ω cos(ω L t), where Ω is the driving amplitude, i.e. Rabi frequency, and ω L is the driving frequency. H b is Hamiltonian of the bath with finite temperature. And H I is the interaction between the spin and bath without rotating-wave approximation (RWA). Here Δ is the coherence tunneling of the spin, σ x and σ z denote the usual Pauli spin matrices. A spin can experience different types of baths. One of the commonly-used the photon or phonon bath can be described by an Ohmic spectrum with the spectral density: J(ω) = k g 2 k δ(ω − ω k ) = 2αωθ(ω c − ω). Understanding and manipulating thermal transport in low-dimensional materials is of significant importance both for scientific development and practical application [10]. The quantumness and limitation of thermal conductance in small quantum systems has recently attracted considerable attention via the quantum dissipation of systems by thermal fluctuations, which results in the quantized thermal conductance of phonons [11,12] and photons [13]. Consequently, the development of low-dimensional quantum thermal devices has stimulated the study of finite temperature small quantum systems dynamics.
Furthermore, it has been found that the effects of counter-rotating (CR) terms are significant in many different interesting topics, such as the quantum Zeno effect [14], interference between driving and dissipation [15,16], coherent destruction of tunneling and Bloch-Siegert shift [17,18], the Rabi model with frequency modulation [19,20], entanglement evolution [21], and so on. The main purpose of this paper is to show the important role of the CR coupling and temperature on the non-Markovian dynamics of DSBM.
The reduced system dynamics of the DSBM has been studied by various analytical and numerical methods, for example, the polaron or generalized Silbey-Harris transformation approach [22][23][24][25], the time-dependent numerical renormalization-group method [26], the quasi-adiabatic propagator path integral [27], stochastic Schrödinger equation (SSE) approach [28] and hierarchical equations of motion approach (HEOM) [29,30], together with the means of a generalized master equation (GME) and the noninteracting-blip approximation (NIBA). Each method has its own regimes of validity depending on the spin-bath coupling (SBC) strength, the bath temperature, and the form of the bath spectral density. In this paper, using the analytically perturbation method based on the unitary transformation [15,16], we study the non-Markovian dynamics of DSBM at finite temperature. Our method is suitable for any bath spectral density function and is in good agreement with the numerical results in the case of (α < 1) and the temperature (T < Δ). This analytical method provides a tool to help us to investigate the interference between driving and dissipation within the non-Markovian frame.
In this paper, we calculate the transient and steady state non-Markovian dynamics of P(t) and analyze the effects of temperature, SBC strength and driving amplitude on P(t). The results show that beyond the weak SBC (α > 0.01), there are obvious differences between the dynamics curves obtained by Markovian approximation and non-Markovian approach. At zero temperature, the curve of steady-state amplitude P lt with Rabi driving amplitude Ω is N-like shaped. While at finite temperature, the curve becomes Λ-like shaped. These nonlinear behaviors are caused by the interference between driving and dissipation. From this interference, it is possible to optimize the parameters: Rabi driving amplitude, temperature and SBC strength, to obtain large amplitude stable coherent oscillation.
For SBM with finite temperature, the renormalized spin tunneling Δ r depends on the temperature, which has been discussed in detail [31,32]. In short, when the bath temperature T is less than a critical temperature T c , a self-consistent solution of the renormalized spin tunneling Δ r can be obtained. However, if T is larger than T c , the renormalized spin tunneling vanishes, Δ r = 0. The presence of a critical temperature T c corresponds to a transition from Δ r = 0 to Δ r = 0. Our discussion in this paper is limited to T < T c .

The model and canonical transformation
Taking into account of the correlation between the spin and bath, we apply a canonical transformation, with the generator S [23,24], and a k-dependent parameter The transformation can be done to the end and the transformed Hamiltonian H can be written as where the unperturbed Hamiltonian reads the first-order coupling term is given by and the higher-order term is given by and Until here, there is no approximation, and the transformed Hamiltonian H = H 0 (t) + H 1 + H 2 is equivalent to the original H(t). In the following, we ignored the terms H 2 , which are related to the doubleand multi-boson off-diagonal transition (like a k a k and a † k a † k ), and their contribution to the physical quantities are O(g 4 k ). Up to the order of g k , the transformed Hamiltonian H can be approximated as whereg Here ηΔ is the renormalized tunneling Δ r , and we define the eigenstates of σ x as |s 1 and |s 2 with σ x |s 1(2) = − (+) |s 1(2) , and the corresponding creation and annihilation operators s + = |s 2 s 1 | and s − = |s 1 s 2 |. Using J(ω) = g 2 k δ(ω − ω k ), one can derive from equation (13) that η is determined self-consistently by the equation The dependence of the renormalization factor η on temperature has been discussed [31]. The interaction terms after the unitary transformation is H 1 = kg k a † k s − + a k s + and the interaction strengths are effectively weakened. Now we solve the equation of motion of the system in the rotating frame with the rotating operator R(t) = exp[i(− 1 2 ηΔσ x + ηΔa † k a k )t] and treat the resonant case ω L = ηΔ. Then the Hamiltonian (14) becomes which has the same form as the Hamiltonian in RWA, but its parameters are renormalized to include the effects of the CR terms of SBC. From the transformed Hamiltonian H , one can see that the ground state of the transformed Hamiltonian H is and the corresponding ground-state energy is −ηΔ/2. Therefore, the ground state of the original Hamiltonian H is given by which is the dressed state of the spin and the bath due to the CR terms.
Here we discuss the property of our unitary transformation. Compared with the well-known variational polaron transformation [23,24], the form of generator (6) in unitary transformation is the same. However, how to determine the parameters ξ k is different. The variational method [23,24] determines the transformation parameters by minimizing the Gibbs-Bogoliubov-Feynman upper bound on the free energy. Our transformation parameters are determined by decoupling the lowest energy state |g from other excited states after the transformation. That is to say, H 1 |g = 0. Thus, we get rid of the variational condition. Although the variational method can predict correctly the delocalized-localized transition point at the scaling limit, it is not suitable for calculating the dynamical properties of SBM. Our approach starts from the unitary transformation equation (5) and then make perturbative calculation based on the transformed Hamiltonian.
The most important physical meaning of our unitary transformation is it effectively weakens the strength of SBC. Specifically, the interaction term after the unitary transformation, H 1 is of the order of bath fluctuation and its thermal average is zero, hence H 1 is a reliable perturbative parameter. Based on this consideration, the Born approximation after our unitary transformation is applied to derive the master equation. If the transformation parameter ξ k = 1, our unitary transformation is equivalent to the full polaron transformation [25], which is accurate for the fast bath case.
To sum up, our transformation has the following advantages: (1) the transformations are not restricted by the spectral density of bath; (2) the interaction term has the simple form of energy conservation, which is important for the convolution operation in solving the equation of motion and leads to an analytic solution; (3) the transformations avoid the infrared divergence in the Ohmic bath [31]. However, the higher-order terms H 2 are ignored in this method, so the calculation results deviate from the numerical exact solutions in the case of high temperature (T > Δ) and strong coupling (α > 1). Moreover, the method in this paper is only applicable to the case of resonance (ω L = ηΔ). If consider the non-resonant case, the time-dependent unitary transformation will be introduced [15], which will be discussed in our future work. In total, as an analytical method concerning the non-Markovian dynamics, it is helpful for us to have a deep understanding for physical problems within the non-Markovian frame.

Equation of motion for the density matrix
The equation of motion of the density matrix ρ SB (t) for the whole system, i.e. driven spin (S) and bath (B) is given by dρ SB (t)/dt = −i[H, ρ SB (t)]. After the unitary transformation in equation (5), we have where ρ SB (t) = exp(S)ρ SB (t)exp(−S) is the density matrix of the total spin-bath system in the Schrödinger picture with the transformed Hamiltonian H (i.e. equation (14)). In the interaction picture with

the interaction Hamiltonian becomes,
Tracing over the degrees of freedom of the bath and returning back to the Schrödinger picture, the equation of motion of the density matrix is given by The derivation details of the master equation are listed in the appendix A.
Note that in the derivation of the master equation, we use Born approximation. Different from the traditional Born approximation, our Born approximation is done after the unitary transformation equation (5), and the unitary transformation effectively weakens the strength of the SBC. Specifically, H 1 is the order of bath fluctuation and its thermal average is zero [23]. Hence, small H 1 guarantees the reliability of the Born approximation treatment. In the case of finite temperature, there are eight summation terms for the bath indicators, and their coefficients depend on n k or n k + 1 respectively, with When the bath reduces to zero temperature, n k = 0, and the eight terms decrease to four ones. In equation (22), the time-dependent variables are all in the form of (t − t) and t . Then, the Laplace transformation and convolution theorem can be used. According to the Kronecker product property and the technique to Lyapunov matrix equation, we expand the density matrix into a vector Note |1 and |2 are the eigenstates of σ x . Based on Laplace transform and convolution theorem, the equation of motion can be obtained whereL(s) is the Liouvillion super-operator in the Laplace space. And considering the symmetry of the super-operatorL(s) and the linear property of Laplace transformation, we apply the following transformation for the equation of motion, Therefore the recombined density matrix becomes The equation of motion are rearranged as, and the regrouped super-operator turns into, The concrete expressions of the matrix elements in the super-operator are shown in appendix B. Finally, the non-Markovian result in the Laplace space is obtained from the equation (27), The initial density matrix |ρ S (0) = |1, 1, 1, 1 /2 is adopted in above derivation. In particular, if the summation of the bath's degrees of freedom is approximated as solution of DSBM naturally reduces to Markovian results [15][16][17][18]. Therefore, the Markovian approximation gives L 23 = L 24 = L 12 = 0, and the solution is simplified as i |ρ 21 From the analytical derivation, we know that the Markovian approximation not only ignores the process with the terms such as (ft 2+ − ft 2− ) and (ft 3+ − ft 3− ), but also, even in the process considered, the frequency shifts caused by the integral principal value are ignored. From the residue theorem, we know that the denominator of the solution determines the resonance frequency and line width of the spin.  non-Markovian dynamics is different from the Markovian ones. At last, the non-Markovian dynamics is obtained by inverse Laplace transformation, which are given in appendix C. To explore the dynamical evolution, the population difference is given Note this method is not limited by spectral function.

The result and discussion
In this paper, taking Δ as the energy unit, we discuss only the case that the driving frequency resonates with the renormalized tunneling, that is, ω L = ηΔ. Before addressing the effect of driving, we show that our results are indeed non-Markovian dynamics. In figure 1, we compare with the numerical exact results, which is the non-Markovian dynamics of different methods: SSE approach [28] and HEOM [29], together with the means of a GME and the well known NIBA. Our analytical results are reasonable agreement with the exact non-Markovian ones, even in not very weak dissipation, α = 0.1 and finite temperature case.
The Markovian dynamics of the DSBM in zero temperature has been investigated [15][16][17][18], whose purpose is to understand the effect of CR terms of both the spin-bath interaction and the driving, on the dynamics. This paper extends this analytical method to the non-Markovian dynamics in finite temperature, discuss the influence of temperature and SBC strength on the non-Markovian dynamics, and show the interference between the driving and dissipation in the non-Markovian dynamics.
Firstly, we consider the non-Markovian dynamics in zero temperature T = 0 and ω c = 100Δ. Figure [15]. From figures 2(a) and 3(a), no matter the SBC strength is α = 0.01 or α = 0.1, the period given by NM-NRWA is the longest and that given by M-RWA is the shortest. The amplitude of steady state P(t) is 0. In figures 2(b) and (c) with weak SBC α = 0.01, compared with the significant effect of CR terms on P(t) dynamics, the effect of Markovian approximation is smaller. However, the effect from the non-Markovian dynamics can be accumulated with the evolution time, which will be discussed in the following steady-state behavior.
When the SBC strength increases to α = 0.1, it can be seen from figures 3(b) and (c) that there is a significant difference between the results of non-Markovian and Markovian. And, the oscillation period and amplitude of non-Markovian are significantly greater than those of Markovian. Recently, the trace distance is proposed to measure the non-Markovian effect based on the information exchange between an open system and its environment [33][34][35]. Here we need to clearly distinguish the non-Markovian effect and the non-Markovian approach. In this paper, the non-Markovian dynamics means that we use the non-Markovian method to solve the quantum dynamics of DSBM. At the same time, there are some important theoretical works recently, to discuss the rich non-Markovian quantum dynamics of open systems, paying particular attention to the rigorous mathematical definition of quantum non-Markovianity (quantum memory effect), to the physical interpretation and classification, as well as to the quantification of quantum non-Markovianity [33][34][35]. Quantifying the quantum non-Markovianity quantum memory effect of the dynamics is a very meaningful topic. Here, we calculate the trace distance by the three approaches of NM-NRWA (solid line), M-NRWA (short-dot lines) and M-RWA (short-dash-dotted line) to quantify the memory effect. Figure 6 shows the time-dependent trace distance for different values of the SBC strength α and temperature T. For weak and intermediate SBM coupling α < 0.1, the trace distance exhibits an overall decay to zero with periodic modulation as time increases, which indicate the presence of memory effect. The results shows that the oscillation period of NM-NRWA is the longest at zero temperature.
In order to characterize the accumulation of non-Markovian dynamics over time, we study the amplitude and phase of the long-time driving oscillation of P(t), which are not easily obtained by numerical methods. The steady-state dynamics is defined as [15] P(t → ∞) = −P lt cos (ω L − Φ), with the amplitude P lt and phase shift Φ. Figures 7(a1) and (b1) shows the amplitude P lt and phase shift Φ versus Rabi driving Ω/Δ at zero temperature with different SBC: α = 0.01, α = 0.05 and α = 0.1. For comparison, the corresponding curves of M-RWA are given in figures 7(a2) and (b2). The results from the three methods (NM-NRWA, M-NRWA and M-RWA) are obviously different. And these distinctions is growing with the increase of the SBC strength α.
In figure 7(a1), at zero temperature, the curve of steady-state amplitude P lt with Rabi driving amplitude Ω is N-like shaped, which is a nonlinear behavior and reflects the interference between driving and dissipation within the non-Markovian frame. That is to say, in the case of the same SBC strength, the steady-state amplitude first increases, then decreases and then increase with the Rabi driving amplitude. In each line of figure 7(a1), we mark two local extreme points as point 1 and point 2, and let their coordinates are (Ω 1 , P 1 lt ) and (Ω 2 , P 2 lt ), which represent the maximum amplitude of coherent oscillation P 1 lt at lower driving Ω 1 and the minimum amplitude P 2 lt at higher driving Ω 2 , respectively. The maximum amplitude P 1 lt increases with the increase of SBC strength α, and the driving strength Ω 1 required to get P lt also increases.  In figure 8(a1) at finite temperature, the steady-state amplitude P lt with Rabi driving amplitude Ω is Λ-like shaped. That is to say, P lt first increases and then decreases with the Rabi driving amplitude. So only the local extreme point 1 exists. As we all know, the driving maintain the coherence of the spin, while the dissipation of bath destroy the coherence. Under the interference of driving and dissipation, the renormalized tunneling ηΔ in the case of zero temperature is nonzero. But ηΔ in the case of finite temperature is suppressed from ηΔ = 0 to ηΔ = 0 with the increase of temperature [31]. The presence of a critical temperature T c corresponds to a transition from ηΔ = 0 to ηΔ = 0. Similar results can be found via adiabatic renormalization [32] from Monte Carlo simulations. The phase shift Φ of steady state of finite temperature is shown in figure 8(b1). For comparison, the corresponding amplitude P lt and phase shift Φ of M-RWA in finite temperature are given in figures 8(a2) and (b2). The results from the three methods (NM-NRWA, M-NRWA and M-RWA) are obviously different.
In figure 7(a1), we also mark another two Rabi driving amplitudes Ω A and Ω B . When the driving amplitude Ω is smaller than Ω A , P lt of weak SBC α = 0.01 is the largest, and P lt with α = 0.1 is the smallest. However, when the driving strength becomes larger than Ω B , the large SBC α produce large P lt . At finite temperature in figure 8(a1), we find that P lt is much smaller than that of zero temperature. Moreover, in the case of intermediate driving, P lt of α = 0.05 is higher than that of α = 0.1 and α = 0.01, which is different from the result of zero temperature. Therefore, in order to obtain large amplitude stable quantum coherent oscillation, we need to optimize the SBC strength, temperature and Rabi driving amplitude in DSBM.

Conclusion
The non-Markovian dynamics of DSBM is investigated by the method based on unitary transformation. After the unitary transformation, the spin is dressed by environmental boson and the strength of the SBC is effectively weakened and the thermal average of the interaction term is zero, then the master equation without Markovian approximation is obtained by applying perturbation theory to the transformed Hamiltonian. This approach extends the regime of validity of the quantum master equation to stronger SBC. The influence of the temperature, SBC strength and Rabi driving amplitude on the transient and steady state dynamics are discussed. From the transient dynamics, when the SBC interaction goes beyond weak regime (α > 0.01), the results of NM-NRWA are quite different from those of RWA and Markovian approximation. The curve of steady-state amplitude versus Rabi driving amplitude is nonlinear and reflects the interference between driving and dissipation. Therefore, by optimizing the driving strength, temperature and SBC strength, a large amplitude of stable coherent oscillation can be obtained. Finally, we hope that this work will stimulate more experimental and theoretical works to investigate the non-Markovian dynamics at finite temperature for quantum information and quantum computation.
We get the equation of motion of the DSBM in the finite temperature,

Appendix B. The expressions of the matrix element in the super-operator
This appendix gives the concrete expression of the matrix element in the super-operatorL (s),

Appendix C. Inverse Laplace transformation of the matrix element
In this appendix, we give details of how to inverse the Laplace transformation in equations (29)- (34). As an example, ds.