Non-Markovian correlation functions for open quantum systems

Beyond the conventional quantum regression theorem, a general formula for non-Markovian correlation functions of arbitrary system operators both in the time- and frequency-domain is given. We approach the problem by transforming the conventional time-nonlocal master equation into dispersed time-local equations-of-motion. The validity of our approximations is discussed and we find that the non-Markovian terms have to be included for short times. While calculations of the density matrix at short times suffer from the initial value problem, a correlation function has a well defined initial state. The resulting formula for the non-Markovian correlation function has a simple structure and is as convenient in its application as the conventional quantum regression theorem for the Markovian case. For illustrations, we apply our method to investigate the spectrum of the current fluctuations of interacting quantum dots contacted with two electrodes. The corresponding non-Markovian characteristics are demonstrated.


I. INTRODUCTION
Open quantum systems, which are of great importance in many fields of physics, refer to a quantum system of primary interest coupled to an environment often called reservoir or bath. [1][2][3][4][5][6] The composite Hamiltonian (H tot ), in general, includes the system (H S ), the bath (H B ), and the coupling between the system and the bath (H ′ ), i.e., H tot = H S + H B + H ′ . It is well-known that the system is described by the reduced density operator, ρ(t) ≡ tr B [ρ tot (t)], i.e., the partial trace of the total density operator ρ tot over the bath space. The corresponding dynamics is determined by the master equation, where the effect of the bath is described by the second term with the self-energy Σ(t − τ ). It contains the memory effect in principle even for weak system-reservoir coupling. Eq. (1) is thus the so-called time-nonlocal (TNL) master equation describing non-Markovian dynamics. As long as one knows the reduced density operator ρ(t), a single-time expectation value of an arbitrary physical observable of the system, e.g.,Ô, is simply obtained via O(t) = Tr[Ôρ(t)]. However, it is not as easy to calculate two-or multiple-time correlation functions. Except for a small number of exactly solvable systems, convenient calculation of the correlation function is possible using the well-known quantum regression theorem which is valid in the Born-Markovian approximation. [1][2][3][4][5]7,8 It is not valid any more for non-Markovian cases, 5,[9][10][11][12] due to memory effects which break the time translation invariance of the non-Markovian propagator.
The calculation of non-Markovian correlation functions for arbitrary system operators is a challenge and long-standing problem. Stimulated not only by fundamental interest but also by great demand because of the rapid progress in experiments which are able to access non-Markovian effects, [13][14][15] there are many efforts to address this issue. For example, by using the stochastic Schrödinger equation approach and the Heisenberg equation, [16][17][18] rather than master equation, the method derived by Alonso and de Vega [16][17][18] is valid for zero-temperature environments and/or Hermitian system-environment coupling operators. Based on a generalized Born-Markov approximation, Budini 19 derived a quantum regression theorem which is applicable for non-Markovian Lindblad equations. Recently, Goan et al., 20,21 developed a scheme for the calculation of twotime correlation functions of the system operators with memory effects in terms of the "time-convolutionless" master equation. Additionally some specific systems have been analyzed [22][23][24] .
In this work, we aim to derive a general formula for non-Markovian correlation functions of arbitrary system operators in terms of the TNL-ME Eq. (1). We will consider weak system-reservoir coupling but short time-scales where non-Markovian effects should dominate. Later we will also analyze the relevant time scales in more detail. By using the fluctuation dissipation theorem of the correlator and introducing the auxiliary density operator in the frequency domain denoted by φ ± (ω, t), it is easy to transform the TNL-ME Eq. (1) into an equivalent set of coupled time-local equationsof-motion (for short TL-EOMs, expressed in Eq. (12)), 25 i.e.,˙ ρ(t) = Λ ρ(t), in the enlarged vector space with T . The corresponding propagator in this enlarged vector space satisfies timetranslation invariance and accordingly the correlator can be treated similar to the Markovian case. The resulting equation takes a form very similar to the quantum regression theorem in a larger space, The most important feature of this approach is the initial condition ρ B (0) = B ρ(0), where ρ(0) is the density matrix which has been time evolved from a initial time t 0 → −∞. We assumed that system and reservoir decouple at the initial time t 0 → −∞ which is the standard assumption for the TNL-ME. In the following we will derive equation (2) and we will discuss in detail the range of validity using the diagrammatic expansion on Keldysh contour 26 . This full non-Markovian description is applicable for both, fermionic and bosonic systems. As an example, we will discuss non-Markovian effects of the electronic reservoir on the current-fluctuation spectrum in quantum dots. The paper is organized as follows. In Sec. II, we first introduce the conventional TNL-ME for weak systemreservoir coupling and then outline the equivalent TL-EOMs. Based on TL-EOMs, we get the formulas for the two-time non-Markovian correlation functions both in the time-domain and in the frequency-domain in Sec. III. We then implement the proposed scheme to study the current-fluctuation spectra of the electron tunneling through quantum dots in Sec. IV. Finally, we conclude in Sec. V.

A. Time non-local master equation
The reservoir with infinite degrees of freedom is described by the non-interacting Hamiltonian H B = k ǫ k c † k c k with the creation (annihilation) operator c † k (c k ). The coupling Hamiltonian between the system and the bath, in general, is given by where (Q + ) † = Q − and (F − ) † = F + , with the operator of the central system Q ± and the operator of the bath F ± . The coupling operator of the bath is defined as F − = k t k c k and contains the coupling coefficients t k . The result will be generalized for the case of coupling to multiple reservoirs in the Appendix B. The Hamiltonian of the small system is composed of the corresponding creation (a † µ ) and annihilation (a µ ) operators, i.e., H S ≡ H S (a µ , a † µ ) which could include many-body interaction terms.
Assuming weak system-bath coupling and performing Born but without Markovian approximation, the self-energy for the expansion up to second-order of the coupling Hamiltonian is expressed as The corresponding diagram is schematically shown in Fig. 1(b). Here, · · · B stands for the statistical average over the bath in thermal equilibrium. The explicit formalism for the self-energy in Eq. (1) Figure 1: The diagrams of (a) Dyson equation and (b) the self-energy with the lowest-order contributions. 26 In (a), Π0 is the free propagator with Π0(t, t0) = e −iL S (t−t 0 ) , and the dashed line in (b) is the correlation function C (±) (t − τ ) of the bath expressed in Eq. (6). thus reads, where we introduced the free propagator with the correlation function of the bath, Consequently, the TNL-ME Eq. (1) is explicitly given by ( = 1), [27][28][29] (7) Note that Eq. (7) is derived assuming initial decoupling at t 0 → −∞, ρ tot (t 0 ) = ρ(t 0 )ρ B and ρ B is the equilibrium state of the bath. The corresponding propagator Π(t, t 0 ) for Eq. (7)( Eq. (1)) defined by satisfies the Dyson equation 26 as shown in Fig. 1 (a) and does not satisfies time-translation invariance, i.e. Π(t, t 0 ) = Π(t, t 1 )Π(t 1 , t 0 ). The conventional quantum regression theorem is thus broken.

B. Time-local equations-of-motion
The key to the calculation of the non-Markovian correlation function, e.g., Â (t)B(τ ) , is how to expand the master equation to an extend space ρ → ρ which preserves again time translation symmetry. This will allow us to case the correlator into the form of the regression theorem Eq. (2).
We adopt the multi-frequency-dispersed scheme 25,30 and define the bath correlation function Eq. (6) in the frequency-domain (C (±) (ω)) as where C (±) (ω) is directly related to the spectral density of the bath depending on the specific operator (F ± ). Correspondingly, the Liouville operator of Eq. (5) in the frequency-domain is [c.f.Eq. (A4)] Furthermore, we introduce the auxiliary density operators in the frequency-domain defined by (11) which means φ ± (ω, t 0 ) = 0 is applicable for the initially decoupled system-reservoir with t 0 → −∞ as we mentioned above. Taking the time derivatives of the auxiliary density operators, it is easy to recast TNL-ME of Eq. (7)( Eq. (1)) in the forṁ which are the so-called time-local equations-of-motion (TL-EOMs) due to the involved time-independent dissipative coefficients. TL-EOMs Eq. (12) is the lowest-tier truncation of hierarchical equation of motion 25 which has the linearity of the hierarchical Liouville space as demonstrated in Ref. 31. Introducing a vector composed of the reduced density operator and auxiliary density operators, i.e., the TL-EOMs Eq. (12) can then be further compacted with˙ Here according to Eq. (12), Λ can be formally written as where we introduced Apparently, Eq. (14) leads to ρ(t) = Π(t, t 0 ) ρ(t 0 ) with Π(t, t 0 ) = e Λ(t−t0) . In this vector space, the propagator satisfies the time-translation invariance, i.e., Π(t, t 0 ) = Π(t, τ ) Π(τ, t 0 ) and the correlation function 2 can be calculated straightforwardly in a form similar to the Markovian case based on the quantum regression theorem.

III. NON-MARKOVIAN CORRELATION FUNCTION
A. Two-time correlation function Using the vector of ρ(t) defined in Eq. (13), a singletime expectation value of the system operatorÂ can be expressed as with the initial condition being that system and reservoir Since time-translation invariance of the propagator has been restored, we can now follow the derivation of the Markovian correlation function based on the quantum regression theorem. 4 Therefore the non-Markovian twotime correlation function in the vector space is given by, where the components of the density operators in the vector are noŵ The initial state at t = 0 is given bŷ where ρ(0) is the density matrix which has been time evolved from a initial time t 0 → −∞.
A similar equation has been derived using linear response theory in Ref. 31. In this case all high-order contributions in the self-energy have been considered and one should keep all hierarchical EOMs for the numerical calculation of the correlation function. 25,31 In this paper we consider the lowest-order contribution for weak system-reservoir coupling. In terms of eqs. Eqs. (17)- (19) together with the TL-EOMs Eq. (12), (see Appendix A for the detail), we finally get with the steady-state ρ(0) = Π(0, t 0 )ρ(t 0 ) =ρ and ρ(t 1 ) = Π(t 1 , t 0 )ρ(t 0 ) =ρ for t 0 → −∞. Compared to the Markovian correlation function, the second modification in Eq. (20) arises from the memory effect and it is the vertex corrections which will be further illustrated in the coming subsection based on a diagrammatic representation. Similarly, it is easy to get the NMK-CF of Â (0)B(t) as expressed in Eq. (A7). In this vector space, Non-Markovian multiple-time correlation functions can be calculated by

B. Spectrum of the two-time correlation function
Since the widely measured quantity in experiments is the spectrum of the correlation function (in Fourier , we will now calculate the spectrum for the stationary two-time correlation function, (21) where the last identity is assumed to beÂ † =B, and Based on either Eq. (20) directly or Eq. (17), we finally obtain (for the detail derivation see the Appendix A), where Π(ω) and C and with the frequency-domain of the self-energy [Eq. (4)] and the frequency-domain of the bath correlation [c.f.Eq. (6) and Eq. (9)] It can be further give by where the real part C (±) (ω) is described by Eq. (9) and the imaginary part is Similarly, the spectrum of the correlation function of Â (0)B(t) , is given by Eq. (A8). The formulas Eq. (20) and Eq. (22) are the main contributions of the present work for the calculation of the non-Markovian two-time correlation function of arbitrary system operators (denoted byÂ andB) in time-domain and the frequencydomain, respectively.
Note that for often relevant Hermitian coupling operator, i.e., (20) and Eq. (22) is simplified to, and

C. Discussion and comments
Now we are in the position to discuss the applicability and the range of validity of the present NMK-CF formula Eq. (20) (or Eq. (22)). For convenience, we consider the coupling Hamiltonian H ′ = QX, with the operator of the bath X containing the coupling coefficient g, as an example. The resulting formula for non-Markovian correlation function is given by Eq. (26). The bath correlator in Eq. (6) is recast to where the second identity is written as a parametric decomposition, 32,33 with η m ∝ g 2 and 1/γ m representing the correlation time of the bath or the so-called the memory time. For the present considered weak systemreservoir coupling, the self-energy expressed in Eq. (4) contains the lowest-order contribution (or the first-order contraction by Wick theorem), e.g., Σ(t) ∝ C(t). Let the index "l" represent the number of contractions. The lowest-order contribution considered in the self-energy means l = 1.
In the following discussion we will not limit us to l = 1, but consider arbitrary number of the contractions such that the self-energy is exact. Based on the diagrammatic technique, 26 we can obtain (see Appendix C 1 for details), which is formally exact. The relevant diagrams are shown schematically in Fig. 2. Here, Σ B is the self-energy of vertex corrections including the operatorB and contain all the unseparable diagrams as well as the self-energy Σ in the propagator Π. It is worth noting that Eq. (28) has the same form as Eq. (26a)(or Eq. (20)) where Σ B can be extracted and Σ is expressed explicitly in Eq. (4) which has the same order of the magnitude to Σ B .
We will now use this formalism to discuss the range of validity of the non-Markovian correlation function. As has been shown in Ref. 34, in general non-Markovian effects are of the same order has higher order contractions in the self-energy. However, for the correlation function we have a well defined time scale and for short times scales, the combination of non-Markovian master equation and lowest order self-energy can be valid.
Let us make the Taylor expansion of the timederivative of the correlation function, i.e, G I (t) ≡ d Â (t)B(0) I dt in H S -interaction, for small t = 0 + , Following the estimation of the order of magnitude of the self-energy in Ref. 34, we roughly get (see Appendix C 2 for the detail), where γ is the smallest decay rate γ m of C(t) in Eq. (27), f (Q,B,ρ) and f (Q) are just the formal expressions arising from Σ B and Σ, respectively. Note that the magnitudes of f (Q)Bρ and f (Q,B,ρ) are nearly equal. We conclude that the expansion is valid for time t 1/γ, even if only the lowest order contributions l = 1 are considered.

IV. CURRENT FLUCTUATIONS OF THE ELECTRON TRANSPORT IN QUANTUM DOTS
For the demonstration of characteristic non-Markovian effects, we consider electron transport through quantum dots (QDs) contacted with two electrodes (left L and right R). This is a typical fermionic open quantum system where we will consider the non-Markovian effects in the spectrum in high-frequency regime. 23,27,35 The coupling Hamiltonian of Eq. (3) is specified by with the system operator Q ± µ being a † µ (a µ ) the creation (annihilation) operator of an electron of the µthlevel of the dot, and the operator of the reservoirs The correlator of lead α in the frequency-domain of Eq. (9) is thus , and β = 1/k B T the inverse of the temperature. For the studied model, we consider C with Ψ(x) the Digamma function.
Here we focus on the investigation of the non-Markovian current fluctuations in the central dots, i.e., Q (t)Q(0) withQ = e dN(t) dt andN = µ a † µ a µ . Since it satisfies the charge conservation ofQ(t) = − I L (t) + I R (t) , the corresponding spectrum is expected to be closely related to the noise spectrum of the transport current I α (t) in the reservoir α.

A. Single quantum dot
Let us first study the simplest model of electron transport through a spin-less one-level QD in the sequential tunneling regime, µ L > ε > µ R , described by the Hamiltonian H S = εa † a. Two states are considered in the dot, the empty state (|0 ) and the single-electron occupied state (|1 ), which leads to a = |0 1|. The numerical result is shown in Fig. 3. It depicts non-Markovian features (solid-line) compared to the Markovian results (dashed-line) (a) for charge fluctuation and (b) for current fluctuation, respectively. For low frequency regime at ω < min{|µ α − ε|} corresponding to the long time limit, the results based on both non-Markovian and Markovian treatments are consistent due to the disappearance of non-Markovian effect. With increasing the frequency higher than the energy-resonance, i.e., ω ω α0 ≡ |ε − µ α |, it enters the non-Markovian regime where the non-Markovian characteristic occurs in the spectra, showing steps at ω ≈ ω α0 in the current fluctuation spectrum (see Fig. 3 (b)). This is consistent with the noise spectrum of the transport current through the reservoirs studied in Refs. 27,35,41 The non-Markovian feature showing steps at ω α0 in S c (ω) provides the information of the energy structure in the central dot. The heights of the steps contain the information of the tunneling rate as demonstrated in the following. For the studied single-level dot in the regime of µ L > ε > µ R , the stationary population of the empty and single-electron occupied states are,ρ 00 = ΓR Γ and ρ 11 = ΓL Γ , respectively. Considering the spectrum in  Fig. 3 the positive frequency regime (ω > 0), it corresponds to the energy absorption processes. Accordingly, when the dot is in the empty state |0 with probabilityρ 00 , the electrons in the right reservoir absorb the energy, i.e, µ R + ω = ε and tunnel to the dot, which leads to the step When the dot is in the occupied state |1 with probabilityρ 11 , the electrons in the dot absorb the energy, i.e, ε + ω = µ L and tunnel to the left reservoir, which leads to the step at ω L0 = |µ L − ε| with the height of h L ≈ρ 11 Γ L = Γ 2 L Γ . The ratio of the heights in the positive frequency regime thus is, Similarly, in the reversed regime (ω < 0) corresponding to the energy emission processes, we get h ′ L ≈ρ 00 Γ L = ΓLΓR Γ and h ′ R ≈ρ 11 Γ R = ΓLΓR Γ , which leads to which is insensitive to the tunneling rate. We further consider the single level in the dot with spin-dependence as described by the Hamiltonian, H S = µ=↑,↓ ε µ a † µ a µ + Un ↑n↓ , wheren µ = a † µ a µ andN = µ=↑,↓n µ . The involved states in the dot are |0 , |↑ , |↓ , and |2 ≡ |↑↓ denoting the empty, two single-occupation spin states, and the double-occupation spin-pair state, respectively. In this state-basis, we have a ↓ = |0 ↓ | + | ↑ 2| and a ↑ = |0 ↑ | − | ↓ 2|. To demonstrate the Coulomb interaction effect more transparently, we fix the dot level with spin-degeneracy, ε ↑ = ε ↓ = ε, and consider spin-independent coupling strength, Γ α↑ = Γ α↓ = Γ α , (α = L, R). The corresponding spectrum of the current fluctuations with different Coulomb interaction is shown in Fig. 4. Besides the steps at the energy-resonance ω α0 , we also find the steps induced by Coulomb interaction at ω αU ≡ |ε + U − µ α |. The different Coulomb interaction modifies the positions and the heights of the steps in the spectrum. In the positive high frequency limit at ω > max{ω LU , ω RU }, the current fluctuation spectra nearly approach the same value due to the absorption of enough energy to open all the tunneling channels.
We identify the regimes as (i) weak U with µ L > ε, ε + U > µ R and (ii) strong U with ε + U > µ L > ε > µ R . In Fig. 4, for U = 0, the ratios of the heights of the steps both in the positive and negative parts are the same as spinless single level. However, compared to Fig. 3, the magnitude of the spectrum is doubled due to the twoenergy levels (spin-up | ↑ and spin-down | ↓ ) involved in the transport. After the similar derivation in spinless single-level dot as demonstrated above, the ratios of the heights of the steps for the positive (denoted by h) and negative (denoted by h ′ ) frequencies in the spectrum are given by, respectively, for (i) weak U (the short-dashed-line in Fig. 4), and for (ii) strong U (the dot-dashed-line in Fig. 4). Since the stationary double occupation is not allowed for strong U , there are no Coulomb-induced steps in the negative part of the spectrum.

B. Coupled double quantum dots
Now let us consider the electron transport through the system of two coupled quantum dots described by the Hamiltonian H S = ε l a † l a l + ε r a † r a r + Un lnr + Ω a † l a r + a † r a l , where U is the interdot Coulomb interaction, n µ = a † µ a µ andN = µ=l,rn µ . The involved states of the double dot are |0 for the empty double dot, |L for the left dot occupied, |R for the right dot occupied, and |2 ≡ |LR for the two dots occupied. Here, we assume at most one electron in each dot. In this space, we have a l = |0 L| + |R 2| and a r = |0 R| − |L 2|. The description of the involved states in this double dots is quite similar to that in the single dot with spin-dependence studied above. However, the essential difference is that the states of |L and |R are not the eigenstates of the system Hamiltonian H S which have the intrinsic coherent Rabi oscillation demonstrated by the coherent coupling strength Ω. The corresponding Rabi frequency denoted by ∆ is the energy difference between the eigenstates (ε ± ), e.g., ∆ = ε + − ε − = 2Ω for the degenerate doubledots system considered here.
The current-fluctuation spectrum for the coupled double dots is numerically displayed in Fig. 5. Similar to the single quantum dot, the spectrum has the feature of a energy-resonance step at ω α0 = |ε ± − µ α | and the Coulomb interaction induced step at ω αU = |ε ± +U −µ α | as shown in Fig. 5 (a). Here, we only illustrate the positive-frequency part of the spectrum due to the similar or less information involved in the negative-frequency regime. It is worth noting that the step behavior reflects the eigenstate energy structure of the dots, say, ε ± rather than ε l,r . Besides the same step behavior as in the single quantum dot, of particular interest is the signal of the coherent Rabi oscillation of the coupled double dots in the current-fluctuation spectrum.
Interestingly, the emergence of the coherent signal of the Rabi oscillation is nearly determined by the strength of the Coulomb interaction U . For weak Coulomb interaction in the regime of µ L > ε ± , ε ± + U > µ R , such as U = 0 (solid-line) and U = 2Γ (dashed line) in Fig 5  (a), no signal at the Rabi frequency ω = ∆ occurs in the spectrum. While in the double-dot Coulomb blockade regime, either ε + + U > µ L > ε ± , ε − + U > µ R or ε ± + U > µ L > ε ± > µ R , the current-fluctuation spectrum always shows a peak at the Rabi frequency ω = ∆ as shown in Fig 5 (a) with U = 6Γ (short-dasheddot-line) and U = 15Γ (dashed-dot-line). Although the coherent signal peak appears for asymmetrical coupling with Γ L = 4Γ R shown in the inset of Fig 5 (a), it is quite weak compared to that induced by Coulomb interaction. This means the coherent Rabi oscillation information in the current-fluctuation spectrum is sensitive to the dynamical blockade channel. This characteristic of Rabi coherence signal is also consistent with Markovian treatment studied in Ref. 42, where the symmetrical current fluctuation spectrum was considered.
Furthermore, by increasing the coherent coupling strength Ω, we find that the coherent signal of the Rabi oscillation is moved to the high-frequency regime with strong non-Markovian effect as shown in Fig 5 (b). Simultaneously, the peak of the coherent signal in the spectrum gradual increases monotonically and sharply increases at the resonance regime where the Rabi frequency approaches the bias voltage, i.e., ∆ = 2Ω = eV , as shown in Fig. 5 (c). This arises from the interplay between the Rabi resonance and the lead-dot tunneling resonance, i.e., ε ± (= ±Ω) = µ L,R (= ±eV /2), combined with strong non-Markovian effect. It may suggest that this resonant regime is good for the observation of the coherent signal in the current-fluctuation spectrum experimentally. Beyond the resonance regime, the system will enter into the cotunneling regime which is beyond the present approach and we have to recur to more advanced approaches such as the hierarchical equations of motion 25,31,43 and the real-time diagrammatic technique 38 for the consideration of higher-order contributions in the self-energy.

V. SUMMARY
In summary, using the frequency-dispersed technique by transforming the typical time-nonlocal master equation into equivalent time-local equations-of-motion, we established an efficient formula for the two-time non-Markovian correlation function of arbitrary system operators in open quantum systems. The key to the calculation of the non-Markovian correlation function is how to restor an effective time-translation symmetry to the propagator. We find that this corresponds to the vertex corrections as further demonstrated by the real-time diagrammatic technique. The final result has an elegant structure and is as convinent to apply as the widely used quantum regression theorem for the Markovian case.
We applied the present method to study the currentfluctuation spectra in the interacting single quantum dot and coupled double dots, respectively, contacted by two electrodes. The typical non-Markovian effect have been demonstrated. We found that the non-Markovian step behavior in the current-fluctuation spectrum of the single quantum dot is consistent with that in the noise spectrum of the transport current through the leads. The sharp peak of the coherent Rabi signal in the double dots occurs at the resonance regime where the eigenenergy levels are comparable to the chemical potential in the leads under the applied bias voltage. From this current-fluctuation spectrum covering the full-frequency regime, the information of the energy structure of the quantum dots, the tunneling rate as well as the Coulomb interaction and even the coherent Rabi signal can be extracted directly.
In this appendix, we give a detailed derivation of the non-Markovian two-time correlation function given by the Eq. (20). For the calculation of the steady state as the initial condition of correlation function, the stationary solution of auxiliary density operatorsφ ± (ω) ≡ φ ± (ω, 0) are given by Eq. (11) Updated with the initial condition of (c.f. Eq. (19)), we get the formula of the auxiliary density operatorŝ Inserting it into Eq. (12a) yieldṡ which has the solution Here, we used the initial conditionρ(0) =Bρ(0) =Bρ. With the use of Eq. (A1) and the relation (c.f. Eq. (9) and Eq. (10)) based on Eq. (17), we finally get the NMK-CF expressed by Eq. (20). We could directly get the spectrum of the correlation function in the frequency domain by a Laplace transformation based on Eq. (20). An alternative way is that we first obtain the stationary solution of auxiliary density operators in Eq. (A1) which reads Then, we transforme Eq. (A3) into the frequency-domain with the relation The result is Finally, according to the first identity of Eq. (17) which suggests L[ A(t)B(0) ] = Tr[Âρ(ω)], we get the formula of the NMK-CF in the frequency-domain as expressed in Eq. (22). Similarly, the correlation function of Â (0)B(t) is given by in the time-domain and in the frequency-domain. For simplicity so far, we showed the derivation of the NMK-CF for a single-operator coupling formalism as expressed in Eq. (3). Realistically, the coupling Hamiltonian is given by the multiple-operator coupling formalism, i.e., H ′ = Q + µ F − µ + H.c. as we illustrated in Sec. IV. The final formulas of Eq. (20) and Eq. (22) can be generalized by simply replacing the operators Q ± and C  dt 1Ã (t) Π(t, t 2 ) Σ B (t 2 , t 1 ) Π(t 1 , t 0 )ρ(t 0 ) .

The derivation for Eq. (29)
We introduce ̺ B (t) ≡ Π(t, 0)ρ B (0) with ̺ B (0) =Bρ, and the auxiliary density operator describing the vertex contribution, with ̺ ΣB (0) = 0 apparently. Since we are interested in the discussion of the order of the magnitude roughly, it is more convenient in the H S -interaction picture. The corresponding time-derivation equations arė withf B (t) = ∞ t dt 1 Σ B (t 1 )ρ. Then the time-derivation of the two-time correlation function (c.f. Eq. (C7)) reads The Taylor expansion of the time-derivative of the correlation function for small t = 0 + is G I (t) = G I (0) + dG I (t) dt t=0 t + 1 2 With Eq. (C9), we get where γ is the minimum decay rate γ of C(t) in Eq. (27), f (Q,B,ρ) and f (Q) are just the formal expression arising from Σ B and Σ, respectively. Inserting Eq. (C12) into Eq. (C11), we finally get Eq. (29).