Decoherence in open quantum systems: influence of the intrinsic bath dynamics

The non-Markovian master equation for open quantum systems is obtained by generalization of the standard Zwanzig-Nakajima (ZN) projection technique. To this end, a coupled chain of equations for the reduced density matrices of the bath $\varrho_{B}(t)$ and of the system $\varrho_{S}(t)$ are written. Formal solution of the equation for $\varrho_{B}(t)$ in the 2-nd approximation in interaction yields a specific extra term, related to the intrinsic bath dynamics. This term is nonlinear in the reduced density matrix $\varrho_{S}(t)$, and vanishes in the Markovian limit. To verify the consistence and robustness of our approach, we apply the generalized ZN projection scheme to a simple dephasing model. We study the obtained kinetic equation both in the Markovian approximation and beyond it (for the term related to the intrinsic bath dynamics) and compare the results with the exact ones.


Introduction
When studying dynamical processes such as relaxation, decoherence, buildup of correlations due to the interaction of an open quantum system with its environment, at a certain stage of investigation one inevitably faces the question: do these phenomena exhibit the Markovian behaviour or not? Many physical systems are believed to be described within Markovian approximation, since the coupling to the environment is weak (Born approximation), and/or the correlations in the bath decay fast on the characteristic time scale. However, there are situations, when the memory effects in the bath cannot be neglected, and the Markov assumption is not applicable anymore. This can be due to strong systemenvironment couplings [2], correlations and entanglement in the initial state [3], at the heat transport in nanostructures [4], or due to the specific character of the finite reservoirs [5]. The latter cases are of a particular interest, because the environment of an open quantum system due to its compactness frequently cannot be regarded as a thermal bath [6]. In such a case, the dynamics of the reservoir should be treated (at least at the initial stage of evolution) on equal footing with the dynamics of the -subsystem.
A powerful tool for dealing with such systems is provided by the projection operator techniques [7,8], which was developed by Nakajima [9], Zwanzig [10], and Mori [11]. This approach has manifested its efficiency in the derivation of generalized master equations and in investigation of the non-Markovian dynamics in initially correlated open quantum systems [12], in the spin star systems [13] and many others.
However, the ZN scheme has some disadvantages. Though the generic master equation is usually written down up to the 2nd order in the interaction, the time convolution in the kinetic kernels is determined by the full evolution operators (including the interaction part of the total Hamiltonian). Thus, to provide a consistent analysis of the system dynamics, one should take into account as many terms in the coupling constant as possible when expanding the master equation in the series in . This is a cornerstone of the time-convolutionless equation [14] and of the related projection technique, when one moves from the retarded dynamics to the equations local in time with the time-dependent generators and the subsequent systematic perturbative expansion scheme.
In the decades that passed since the pioneering works of Nakajima and Zwanzig, the projection methods technique has been significantly improved [15]. The main efforts are usually directed toward the construction of correlated projection superoperators [16,17], which consider the relevant part of the dynamics as a correlated system-environment state, rather than a tensor product state ( ) ⊗ , which allows one to reproduce the exact solutions already in the lowest order in the interaction [18]. In the very recent paper [19] the correlated projector technique has been generalized, yielding the so-called adapted projection operators (APO). The APO master equation encloses the full dependence on the initial correlations in a homogeneous term, which is essential to avoid the problems [14] dealing with an inhomogeneous time contribution arising from the unfactorizable initial states. An application of the APO methods manifested its efficiency in several cases, e. g., in investigating the polarization degrees of freedom of a photon going through a quartz plate -the system, which can be described by a similar dephasing model. Another kind of projections, namely, the projections to the environment subspaces corresponding to a given energy, was proposed in [20], which allowed the authors to derive an equation for the time-dependent reduced density matrix of the environment ( ) and to closely approach the problem of taking account of the intrinsic bath dynamics (BD).
However, none of the above methods considers the dynamics of the environment from the very beginning, on equal footing with that of an open quantum system . Usually, the evolution equation for the reduced density matrix of the environment ( ) is neglected too. Though for the infinite thermal baths this simplification seems to be reasonable, its validity could be questioned, when one deals with the finite size reservoirs [5] or tries to investigate the impact of the running correlations [21]. In other words, the open system dynamics should be definitely in accord with the concept [22,23], claiming that an evolution towards local equilibrium is always accompanied by an increase of the system-bath correlations.
In this paper, we generalize the standard ZN projection technique on our way to derive the non-Markovian master equation for open quantum systems. In section 2, we start from the coupled set of equations for the reduced density matrices of the system ( ) and of the environment ( ). Then, in section 3, we insert a formal solution of the dynamic equation for ( ) into the equation for the reduced density matrix of the -subsystem and obtain a generalized master equation in the second order in interaction. The above generalization gives rise to the extra term of a very specific structure: i) it is nonlinear in ( ) and ii) vanishes in the Markovian limit. These two points make the situation very similar to what one faces when studying the inset of the running correlations: the generic kinetic equations are strongly nonlinear, and the correlational contribution to the collision integral tends to zero in the Markovian limit [22,24]. We also briefly compare our results with those following from the standard ZN scheme.
To verify the consistency of our approach, in section 4 we apply the generalized projection scheme to a simple dephasing model [25] and obtain the kinetic equation for the system coherence. It is pointed out that the application of the Markovian approximation (MA) allows us to reproduce the exact results reported in [25] in the lowest orders in the interaction. On the other hand, going beyond the MA in the term dealt with the BD leads to the renormalization of the phase shift and the generalized decoherence. We demonstrate that in a certain domain of the state parameters, the above renormalization allows one to describe the correlational contribution to the generalized decoherence more precisely than the ZN scheme.
In the last section, we make conclusions and discuss the obtained results in their relation to the concept of the running correlations build-up.

Basic equations
Suppose that the composed system ( + ) consists of the open quantum system (subsystem ) and its surrounding , which usually can be considered as a thermal bath. The total Hamiltonian of this system consists of the term ( ) corresponding to the open quantum system (which in the general case is allowed to depend on time due to the action of the external fields), the summand related to the bath, and the interaction term .
The density matrix ( ) of the composed system obeys the quantum Liouville equation (hereafter we put ℏ = 1), For further convenience, let us pass to the interaction picture ( ) for the total density matrix, which is defined as where we have introduced the unitary evolution operators wherefrom it is easy to present the quantum Liouville equation (2.2) in the interaction picture, Hereafter, operators in the interaction picture are defined as If ( ) is an operator acting in the Hilbert space of the subsystem (bath ), then the corresponding interaction representation can be written down as follows: The reduced density matrices for the system and the environment are introduced in a usual way, by taking a trace over the environment (the system) variables. The corresponding time averages for operators ( ) can be introduced as follows: The reduced density matrix ( ) of the open quantum system is defined as (2.11)

13302-3
Analogously, in the interaction picture, the reduced density matrix ( ) of the environment reads It is obvious that ( ) = Tr ( ) and ( ) = Tr ( ). To derive the evolution equations for the operators ( ) and ( ), we apply the following decomposition for the total density matrix:  Using the decomposition (2.15) and taking into account equations (2.22)-(2.24), one can obtain the equation of motion for the correlation contribution to the total density matrix, The above equation can be rewritten in a more compact form, using the superoperators It can be shown that, when acting on operators with zero trace, Tr = 0, the superoperator P ( ) satisfies the relation P 2 ( ) = P ( ), i.e., it is a projector. A formal solution of equation (2.27) is where the superoperator U ( , ) satisfies the equation of motion and has an explicit form Substituting expression (2.29) into equations (2.23) and (2.24), we arrive at the closed system of equations for the reduced density matrices, However, equation (2.32) cannot be considered yet as a master equation for the open quantum system since it depends on the reduced density matrix of the environment. Solving equation (2.33) with respect to the density matrix of the environment ( ) seems unrealistic. Thus, the only way is to use some approximations for ( ) obtainable from equation (2.33). This is a subject of the next section.

( ): weak coupling approximation
Let us start this section with two assumptions: • If the subsystem is small compared to the environment , it is reasonable to suppose that, for sufficiently small time , the state of the environment is close to (0). Based on the above assumption, we write the solution of equation (2.33) as The above assumption means that we restrict the equation for ( ) by the first order in the interaction.

13302-5
• On the right-hand side of equation (2.32), we set U ( , ) = 1 and ( ) = (0), ensuring the integral terms to be precisely of the 2nd order in interaction.
Based on these approximations, we get where After some manipulations with the r.h.s. of equation (3.2), which are described in detail in the Appendix A of our recent paper [26], we obtain the final master equation in the interaction picture: where we have used the denotation for the kernel of the first integral in equation (3.4), while the operator L (0) ( ) is defined as Some comments on the structure of master equation (3.4) are quite pertinent at this stage: • The first term on the r.h.s. vanishes if (0) is the equilibrium density matrix of the bath and Tr (0) = 0.
• In the above case, the term L (0) ( )L (0) ( ) ( ) also vanishes. It is shown in section 4 that the terms with L (0) describe an influence of the initial correlations in the open quantum systems. In particular, the non-integral term contributes to the non-dissipative properties, whereas the integral term represents a new channel of dissipation related to the initial preparation of the system.
• The second term in the r.h.s. has some interesting properties: 1) It vanishes in the Markovian limit , when the dynamics of the kinetic kernel ( , ) is considered fast in comparison with the time evolution of the density matrix; then we set ( ) ≈ ( ).
2) This term is nonlinear in since the operator ( ) appearing in equation (3.5) depends on ( ) [see also equation (2.21)]. A nonlinearity appears due to consideration of the intrinsic dynamics of the environment by means of the time evolution of the reduced density matrix ( ). There is some resemblance with the basic points following from the quantum kinetic theory [22,24], where the running correlations yield nonlinear terms in the kinetic equations, vanishing in the Markovian limit.
One should distinguish two kinds of approximation for the integrand. The first one looks as This approximation is often called the MA of the 1st kind [15] or the Redfield approximation (dealt with the corresponding Redfield equation [14] for the reduced density matrix of the subsystem ). The other approximation, is known as the MA of the 2nd kind [15], very often [22] referred to just as the Markovian approximation. In this paper, we use the MA of the 1st kind only.
It was mentioned in the end of the previous section that the system of equations (2.32)-(2.33) for the reduced density matrices is unmanageable, and obtaining of the master equation for ( ) in higher orders in interaction is possible only via some iterative scheme. The eventual master equation would be too cumbersome to be manipulated with. Nevertheless, some possible conclusions can be anticipated by remembering the results which follow from the quantum kinetic theory. For instance, a contribution of the dynamical correlations is known to remain in the Markovian limit [22] if high order approximation in the coupling strength is performed. The origin of the above dynamical correlations [27] can be mapped to our concept of the intrinsic bath dynamics. Thus, we could expect the presence of the BD term even in the Markovian limit, as opposed to the above mentioned zeroth contribution following in MA from equation (3.4).
On the other hand, as it is shown in section 4, the obtained non-Markovian master equation (3.4), though it is formally of the second order in L ( ), contains explicitly all the orders in the interaction, if one is going to pass to its convolutionless counterpart. Thus, the latter approach seems to be much easier but not less effective, and generally there is no special need to construct the non-Markovian master equations of higher order in the coupling strength.
Now, it would be instructive to compare the obtained results with those that follow from the traditional ZN projection scheme. In this scheme, one usually uses a different decomposition for the total density matrix instead of (2.15), It is important to note that [c.f. equations (2.16)] Tr In fact, decomposition (3.7) means that one neglects the intrinsic dynamics of the environment by assuming the reduced density matrix of the bath to be equal to its initial value at any time, ( ) ≡ (0) . Note that the initial state (0) of the environment should not necessarily be equilibrium: it may contain information about the initial system-bath correlations arising due to the selective [25] quantum measurements, as it occurs in the case under study (see section 4) or non-selective measurements [14,31,32].
It is pointed out in [28] that the ZN projection technique can be viewed as a particular version of the non-equilibrium statistical operator method (NSO), when the interaction term is not taken into account in the relevant distribution rel ( ). The relevant statistical operator is usually considered as the distribution, relating the initial state of the system to the NSO ( ) through the decomposition [c.f. equations (2.15) and (3.7)]. Since in the limit → ∞ the system attains an equilibrium state [22], the correlational part Δ ( ) of the total density matrix approaches zero only if the interaction is included into the relevant distribution. Otherwise, like in the case (3.7) of the ZN projection method and its generalization (2.15), lim →∞ ( ) ≠ 0 and lim →∞ Δ ( ) ≠ 0, since the equilibrium distribution of the composite system is not factorisable in the operators, belonging to the different subsystems and .
Here, we leave aside all the peculiarities dealing with the system tending to equilibrium, and refer our readers to the papers [27,28]. Secondly, the expression (3.8) renders the main disadvantage of the ZN projection method quite obvious: the state of the environment is postulated to be described by an equilibrium distribution plus some correction Tr ( ), which seems reasonable only for the case of small interactions. On the contrary, the generalization of the ZN method implies that the reduced density matrix of the environment has its own dynamics, completely incorporated in the distribution ( ), see equations (2.13) and (2.17) in [29], whereas approximation (3.1) can be treated as its linearized version.
It should be noted that in the original ZN scheme [14] it is supposed that [ , In the weak coupling limit (see Appendix B in [26] for details), the ZN master equation looks as follows: . As a consequence, in contrast to equation (3.4), the master equation (3.10) does not contain a nonlinear term with ( , ) due to the intrinsic dynamics of the environment being excluded from consideration. Nevertheless, in the Markovian limit, both master equations are identical. This is quite obvious, since the effect of BD [considered in terms of ( )] is expected to manifest itself only in the initial stage of the system evolution.
To check out the consistency and robustness of the generalized projection scheme, one should use the master equation (3.4) to derive the quantum kinetic equations for the observables in some simple models (preferably, exactly solvable ones). This is done in the next section, where we obtain a dynamic equation for coherence, considering the purely dephasing model [25].

Non-Markovian quantum kinetic equation for a dephasing model
We consider a simple version of a spin-boson model describing the two-state system ( ) coupled to the bath ( ) of harmonic oscillators [25,30]. In the spin representation for a qubit, the total Hamiltonian of the model is written as where 0 is the energy difference between the excited state |1 and the ground state |0 of the qubit, and 3 is the 3rd Pauli matrix , = {1, 2, 3}. Bosonic creation and annihilation operators † and correspond to the -th bath mode with the frequency . Since 3 commutes with Hamiltonian (4.1), it does not evolve in time, 3 ( ) = 3 . Hence, the interaction operator can be easily expressed as Let us introduce in a usual fashion the spin inversion operators ± = ( 1 ± 2 )/2, which obey the permutation relations [ 3 , ± ] = ±2 3 . Our task is to obtain the quantum kinetic equations for the mean values ± = Tr { ( ) ± } dealing with the system coherence. For simplicity, we consider the case with + , since the equation for − can be easily obtained in a similar manner.
Using the basic equation (3.4) for the reduced density matrix ( ) and taking a trace over the system variables, one can derive the following kinetic equation: Having calculated all commutators in equation (4.3) and having performed thermal averaging (see Appendixes C and D in [26] for details), we can write down the final kinetic equation for the generalized coherence: where the quantity is related to the initial correlations (in other words, it deals with the initial preparation of the system due to some kind of the quantum measurement [25,31]). In the kinetic equation (4.4), the qubit-bath coupling is modelled by the spectral weight function ( ) taken in its standard form [25] ( ) = Ω 1− exp(− /Ω). The parameter ∼ | | 2 in equation (4.6) denotes a dimensionless coupling constant, while Ω is the cut-off frequency. Depending on the exponent , we can distinguish three coupling regimes: the sub-Ohmic at 0 < < 1, the Ohmic at = 1, and the super-Ohmic at > 1.
Let us discuss some peculiarities of the quantum kinetic equation (4.4). First of all, we recall that the second term in the r.h.s. of (4.4) vanishes in the Markovian limit. On the other hand, we can always pass from the mean values in the interaction picture to the averages taken with the statistical operator ( ) according to the rule + → exp(−i 0 ) + . Then, the second term in the r.h.s. of (4.4), being an imaginary one, along with the first (non-dissipative or "quasi-free") term should add up to 0 , yielding the corresponding phase shift [25].
Secondly, the initial state of the qubit contributes to the quasi-free term as well as to the last summand in the r.h.s. of equation (4.4) via init . Besides, while the kinetic kernels in the 2nd and the 3rd terms have a usual form K ( − ), this is not the case for the last term dealing with the initial correlations in the "qubit+bath" system. This is not surprising: the dynamics of initial correlations is not a stationary process with a typical convolution dependence; the initial correlations decay in time due to the aging effect [33]. Thirdly, the contribution of the initial correlations is of the 4th order in interaction.
In our recent paper [26] we solved the generic non-Markovian equation (4.4) numerically. We considered the Ohmic and super-Ohmic (with = 3/2) qubit-bath coupling and studied the qubit dynamics at both low and high temperatures. The results were not quite satisfactory from the physical point of view: the BD being taken into account led to an uncontrollable increase of Im + ( ) at low temperatures. The picture was not improved much even at high temperatures, since the imaginary part of the system coherence tended to saturation. This contradicts the conclusions made in [25] that the partial decoherence is admitted only at the strong super-Ohmic regime with > 2, while at smaller values of the ohmicity index the system coherence always tends to zero at large times. Since no systematic analysis can be proposed in the non-Markovian regime but only the possibility to solve the kinetic equations numerically, and since  the reasons leading to the above mentioned unphysical behaviour of the coherence were formulated only qualitatively, we perform the MA of equation (4.4) to obtain exact analytical expressions and to study the system dynamics more carefully.
The first term in the the r.h.s. of equation (4.4) is proportional to the time derivative of the phase shift, which is nothing else but the linearized over Φ( ) version of the expression (32) from [25]. Furthermore, taking the generalized coherence out of the integral in the third term and integrating the kinetic kernel over , we come to the conclusion that it deals with the time derivative of where vac ( ) ( th ( )) denotes the vacuum (thermal) contribution to the generalized decoherence.
In a similar way, performing the MA, we verify that the last term of equation (4.4) is related to the time derivative of the correlational contribution to the generalized decoherence, (4.9) Again, the expression (4.9) is nothing else but the series expansion of equation (31) from [25] up to the fourth order in interaction. Now, let us examine more thoroughly the role of the intrinsic BD, which is described by the second term in the r.h.s. of the kinetic equation (4.4). This term has a very specific structure: it does not depend on the temperature in spite of having been generated by the qubit environment. In some sense, it is similar to the vacuum contribution vac ( ) to the decoherence function. It is obvious that the above similarity is rather coincidental, being a cumulative effect of the applied approximation (3.1) for˜( ) and a specific feature (non-ergodicity) of the dephasing model. To get some more information about the role of BD, let us go beyond the MA and expand the difference of the generalized coherences up to the first order in = − : Since the time derivatives of both the phase shift and the generalized coherence are of the 2nd order in the coupling constant, the second term in the r.h.s. of (4.4) is proportional to | | 4 . Thus, taking (4.7)-(4.10) into account, it is straightforward to solve equation (4.4) and to obtain the generalized coherence in its usual form Here, the renormalized values of the phase shift¯( ) and the generalized decoherence¯( ) up to the 4th order in interaction read as follows: where the function can be calculated explicitly, but is too cumbersome to be presented here. In the integrand of (4.12) we omitted cor ( ), since it is already of the 4th order in the coupling constant. At the chosen order in the coupling constant, the BD contribution to the decoherence function vanishes for the initially uncorrelated subsystems [see the equation (4.13)], since the phase shift is zero in such a case [25]. It can be said that the initial correlations and the intrinsic BD are working together in the composite quantum system to enhance or decrease the coherence, depending on the sign of 3 . Since in the dephasing model the exact value of the correlational contribution to the decoherence function is known [25], 15) it is reasonable to treat the last term in (4.13) as a renormalization of (4.9) and to comparē with the above exact result (4.15).
As it is seen in figures 1 and 2, there is a domain of the state parameters (the temperature 1/ and the inversion of the mean levels populations 3 ), where the renormalized values (4.16), obtained within our generalized projection technique, reproduce the exact result (4.15) better than those derived using the ZN approach (4.9). This region is found to be characterized by the values −0.2 < init < 0, meaning a slight qubit frequency downshift. It should be noted that such a behaviour of cor ( ) is typical both in the high (figure 1) and low (figure 2) temperature limits. It could be also shown that in the moderate super-Ohmic regime, the obtained dynamics shows tendencies similar to those presented at = 1. However, at the sub-Ohmic coupling, the behaviour of exact cor ( ) is known to be non-monotonous in time [25], and none of the approximations (4.9) or (4.16) can reproduce the true dynamics of the correlational contribution to the decoherence function at large times due to a specific kind of coupling [15].
On the other hand, the lowest series expansion in (4.10) does not allow us to reliably estimate the correction to the phase shift according to equation (4.12), it arises even at ( ) = 0, contradicting the well-known statement [25], which attributes the qubit frequency renormalization solely to the appearance of the initial correlations in the composite + system. Going beyond the 4th order in the interaction in calculations of the renormalized values of ( ) and cor ( ) could eventually improve the situation. However, it also implies taking into account the next terms in the series expansion (4.10) which automatically would deviate the expected result from the exact exponential form (4.11). Nevertheless, the results presented in figures 1 and 2 show that the proposed generalization of the projection technique can be considered as a good alternative to other perturbative methods of solving the quantum master equations [14], which perceptibly improves the ZN results already in the lowest order in the coupling strength.
To conclude this section, we mention that all the results are obtained for the system that was initially prepared by the selective quantum measurement. The quantity (4.5), entering all the expressions dealing with initial correlations, is uniquely defined by the above preparation measurement. However, there is a much more interesting situation, when the system is prepared due to a non-selective measurement scheme [14]. In such a case, not only the decoherence induced by the initial correlations is observed but also the inverse process (recoherence, or the coherence enhancement) is possible. In [31], we proposed some kinds of measurements leading to the permanent coherence growth at small times ∼ 1/Ω. Obviously, at times 1/Ω, the other dissipation channels (thermal and vacuum ones) prevail, yielding a decrease of the system coherence. There is also a possibility to prepare the system in such an initial stage, when the periods of recoherence and decoherence alternate. In the very recent papers [34,35], the evolution of quantum coherence for the same dephasing model is studied in detail over short and long times.
Thus, we are in a position to perform a similar analysis of the BD contribution to the system, which is initially prepared after the non-selective quantum measurements, since the general ideas leading to the "renormalized" functions (4.12)-(4.13) are valid in this case as well, and the general expression (4.11) for the system coherence remains unaltered. However, the expression for the "renormalized" correlational contribution¯c or ( ) to the decoherence/recoherence rate would consist of too many parameters to be analyzed reliably. As far as our case of selective preparation measurements is concerned, the short or long time analysis of the renormalized function (4.16) can hardly be informative or useful, because of its small contribution to the (permanent!) system decoherence. Therefore, we prefer to focus our efforts on the investigation of a measure of its deviation from the exact result (4.15) in order to verify the validity and robustness of our generalized projection scheme beyond the Zwanzig-Nakajima framework.

Conclusions and outlook
In this paper, we generalize the ordinary ZN projection scheme in order to derive a master equation for the open quantum system weakly coupled with its surrounding. We start from the chain of equations for the density matrices of the -and -subsystems and obtain a non-Markovian master equation for ( ).
We simplify this equation, restricting ourselves by the lowest approximation in the coupling constant. The consideration of the intrinsic BD yields an extra term in the master equation, which is found to be nonlinear in ( ) and vanish in the Markovian limit. When the dynamic equation for ( ) is neglected, we come to the standard ZN projecting scheme [14].
In order to test the developed method we applied it to the open quantum system described by a very simple dephasing model. The non-Markovian quantum kinetic equation for the generalized coherence has been derived. In the MA, which means nothing else but the time convolutionless ZN result, its solution up to the 2nd order in interaction coincides with the exact one obtained in [25]. We also analyzed the kinetic kernel related to the BD, which is found to be temperature independent. When going beyond the MA in the above mentioned term, we obtain the renormalized values of the phase shift and the generalized decoherence. We compared the renormalized correlational contribution to the decoherence function with the exact result [25] and with that of the ZN approach. It is found that consideration of the intrinsic BD brings the corresponding curves closer to the exact data compared to the results obtained within the ordinary ZN projection scheme. However, some model peculiarities, e.g., its non-ergodicity do not allow us to take the full advantages of our method: in the lowest order in interaction, the kernel dealing with the BD, turns out to be temperature independent, which seems a bit strange if one speaks about the effect of the qubit surrounding. Other difficulties arise from the necessity of taking the higher order terms into consideration which inevitably causes significant computational complications and, from the mathematical point of view, the solution for the generalized coherence can even deviate from its exact exponential form [25].
There is also a different approach to treat the influence of BD, which is based on the consideration of dynamical correlations in a composite system [21]. These long-living dynamical correlations, which are associated with the total energy conservation, play an important role in the transition to the Markovian regime and the subsequent equilibration of the system. In our recent paper [27], we applied this approach to the dephasing model and obtained nonlinear kinetic equations, involving the so-called correlational temperature. Though the above method does not admit consideration of the intrinsic BD explicitly, it allows one to study the build-up of the dynamic correlations between the qubit and the bath, as well as its time evolution from the initial stage up to the system equilibration. It is interesting that the dynamic correlations act in a similar way, renormalizing the qubit frequency. Thus, we dare to say that the studies of the BD influence on the state of open quantum systems, the build-up of the system-bath correlations, as well as a choice of the most suitable models for verification of these approaches, are still far from completeness and should be continued in the future.

Note added in proof
My attention (V. V. I.) has been drawn to a very recent paper [36] where the problems of intrinsic bath dynamics were studied in the framework of the so-called correlation picture, which connects a correlated bipartite state to its uncorrelated counterpart. Some of the results and conclusions are very close to ours, obtained in the current paper and previously in [26], namely: i) generalization of the projecting operators; ii) the chain of equations for the reduced density matrices of the system and the bath. However, there are also some distinctive points: i) the models utilized in [36] in order to verify the robustness of their approach differ from ours; ii) we have separated the initial correlations and the BD contributions, yielding different forms of the corresponding time-convolutionless master equations. Similar results, obtained independently by different approaches, undoubtedly confirm the importance of dynamical correlations in the open quantum systems and the necessity to take them into account.
In this context it should be noted that a formal application of the perturbation theory in the interaction can sometimes be misleading: in the case of the dephasing model under study, any eventual higher-order corrections to the vacuum vac ( ) and thermal th ( ) terms should be rejected as contradicting the exact solution.