Operator thermalisation in $d>2$: Huygens or resurgence

Correlation functions of most composite operators decay exponentially with time at non-zero temperature, even in free field theories. This insight was recently codified in an OTH (operator thermalisation hypothesis). We reconsider an early example, with large $N$ free fields subjected to a singlet constraint. This study in dimensions $d>2$ motivates technical modifications of the original OTH to allow for generalised free fields. Furthermore, Huygens' principle, valid for wave equations only in even dimensions, leads to differences in thermalisation. It works straightforwardly when Huygens' principle applies, but thermalisation is more elusive if it does not apply. Instead, in odd dimensions we find a link to resurgence theory by noting that exponential relaxation is analogous to non-perturbative corrections to an asymptotic perturbation expansion. Without applying the power of resurgence technology we still find support for thermalisation in odd dimensions, although these arguments are incomplete.


Introduction
Free field theories are the simplest and most prominent examples of (super-)integrable quantum field theories (QFTs), rendered exactly solvable by the existence of an infinite set of conserved charges. A direct consequence of the presence of such charges is a severely constrained time evolution even in thermal backgrounds. In particular, simple operators in free QFTs fail to satisfy the requirements of the eigenstate thermalisation hypothesis [1,2] and their late time behaviour is therefore unlikely to approach ensemble averages, tantamount to the absence of thermalisation.
Nonetheless, it is known that nontrivial interference effects can effectively mimic equilibration. For example, after quantum quenches [3][4][5][6], correlation functions in free QFTs approach those of a generalised Gibbs ensemble [4,[7][8][9], characterised by chemical potentials for all conserved charges which in free QFTs is equivalent to a momentum-dependent temperature. Similarly, nontrivial time dependence arises when considering composite operators. Such operators can in fact interact with the thermal bath and as such exhibit a range of phenomena that are usually attributed to their interacting counterparts. For instance, their correlation functions can exhibit exponential decay at late times [10] and their spectral densities have support in the deeply off-shell regime [10,11], reminiscent of collision-less Landau damping [12].
Clearly, the composite nature of an operator is a necessary condition for its effective thermalisation, since only then does it couple to a thermal bath, indicated by a temperature dependence of its response functions. On the other hand, to which extent it is also a sufficient condition is less understood. In recent work [13,14], a simple criterion has been formulated that guarantees the absence of thermalisation of a given operator, characterised by a lack of exponentially decaying contributions to its linear response function. In addition, it was conjectured that a converse statement can be made and any operator that fails this non-thermalisation condition in fact thermalises. This conjecture was introduced as the Operator Thermalisation Hypothesis (OTH).
This note aims to shed light on several remaining puzzles. First, in singlet models [10], the calculated correlation functions were observed to display exponential decay in even dimensions d > 2. This decay is directly related to the thermalisation later extracted in [13] by arguments which however do not resolve a difference between even and odd dimensions. Since the odd-dimensional singlet model correlation functions do not decay exponentially, the results appear to be in tension with each other. There is something to learn about thermalisation or singlet models, in fact both, from a closer study. The requisite developments of concepts indeed leads to a more precise formulation of the OTH. Second, the singlet model study demonstrated how phases below or above a critical temperature exhibit different relaxation properties, most plainly for auto-correlators. Since response functions diagnose thermalisation, previous singlet model studies should be extended with results on response functions in different phases.
Thus, we will study OTH in a particular class of free field theories, namely those with a large-N singlet constraint. These theories have received widespread attention in the context of gauge/gravity duality as the holographic duals of gravitational theories with an infinite tower of massless fields of higher spin. They exhibit an interesting thermal structure on compact spaces with a large N confinement/deconfinement phase transition. In ordinary AdS/CFT, this transition is also present and can be mapped to the Hawking-Page transition from thermal AdS to the large AdS black hole in the bulk. In the deconfined phase, thermalisation in holographic gauge theories is in direct correspondence with black hole formation and equilibration in the bulk. Understanding thermal properties of free singlet models thus provides insight on putative black holes in higher spin gravity. More generally, however, they allow one to disentangle generic properties of composite operators from those particular to strong coupling, thereby teaching valuable lessons on the inner workings of gauge/gravity duality.
In the low temperature phase we observe the absence of thermalisation to leading order in 1/N , in complete accordance with the OTH. Below the phase transition, a composite operator tr(Φ(x)Φ(x)), built of N adjoint scalars, plays the role of a generalised free field with interaction strength of order 1/N . To leading order, it obeys our generalisation of the non-thermalisation condition, confirmed by the absence of temperature dependent contributions to its response functions and in particular the lack of exponential damping. This is generic to all QFTs that admit a description in terms of generalised free fields and the thermal version of this concept will be presented below. At high temperatures, the time dependence becomes significantly richer. Response functions become temperature dependent and are characterised by non-analyticities off the real axis in the complex frequency plane. They describe a damped response to sources, with a power law tail and sub-leading exponentially decaying contributions. The latter contribution is the exponential damping predicted by the OTH. The presence of the power law tail implies that information about the source is retained to a larger degree than in standard thermalisation, although parts are effectively lost in exponentially damped terms.
The general lessons from our study concern details of the formulation of OTH, and the difference between even and odd dimensions. Indeed, it is well known that the interior of the light cone plays a fundamentally different role in wave propagation in even and odd dimensions (cf Huygen's principle and Hadamard's problem [15]). By explicitly focusing on evaluating correlators close to the light cone we reduce the difference between odd and even dimensions, and identify the damped quantities that continue analytically between different dimensions, to put d > 2 OTH on a firmer footing. This light cone limit notwithstanding, crucial differences between even and odd dimensions remain. While OTH can be confirmed straightforwardly in even dimensions, we observe that subtleties involved in isolating exponentially decaying terms in the response functions become critical in odd dimensions. We describe the difficulties and find some support for thermalisation, but also indications that the resolution requires more powerful tools from the theory of resurgence [16][17][18][19]. That cautionary observation aside, our scrutiny of OTH permits us to give a more precise formulation of both the hypothesis and the converse non-thermalisation condition in all d > 2.
Our paper is organised as follows. In section 2.1, we introduce the concept of operator thermalisation, non-thermalisation and the role stable thermal quasi-particles and generalised free fields. Before going into basics of singlet models in 2.2 we also introduce the potential relation of thermalisation to resurgence. In section 3.1, we then deduce and discuss absence of exponential relaxation in singlet model response functions in the low temperature phase. The high temperature phase, which displays relaxation in even dimensions and appears to allow for it in odd dimensions, is analysed in section 3.2 and a discussion in section 4 leads up to our conclusions 5.

Operator thermalisation
Sabella-Garnier et al formulated the operator thermalisation hypothesis in [13] and considered the thermalisation properties of operator correlation functions in a fixed background rather than operator expectation values in the presence of of assumptions on the energy spectrum, as done by the eigenstate thermalisation hypothesis, ETH [1,2]. They discuss thermalisation in terms of an exponentially fast return to equilibrium of operator expectation values in response to a perturbation by the operator in question. More precisely, in [13] the retarded Green's function of the operator in question is taken to define thermalisation of a perturbation, when it decays exponentially, in line with the retarded Green's function encoding the linearised response of the operator O induced by a perturbation by the same operator O. The requirement of exponential decay for a perturbation to thermalise corresponds to the intuition that a thermalising perturbation is "forgotten" by the system at late times. The latter means that exponential precision would be required in order to fully reconstruct the source from the response of the medium.
A motivation behind the operator thermalisation hypothesis, and one of its strengths, is that it can be used to study surprising similarities between ordinary interacting systems and free or integrable systems [6,10,11,13]. While pure exponential decay occurs in free systems in contrast to naive expectations, it is generally masked by leading power law decay for d > 2, as well as multiplied by inverse powers of time. We will provide such examples below. In reviewing the operator thermalisation hypothesis, we will therefore introduce new terminology which precisely captures these features. In effect, we demonstrate an operator non-thermalisation condition which excludes this kind of partial thermalisation, and state a converse partial operator thermalisation hypothesis. Our arguments are essentially copied from [13], and the "partial" qualifier only indicates a slight shift of definitions. The new definitions are important for consistency with the examples we discuss, but the idea is approximately the same.

Partial operator thermalisation
We define partial thermalisation of an operator O to mean that: The retarded Green's function of O contains terms with exponentially damped factors at late times. This definition allows for leading power-law decay, and exponential terms which are only sub-leading 1 . In such cases, time evolution still "forgets" part of the initial perturbation, but not all of it. Clearly, partial thermalisation includes the thermalisation notion discussed in [13] and the more conventional notion of approach to a thermal ensemble, but it is a broader concept 2 .
Crucially, partial operator thermalisation captures the observation that conservation laws prevent some operators in free or integrable theories to thermalise, but that almost all other operators thermalise partially. A special class of non-thermalising operators was characterised by Sabella-Garnier et al [13]. We will see that these operators do not even thermalise partially. In essence, these non-thermalising operators are generalisations of free fields which satisfy a sharp dispersion relation relating energy to momentum. Formally, the conditions on the operators are given by the mathematical descriptions below. Physically, they correspond to stable thermal quasi-particle fields having clear-cut dispersion relations, which are permitted to differ from those of free relativistic particles. One may invoke the Narnhofer-Requardt-Thirring theorem [21] to argue that they describe a sector of the thermal system which is completely free from interactions, except for modified dispersion relations. The theorem permits other sectors, but they are completely decoupled from the quasi-particles.
We interpret the operator thermalisation hypothesis proposed in [13] to state that any other local operator, not representing a stable quasi-particle field, will thermalise. This is the converse of the above non-thermalisation condition. For it to hold, the notion of thermalisation has to be weakened to partial thermalisation. Thus, we propose a more precise partial operator thermalisation hypothesis: Any local operator not representing what we call a thermal generalised free field 3 , or a generalised quasi-particle field, thermalises partially. Note that we still have not proven this hypothesis, though we find it reasonable. All the plausibility arguments in [13] still apply.

Non-thermalisation and the thermalisation hypothesis
We now proceed to essentially repeat the arguments of [13], expressed in our terminology.
Consider a stable thermal quasi-particle operator, which we denote Q(t, x) to distinguish it from more general local operators O(t, x). By definition it has a definite dispersion relation. In finite volume and in a basis which simultaneously diagonalises energy and momentum, this means that the transitions Q can mediate between momentum states determine the simultaneous transitions between energy eigenvalues. We do not need to know if there is a single functional relation between momentum and energy for the operator Q or if there are several branches of solutions to the dispersion relations coupling to Q. To reproduce branch cuts which can be found, for example, in singlet models, it will turn out to be important to allow for a growth of the number of solutions to dispersion relations with volume.
The retarded thermal Green's function is which can be expanded in a sum of expectation values 2) where s is the spin of the operator Q. Making use of translations and inserting a complete set of states In Fourier space where ǫ > 0 is infinitesimal. Now, the special properties of the quasi-particle operator Q lead to a proof of nonthermalisation. Denoting by M the number of different branches of solutions labeled j = 1, . . . , M to the dispersion relations for Q and defining the residue functions we find (2.8) The frequencies and wave numbers above are related to the matrix elements m| Q |n by signifying that all contributions from the operator Q are due to transitions between states whose energies and momenta differ by amounts related by the allowed dispersion relations in eq. (2.6). For a more detailed analysis of thermodynamic and large N limits, it may become useful to allow an effective temperature dependence in the dispersion relations contributing to eq. (2.8). Noting that Ω Q j (k) has to be real by definition, the retarded Green's function only has singularities on the real axis, which is tantamount to non-thermalisation of the operator Q. For finite M , the singularities are manifestly poles. If M grows without bound in the thermodynamic or large N limit, branch cuts may also arise, but they will be on the real axis. There will still not even be partial thermalisation, since only singularities off the real axis can produce exponentially decaying terms.
The converse of the original non-thermalisation result would be that only stable quasiparticle operators are non-thermalising. Allowing for partial thermalisation, which includes power law fall-offs related to branch cuts on the real axis, it seems judicious to consider branch cuts also in the non-thermalisation results. Thus we are led to allow unbounded M . This generalisation replaces quasi-particles with generalised quasi-particles or thermal generalised free fields.
The OTH in our version becomes: All local operators which are not generalised quasiparticles thermalise partially. The original plausibility arguments of [13] remain, and this adjusted version survives all tests we have considered.

Thermalisation and resurgence
Below we will introduce examples of retarded Green's functions with asymptotic late time expansions containing both inverse powers and damped exponentials of time. In free systems, they force us to consider the partial, and more general, version of operator thermalisation, which allows for the possibility that exponential damping terms are sub-leading. Unfortunately, the price for the generalisation is another level of mathematical sophistication. It is required for a physical reason: Only under very special circumstances, e.g. when an asymptotic series of inverse powers terminates, is it possible to operationally separate sub-leading exponentials from more important inverse powers. Only under these special circumstances can we have a chance to resolve and observe the damped exponentials, even in principle.
This discussion is parallel to the potentially more familiar discussion about prescription dependence of non-perturbative terms in quantum mechanics and in quantum field theory. There, one encounters non-perturbative exponentials e −1/g 2 complementing power series in a coupling g. Substituting where t is time and β is inverse temperature, we are alerted to the possibility that thermalisation, signalled by exponential damping at late times, can be analogous to nonperturbative effects. The analogy indeed holds for standard Green's functions: Their late time expansion in inverse powers of t/β is typically asymptotic rather than convergent, and exponential terms can sometimes be extracted from integral representations of the Green's functions. Cases with terminating or at least a convergent (inverse) power series would be useful in practice, and would allow unambiguous identification of exponentials, but are exceptional.
The beautiful idea that there is a relation between the form of non-perturbative terms and the divergence of asymptotic series [22] can be systematised in non-perturbative techniques like Borel resummation, but does not always yield a unique answer for the series. To be clear, for thermal Green's functions in free field theory, the integral representations are unambiguous. A series representation does not improve the already complete encoding of a response function. However, a well-defined representation of the result of the integral in a double series expansion with inverse powers and exponentials as above, a trans-series in the framework of resurgence theory [17][18][19], would lend itself nicely to an extended definition of partial thermalisation. The response function would be said to thermalise partially if the series contained exponentials 4 .

Thermal singlet models
In order to distinguish low and high temperatures, we consider free field theories on R×S d−1 leading to a characteristic temperature scaling as 1/R, the inverse of the radius R of the sphere S d−1 . To make the distinction sharper we consider a large number N of fields. A large N will then allow for qualitatively different limits for physics below and above the characteristic temperature. We consider a scalar field transforming in a representation, usually fundamental or adjoint, of some large N symmetry group, for example U(N ) or O(N ). Projection onto the singlet sector is achieved by weakly gauging the symmetry, i.e. introducing a gauge field A µ in the limit of vanishing gauge coupling, where only the zero mode α ∼ S d−1 A 0 that imposes the Gauss' law constraint remains. We will focus attention on correlation functions on scales much smaller than R, corresponding to times and distances t ≪ R, |x| ≈ Rθ ≪ R, where θ is the polar angle on the sphere. The entire difference between low and high temperature physics in effectively flat space can then be encoded completely in functions ρ(λ), which appear as eigenvalue densities in the more detailed description in the next two paragraphs.
At finite temperature, the integral over the gauge field can be recast into a unitary matrix model, where the projection onto singlets results from the integral over the gauge group over unitary matrices [24] corresponding to the Polyakov loop operator, P ∼ e i S 1 dτ α , or gauge holonomy around the thermal circle [25]. The distribution of the large N number of matrix eigenvalues then controls the thermal behaviour.
At large N , the model can be solved in a saddle point approximation [24]. This is readily achieved by introducing the eigenvalue density ρ(λ). At low temperatures, T < O(1), the dominant saddle corresponds to a constant eigenvalue distribution, ρ(λ) = 1 2π . This is the confined phase, with a free energy of order N 0 . At intermediate temperatures whose N scaling depends on the representation under consideration, there is a transition to a deconfined phase, characterised by a free energy that is extensive in N . At very high temperatures, the eigenvalue distribution becomes a delta-function 5 , ρ(λ) → δ(λ).
Correlation functions of singlet operators can be constructed through finite temperature Wick contractions. For simplicity, we focus here on the scalar singlet primary, x)) for scalars in the adjoint representation, whose time ordered twopoint function is given by [10] (2.10) where we have used rotational and time translational invariance to set one of the insertion points to zero. The pre-factor has been chosen to simplify the expression, while the operator is normalised such that its two-point function is of order N 0 . Eq.(2.10) can formally be derived using the aforementioned Wick contraction, as well as the fact that the unitary matrix is represented in the scalar kinetic term like a temporal gauge field. The retarded Green's function can be extracted using its definition, G R (t, x) = Θ(t)Im G(t, x). It is simple to see that the purely thermal contributions to eq. (2.10) are real. An imaginary part can thus arise only from the vacuum piece, and the mixed thermal-vacuum term.
Explicitly, one finds [10] G R (t, x) = Θ(t) Im 1 (cos t − cos θ) d−2 where in the last line we have introduced the k-th Fourier cosine coefficient of the eigenvalue distribution, ρ k = dλ ρ(λ) cos(kλ). We note that the infinite series in the second term captures all temperature dependence, and in fact is precisely that of the thermal Feynman propagator of the fundamental scalar field, when all ρ m become equal, which is the case in the high temperature limit.

Thermalisation in singlet models
A number of challenges to the OTH may be tested in thermal singlet modes in d > 2.
In this section, we describe our technical results, which support the hypothesis in even dimensions, given the adjustments we have introduced in section 2.1. In odd dimensions the interpretation of results is intricate, and will be deferred to the discussion 4. The low and high temperature phases of the singlet models are qualitatively different and are discussed separately below, with equations specialised to scalars in the adjoint representation. In both cases, the concrete operator under study is the lowest dimension singlet operator O(t, x) = 1 N tr(Φ 2 (t, x)).

Low temperatures: T < T H
As noted in the thermal singlet model section 2.2, the eigenvalue distribution at low temperatures is constant, ρ(λ) = 1 2π , and thus ρ 0 = 1 and ρ k =0 = 0. For the retarded Green's function (2.11), this implies in the large N limit. From the explicit lack of exponentials we see that O(t, x) fails to thermalise at low temperatures. For completeness, let us take the "thermodynamic limit" of large R corresponding to t, θ ≪ 1, and Fourier transform, thus for example obtaining where µ is a renormalisation scale. In this Lorentz invariant expression, there is only one branch cut located at ω 2 > k 2 on the real line, representing a continuum of physical excitations on top of the vacuum state. That (3.2) is analytic everywhere off the real line corresponds one-to-one with the fact that the corresponding expression in configuration space lacks exponentially decaying contributions. We see explicitly that it is useful to extend the notions of the non-thermalisation condition beyond poles in the frequency plane to cuts, as long as they are on the real axis. In position space, the corresponding thermodynamic limit of (3.1) involves power-law fall-off, and we will find similar fall-offs to be general consequences of free field dynamics in d > 2 below, even for response functions of operators that display relaxation after long time in exponentially decaying terms. The physical origin of non-thermalisation is clear from large N considerations. At low T , one finds for the connected components of n-point functions Here, the conservation of individual momentum modes is explicit to zeroth order in 1/N . In consequence, quasi-particles remain intact to this order. Of course, this argument is rather superficial, but can be made more precise by properly constructing the effective action, for example using collective field theory [27]. By the above argument, taking into account 1/N corrections will reveal nontrivial features in the response functions even below the phase transition. While this requires a finite N analysis, and is therefore beyond the scope of this work, even the leading order behaviour can change drastically once occupation numbers in the thermal background are of order of the inverse coupling. Indeed, as we will show now, this is what happens in the high temperature phase.

High temperatures: T ≫ T H
At very high temperature, the eigenvalue distribution can be approximated by a deltafunction. One thus obtains for the Fourier cosine coefficients Note that one should only really expect effective thermalisation in the "thermodynamic limit" of large R, here corresponding to t ≪ 1, θ ≪ 1 and β ≪ 1. In this regime the retarded Green's function (2.11) becomes upon insertion of (3.5) Clearly, the operator now responds to the thermal bath, which may induce thermalisation. In fact, the second term represents the cross term between vacuum and thermal propagation contributing to the response function of the quadratic composite operator O(t, x).

d = 4
To get a better understanding of the precise dynamics, we will confine ourselves to d = 4, since generalisation to higher even dimensions is simple once the basic ingredients are understood. There, which can be simplified to Evidently, the Green's function falls off as a power law, with a power that is smaller than in vacuum. This is in fact a manifestation of the effective dimensional reduction that is prevalent in generic thermal systems in the high temperature limit (see e.g. [28]). However, judging by the coth term there are sub-leading exponentially decaying contributions. This may be further illuminated by Fourier transforming (3.8), yielding [10], This expression allows us to map the late-time dominant behaviour of (3.8) to the branch cut in the ω plane located between −k and k on the real line and the subdominant exponential decay to the branch cuts located off the real line. Similar analytic structures are discussed in [29]. It can be contrasted with that of eq. (3.2). Let us now return to how thermalisation could be consistent with the effective action arguments presented in the low temperature discussion 3.1. Only large N counting, which is the same at high temperature, seemed to be important. The large N suppression of interactions is indeed the same as at low temperature, but the action (3.4) assumes the vanishing of thermal one-point functions O β . Above the phase transition, the equilibrium background expectation value is non-zero and of order N , which invalidates the argument that individual O momentum modes are conserved in the large N limit, due to order N 0 interactions with the background. Generalised quasi-particles are then not intact in the large N expansion, although their response functions are well-defined. The thermalisation of O ensures that O does not represent a generalised quasi-particle, by the arguments of subsection 2.1.2. As explained above, this is consistent with large N counting, thanks to the thermal condensate of O above the critical temperature, which eq. (3.6) thus probes indirectly.

General d > 2.
The above retarded Green's function of an operator quadratic in free fields clearly separates into a vacuum-vacuum term and a mixed vacuum-thermal term. Higher powers of free fields also decompose analogously. (Purely thermal terms will not contribute to the retarded propagator.) Now, the vacuum factors differ significantly in behaviour between odd and even dimensions. The imaginary part of (x 2 − t 2 ) −(n+1) for integer n ≥ 0 is given by which demonstrates that the support of the retarded Green's function is confined to the light cone for even d, while square root branch cuts ensures support also inside the light cone for odd d. This is a known property of the wave equation, which evidently is inherited by thermal systems probed by composite operators built of powers of free fields. While the behaviour in the interior of the light cone is interesting, both for other correlation functions than the retarded Green's function, and for odd d, the temperature dependent term of the response function is entirely determined by the factor which multiplies a simple light cone divergence. We thus factor out its singular light cone behaviour and study the behaviour of what might be called the position space "residue" of the singularity, by abuse of terminology. The light-cone factor isolated from eq. (3.6) is then This expression lends itself to a comparatively uniform treatment independently of dimension, and it measures the effect of the heat bath on the light cone in position space. The functionsS d (t, x) and its light cone limit S d (t) are discussed in the appendix A.
In even dimensions the calculation confirms thermalisation on the light cone, essentially by expressing S d (t) as an expansion in modified Bessel functions each of which equals a decaying exponential times a terminating sum of inverse powers for even dimensions (i.e. half-integer orders of the Bessel function). Details concerning finiteness of the expressions are also given in the appendix A. For all even dimensions the positive m terms in the series above explicitly yield exponentially decaying terms, which are of a form that cannot cancel with other exponentials. The negative m terms produce power law fall-off. Thus, the response functions signal partial thermalisation on the light cone in even dimensions. The odd-dimensional case is significantly more subtle. A similar treatment of the Bessel function terms in series (3.12) now leads to an asymptotic expansion in inverse powers for each m, which does not terminate. Hence, it is far from clear what significance to attach to exponentially small terms. If the asymptotic series is truncated, the error terms will be larger than the exponentials we have extracted, even if exponentially improved expansions are used [30].

Discussion
Our description of an important class of non-thermalising operators as generalised free field operators in section 3 connects to the intuition that such operators should not thermalise. Generally, however, naive intuition is treacherous, and our study is founded on the observation that free field equations of motion do not generally guarantee absence of relaxation for operators which are non-linear in free fields. Composite operators are regularised operators belonging to this class, and as described in the introduction 1, they have in many instances been shown to display the decay which we take to define thermalisation. The idea of the OTH is that the implication could go in the other direction: Operators that do not thermalise partially would have to be generalised quasi-particle operators. Or equivalently, any other operators thermalise partially.
Known thermal behaviour of singlet models motivated a closer study of response functions in order to compare with the OTH in dimensions d > 2. The d > 2 treatment of [13] is somewhat less detailed than the d = 2 discussion, and we were able to resolve new even/odd dimension differences in eqs. (3.11-3.12) from the high temperature phase of the singlet model response functions. To get expressions which depend analytically on d, it proved important to focus on the light cone. Thereby, the qualitative difference between the support for the Green's functions in the light cone, related to the absence or not of Huygens' principle were factored out.
In the large N limit, where there is a phase transition, and below the critical temperature, the response function (3.1) lacks exponentially decaying terms and the individual momentum modes are independently conserved. This means non-thermalisation and also the presence of generalised quasi-particles of definite momenta, as expected from the general non-thermalisation results. Clearly, we can only expect a precise match to the general operator thermalisation theory, described in section 2.1.2, to be valid at leading order in small 1/N . To the extent that the generalised quasi-particle picture holds, we can rely on non-thermalisation. Indeed, the idea that the general theory applies parametrically close to ideal cases makes the results much more powerful. In this example, we see how it works.
Above the critical temperature, the response functions develop exponentially decaying terms, as for example in the d = 4 expressions (3.8) and (A.3). These examples clearly show how power law tails and damped exponentials combine non-trivially. Indeed, such terms which generally appear in d > 2 motivate us to consider partial operator thermalisation. This partial thermalisation concept also simplifies the non-thermalisation results for stable thermal quasi-particles, by allowing branch cuts in the thermodynamic limit as in the paragraph after eq. (2.9).
To conclude the match of the OTH and singlet model response functions we should now argue that the operators we consider fail to be generalised quasi-particle operators above the critical temperature. Without going deeply into the physics of singlet models, we have found a suitable mechanism, namely that the background condensate in the high temperature phase modifies the propagation of perturbations at order N 0 which is too much for a generalised free field, unless there is extreme fine-tuning.
Partial thermalisation diffuses the dichotomy between thermalising and non-thermalising operators to some degree, but in even dimensions calculations like (A.9) and (A.12) demonstrate the general structure from modified Bessel functions of order d−3 2 . At half integer order the resulting functions are simple polynomials of exponentials exp (−4πt/β) and powers of β/t. The further sums in the thermal response functions primarily gives rise to an infinite series of higher order terms exp (−4πnt/β), where n are integers, but the expansion in β/t terminates and the damped exponential terms can be distinguished from the resulting polynomial. Thermalisation can be confirmed although with a bit more work than if power law tails had not been present. This comparatively simple procedure works in even dimensions, when the whole effect of the induced thermalisation is confined to the light cone by Huygens' principle.
In odd dimensions Huygens' principle does not apply and some of the induced thermalisation diffuses into the interior of the light cone. The series encoding the thermalisation on the light cone, which corresponds to the polynomial in β/t, now fails to terminate. Instead it produces an infinite asymptotic expansion controlled by the asymptotic expansion of modified Bessel functions of integer order. The resulting series in β/t is divergent and the task to identify sub-leading exponentials becomes quite subtle. The potential meaning of exponential terms can only be ascertained within a larger framework, such as the study of resurgence of asymptotic series. In such a framework one should be able to assign a meaning to partial thermalisation of composite operators in odd-dimensional free field theories, but a firm conclusion is beyond the scope of the present work.
Tentatively, the Borel summability of modified Bessel asymptotic expansion indicates that there are no exponential correction terms in its asymptotic expansions, which would suggest that the exponential terms we actually find in the appendix are not masked, but on the other hand error terms of even doubly improved asymptotic series are of the same order as the sub-leading exponential terms.

Conclusions
We have refined the operator thermalisation concept and the OTH, and related it to generalised free fields. Except in the special case d = 2, operator thermalisation is generally incomplete and partial, since there are power law tails that dominate exponentially decaying terms at late times. This finding establishes the intermediate nature of thermalisation in free field theories: while exponential relaxation is ubiquitous, it typically coexists with the more unyielding time dependence expected from the presence of conservation laws. In our model system, large N singlet models, we have found both non-thermalising and thermalising behaviour of the same operator: generalised quasi-particle behaviour without exponential relaxation below the critical temperature, and thermalising exponential behaviour above the critical temperature. Importantly, the operator thermalisation concepts turn out to be applicable to operators which only satisfy the theoretical conditions in a limit, in this case when 1/N vanishes. This enhances the scope of our analysis.
The analysis is comparatively straightforward in even dimensions, where Huygens' principle holds and ensures that the thermalised responses induced by a heat bath are localised to the light cone. In contrast, the thermalised responses in odd dimensions are quite intricate due to their distribution over the forward light cone and the whole of its interior. We refrain from formulating a definite conclusion in odd dimensions, since we believe in a deeper conceptual analysis. The importance of simultaneous infinite expansions in inverse powers, and decaying exponentials, of time, suggests resurgent analysis. A connection between thermalisation of integrable systems and resurgence may find further applications.
Some properties of singlet models that are highlighted by our study generate further questions. For example, an efficient description of the high temperature phase remains elusive. We expect that all standard composite operators will thermalise and no longer represent generalised quasi-particles. The fundamental free fields Φ describe the thermodynamics of the high temperature "deconfined" limit well, but they do not represent physical singlet states. Do they provide the best description, or are there better alternatives? There are also holographic gravity duals to these questions, since singlet models are limits of large N gauge theories, some of which are conformal.
Finally, we find it inspiring to contemplate other conformal or integrable systems, in particular in odd dimensions, where resurgence appears to be fundamental.
Since in even d the retarded Green's function only has support on the light cone we will only evaluate the sum there. We have The sum (A.1) can be rewritten using the formula 1 y α = 1 Γ(α) where ϑ(z, q) is the third Jacobi theta function andt = t/β. Using the modular transformation property of ϑ, we obtain ϑ(rt, e −r ) = π r e −rt 2 ϑ(iπt, e − π 2 r ) = π r e −rt 2  The last term is divergent but will be canceled by a divergence stemming from the sum. In order to allow for such a cancellation we regularise the integral by introducing an exponential suppression e −r/R with R taken to infinity after performing the integral. We have which is valid for d > 3. As we are interested in large t-behaviour we employ the asymptotic expansion As a check we set d = 4 to compare with (A.3). We obtain which agrees with (A.3).
For odd dimensions we investigate (A.11) for d = 3 even though the expression is technically only valid for d > 3 because of the logarithmic divergence in the first term. As can be seen, the asymptotic power series does not terminate and the sub-leading exponentials are masked, suggesting that more powerful asymptotic or resurgence method should be employed.