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 late-time 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.
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 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 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] 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). 1 We use Weinberg's curvature conventions [41], that differ from those of Misner, Thorne and Wheeler [42] only by an overall sign in the definition of the Riemann tensor.
2 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.
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 (axion-like) 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 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 (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 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 timeevolution problem is completely solved if the time-dependence of the reduced density matrix, is known, with initial condition (0) = 0 . (2.21) 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 .
When the dust settles, a calculation identical to that in [21] shows the qubit's interaction-picture reduced density matrix I (τ ) := e +ihτ (τ ) e −ihτ (2.22) evolves -at O(g 2 ) -according to the following Nakajima-Zwanzig equation: Here the function W BD denotes the Bunch-Davies Wightman function BD|φ(x)φ(x )|BD , evaluated along the qubit's trajectory 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 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τ Table 1. Leading asymptotic behaviour of Re[WBD(τ )] 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.
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)  [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 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. Subleading 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 Fig. 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 comparatively 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 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 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 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 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: S = ∪ n S n .
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 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 : 2ξ 11 = ξ 12 = ξ M , where (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 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.
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 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 3. We note that satisfying all conditions in (3.24) and (3.25) can be done in all three cases by making g and ω 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 Table 2. The third column indicates whether parameters exist that make the previous two terms small, without checking whether or not (3.25) is also satisfied. Two checks indicate cases where ω/H is not required to be small in an M eff /H-dependent way. (Notice in particular that checks only appear in regimes where ω/H is small).
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].
The simplest way to establish preservation of hermiticity and normalization is to rewrite the Schrödinger-picture version of eqs.  Table 4. Leading-order behaviour for the quantities appearing in conditions (3.24) and (3.25) for the three parameter regimes that receive checks in Table 3. Although the first four rows here duplicate information in Table 3 about the size of terms appearin in (3.24), the last two rows compare this information with the asymptotic form for the combination appearing in condition (3.25).
the form for a basis of 2×2 matrices, F j = 1 2 σ j , and a collection of coefficients c = [c jk ] (called the Kossakowski matrix). The utility of writing the evolution in this form is that it is known to preserve hermiticity and normalization so long as the Kossakowski matrix is positive semi-definite.
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 BDsee 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.

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 criticalslowing-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 timescale 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 defineC 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).C log(e γ+1 ω ) 2π 2 Table 5. Leading-order behaviour for the functionsCBD,C BD ,∆BD and∆ BD in various regimes of relative sizes of ω, M eff and H, where M eff /H 1. These functions are markedly better behaved near (M eff , ω) = (0, 0) than the functions listed in Table 2.  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 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 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 (τ )  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 Taking the complex conjugate of this gives a second evolution equation Using the property 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 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 of C 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 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 where the ellipsis denote negligibly small O(g 4 ) combinations of the dimensionless variables in (4.28)-(4.30). Clearly the real part of each of the three eigenvalues are negative which means that the solution sinks towards the steady state x static = 0. The first two eigenvalues clearly capture information relating to the Markovian approximation (with the Markovian time-scale ξ M ), while the latter eigenvalue corresponds to a novel time-scale 1/κ + O(g 2 ).
The right eigenvectors r O corresponding to the eigenvalues (4.33) are which can be used to solve for the solution. An ansatz of the form x(τ ) = j c j r O j e λ O j τ 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 c.f. (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).
Validity relations for the off-diagonal equation which are almost the same as (C.14) derived in the Markovian limit in §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 10 Given |a| 1 and |a − 2b| 1 then |b| 1 has been used here.
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 §3. 3). Table 6. Leading-order behaviour of the functions in the validity bounds (4.38) and (4.39) (the sizes of all which must be 1), using information from Table 5 (and also Table 2). The bounds (4.28)-(4.30) are also included, keeping in mind that g 2 CBD/κ 2g 2 Re[W0]/(κ 2 + ω 2 ) in the regimes quoted. Note that some of the bounds are trivially satisfied in the quoted regimes (for example |κC BD /(2CBD)| 1).

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 §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 (ie. e −κτ 0), such that Furthermore when the Markovian validity relations (given in Table 4) are satisfied 11 then the combination 2g 2 Re[W0] κ 2 +ω 2 is negligibly small. Using this fact it follows 12 that which is precisely the solution for the Markovian off-diagonal equation (see (C.6) in the interactionpicture).

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 §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ȯ 11 This follows for example by noting ω/κ O(g 2 ) when the Markovian condition |ω∆ BD /C BD | 1 is satisfied. In this case × ω κ 0 within the approximations used. 12 In the ω κ regime, the factor 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 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. The validity conditions are again derived with an exponential ∼ exp ([−1/ξ + iΦ] τ ) ansatz, which leads to linear combinations The bounds in Table 6 (that were set when solving the off-diagonal non-Markovian solution) are sufficient to satisfy the above bounds also. Finally, as before it is also of no surprise that the above solution simplifies to the Markovian solution when the limit κτ 1 is taken, and factors of g 2 C BD /κ 2g 2 Re[W 0 ]/(κ 2 + ω 2 ) are neglected (as demanded in the ω κ Markovian regime).

Interacting field theories and secular growth
Issues of secular growth can be explored as in our study of the accelerated qubit in Minkowski space [21,71] by adding a λφ 4 interaction to the massless theory, which we briefly discuss here. Including For the case of ν = iµ where µ = M 2 eff /H 2 − 9/4 (when M eff /H > 3/2), we write the relevant factors in polar form such that .

B
-dependence of divergences in ∆ BD and ∆ BD Here we derive the asymptotics for the divergent function ∆ BD . Using the Wightman function (A.3), but with a small-distance regulator . The integral (3.8) for ∆ BD is explicitly 14 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 (B.2) 14 We replace sinh( as 2 ) − i H 2 → sinh( We find no way to analytically compute (B.10) here, although this form is easier to use for numerical integration with computer algebra systems. Numerical investigation reveal the logarithmic -dependences (in terms of H 1 or ω 1) given in Table 2 (differentiating these forms with respect to ω leads to the -divergences for ∆ BD ). using 1/ξ M = g 2 C BD from (3.23),and defining the phases Φ 1 = 0 and Φ 2 = 2ω for the off-diagonal ansatz (A j , B j ∈ C are placeholders for the time-independent amplitudes in the corresponding solutions). These ansatzes will be used to bound the derivative terms in (C. 3

) and (C.4).
It is simpler to begin with the diagonal equation. Using the ansatz I 11 (τ ) ∝ e −2τ /ξ M results in the RHS of (C.3) containing terms To neglect the derivatives here such that I 11 (τ − s) I 11 (τ ) therefore means to satisfy the bound To derive validity relations for the off-diagonal equation, a totally analogous procedure is followed. Using the earlier ansatz (C.8) which is a linear combination of terms ∝ B j e −τ /ξ M +iΦj τ , the RHS of (C.4) contains terms Dropping the derivatives in the Taylor series therefore means to satisfy the bounds Since |C BD −i∆ BD | ≥ |C BD |, bounding the first condition automatically bounds the second. In addition to this, |C BD | is larger than each of the real and imaginary parts of the complex number on the inside of the modulus in (C.12). This implies that are sufficient to bound the derivatives and end up with the desired Markovian equation of motion. For Φ 1 = 0, the above bounds imply the first of which is essentially the same as (C.10) (up to a factor of 2). These can be put in the more useful form 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 , (C. 18) which are precisely the validition relations given in (3.24) in §3.3.
D Singularities and ordering of the limits ω → 0 and M eff → 0 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. More precisely, using the asymptotic form W BD (τ ) W 0 e −κτ in the M eff /H 1 regime (see (2.32)) the relevant functions have the following behaviour C BD 2Re[W 0 ]κ κ 2 + ω 2 and ∆ BD 2Re[W 0 ]ω κ 2 + ω 2 , (D.1) and so where κ M 2 eff /(3H). These functions therefore have different leading-order asymptotic behaviours near (ω, M eff ) = (0, 0) depending on which path is chosen to explore it in the M eff − ω plane.
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 subleading 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).