Hot cosmic qubits: late-time de Sitter evolution and critical slowing down

Temporal evolution of a comoving qubit coupled to a scalar field in de Sitter space is studied with an emphasis on reliable extraction of late-time behaviour. The phenomenon of critical slowing down is observed if the effective mass is chosen to be sufficiently close to zero, which narrows the window of parameter space in which the Markovian approximation is valid. The dynamics of the system in this case are solved in a more general setting by accounting for non-Markovian effects in the evolution of the qubit state. Self-interactions for the scalar field are also incorporated, and reveal a breakdown of late-time perturbative predictions due to the presence of secular growth.


Introduction
Many of the foundational paradoxes of gravitating quantum systems -e.g. black-hole information loss, eternal-inflation and multiverse issues -arise due to puzzling behaviour displayed by simple systems at very late times. The simple systems used are often free quantum fields evolving in gravitational backgrounds, often spacetimes with horizons, that are chosen because explicit calculations can be made. When using these systems to make latetime inferences an implicit assumption is that it is the interaction with the background that always dominates, and any other neglected interactions can be treated as perturbations.

JHEP02(2020)053
Similarities between the physics of quantum systems in gravitational spacetimes (especially with horizons) and open systems [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] show that this assumption is actually unlikely to be true. The problem is that open systems (by definition) always have an 'environment' whose properties are not measured (in this case, perhaps, the degrees of freedom behind the horizon). But a perturbative treatment of the interactions with such an environment essentially always fails at sufficiently late times. Ultimately it fails because the environment never goes away. Given enough time the effects of any interaction -regardless of how weak it might be -eventually accumulate to become large. Concretely, no matter how small an interaction Hamiltonian, H int , might be, there is eventually a time t for which the evolution operator, U (t) := exp[−i(H 0 + H int) t], is not well-described by a finite number of powers of H int .
In practice this problem usually manifests itself through the phenomenon of 'secular growth', where the coefficients in a perturbative expansion contain growing powers of time. For example if an observable O(t) is computed in powers of a small coupling g 1 then the coefficients c n (t) typically grow without bound as t gets large. Part of the motivation for using open-system tools [4][5][6][7][8][9][10][11] is that they provide systematic ways to resum this late-time growth, often allowing perturbative results at short times to be converted into expressions that for large t include all orders in g 2 t but drop contributions of order g n t for n > 2. They thereby promise controlled and reliable approximations for late-time behaviour that straight-up perturbative methods cannot. This work accompanies a companion paper [21], which uses these tools to track the late-time evolution of an Unruh-DeWitt detector [22,23]: a simple two-level system (or qubit) as it uniformly accelerates in flat space while coupled to a simple quantum scalar field (prepared in its Minkowski vacuum). Such a simple system allows these open-system tools to be explored in a very concrete and explicit way (see also [24][25][26][27][28][29][30][31][32][33][34][35]). Ref. [21] treats the field as an environment and integrates it out to set up the Nakajima-Zwanzig equation [36,37] describing its perturbative effects for the qubit. This is an integro-differential equation that is difficult to solve, but which simplifies at late times under certain assumptions to give an approximate late-time Markovian evolution. In particular much attention is given to the precise parameter range that controls this approximation. Not surprisingly the classic transition rates computed for Unruh-DeWitt detectors decades ago [22,23,38] prove to break down at very late times. Evolution at much later times instead describes thermalization and decoherence as the qubit gets heated to the Unruh temperature.
We here apply open-system tools to a similarly simple late-time question: what happens to such a qubit (again coupled to a scalar field) moving for very long times along a co-moving trajectory in de Sitter space. We again identify the relevant master equation for differential qubit evolution once the field is integrated out, and again find that a Markovian approximation works at sufficiently late times, asymptotically approaching a thermal state in much the same manner as in [21]. de Sitter space brings an important complication, however: the length of time required for this Markovian limit to apply grows like an inverse JHEP02(2020)053 power of the scalar mass if this mass is sufficiently small. That is, the qubit responds to field self-correlations that become increasingly persistent and extraordinarily long-lived. As a direct consequence of this, the qubit's approach to equilibrium becomes extremely slow in this regime (a phenomenon reminiscent of critical slowing down [39,40]).
Another new feature of the de Sitter example is the emergence of a non-Markovian regime for which an approximate form of the late-time evolution can nevertheless be explicitly integrated to give a closed-form solution. This allows us to track a portion of the memory effects that work to foil the Markovian limit, and find a more general solution that settles to the expected late-time thermal state. We again develop precise conditions for when this solution is a good approximation, and recover the earlier Markovian solution as a limit.

Co-moving qubits and fields in de Sitter space
This section reviews for later use some basic properties of de Sitter space, with details of the qubit/field system to be studied. The section then closes with a statement -following [21] -of the Nakajima-Zwanzig equation that governs qubit evolution once the scalar field is integrated out.
The time evolution to be followed in later sections is along the coordinate direction within the flat slicing of de Sitter space, with line-element where H > 0 is the Hubble constant and the conformal time η is related to comoving time τ by η = −H −1 e −Hτ [41], with the negative sign chosen to ensure dη/dτ > 0. The range −∞ < η < 0 corresponds to −∞ < τ < ∞. This metric is maximally symmetric and so has constant Ricci curvature 1

Scalar fields in de Sitter space
We consider a real scalar field with Lagrangian density which includes a nonminimal interaction with the metric's Ricci scalar, R, with coupling parameter ξ. Because the Ricci scalar for de Sitter space is constant, from the point of view of the scalar field the nonminimal coupling effectively shifts the scalar field's mass from m to

JHEP02(2020)053
The canonical Hamiltonian that generates the scalar field's evolution 2 in τ is therefore where Σ τ is a sheet of fixed comoving time τ . In what follows we take the scalar field to be prepared in the Bunch-Davies vacuum |BD . Of particular interest for later evaluation of qubit evolution is the scalar field's Wightman function evaluated in this vacuum, which turns out to be given by [43][44][45][46] (2.6) where 2 F 1 is the hypergeometric function and → 0 + is an infinitesimal whose presence is determined by the Wightman boundary conditions and which determines how integrations should navigate around the singularity in the coincidence limit (where the points (η, x) and (η , x ) are lightlike separated).
The parameter ν in (2.6) is defined as the following function of ξ, m and H: In our conventions the special case of a conformally coupled scalar field is the choice m = 0 and ξ = − 1 6 , in which case M 2 eff = 2H 2 and ν = 1 2 . A minimally-coupled massless (axionlike) field by contrast satisfies m = ξ = M eff = 0 and so ν = 3 2 . The Wightman function in the conformally coupled case is particularly simple, reducing to (2.8)

Qubit/field couplings
To this field we couple a qubit, following the construction used in [21]. The result is an Unruh-DeWitt detector coupled to a scalar field, along the lines of that first introduced in [22,23]. We take the qubit's free Hamiltonian to be where ω 0 > 0 denotes the splitting between the two qubit energies. We suppose the qubit sits indefinitely at a fixed position y ∈ R 3 , and so follows the (geodesic) trajectory of a co-moving observer This is a special instance of the Klein-Gordon Hamiltonian for a static spacetime which admits a foliation with metric ds 2 = −dt 2 + γij dx i dx j , with Σt a sheet of fixed x 0 = t.

JHEP02(2020)053
along which the coordinate τ is the proper time as measured with the spacetime metric of (2.1). The Hilbert space of states for the combined qubit/field system is the product of the Fock space for the field with the qubit's 2 × 2 space of states. The free Hamiltonian acting on this Hilbert space is then where I and I are identity operators. 3 The total Hamiltonian is H = H 0 + H int 0 where the qubit/field coupling is described by the interaction Hamiltonian 12) and the dimensionless coupling 0 < g 1 is small enough to justify a perturbative treatment. We follow common convention and choose m = σ 1 , but all that really counts is that m and h do not commute with one another so that H int 0 drives transitions between the zeroeth-order qubit energy eigenstates.
Naively one would perform perturbation theory simply by expanding in powers of H int 0 . It happens, however, that at O(g 2 ) the qubit/field interaction shifts the qubit energy splitting so that ω := E ↑ −E ↓ = ω 0 +g 2 ω 1 for calculable ω 1 . A better choice for perturbation theory uses the physical value ω in the unperturbed Hamiltonian, and so writes With this choice H = H 0 +H int 0 = H free +H int and so the interaction Hamiltonian acquires a new counter-term, with (2.14) Although not motivated by divergences, this counter-term interaction has an added benefit inasmuch as the quantity ω 1 happens to cancel an ultraviolet divergence that arises at second order in g.

Time evolution and the Nakajima-Zwanzig equation
In principle the question of calculating the system's time evolution perturbatively in powers of the coupling g is a solved problem. One converts to the interaction picture, by performing a unitary transformation, As applied to the system's state (described by its density matrix, ρ) this transformation removes the 'free' part of the evolution, leaving ρ to be evolved by the interaction-picture Liouville equation Here φ I (τ, x) := e +iHτ φ(x)e −iHτ is the interaction-picture field and the interaction-picture qubit coupling matrix is m I (τ ) := e +ihτ m e −ihτ . (2.18) Standard arguments then solve eq. (2.16) interatively to any required power of V (τ ), subject to an initial condition, ρ(0), at τ = 0. In what follows we take the field to be initially prepared in the Bunch-Davies vacuum |BD , and assume the initial qubit state to be uncorrelated with the field degrees of freedom: where 0 is the qubit's initial 2 × 2 hermitian density matrix, satisfying tr 0 = 1. But a complete solution for ρ(τ ) is overkill if the goal is simply to predict the behaviour of qubit observables, with no measurements made in the scalar-field sector. For such observables the time-evolution problem is completely solved if the time-dependence of the reduced density matrix, is known, with initial condition In (2.20) the partial trace is only over the field-theory sector (and not also the qubit). The Nakajima-Zwanzig equation provides a formal solution to the problem of identifying the evolution equation for (τ ) that should replace (2.16). It is obtained by perturbatively solving for the evolution of the scalar-field state and using this to eliminate the scalar completely 4 from (2.16), leaving an evolution equation that involves only .

JHEP02(2020)053
It is the symmetry of the spacetime and the choice of the trajectory (2.10) that ensures that the Wightman function depends only on the difference τ 2 − τ 1 . In components, after some matrix algebra (2.23) yields the equations of motion where the integration variable are also switched, s → τ − s. The properties tr (τ ) = 1 and † (τ ) = (τ ) have also been used to eliminate 22 (τ ) = 1 − 11 (τ ) and 21 (τ ) = * 12 (τ ). Eqs. (2.25) and (2.26) make it clear that the diagonal and off-diagonal components of evolve independent of one other.

The Wightman function along the qubit trajectory
The Nakajima-Zwanzig equation boils all effects of the scalar field to the correlations encoded in the Wightman function W BD (τ ) evaluated along the qubit trajectory (2.10). Specializing (2.6) to this trajectory gives the explicit form where we remind the reader that → 0 + is taken at the end of any calculation. As is easy to verify this function satisfies the identities and This last condition -the Kubo-Martin-Schwinger (KMS) [55,56] condition -expresses detailed balance and ultimately ensures that (τ ) eventually asymptotes to a thermal state with Gibbons-Hawking temperature 5 [57] In the special case of a conformal scalar case -i.e. where M 2 eff = 2H 2 -the Wightman function simplifies to Table 1. Leading asymptotic behaviour of Re[W BD (τ )] in the limit Hτ 1 for several possible values of M eff /H (details given in appendix A). In the first row ν = iµ lies on the positive imaginary axis with µ ≡ M 2 eff /H 2 − 9/4 > 0 ensured because M eff /H > 3 2 . The leading coefficient in the asymptotic expansion for ν ∈ [0, 1) ∪ (1, 3 2 ) vanishes when ν → 1 ± which is why there is a separate row for the special case of ν = 1. which (after replacing H with the acceleration parameter) agrees with the Wightman function for a massless field in Minkowski space evaluated along a uniformly accelerated trajectory [21,58,59].
As discussed in [21], for the present purposes (late-time evolution) what is important about expressions (2.27) and (2.31) is their asymptotic form in the limit Hτ 1. As is shown (for various choices of parameters) in table 1 the function Re[W BD (τ )] always falls off exponentially, W BD ∝ e −τ /τc , for large enough τ . This falloff is important because it makes possible the existence of a simpler Markovian limit, at least for times, τ τ c . We see from this table that -with one exception -the time-scale τ c over which Re[W BD (τ )] is exponentially suppressed is generically of order the Hubble time. The exception is close to the massless, minimally coupled case, for which M eff → 0 and ν → 3 2 . Writing the asymptotic form as provided M 2 eff = 5 2 H 2 the parameters W 0 and κ are given by (see appendix A for details) For the limit M eff → 0 both W 0 and κ −1 blow up, with asymptotic forms (2.34) Evidently in the regime 0 < M 2 eff 2H 2 the function Re[W BD (τ )] falls off much more slowly, with the longer-range correlations noted in [2,30,60].
For later sections it is also useful to know the subdominant terms in the expansions about the Hτ → ∞ and M eff /H → 0 limits. The next-to-leading terms in the M eff /H JHEP02(2020)053 expansion are For Hτ 1, if the sub-leading corrections to the late-time asymptotic series for W BD (τ ) are defined by then f (τ ) falls much faster with τ than does W BD (τ ) in the limit of vanishing effective mass. Sub-leading corrections to f (τ ) in (2.37) are O(e −Hτ ) even in the M eff /H 1 limit. 6 These asymptotic properties are also visible in the numerical plots given in figure 1.

Late-time Markovian behaviour
This section exploits the exponential falloff of the Wightman function, W BD (τ ) W 0 e −κτ , to derive an approximate form for eqs. (2.25) and (2.26) that captures well the late-time evolution.

Markovian approximation
The idea behind the approximation is simple: to the extent that one's interest is in slow evolution at very late times (i.e. τ τ c ) then (τ ) does not vary significantly over the com-

JHEP02(2020)053
paratively short time-scales over which the integrals in (2.25) and (2.26) have appreciable support. It should be a good approximation in this regime to Taylor expand within the integrand [61,62], in which case the dependence on (τ ) comes out of the integral. Keeping only the leading term in this expansion allows (2.25) and (2.26) to be rewritten in the form As quantified below, the error in dropping subdominant terms in the integrand should be of order (τ c ∂ τ ) n (τ ), and therefore should be small when only varies slowly over times of order τ c . Specialization to times τ τ c allows further simplification because in this regime the limits of integration in (3.2) and (3.3) can also be taken to infinity up to exponentially small corrections. Switching back to the Schrödinger picture, this allows the above equations to be written as Here the coefficient functions are defined by the integrals

and [63, 64]
The function ∆ BD diverges at the s → 0 end of the integration, and its form is evaluated more explicitly in appendix B. We choose to regulate this divergence by making the JHEP02(2020)053 small-distance regulator in the Wightman function (2.27) small but finite, leading to a divergence that is logarithmic in whose form is derived more explicitly in appendix B. Notice also that the disappearance of the quantity ω 1 when passing from (3.3) to (3.5) occurs because we choose ω 1 = −∆ BD to cancel a term i∆ BD 12 in (3.5) in order to ensure that ω appears in the same way in the evolution of as should the physical qubit energy splitting.
The identities (2.28) and (2.29) impose some relations amongst the above integrals. First, (2.28) implies that Re[W BD (τ )] and Im[W BD (τ )] are even and odd functions in τ , respectively, and so which is the detailed-balance relation [58] (and so underlies the thermal nature of many of the late-time equilibrium properties). Because C BD and S BD are even and odd in ω respectively, the detailed-balance relation also implies from which we see

Late-time evolution
The quantity R BD appears throughout the literature on Unruh-DeWitt detectors [23,38,58,59] because it governs the perturbative excitation rate of a qubit that is initially prepared in its ground state: 0 = vac , where With this choice 11 (0) = 12 (0) = 0 and so both 11 (τ ) and 12 (τ ) are at most of order g at later times. Consequently any appearance of these quantities on the right-hand-side of eqs. (3.4) and (3.5) can be dropped if we wish to compute ∂ τ (τ ) only to order g 2 . This leads to the standard lowest-order prediction for the rate of accumulation of probability into the excited qubit state: As is clear from the above derivation, this rate strictly only applies at times τ τ c and in this sense is usually interpreted as the qubit's late-time transition rate. However it JHEP02(2020)053 is also clear that the prediction (3.15) cannot apply at asymptotically late times. Eq. (3.15) must break down at sufficiently late times because unitarity requires 0 ≤ 11 (τ ) ≤ 1 for all τ and so forbids eternally accumulating probability into the state |↑ with constant rate. 7 What happens at these later times? This is where eqs. (3.4) and (3.5) prove their worth, because their domain of validity can extend to much later times than their derivation naively would indicate [5,21]. Their extended domain of validity arises because -unlike for (say) eqs. (2.25) and (2.26) or eqs. (3.2) and (3.3) -neither of (3.4) or (3.5) make explicit reference to the initial time at which the evolution starts. This means that although they were derived starting from an assumed initial state at τ = 0, they could have equally well been derived starting at some other initial time, τ 1 or τ 2 . Although the derivation that starts from τ = τ n is only valid over a limited window of time, τ ∈ S n , to the future of τ n , since it is the same differential equation that is derived for all n, the differential equations themselves -i.e. eqs. (3.4) and (3.5) -remain valid for the union of all of these intervals: This argues that at much later times it is eqs. (3.4) and (3.5) that govern the evolution of (τ ), and it is the presence of the other terms on their right-hand-side beyond the g 2 R BD term that build in the constraints of unitarity that are missing in (3.15). Furthermore, keeping only terms to O(g 2 ) in (3.4) and (3.5) implies that their solutions can be trusted to all orders in g 2 τ as g → 0 and τ → ∞, though they do not properly capture terms of order g n τ with n > 2. This is most easily seen by scaling τ →τ := g 2 τ in eqs. (3.4) and (3.5). We return to a more careful determination of the domain of validity of eqs. (3.4) and (3.5) below.
What do these equations imply for late-time evolution? Eq. (3.4) has the solution where the relation (3.13) has been used to eliminate R BD . Eq. (3.5) similarly integrates to give As later sections show, the domain of validity of (3.17) is slightly more complicated to justify in detail because of the potential disparity in scale between the qubit gap ω and the O(g 2 ) factors. The discussion is simplest in the 'non-degenerate' regime, for which however, because then O(ω) and O(g 2 ) effects do not interfere with one another. We therefore assume (3.19) to be true in what follows (though the limit of smaller ω goes 7 In this qubit evolution differs from standard discussions of exponential decays of unstable states. Decays can proceed indefinitely with constant rate because decay daughters can escape to infinity and so do not accumulate probability in the same way as does a qubit.

JHEP02(2020)053
through similarly, as described in [21]). With this choice we may use Σ ω and the solution (3.17) becomes which is valid for times as late as ωτ ∼ O(ω/(g 2 C BD )). Note in particular that (3.19) implies that both g 2 C BD /ω 1 and g 2 ∆ BD /ω 1, which are small but not necessarily as small as O(g 2 ). The assumption (3.19) also implies that the oscillations of (3.20) are underdamped, in that they are very fast compared to the solution's decay time.
Because C BD > 0 (see (3.6)) these solutions describe relaxation towards a late-time static state which is clearly a thermal state -of the form e −βh /Tr[e −βh ] -with temperature given by the Gibbons-Hawking formula The relaxation is predicted to be exponential -proportional to e −τ /ξ -with a characteristic time-scale that differs for the diagonal and off-diagonal components of : (3.23) Asymptotic forms for this expression are given in subsequent sections.

Validity of the Markovian limit
We next quantify the domain of validity of the Markovian approximation as a function of the problem's parameters: H, M eff (or m and coupling ξ), ω, and g. As emphasized above, the Markovian evolution eqs. (3.4) and (3.5) rely on two conditions. First, the focus must be on late times, τ τ c , compared to the correlation (or fall-off) time defined by W BD (τ ); second, higher terms in the Taylor expansion of I (τ − s) in powers of s within the integrand should be negligible. Since Re[W BD (s)] falls off exponentially with s for s τ c , after integration the neglect of derivatives relies on τ c ∂ τ being small compared with .
It is this second condition whose validity imposes constraints on the system parameters. To see how, recall that the Markovian solutions found above for I 11 (τ ) and I 12 (τ ) are linear combinations of exponentials of the form exp [(−1/ξ + iΦ) τ ], where 1/ξ ∼ g 2 C BD and the phase Φ is either zero or ω. Consistency therefore requires | − 1/ξ + iΦ| 1/τ c . Chasing through the definitions this turns out to imply (see appendix C or [21] for details) the conditions

24)
JHEP02(2020)053   where primes denote differentiation with respect to ω. We explore the implications of these conditions in the non-degenerate special case when (3.19) is satisfied, and so Similar conditions for validity of the Markovian approximation can also be made when (3.19) is not satisfied but for simplicity we do not pursue these further here (see, however, [21]).
To explore the implications of (3.24) for the parameters ω, H, M eff and g it helps to have asymptotic expressions for the integrals C BD and ∆ BD in various limits. These are summarized for convenience 8 in table 2. Since the first two conditions in (3.24) are proportional to g they tend to be satisfied once g is chosen deep enough in the perturbative regime. How deep depends on the relative sizes of the other parameters, as can be seen from Columns 2 and 4 of table 2.
JHEP02(2020)053 Table 3. Asymptotic forms in different parameter regimes for the two g-independent conditions for the validity of the Markov approximation, as computed using the asymptotic forms provided in It is the second two conditions of (3.24) that cannot be ensured simply by making g small. The implications of the estimates in table 2 for these two quantities are summarized in table 3. The final column indicates with checks or crosses whether these quantities can be small enough to allow a Markovian approximation. The first observation emerging from these tables is that a Markovian approximation necessarily requires ω/H 1. (That is, the Markovian approximation can be satisfied only if ω is much smaller than the Gibbons-Hawking temperature. This result is intuitive because, if ω were much bigger than the temperature, interactions with the field (which is for the qubit effectively a thermal bath) become very inefficient at erasing the correlations whose absence underlies the qubit's thermalization and Markovian evolution.) Inspection of the last two rows of table 3 shows that it is the condition |ωC BD /C BD | 1 that generically fails if ω H. More generally, once ω H is satisfied conditions (3.24) can also be satisfied at sufficiently late times for any value of M eff /H by choosing ω and g to be sufficiently small. In practice, for some choices for M eff the allowed value of ω can be small enough to push JHEP02(2020)053 it to the point where (3.25) is no longer true. In such a case the above formulae need not hold and a separate analysis must be done (along the lines given in [21]). Two check marks in table 3 indicate that the required value for ω/H is not parametrically small in a way that depends sensitively on M eff /H. Table 4 displays the asymptotic behaviour of all of the conditions in both (3.24) and (3.25) to see how restrictive condition (3.25) is for situations that have only one check mark in table 3. We note that satisfying all conditions in (3.24) and (3.25) can be done in all three cases by making g and ω small enough. The allowed range for ω/H for which the Markovian limit applies becomes smaller and smaller for M eff /H 1 and ω/H e −πM eff /H is exponentially small when M eff /H 1. By contrast no obstruction seems to arise to a Markovian limit when M eff /H is O(1). We remark that in the regimes of both large and small M eff /H inspection of table 2 shows that the relaxation time-scale ξ = (g 2 C BD ) −1 also becomes parametrically long. For large M eff /H this seems due to the inefficiency of thermal interactions due to the Boltzmann suppression of excited field states. For small M eff /H we think the system instead displays critical slowing down due to the large fluctuations that become possible as M eff → 0.
Positivity of the Markovian limit. A reality check for these approximation schemes is to verify that the Markovian limit preserves the hermiticity and normalization of the qubit density matrix. The argument given here to this end closely parallels the one given in [21].
Eqs. (3.4) and (3.5) can indeed be written in the form of (3.26) provided that the Kossakowski matrix is given by where the identity R BD = C BD + S BD from (3.13) has been used. The eigenvalues of the Kossakowski matrix (3.27) are and λ c 3 = 2g 2 (C BD − C 2 BD + S 2 BD + ∆ 2 BD ) .
At first sight the third eigenvalue is a problem because it is negative. We now briefly argue why this is not a problem: we show that within the domain of validity of the Markovian limit the eigenvalue λ c 3 is actually consistent with zero. To see why this is so, consider for illustrative purposes the simplest regime, for which ω/H 1 and M eff /H √ 2. Table 2 shows that in this regime C BD H/(4π 2 ) is systematically larger than is ∆ BD ω log(H )/(2π 2 ). The same is also true for S BD = −C BD tanh(πω/H) −ω/(4π) C BD -see eq. (3.12). As a consequence the two nonzero eigenvalues can be approximated by The negative eigenvalue is therefore comparable to terms that were neglected when deriving the Markovian limit, and is consistent with zero at the order being studied. This also makes sense: unitarity and hermiticity are properties of the full theory's density matrix, and cannot be ruined by any approximation that accurately captures this underlying physics.

JHEP02(2020)053 4 Solved non-Markovian evolution when M eff H
The previous section argues that the Markovian limit only exists when ω H, and is most robust in this regime when M eff ∼ H. The Markovian approximation got harder to achieve in the critical-slowing-down regime where M eff H. Fundamentally, Markovian evolution becomes harder to achieve for small M eff because the Markovian approximation requires a hierarchy between the correlation time τ c → κ −1 of the Wightman function and the longer time ξ over which the qubit evolves. As M eff is decreased the correlation time κ −1 ∝ H/M 2 eff grows larger and for generic choices of the parameters the qubit does not evolve slowly enough for the derivative expansion in the Nakajima-Zwanzig equation to be a good approximation.
In this section we focus on this M eff H regime, and show how to solve for the system's evolution without requiring access to the Markovian limit. This allows us to track explicitly the memory effects coming from the very long tail of the Wightman function. By doing so, we widen the regime of parameter space for which solutions for the late-time evolution of (τ ) are explicitly known.

Non-Markovian evolution
We first recall the form of the subleading tail f (τ ) defined in (2.37), In particular the sub-leading tail f (τ ) ∼ O(e −Hτ ) falls off exponentially on the much faster time-scale 1/H than does W BD (τ ). Our strategy for solving the Nakajima-Zwanzig equation is to include the memory effects contained in W 0 e −κτ explicitly and use the Taylor-expansion argument only for integrals containing the much more quickly-falling function f (τ ).
In doing so, it will be convenient to define two integral transforms (closely related to C BD and ∆ BD ) which will prove useful in the critical slowing-down regime. In analogy to (3.6) and (3.8) we definẽ As can be seen in table 5 these functions along with their ω-derivatives,C BD and∆ BD , are notably better behaved in the M eff /H 1 regime (as opposed to their counterparts in table 2 which have singularities at M eff = ω = 0).
The first step is to rewrite the Nakajima-Zwanzig equation in terms of f (τ ) by explicitly expanding W BD (τ ) = W 0 e −κτ + f (τ ) in (2.23), so that where we (as before) set the counter-term ω 1 = −∆ BD . In components (again using

JHEP02(2020)053
The new step is to perform the same Taylor expansion of (τ − s) as before, but only for those integrals involving f (s). For these integrals this Taylor expansion is a much better approximation because the function f (τ ) is correlated over much shorter times, of order H −1 . Keeping only the leading term in these Taylor expansions then implies where higher derivatives in the Taylor series are now order (H −1 ∂ τ ) n ij and have been dropped. Note that the above equations explicitly track convolutions of I with the slow correlations ∼ e −κs .
Again replacing the upper limits of integration by ∞, but only for the integrals involving f (s) then leads to ds e −κs cos(ωs) I 11 (τ −s), (4.11) which uses the definitions (4.4) and (4.5) and, as before, R BD is given by (3.7). We next solve these equations explicitly without the need for further approximations. We return at the end to quantify the domain of validity of dropping the subdominant terms in the Taylor expansion.

The non-Markovian off-diagonal equation
It turns out that it is simpler to first solve the off-diagonal equation of motion. Using the relation 12 (τ ) = e −iωτ I 12 (τ ), we transform (4.12) to the Schrödinger-picture, which is JHEP02(2020)053 easier to solve: We solve this equation by converting it to a system of first-order linear differential equations. To this end define the three functions x 1 (τ ) := 12 (τ ) (4.14) x 2 (τ ) := * 12 (τ ) (4.15) We regard x 2 (τ ) and x 3 (τ ) only as auxiliary functons, and our real interest is in the solution for x 1 (τ ) = 12 (τ ). We now show that the functions x i (τ ) satisfy a set of coupled first-order differential equations that are easily solved. In terms of these new variables the equation of motion (4.13) can be recast as (4.17) Taking the complex conjugate of this gives a second evolution equation In vector notation, defining x(τ ) = x 1 (τ ), x 2 (τ ), x 3 (τ ) , the above three equations close as a system of differential equations of the form and initial condition

JHEP02(2020)053
Using (4.5) to relate∆ BD to ∆ BD , the determinant of the matrix O is which is non-zero in the regime of interest that follows (for which to good approximation det O −κω 2 ). Because O is invertible the unique static solution for the system is The general solution to (4.21) is given by where specifically solving the first component y 1 (τ ) = 12 (τ ) is of interest. Although the above matrix exponential can be exactly computed, the resulting solution is extremely unwieldy -we will instead compute the above solution perturbatively using standard methods [69,70]. The above equation of motion can be used to probe various regimes of small g, ω/H, M eff /H: we demonstrate it's utility here by perturbing the solution in g 1. The dynamics of the matrix exponential are governed by the eigenvalues λ O of O, which can be computed from its characteristic polynomial whereC BD and∆ BD have been related to C BD and ∆ BD using (4.4)-(4.5) (although a single factor ofC BD has been left in the first line of the above). If the qubit-field coupling were to vanish, the above equation becomes which yields 'free' eigenvalues ±iω, −κ. We here seek a solution slightly deviated away from the above free equation, as a perturbation in g 1. In doing so, several dimensionless quantities in (4.26) need to be bounded in order to control the proceeding expansion for the eigenvalues. It turns out that sufficient conditions to bound the perturbation series are In particular in the third bound, to leading order uses the fact that

JHEP02(2020)053
which may be easily verified in the regimes of interest (either when κ ω or ω κ). In particular, the bounds g 2 ∆ BD /ω 1 and 2g 2 Re[W 0 ]/(κ 2 + ω 2 ) 1 ensure that the O(g 4 ) terms in (4.26) are negligibly small. The remaining O(g 2 ) terms are therefore small given that g 2 κC BD /ω 2 1 also holds. 9 Satisfying these bounds and perturbing the eigenvalues λ O in g yields The right eigenvectors r O corresponding to the eigenvalues (4.33) are explicitly satisfies (4.20), with c j being coefficients that can be perturbatively solved for using the initial condition (4.25). After a straightforward computation, the solution for x 1 (τ ) = 12 (τ ), which then transformed to the interaction-picture is found to be This solution contains much more information than the earlier Markovian solution cf. (C.6) (or (3.20) in the Schrödinger-picture). Note that the bounds (4.28)-(4.30) ensure all the O(g 2 ) terms are small, as expected (in particular g 2 C BD /ω 1 follows from (4.28) and (4.30) in both regimes ω κ and κ ω of interest).

JHEP02(2020)053
Validity relations for the off-diagonal equation. Having established the above non-Markovian solution, the validity conditions for dropping derivatives in the expansion I 12 (τ − s) which are almost the same as (C.14) derived in the Markovian limit in section 3.3 (see also [21]). The derivation of the Markovian validity relations relies only on the exponential form of the ansatz exp [−1/ξ + iΦ]τ , which is why it applies here as well. The main difference in the above is that functions involved are the integral transforms associated with the correlation function f (rather than W BD as in the Markovian limit). All that remains is to insert (ξ, Φ) ∈ {(ξ M , 0), (ξ M , 2ω), (1/κ, ω)} into the relation (4.37) (where we recall that the Markovian time-scale is defined by 1/ξ M = g 2 C BD ). Using (ξ M , 0) and (ξ M , 2ω) leads 10 to the bounds And using (1/κ, ω) (along with the latter two bounds in (4.38) above) leads to the bounds All six of these bounds need to be satisfied in order for the above non-Markovian solution to be valid. Note that once again these relations only allow ω/H 1 (as in the Markovian case). In table 6, the functions in the above six bounds are provided in the allowed regimes (where ω/H 1 and of course M eff /H 1). We also list the bounds (4.28)-(4.30) for completeness. The bound in table 6 require less extreme values of g, ω/H and M eff /H to be satisfied (than that required by the Markovian limit in section 3.3).
Recovery of the Markovian limit. The Markovian solution for 12 (τ ), which is valid in the ω/H κ/H M 2 eff /3H 2 1 regime only (as derived in section 3), is a limit of the more general solution (4.36). The non-Markovian solution is valid for any Hτ 1, while the Markovian approximation is valid for times κτ 1. Considering times κτ 1 the exponentials in (4.36) with time-scales 1/κ become negligible (i.e. e −κτ 0), such that 10 Given |a| 1 and |a − 2b| 1 then |b| 1 has been used here.

The non-Markovian diagonal equation
The solution to the non-Markovian diagonal equation (4.11) follows through in much the same manner, so we provide a more concise overview of its derivation. As in section 4.2, the equation can be converted to a system of first-order linear differential equations. Defining the vector of variables y(τ ) := y 1 (τ ), y 2 (τ ), y 3 (τ ), y 4 (τ ), y 5 (τ ) where and differentiating these functions makes the system of differential equations close again tȯ The determinant det D = −2g 2 C BD (κ 2 + ω 2 ) 2 = 0 implies that D is invertible. The unique steady-state solution is then given by y static = −D −1 b. Using (D −1 ) 11 = −1/(2g 2 C BD ) means that the steady-state solution for static 11 = y static 1 is static 11 reproducing the expected thermal result. The general solution to (4.46) is which shows that the matrix exponential governs the approach to the thermal steady state. Again, the eigenvalues λ D of the matrix D govern the dynamics of the matrix JHEP02(2020)053 exponential, which are the roots of the characteristic equation As a perturbation in g 1 the above yields eigenvalues, where the last two eigenvalues are exact. In performing the above expansion the bounds (4.28) and (4.30), repeated for convenience, are sufficient to control the expansion of the eigenvalues (4.52). The reason for this is that 13 which bounds the O(g 2 ) terms in the second line of the characteristic equation (4.51). As before, g 2C BD /κ 1 on account of (4.54). The solution relaxes towards the thermal steady state here too with time-scales ξ M /2 and 1/κ. A tedious calculation results in the perturbative solution for x 1 (τ ) = I 11 (τ ) where which is real-valued and again valid for Hτ 1.
In the interacting theory the effective mass M λ therefore lies in precisely in the regime of parameter space associated with critical slowing down. At lowest-order in the expansion of the kernel in the Nakajima-Zwanzig equation (see [21]), the equations of motion are unaltered by the introduction of H λ . For this reason, the entire earlier analysis applies here with the mass replaced with M λ , and in particular the non-Markovian analysis in the critical slowing down regime. The qubit relaxes with time-scales where of course the time-scales must satisfy the validity relations in table 6.

JHEP02(2020)053
which is valid for a − b / ∈ Z and z / ∈ (0, 1). Using the first term in the series representation of 2 F 1 we get the |z| → ∞ asymptotics in the case that a − b / ∈ Z, where B -dependence of divergences in ∆ BD and ∆ BD We cannot take the limit → 0 + here, so we keep small (compared to the other scales) but finite. For s approaching the coincident limit, the Wightman function has the behaviour We proceed in the same manner as the analogous calculation in [21] by subtracting and adding (B.2) in the expression for ∆ BD giving ∆ BD = ∆ On the other hand, for Φ 2 = 2ω the earlier bounds imply (C.17) By explicitly using 1/ξ M = g 2 C BD in the validity relations (C. 15) and (C.17), these can be put into the simpler forms g 2 |∆ BD | 1 , g 2 |C BD | 1 , g 2 ∆ BD − 2ωC BD C BD 1 and g 2 C BD + 2ω∆ BD C BD 1 , There is a subtlety that arises when computing the asymptotic forms for the functions C BD and ∆ BD used in the main text. The subtlety arises because these functions diverge in the limit that M eff and ω both vanish together, and this makes the asymptotic form found in the regime where both M eff H and ω H depend somewhat on the order in which these parameters are made small.
Because of this table 2 focusses on the leading-order behaviour for these functions in the ω/H M eff /H 1 regime and vice-versa. For the ω/H M eff /H 1 regime, it turns out that the sub-leading terms in the asymptotic series for C BD , ∆ BD and their derivatives are O(ω/κ) which means that in this case we probe the ω κ 1 regime (see table 2). For the opposite case of M eff /H ω/H 1 the sub-leading terms are O(κ/ω) which means we actually probe the κ ω 1 regime (also given in table 2).

JHEP02(2020)053
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.