Exact quantum dissipative dynamics under external time-dependent fields driving

Exact and nonperturbative quantum master equation can be constructed via the calculus on path integral. It results in hierarchical equations of motion for the reduced density operator. Involved are also a set of well--defined auxiliary density operators that resolve not just system--bath coupling strength but also memory. In this work, we scale these auxiliary operators individually to achieve a uniform error tolerance, as set by the reduced density operator. An efficient propagator is then proposed to the hierarchical Liouville--space dynamics of quantum dissipation. Numerically exact studies are carried out on the dephasing effect on population transfer in the simple stimulated Raman adiabatic passage scheme. We also make assessments on several perturbative theories for their applicabilities in the present system of study.


Introduction
The central problem of quantum dissipation theory is to study the dynamics of quantum system embedded in quantum thermal bath. The primary quantity of interest here is the reduced density operator, ρ(t) ≡ tr B ρ T (t), after the bath degrees of freedom are all traced out from the total composite density operator. Due to its fundamental importance, quantum dissipation theory has remained as an active topic in diversified fields [1][2][3][4][5][6][7][8][9][10]. The challenge here, from both formulation and numerical aspects, is nonperturbative dissipation, with multiple time scales of memory, under time-dependent external field driving.
For Gaussian stochastic force, the influence of bath on system can be characterized by force-force correlation functions. Exact formalism had then been established via the Feynman-Vernon influence functional approach [1][2][3][4][5]. Direct numerical integration methods, based on discretization of the path integral and summation up of the memory correlated terms, have been put forward such as the quasi-adiabatic propagator method [11][12][13][14] or the real-time quantum Monte Carlo scheme [15][16][17][18][19]. ‡ Author to whom the correspondence should be addressed.
The alternative is the differential approach, especially in a linear form to maximize the numerical advantage. It has also the advantage in the study of various dynamics such as the spectroscopic or control problems [20]. The calculus-on-path-integral (COPI) method is hence proposed to construct the differential counterpart of the path integral theory, reported as the hierarchical equations of motion (HEOM) formalism [21][22][23][24][25][26][27][28][29][30]. This formalism can also be derived via the stochastic description of quantum dissipation [31][32][33][34][35]. The COPI algorithm provides a unified approach to the influence of quantum environment ensembles, either canonical or grand canonical, and either bosonic or fermionic [28,29]. The COPI algorithm also takes into account the combined effects of multiple memory time scales, system-bath coupling strengths, and system anharmonicity. The resulting HEOM formalism is therefore nonperturbative in nature, and always converges in principle. Moreover, the HEOM formalism is exact, not just for its propagation equivalent to the path integral theory, but also for the fact that the initial correlations between system and bath can now be incorporated by the steady-state solutions to the HEOM, before external time-dependent fields taking effect. Recently, we have further developed a numerical efficient filtering method for the propagation of the HEOM [36,37].
In this work, we report a HEOM-based study on population transfer with dephasing in the scheme of stimulated Raman adiabatic passage (STIRAP) [38]. The laser control of dissipative systems has been addressed extensively [39][40][41][42][43][44][45][46][47][48][49], but mostly on the basis of weak dissipation treatment. The correlated influence of driving and dissipation is often important, as demonstrated previously [50,51]. With the aid of the numerically exact results, we analyze the dephasing effects on transfer dynamics in relation to the STIRAP mechanism and examine some second-order quantum dissipation theories for their applicabilities in the systems of study.
The remainder of paper is organized as follows. We present the HEOM formalism together with comments on its numerical implementation in Sec. 2, and the derivations in Appendix. In Sec. 3, we study the dephasing effect on population transfer dynamics in the STIRAP scheme. We report the numerically exact results via the HEOM formalism, followed by discussions in relation to the STIRAP mechanism. In Sec. 4, we present the details of numerical performance of the HEOM results, and make concrete assessments on several approximated quantum dissipation theories. Finally we conclude the paper.

Description of stochastic bath coupling
The total system-plus-bath Hamiltonian can be written in general as The last term denotes the multi-mode system-bath interactions. The involving system operators {Q a } are called the dissipative modes, through which the generalized Langevin forces {F a (t) = e ihBtF a e −ihBt } from the bath (h B ) act on the system. For convenience, let the dissipative modes be dimensionless. The time dependence in the system H(t) arises from external driving fields. Throughout this paper, we denote the inverse temperature β ≡ 1/(k B T ) and seth ≡ 1.
We treat the Langevin forces as Gaussian stochastic processes. Therefore their effects on the system are completely characterized by the correlation functions, Here, Ô B ≡ tr B (Ôρ eq B ) denotes the thermodynamics average over the canonical ensembles of the bosonic bath. The correlation functions satisfy the symmetry and detailed-balance relations, or equivalently the fluctuation-dissipation theorem [3,20]: with J ab (ω) = −J ba (−ω) = J * ba (ω) being the bath spectral density functions. The HEOM formalism requires C ab (t) be expanded in certain series form, so that the hierarchy can be constructed via consecutive time derivatives on path integral. Various schemes [22,28,52,53] have been proposed to expand C ab (t) in exponential series, on the basis of analytical continuation evaluation of Eq. (3). In particular, the hybrid scheme that also exploits quadrature integration method is applicable for arbitrary spectral density functions [30].
For simplicity we set C ab (t) = C aa (t)δ ab . In this case the contributions from different dissipative modes {Q a } are additive. Without loss of generality, we present the formalism explicitly only for the single-dissipative-mode case, Q a = Q. We thus omit the index a for clarity of formulation. We also adopt the super-Drude model, The corresponding correlation function can be analytically evaluated as [20,28,53] All coefficients here are real and given in Appendix [cf. Eqs. (A.9) and (A. 10)]. The first term arises from pole of the spectral density function, which is of rank two. The second term is from the Matsubara poles, withγ m ≡ 2πm/β being the Matsubara frequency. The last term is the Matsubara residue, which would approach to zero if M → ∞. In this work, we adopt the Markovian residue ansatz [25,35], i.e.,γ m e −γmt m>M ≈ δ(t); thus,

The HEOM formalism
The dynamics quantities in the HEOM formalism are the reduced density operator ρ(t) and a set of auxiliary density operators (ADOs), {ρ n (t)}, that hierarchically resolve the memory contents of the bath correlation functions in the exponential series expansion of Eq. (5). The index n that specifies an N th -tier ADO ρ n consists of a series of nonnegative integers, n ≡ {n, n ′ ,n,n ′ ,ň 1 , · · · ,ň M } , with n + n ′ +n +n ′ +ň 1 + · · · +ň M = N. (7) Comparing to the reduced density operator ρ(t) ≡ ρ 0 (t) of primary interest, the specified ρ n would have the order of |ν| n+n ′ · |ν r + iν i |n +n ′ · M m=1 |ν m |ň m , for its dependence on the individual components of interaction bath correlation functions in the series expansion of Eq. (5). These scaling factors will be incorporated properly in the final dimensionless ρ n , in order to validate a filtering algorithm for the numerical efficiency of the HEOM formalism. On the other hand, the indices in the set n of Eq. (7) cover all accessible derivatives of the Feynman-Vernon influence functional; see Appendix for the details.
The final HEOM formalism is summarized as follows. It has the generic form oḟ Here, L(t)ρ n ≡ [H(t), ρ n ], which depends in general on the external driving fields; Γ n is the damping parameter that collects all related exponents, and δR is the residue dissipation superoperator due to δC(t). For the bath correlation function in the series expansion, Eq. (5) with Eq. (6), they are given respectively by The last three terms in Eq. (8) denote how the specified N th -tier ADO ρ n depends on other ADOs of the same tier, the (N −1) th -tier, and the (N +1) th -tier, respectively. For the bath correlation function in Eq. (5), they are given explicitly by Here, λ n = n|ν| ,λň m = ň m |ν m | ,λ r/i n = n|ν r/i | , and λ r/i n = γλ r/i n+1 n/|ν| , with the italic-font indices being from those in n of Eq. (7). The indexes variations in Eq. (10) that specify those ADOs participating in the equation ofρ n are exemplified as follows: Similarly, n ′ differs from n of Eq. (7) only by changing (n ′ ,n ′ ) to (n ′ + 1,n ′ − 1), whilě n ± m by changingň m toň m ± 1, and so on. Also note that ρ n is an N th -tier ADO, while ρ n ± is of an (N ± 1) th tier, as inferred from the second identity of Eq. (7).
The initial conditions to the HEOM in the study of driven dissipative dynamics are obtained via the steady-state solutions to Eq. (8), before the time-dependent external fields interactions. For the steady-state solutions satisfyingρ st n = 0, Eq. (8) reduces to a set of linear equations, under the constraint of Trρ 0 = 1. The resulting ρ st n is used as the initial ρ n (t 0 ) to the HEOM. The initial system-bath correlations are accounted for by those nonzero initial ADOs.

Comments on numerical implementation
For the numerical HEOM propagation, we would like to have certain convenient working index scheme to track the multiple indices, denoted now as an ordered set of n = {n 1 , · · · , n K }, that specifies ρ n . Here we will provide two such schemes. The number of the N th -tier ADOs, with n 1 + · · · + n K = N , is (N +K−1)! N ! (K−1)! ≡ N K . In one scheme, the ADOs are arranged as ρ n ≡ ρ jn with j n initialized by j n=0 ≡ 0 and then Let L be the maximum level of the hierarchical tier. The total number of the ADOs In another scheme, ADOs can also be arranged as ρ n ≡ ρ ln ; l n = 0, · · · , N − 1, with Both schemes, Eq. (12) and Eq. (14), allow easy tracking of the coupled ADOs in the HEOM. The former [Eq. (12)] is somewhat more convenient in the filtering propagator described soon since it does not depend on L.
The major difficulty in implementing the HEOM formalism is its numerical tractability. The number of ADOs, N of Eq. (13), itself alone can be huge in the case of strong non-Markovian system-bath coupling and/or low temperature, as large L and/or large K implied. Thus, a brute-force implementation is greatly limited by the memory and central processing unit (CPU) capability of computer facility, even for a two-level system where each ADO is a 2 × 2 matrix.
To facilitate this problem, Shi, Xu, Yan and coworkers have recently proposed an efficient numerical filtering algorithm that often reduces the effective number of ADOs by order of magnitude [36,37]. In reality, there is usually only a very small fraction of total ADOs significant to the reduced system dynamics. To validate the accuracy-controlled numerical filtering algorithm, the present HEOM formalism has been scaled properly so that all ADOs {ρ n (t)} are of a uniform error tolerance. This remarkable feature is suggested by comparing the HEOM theory with the stochastic bath interaction field approach in the case of Gaussian-Markovian dissipation [37]. The involving ADOs are just the expansion coefficients, over the normalized harmonic wave functions that are used as the basis set for resolving the diffusive bath field [37]. Our numerical HEOM propagator exploits the filtering algorithm [36]. It goes simply as follows. If a ρ n (t) whose matrix elements amplitudes become all smaller than the pre-chosen error tolerance, it is set to be zero. Apparently, the filtering algorithm also automatically truncate the required hierarchy level on-the-fly during numerical propagation. By far the truncation for the Matsubara expansion still goes by checking convergency.

Numerical results
The STIRAP is celebrated as an efficient and robust method for population transfer [38]. It is characterized by its counterintuitive field configuration. For a three-level Λ-system as Fig. 1, the Stokes pulse proceeds the pump pulse, and the intermediate state remains effective in dark. The STIRAP mechanism [38] is rooted at the existence of the coherent population trapping state under the two-photon resonance condition, ω S − ω P = ǫ 1 − ǫ 3 , in the Λ-system. Dephasing destroys this condition in terms of resonance and/or the existence of coherent population trapping state. The effect of dephasing on the simple STIRAP scheme has been studied extensively, but with approximations. These include phenomenological/perturbative methods [44][45][46], or classical/stochastical bath treatments [47][48][49].
We revisit the dephasing effect on the simple STIRAP-based population transfer, as the exact dissipative dynamics are now established with the present HEOM formalism. We also examine three schemes of second-order approximation [20,[52][53][54][55]: (i) The Redfield theory, which neglects the correlated driving-and-dissipation effect; (ii) CS-COP, which is the conventional time-nonlocal quantum master equation, including the field-dressed dissipation contribution, and equivalent to the present HEOM truncated at the first tier; (iii) CODDE, in which the driving field-free part of dissipation superoperator is time-local, while the field-dressed part is time-nonlocal. Neglecting the latter leads it to the Redfield theory.
The total Hamiltonian under the rotating wave approximation assumes Here,D P ≡ |1 2| + |2 1| andD S ≡ |2 3| + |3 2|, while Ω P (t) and Ω S (t) denote the Rabi frequencies of the resonant pump and Stokes fields, respectively. The dissipative mode Q a = |a a| is responsible for dephasing. The interaction spectral density function J a (ω) ≡ (π/2) j [c 2 aj /(m j ω j )]δ(ω − ω j ) assumes super-Drude as Eq. (4). The system Hamiltonian is then The Caldeira-Leggett renormalization energy [56,57] is δǫ a = 1 π ∞ 0 dω J a (ω)/ω = We set the pump and Stokes fields be of same Gaussian shape, Ω P (t + t P ) = Ω S (t + t S ) = A exp[− 1 2 (wt) 2 ], but center them at t P = 200 β and t S = −200 β, respectively and counter-intuitively. The driving strength and inverse duration parameters are set to be βA = 0.1 and βw = 0.005. The corresponding dissipationfree transfer dynamics is shown in Fig. 1(b). As here the bath influence is considered to be pure-dephasing in the absence of fields, the initial system is just chosen to be completely on the |1 state and all the ADOs are zero. For the effect of bath, we set the coupling strength η = 0.64 [cf. Eq. (4)], and consider both the Markovian and non-Markovian cases as follows.
The Markovian transfer dynamics, under the influence of single dephasing mode of either Q 1 = |1 1| or Q 2 = |2 2|, is exemplified in Fig. 2, with βγ = 5. We observe: (i) The Q 1 -mode effect shown in Fig. 2(a) leads to all three populations about 1/3 after the driving; (ii) The Q 2 -mode effect shown in Fig. 2(b) is less sensitive than its Q 1 counterpart, achieving a higher transfer efficiency, despite it is only about 0.55.
The non-Markovian transfer dynamics is exemplified in Fig. 3, with βγ = 0.5. In comparison with the Markovian counterparts, we observe: (iii) The Q 1 -mode case behaves about the same; (iv) But the Q 2 -mode results in a higher transfer yield, increasing to about 0.73 via the exact calculation. We have also calculated the influences of Q 3 = |3 3| for both Markovian and non-Markovian cases. The results (not shown here) are similar to those of Q 1 , except for some small oscillations. Two double-modes (Q 1 +Q 2 and Q 1 +Q 3 uncorrelated) non-Markovian dephasing dynamics are shown in Fig. 4(a) and (b), respectively. They are insensitive to the non-Markovian parameter, and both reach at the final equal-populations, based on the numerically exact results. Comments on the approximated schemes, the CODDE, CS-COP, and Redfield theory, presented in Figs. 2-4 will be given later; see Sec. 4.2.

Discussions
The above observations can be understood by the well-established STIRAP mechanism [38]. The Q 1 -mode, which associates with the fluctuation of level |1 , easily destroy the two-photon resonance (TPR) condition, as described at the beginning of Sec. 3. Thus, it ends up with the observed equal populations in all accessible levels by the strong fields, as consistent with the analysis in [44]. The similarity between the Q 3 and Q 1 influences is also explained. The same reason accounts further for the case of uncorrelated two modes (Q 1 + Q 2 or Q 1 + Q 3 ) dephasing, as depicted in Fig. 4. It is anticipated that when γ ≪ w (termed as the linear adiabatic limit below), the equal-population will be broken to be in favor of |3 , due to the marginally partial fulfilment of the TPR condition.
On the other hand, the Q 2 -mode is associated with the fluctuation of the intermediate level |2 . It alone does not affect the TPR condition. However, this condition, based on the numerically exact results shown in this work, is not sufficient to retain the coherent population trapping state, chosen ad hoc earlier for the dephasingfree STIRAP scenario in Fig. 1(b). It is anticipated that the coherent population trapping state may be recovered in the aforementioned linear adiabatic limit. This is in line with the observation-(iv), where the non-Markovian population transfer with single Q 2 -mode dephasing [ Fig. 3(b)] is of higher efficiency than its Markovian counterpart [ Fig. 2(b)]. The previous study based on perturbative dephasing dynamics [44] has also shown that the single Q 2 -mode does not affect the transfer efficiency in the linear adiabatic limit. Nevertheless, STIRAP in the presence of complex dephasing, if 100% transfer ever achievable, would require dynamics feedback control of pump or Stokes laser frequency [59]. This would involve chirp and realize STIRAP in a nonlinear adiabatic condition, rather than the linear simplification considered here.

Numerical performance of the HEOM formalism
The numerical performance of the HEOM formalism with filtering is summarized in Table 1, for the systems reported in the three figures' (b)-panels. The CPU time is for a single Intel(R) Xeon(R) processor@3.00GHz to calculate the exact result in each (b)-panel for the time period −1200 β < t < 2000 β with the time step dt = 0.01β using the fourth-Order Runge-Kutta propagator; N max denotes the largest number of active ADOs and L max the highest tier level, ever survived in the entire time span of the numerical propagation. The filtering error tolerance is chosen to be 10 −6 , following our previous work [36]. We input M = 6 for the number of Matsubara terms being explicitly included, which has been tested to give converged results of ρ(t) = ρ 0 (t) in all calculations. The total number N of mathematical ADOs follows Eq. (13) and is given inside the parentheses. The effect of filtering is clearly seen. The number of active ADOs with filtering is insensitive to the input M , as long as it is large enough.
In the present study, the number of active ADOs reaches N max only during the period about −250 β < t < 500 β and grows up or drops down dramatically outside that period with the fields turning on or getting over. Apparently, N max increases with the number of dissipative modes. At least one (L max ) th -tier ADO actively participates during the HEOM propagation. Its leading contribution to the reduced density operator is of (2L max ) th order in the system-bath interaction. Physically, L max is closely related to the modulation κ-parameter [27], introduced originally by Kubo for motional narrowing problem [6]. This dimensionless parameter is determined via κ ≡ γ/ √ ν, or similar, for each individual exponential component in Eq. (5). The last two columns of Table  1 are the modulation parameters κ andκ m=1 of the leading Matsubara term. The modulation κ-parameter relation to the value of L max [27] can be clearly seen. In both the Markovian and non-Markovian cases of the present study,κ m =γ m / |ν m | monotonically increases with m, cf. Eq. (A.10). Actually, the Matsubara series truncation M in Eq. (5) can be estimated via its reaching the fast modulation condition,κ M ≫ 1. As the temperature decreases,κ m getting smaller, and eventually cause the value of L max be pretty large. The present HEOM construction is based on the Matsubara series expansion, which may no longer be numerically implementable in the extremely low temperature regime. Alternative expansion method such as the hybrid scheme [30] is needed to the required HEOM construction.

Assessments on three second-order approximated theories
With the exact results, we can now make concrete assessments on the three secondorder approximated schemes, the Redfield theory, CS-COP, and CODDE, exploited in the numerical demonstrations. The dissipative modes {Q a = |a a|} considered in Sec. 3 are all of pure dephasing, in the absence of external fields. The Redfield theory would be exact if there were no correlated driving-and-dissipation effect [20,54]. Therefore, the non-Markovian dynamics manifest here as the correlated driving and dissipation. Apparently, the Redfield theory, by its Markovian nature, is independent of the width γ-parameter of bath spectral density. Observed is also the fact that the schemes of approximation are sensitive to the Q 2 -mode rather than the Q 1 -or Q 3mode dephasing. This fact is also easily understood, by considering their relations to the STIRAP mechanism as discussed earlier. Further remarks on the approximated theories for their applicabilities in the systems of study are as follows.
The CS-COP theory [20,53,55] is overall most unsatisfactory, despite it contains formally a description of memory and driving-and-dissipation correlation. Even in the non-Markovian Q 2 -mode case of Fig. 3(b) where it appears to be superior than the Redfield theory, the CS-COP results in a decreased transfer efficiency from its Markovian counterpart shown in Fig. 2(b). This is qualitatively contradictory to the physical anticipation, as discussed earlier.
The CODDE [20,53,55] appears to be the most favorable perturbation theory. It gives the best approximated transfer dynamics in all cases presented in Sec. 3, except the one to be discussed soon. Its overall superiority is also true in the driven Brownian oscillator systems [20,51]. The CODDE is actually a modified Redfield theory, with inclusion of correlated driving-and-dissipation effects. The involving fielddressed dissipation kernel is time-nonlocal but constructed with a partial ordering resummation, rather than the chronological ordering prescription that characterizes the CS-COP [20,53,55].
The only exception is the Markovian Q 1 -mode case shown in Fig. 2(a), where the Redfield dynamics is almost exact. The reason for this exception is also accountable. As we mentioned earlier, the non-Markovian dynamics manifest as the correlated driving and dissipation. This correlated effect diminishes in both the fast-and slow-modulation regimes, as inferred from the exact and analytical results of driven Brownian oscillator systems [58]. This conclusion can be carried over to the present system of study, as suggested here. Apparently, the identical value of βγ = 5, adopted in the two cases of Fig. 2, acquires the fast-modulation limit for the Q 1 -mode, but not yet for the Q 2 -mode. In the latter case, the CODDE resumes its superiority.

Closing remarks
In summary, we have presented a hierarchical Liouville-space approach, which is exact and also quite tractable numerically, to general quantum dissipation systems under external driving fields. The auxiliary density operators are all of a uniform error tolerance, as that of the reduced density operator. We also make comments on the numerical facilitation about the multiple-index assignment and the filtering algorithm. We numerically study the dephasing effects on the population transfer, with a fixed simple STIRAP configuration, and present a concrete assessment on various approximation schemes.

Appendix. Construction of HEOM via the COPI approach
This appendix gives the details of the COPI approach to the HEOM formalism. It starts with the influence functional in path integral. Let {|α } be a basis set in the system subspace and set α ≡ (α, α ′ ) for abbreviation. Denote the evolution of reduced system density operator in the α-representation by Here, the reduced Liouville-space propagator is The effects of bath on the reduced system are contained completely in the influence functional F . For Gaussian stochastic forces {F a (t)} from fluctuating bath, it assumes the Feynman-Vernon form [1], which can be recast as [27][28][29]: and (A.6) Here, C ab (t) is the bath correlation functions, defined by Eq. (2). The functional A a [Eq. (A.4)] depends only on the local time and its operator level form is just the commutator of the dissipative mode Q a . The functional B a [Eq. (A.5)] does however contain memory, which is to be resolved via the COPI algorithm of consecutive time derivatives on all memorycontained functionals. To construct a close set of HEOM via the COPI algebra, it needs a proper expansion of C ab (t) such as the exponential series, while maintaining the fluctuation-dissipation theorem of Eq. (3). A super-Drude parametrization scheme and the resulted HEOM formalism have been presented in our previous work [28]. A hybrid scheme, exploiting the analytical continuation and quadrature integration methods to evaluate the fluctuation-dissipation theorem, has also been proposed [30].
To proceed we denote for every distinct exponent terms in Eq. The scaling factor s n is defined as |ν ab r | n ab |ν ab i | n ′ ab |ν ab r |n ab |ν ab i |n ′ ab n ab ! n ′ ab !n ab !n ′ ab !