Analytic structure of nonhydrodynamic modes in kinetic theory

How physical systems approach hydrodynamic behavior is governed by the decay of nonhydrodynamic modes. Here, we start from a relativistic kinetic theory that encodes relaxation mechanisms governed by different timescales thus sharing essential features of generic weakly coupled nonequilib- rium systems. By analytically solving for the retarded correlation functions, we clarify how branch cuts arise generically from noncollective particle excitations, how they interface with poles arising from collective hydrodynamic excitations, and to what extent the appearance of poles remains at best an ambiguous signature for the onset of fluid dynamic behavior. We observe that processes that are slower than the hydrodynamic relaxation timescale can make a system that has already reached fluid dynamic behavior to fall out of hydrodynamics at late times. In addition, the analytical control over this model allows us to explicitly demonstrate how the hydrodynamic gradient expansion of the correlation functions can be Borel resummed such that the full nonperturbative information is recovered using perturbative input only.

How physical systems approach hydrodynamic behavior is governed by the decay of nonhydrodynamic modes. Here, we start from a relativistic kinetic theory that encodes relaxation mechanisms governed by different timescales thus sharing essential features of generic weakly coupled nonequilibrium systems. By analytically solving for the retarded correlation functions, we clarify how branch cuts arise generically from noncollective particle excitations, how they interface with poles arising from collective hydrodynamic excitations, and to what extent the appearance of poles remains at best an ambiguous signature for the onset of fluid dynamic behavior. We observe that processes that are slower than the hydrodynamic relaxation timescale can make a system that has already reached fluid dynamic behavior to fall out of hydrodynamics at late times. In addition, the analytical control over this model allows us to explicitly demonstrate how the hydrodynamic gradient expansion of the correlation functions can be Borel resummed such that the full nonperturbative information is recovered using perturbative input only.

I. INTRODUCTION
A broad range of physical phenomena is involved in how relativistic nonequilibrium systems reach thermal equilibrium. For near-equilibrium systems, these mechanisms are expected to leave characteristic traces in the analytic structure of the retarded correlation function of conserved quantities G R (ω, k). On the one hand, the prototypic longtime behavior of the correlation functions that describes collective excitations evolving towards global equilibrium is given by hydrodynamic poles, whose locations and residues are dictated by the fluid dynamical gradient expansion. On the other hand the question at which time scales hydrodynamic behavior emerges, and with which confounding mechanisms it may compete, is related to the existence and properties of other nonanalytic structures in the lower complex half plane of the correlators. These nonhydrodynamic modes have been seen to govern the approach to hydrodynamics -or hydrodynamization -not only in static but also in rapidly evolving backgrounds, used in the phenomenological description of heavy-ion collisions [1][2][3][4][5][6]. While much of the recent work on nonhydrodynamic modes has focused on strongly coupled theories [7][8][9][10], the present study will deal with nonhydrodynamic modes in weakly coupled theories.
Additional motivations for studying nonhydrodynamic modes in relativistic equilibrating systems come from the apparent phenomenological need to understand how fluid dynamical behaviour arises in nucleus-nucleus, nucleusnucleon and possibly proton-proton collisions [11,12]. Some phenomenologically successful descriptions of these systems interface hydrodynamics with transport models (see e.g. [13,14]), while others do not invoke hydrodynamics explicitly (see e.g. [15]). This asks for a better understanding of where and how kinetic theory differs from hydrodynamics. The standard way of relating kinetic theory to viscous hydrodynamics is to derive the latter by truncating the former to a finite set of moments of the distribution function [16,17]. However, this truncation is based on the assumption that hydrodynamics works. To understand whether, when, and how it breaks down necessitates investigating kinetic theory beyond the moment expansion. The purpose of the present manuscript is to do so by studying how small deviations from thermal equilibrium relax in a full kinetic theory framework.

Analytic structure at strong and weak coupling
In known examples of strongly coupled systems at large N c , the remarkable simplicity of the microscopic structures of nonabelian plasmas is reflected in a remarkably simple analytic structure of the full field theoretic correlation functions. More specifically, in N = 4 SYM theory in the limit of large number of colors N c → ∞ and strong coupling λ = g 2 N c → ∞, the retarded correlation functions are known to exhibit an infinite set of nonhydrodynamical poles located (asymptotically for large n) at ω ± n = ω ± 0 ± 2πnT (1 ∓ i), with n ∈ [1, 2, 3, . . .], and ω ± 0 /πT = ±1.2139 − 0.7775i [7][8][9]. In addition, in the channels where energy momentum conservation demands, the correlation functions exhibit poles whose locations and residues are dictated for small k by the hydrodynamic gradient expansion.
In weakly coupled theories, the analytic structure of retarded correlation functions is much richer. In these theories, there is a scale separation between the typical size of the wave packets 1/T and the mean free path between the individual scatterings t scat . Therefore for time separations larger than ∆t 1/T , when interference effects can be neglected, the correlation function is determined by Boltzmann transport theory, in which the collision kernels are given by in-medium scattering processes in the field theory [18][19][20][21]. The nonanalytic features of the full field theory that are absent in the transport theory are well known (see Sec. II A). However, the nonanalytic structures appearing in the transport theory are less well understood, and will be the topic of this contribu-tion. As transport theory has a wider regime of validity than hydrodynamics but encompasses it, understanding these structures provides a technically controlled in-road to understanding the onset of hydrodynamic behaviour in weakly coupled theories.

Kinetic theory in the relaxation time approximation
While there have been numerous numerical studies of the full collision kernel in nonabelian gauge theories [22][23][24][25][26][27][28][29], including computations of equilibrium and nonequilibrium retarded correlation functions [30], the question of analytical structures has been addressed only recently [31] in the simplest possible model of the collision kernel -that of simple relaxation time τ R . In this relaxation time approximation (RTA), an ostensibly crisp and simple picture of the onset of fluid dynamic behaviour appears by a migration of a hydrodynamic pole through a nonhydrodynamic cut for a specific value of Knudsen number K = k τ R where k is the wave number of the perturbation [31]. However, this simple model forgoes much of the structures of the collision kernel in favour of a single relaxation time. The question of whether this simple picture survives the inclusion of more realistic collision processes is the starting point of this paper.
The full weak coupling dynamics contains nonhydrodynamic excitations at different energy scales that relax at widely different time scales. A minimal way of incorporating this generic qualitative feature while maintaining an analytically tractable model is to extend the standard RTA to a model with a momentum dependent relaxation time For a power law form of the relaxation time such a model has been used before to gain insight into freeze-out dynamics [32]. By including the scale dependence of τ R (p), we supplement the standard RTA approximation with features that are known to exists in QCD and other field theories of nonabelian plasmas. In particular, for extreme out-ofequilibrium perturbations, a.k.a. jets, the relaxation is related to the famous jet stopping time [33,34] corresponding to the value ξ = 1/2 in our model. Moreover, this generalized model shares features of bottomup thermalization [35] in the sense that decaying particles will heat up the thermal bath locally (see discussion at Sec. III). Both features appear generically for ξ > 0 while they are not realized in the exceptional case 1: Analytic structure of the retarded energy momentum correlation function in the shear channel G 0x,0x (ω, k) in the complex frequency plane ω for the kinetic theory (1). The parts of the cut marked with red crosses correspond to medium constituent particles with lifetimes longer than the hydrodynamical decay time and will eventually dominate the correlation function at late times. The upper complex half plane is analytic by causality whereas for |Re ω| > k the correlation function is analytic by locality of the scattering kernel. The nonanalytic features of the function are confined to the grey area. ξ = 0. Other characteristic features of QCD thermalization processes are not realized in the simple model (1). For instance, according to (1), hard particles decay directly to the thermal bath while this process proceeds in full QCD via a cascade of intermediate quasi-democratic splittings [35,40]. Therefore, we cannot exclude that additional analytical structures of retarded correlations functions might arise in full QCD that cannot be illustrated in an analysis of (1). However, as the analytic structures established in this manuscript for the model (1) arise from generic features of kinetic theory, we expect them to be realized in more complete descriptions, too.
The main result of the present paper is to establish the analytic structure of the retarded correlators of the energy-momentum tensor for the model (1). This result is sketched in Figure 1 for the (analytically continued) shear channel correlation function obtained from the model. Causality and the stability of thermal equilibrium make the correlation function analytic in the upper complex half-plane, while the locality of the collision kernels in the Boltzmann equation makes the correlation function analytic for |Re ω| > k. In addition to the hydrodynamic pole, the model exhibits two nonhydrodynamic cuts whose branch points are located at ω = ±k. For any k, the cuts extend to smaller imaginary parts than the hydrodynamic pole; it is these structures that are responsible for a nontrivial competition between hydrodynamics and nonhydrodynamic modes that we discuss in detail. The paper is organized as follows: In section II, we first provide simple qualitative arguments for the physical mechanisms and corresponding analytic structures arising in full gauge theories. For the class of models (1), section III derives then explicit expressions for the retarded correlation functions. For the case ξ = 1, these correlation functions can be expressed in terms of one single, analytically known generating function H that largely determines the analytic structure of the correlation functions. A detailed discussion of this analytic structure, its physical meaning, and its ambiguities is the focus of section IV, before we turn in section V to a discussion of the physical response on pre-hydrodynamic, hydrodynamic and post-hydrodynamic time scales. As our study provides explicit analytic control over a model of significant physical complexity, it is also an interesting scholarly playground for understanding how Borel resummation techniques can be applied to the asymptotic hydrodynamic gradient expansion. This will be discussed in section VI, before we conclude with a short summary of main results and open questions.

II. GENERIC ANALYTIC PROPERTIES OF RETARDED CORRELATORS AND THEIR PHYSICAL ORIGIN
Before analyzing in detail the model (1) in subsequent sections, we discuss here generic features of the analytic structure of retarded correlation functions of the energy momentum tensor. In particular, we aim at providing physical intuition for the features appearing in kinetic theory.
A. Analyticity properties of retarded correlation functions in gauge theories at finite temperature At weak coupling the analytic structure of the retarded correlation function for ω 1/t scat ∼ g 4 T has been discussed in the context of theories with different field content as well as in terms of different operators [9,18,[36][37][38][39]. Quite generally, the two point function of composite operators constructed from two field operators (such as T µν or the electromagnetic current J µ of a charged field) is given to leading order by the simple one loop diagram depicted in Fig. 2. In the time domain, this diagram is of the generic form where D R stands for the retarded propagator, D rr = D > + D < is the symmetric one, and E p = p 2 + m 2 denotes the energy associated with an excitation of momentum p. The vertices combine to a function V which depends on the theory and the particular channel studied and is a function of momenta p and k. For specific cases, see [9,36] for gauge theories. The correlator (4) can be decomposed naturally into two parts that contain slowly oscillating modes of frequencies ω = E p − E p+q , and rapidly oscillating modes of frequencies ω > E p + E p+q , respectively, 1. Rapidly oscillating part D(t, k) The Fourier transform of the rapidly oscillating part has a cut that extends from m + √ k 2 + m 2 < ±ω < ∞. It will be a recurrent theme in this paper that the analytic structure of retarded correlation functions is ambiguous in the sense that different analytic continuations in the complex frequency plane can account for the same physical response in the time domain. In the present case, this can be illustrated by inserting for the massless theory the Matsubara representation n(p) + 1 2 = ∞ n=−∞ βp (2πn) 2 +(βp) 2 into (8) and integrating over p. It can be seen that by choosing a suitable analytic continuation of D in the lower complex half-plane, the cuts m + √ k 2 + m 2 < ±ω < ∞ along the real axis can be exchanged into a series of cuts that are positioned deep in the negative imaginary region at (for m = 0) Imω = −4πn T and −k < Re ω < k with n ∈ [1, 2, . . .], see figure 3 of Ref. [9]. As the nonanalytic structures in D have a distance O (T ) from the real ω-axis, the contribution D decays on timescale 1/T , and it is insignificant at late times when fluid dynamic behaviour is expected to take place.
2. The slowly oscillating part C(t, k) and kinetic theory As argued in [18], the slowly oscillating part C arises from contributions that can be written in terms of expectation values of number operators. This suggests that for small k, the physics contained in C(t, k) can be captured by kinetic theory. In Fourier space, the slowly oscillating nature of C(ω, k) is reflected in a branch cut that extends along the real axis over the limited range −k < ω < k (for all masses). For small k, this expression can be expanded to give where v = ∂ p E p is the group velocity and the term in square brackets is the ballistic propagator of a free streaming point particle. We shall encounter the same branch-cut −k < ω < k and the same integral (10) when we discuss the free kinetic theory in section II B. The free theory calculation recalled here and presented, e.g., in [9] is insufficient for ω ∼ 1/t scat ∼ g 4 T , where interactions change the dynamics qualitatively. It therefore does not reveal the hydrodynamic pole which is close to the origin at ω ∼ g 4 T . To obtain even at leading order complete results in this region, a class of ladder diagrams needs to be resummed [19]. Such resummation can be dressed in the language of an effective kinetic theory [20,21] of nearly massless quasiparticles, where the resummed diagrams appear in the particular scattering kernels of the kinetic equation. The effective kinetic theory is suitable for the computation of correlation functions of the quantum field theory with external momenta ω, k 1/t scat , and therefore it is suitable for studying the vicinity of the slowly oscillating cut of C in more detail than the unresummed calculation. However, this resummation fails for larger (negative imaginary) values of ω and does not capture the physics of cuts of D.

B. Analytic structure of retarded correlation functions in kinetic theory
In this subsection, we develop an intuitive understanding for the analytic structures accessible via kinetic theory.

Massless kinetic theory without interaction
As sketched on the left hand side of Fig. 3, a sound channel perturbation in an equilibrium system may be viewed as embedding alternating sheets of overdense and underdense regions that are separated in the z-direction by a distance 2π/k. Analogous sketches can be given for perturbations in other channels. Computing the retarded response at time t amounts then to studying the state of the system at some arbitrary point x which initially is on the peak of the overdense region at t = 0 when the perturbation is introduced.
In a massless kinetic theory without interactions, particles move on straight lines at the speed of light. What determines the state at the point x at time t is then the average over a sphere of radius ct. As the overdense regions are spaced 2π/k apart, the particles moving in -z direction will give rise to a signal oscillating with frequency ω = k. This corresponds to a pole at k in the complex ω plane. Particles coming from any other direction with velocity v will result in an oscillating signal with smaller frequency ω = k · v, corresponding to a pole at k · v in the complex ω plane. Integrating over all orientations v from which particles reach the point x, one finds a string of poles between −k < ω < k that assemble to a logarithmic cut This cut is also well known in the physics of hard thermal loops, where it gives rise to Landau damping [41]. We conclude that the simple picture of a homogeneous and isotropic free-streaming dynamics explains the logarithmic branch cut found in interaction-free massless kinetic theory for retarded correlation functions like, e.g., the correlation function in the sound channel calculated in [31] G 00,00 2. Massless kinetic theory in the standard RTA Romatschke [31] has studied the effect of adding interactions to the free kinetic theory in a simplified model of momentum-independent relaxation time approximation with collision kernel where f eq is the local equilibrium distribution function to which the system wants to relax, determined by energy and momentum conservation. The inclusion of these interactions has two qualitative effects. First, trivially, the free particle propagator will be damped at length scales of ∆x ∼ t R , shifting the cut into the negative complex plane by an amount of −i/t R dΩ 4π A more subtle effect arises as a consequence of energymomentum conservation (see eq. (23) for technical details). As the energy and momentum from the lost particles need to go somewhere, a new collective excitation is dynamically created in channels where the conservation demands it (sound G 00,00 and shear G 0x,0x ). For small k, the location and residues of these poles are dictated by the hydrodynamic gradient expansion. We will call this pole in the following hydrodynamic pole. For k ≥ π/2t R , the pole crosses the cut and enters the next Riemann sheet, thus disappearing from the physical plane. Therefore, the model has two distinct kinematic regimes: one where the pole is above the cut and the late time behaviour of the system is dictated by the hydrodynamic pole, and the other where the cut dominates the dynamics at all times. This was called the hydrodynamic onset transition in [31].

Massless kinetic theory with scale-dependent RTA
How does the analytic structure of the retarded correlator indicate that the kinetic theory of a free-streaming gas has been supplemented with the scale-dependent relaxation dynamics of (1)? In close analogy to the angular integrals (11) and (14), we expect that qualitative properties of the analytic structure of retarded correlation functions are captured in this case by the integral .

Hydrodynamic pole
Non-hydrodynamic cut This indicates that relaxing the assumption of a single relaxation time will render the correlation function nonanalytic in the entire strip −k < Re ω < k, Im ω < 0, where poles at different Re ω correspond to different angles of the particles, and different Im ω correspond to different p. We shall establish this picture in an explicit calculation in section IV. It implies that the hydrodynamic pole is always embedded in the nonanalytic structure. The existence of a clear onset transition of hydrodynamics is therefore a spurious feature of the simple assumption of a single relaxation time in (13). As we discuss in the next subsection, emersing the hydrodynamic pole in a nonanalytic strip results in a subtle interplay between hydrodynamic and nonhydrodynamic modes that can lead to a qualitatively novel phenomenon in the long-time behavior.

C. Dehydrodynamization in kinetic theory
With the simple extension to a scale-dependent relaxation time, the notion of a unique Knudsen number is obscured, as for any arbitrarily small wavenumber k, physics of different energy scales enters the transport on different time scales. To illustrate this parametrically, consider a generic small deformation of the thermal equilibrium. As by assumption the deformation does not take the system far from equilibrium, the number of perturbed modes will be, for large p, proportional to e −βp . Each of these modes will then evolve toward equilibrium in a timescale τ R (p), such that the overall magnitude of the nonhydrodynamic part of the perturbation can be estimated at time t by For a given t, the integral is dominated by the decay of modes at a characteristic scale and the perturbation has then an overall magnitude of In channels where conservation laws so demand, the deformation may also excite modes which relax on hydrodynamic time scales where D ∼ t R is the appropriate diffusion coefficient in the channel in question. For ξ = 0, corresponding to a single relaxation time, both contributions turn out to be exponentials and the origin of the well defined hydrodynamization scale discussed in the previous subsection is related to the question which contribution decays faster. However, for general ξ the situation is obviously more intricate. Amusingly, for ξ > 0, the contribution arising from the nonhydrodynamic sector is subexponential, and dominates the signal at late times t t out so that one expects that at some late time a system that was hydrodynamic, will again lose its universal fluid dynamical description and be again described by specific microscopic physics related to the dynamics of the nonhydrodynamic modes. This dehydrodynamization mechanism will be seen at work in the model (1) of scale-dependent relaxation time, where hard particles still decay directly to a thermal bath and hydrodynamic fluctuations of the thermal bath are ignored. In the full QCD collision kernel, however, the same process proceeds via a cascade of intermediate quasi-democratic splittings [27,28,35,40]. Also, due to the fluctuation-dissipation theorem, there are other sources of hydrodynamic perturbations and long-time hydrodynamic tales [18,37,42,43]. Therefore, while the mechanism discussed here is part of full QCD, it may not dominate the late-time behavior of the full theory.

III. THE MODEL: MOMENTUM DEPENDENT RELAXATION TIME
We consider a kinetic theory of the form (1), coupled to an external force F α The particle distribution f ( x, p, t) fulfils the massless onshell condition p α p α = 0, and it is taken to be a function of spatial momenta only, such that the partial derivative ∇ (p) 0 f ≡ 0. We write p = p 0 = | p|, and our metric convention is mostly plus η µν = diag(−1, 1, 1, 1). τ R (p α u α ) is the momentum dependent relaxation time defined in (2). The local target equilibrium distribution function depends on four macroscopic variables, the inverse temperature β = 1/T and the flow field u, with u a u a = −1 that need to be adjusted locally such that the time evolution conserves energy and momentum locally in the absence of the external force F α = 0. According to (21), the condition ∂ µ T µν = 0 implies For the case of a scale-independent relaxation time approximation when ξ = 0, eq. (23) implies that the target thermal system has the same local energy density as the perturbed system. In contrast, for ξ = 1 when τ R (p α u α ) = t R p α u α /T , it is the particle number density that is the same in both systems. For the case ξ = 1/2 it is something in between. Therefore, for ξ > 0, the evolution of the perturbed system to the local target equilibrium will increase the energy density of the local target equilibrium system, i.e., it will heat it up. It is in this sense that the model displays features of bottom-up thermalization for ξ > 0.

A. Solution for linear perturbations induced by an external source
The application of an external force F α reshuffles energy and momentum such that, at a given point in space, the local target thermal distribution f eq is no longer the global equilibrium distribution f g eq but rather the local thermal distribution given by the local energy and momentum densities, f eq = f g eq + δf eq . Here δf eq accounts for the change of the target local equilibrium distribution due to the external force. For a Maxwell distribution -relevant for the high-momentum particles that we are concentrating on -the δf eq can be written as a local perturbation of the global distribution with v ≡ p/p.
In the presence of a small external force F α , the evolution of linear perturbations δf on top of the global thermal equilibrium f = f g eq + δf can be expressed by formulating eq. (21) in Fourier space Our convention for the Fourier transform is Q(ω, k) = dtd 3 ke iωt−i k· x Q(t, x). In eq.(25), we have used the relation f − f eq = δf − δf eq . This relation implies also that up to linear perturbations, eq. (23) translates into constraints for four particular integral moments of δf eq and δf , namely As a consequence, both sides of eq. (25) depend on δf , the left hand side explicitly and the right hand side implicitly through δf eq . The rewritten condition (26) for energymomentum conservation makes this implicit dependence manifest. The task is to solve the four equations (26) selfconsistently for the four local perturbations of the target temperature δT ( x, t) and target flow fields δ u( x, t) that define δf eq . This is done by inserting (25) into (26), thus finding a closed set of four equations for the four variations δT and δu i . The solution of this set of equations is where the integral moments and sources are defined by (27)- (30) for the perturbations of the local target temperature and flow velocity fully define the deviation δf eq of the local target equilibrium distribution from the global equilibrium distribution. This allows one to write explicit expressions for all terms on the right hand side of eq. (25). Therefore, in terms of these solutions, eq. (25) contains now the full microscopic information of the system.
We note as an aside that the following discussion could be easily extended to the case of Bose (or Fermi) statistics, replacing (22) by the corresponding sum over exponentials In particular, the integral moments (31) can be simply calculated for this statistics, resulting in

B. Retarded correlation functions
We follow the standard procedure of sourcing the departure of the energy-momentum tensor from equilibrium, by a perturbation of the metric g µν = η µν + h µν . This amounts to applying an external force where the Γ i αβ denote Christoffel symbols. The retarded correlation functions G µν,αβ R define then the response of the energy momentum tensor to the metric perturbation, and they can be evaluated in terms of functional derivatives The disturbance δT µν of the energy momentum tensor is given explicitly in terms of equation (25), with δf eq defined in terms of eqs. (24) and (27)- (30). The evaluation of the functional derivative δT µν /δh αβ is then straightforward and one finds These retarded correlators describe the response in the spin 2 tensor channel (39) induced by h xy , in the spin 1 shear channel (40) induced by h 0x (or h 0y , h xz , h yz ) and in the spin 0 sound channel (41) induced by h zz (or h 00 , h 03 , h xx , h yy ), respectively. The remaining components of the correlation functions can be obtained from relations imposed by energymomentum conservation, such as We have explicitly checked (up to high orders in the gradient expansion) that the various correlation functions satisfy these nontrivial Ward identities that are not apparent in the above calculation. We have also checked explicitly that for the special case of a momentum-independent relaxation time, ξ = 0, the retarded correlation functions (39), (40), and (41) reduce to the results of Ref. [31]. C. The fluid dynamic limit of GR Up to second order in the gradient expansion in small ω and k, the form of retarded correlation functions is dictated by second order fluid dynamics, namely where dots indicate terms of higher power in k or ω. These fluid dynamic expressions depend on entropy s, temperature T , sound velocity c 2 s , as well as shear viscosity η, the shear viscous relaxation time τ π and the second order transport coefficient κ. To determine these fluid dynamic parameters for the kinetic theory with scaledependent relaxation time, we want to compare the gradient expansion of (39), (40) and (41) to the hydrodynamic expressions (42), (43) and (44). To this end, we expand the integrand of the integral moments (31) to arbitrary order N in ω and k, and we perform the p-integration for each term in this expansion. This leads to where 3F2 is the regularized generalized hypergeometric function. If one of the upper indices of the hypergeometric function is zero or negative integer, the sum truncates to a hypergeometric polynomial, which is the case here when R = 0. For example which, modulo prefactors, determines the gradient expansion of the tensor channel G xy,xy R in (39). We note that this is an asymptotic series. Comparing these gradient expansions to the hydrodynamic limits, one finds see also Ref. [32]. Given that the retarded correlators (39), (40) and (41) are those of a kinetic theory of massless particles, the speed of sound takes of course the value expected for a conformal theory. The expressions for shear viscosity η, the shear viscous relaxation time τ π , and κ are genuine kinetic theory results. As the present evaluation is based on a linearized response to perturbations, it is not sufficient to determine those second order transport coefficients λ 1 , λ 2 , λ 3 which depend nonlinearly on perturbations.
Hydrodynamic poles arise as a consequence of energy momentum conservation. In the kinetic theory calculation of section III B, the structures in the retarded correlators that arise from energy-momentum conservation are related to the term δf eq on the right hand side of (25). Inserting the disturbance (25) into (35) and performing the functional derivative δT µν (ω, k)/δh αβ (ω, k), one finds that it is exactly the nontrivial denominators in (40) and (41) that arise from the terms proportional to δf eq . The hydrodynamic poles in (40) and (41) are therefore given by the zeroes of the nontrivial denominators in these two channels that arise from energy momentum conservation.
To make the pole structure of the fluid dynamic limit of retarded correlation functions more explicit, one can write the fluid dynamic limit of the shear and sound channels as and These expressions agree up to second order in gradient expansion with (43) and (44), respectively. Higher orders in the gradient expansion of the full retarded correlators cannot be expected to be reproduced correctly by (51) and (52). In this sense, the precise location of fluid dynamical poles is beyond the scope of a second order gradient expansion. We shall discuss it in section V without taking recourse to the gradient expansion.

IV. ANALYTIC STRUCTURE OF THE RETARDED CORRELATION FUNCTION IN MOMENTUM DEPENDENT RELAXATION TIME APPROXIMATION
The full retarded correlation functions are defined in terms of the integral moments I a,b,c (ω, k) . To study these correlation functions beyond the simple gradient expansion, one needs to evaluate I a,b,c (ω, k) for nonzero ω and k. A numerical evaluation of I a,b,c (ω, k) in (32) is possible for arbitrary momentum dependencies of the relaxation time approximation (2), i.e., for arbitrary ξ. However, analytical control is advantageous for studying the analytic structure. We therefore focus in the following sections on the case ξ = 1 for which explicit analytical results can be obtained. However, we expect that the qualitative features found for the case ξ = 1 extend to the generic case ξ > 0.
The simplification in the case ξ = 1 arises from the fact that all integral moments can be related explicitly to a single generating function Here, R a,b,c 1 and R a,b,c 2 are simple rational functions of ω and k, and in R 2 the derivative ∂ ρ appears only in the numerator of the rational function. The generating function reads Appendix A provides details of this reduction. According to the procedures presented there, a symbolic computation program for algebraic reduction [45] can be employed to obtain explicit expressions for the rational functions R a,b,c 1 and R a,b,c 2 that enter all moments I a,b,c of interest. We note as an aside that we have attempted to derive expressions similar to (53) for other values of ξ. For other rational values, such as ξ = 1/2, ξ = 1/3 etc, one finds typically expressions in terms of more than one generating function, but we were not able to bring all of them into closed analytical form.
To evaluate H(ω,k, ρ), we start from the representation in terms of the function Here, the integration contour crosses the pole whenω takes negative imaginary values. The incomplete gamma function Γ 0, iρ ω therefore has a logarithmic branch cut for negative imaginaryω.
A. Analytic structure of the generating function H The analytic structure of the retarded correlation functions is determined by the analytic structure of the integral moments I a,b,c which in turn is mainly determined by the analytic structure of the generating function H. We therefore discuss now the properties of H in detail. To perform the integral in (57), we note that ρG(ω, ρ) is a function ofω/ρ only. The derivative with respect to ρ can therefore be replaced by a derivative with respect to x, When integrating this total derivative, one needs to note that forω-values with real part in the range −k < Reω < k, the x-integration crosses between x = Reω− and x = Reω + the branch cut of Γ 0, iρ ω−x for all valuesω with Im (ω) < 0. The corresponding discrete contribution to the integral is proportional to Integrating the total derivative in (59) therefore yields The analytic structure of the full retarded correlation functions inherits the analytic structure of the generating function H in the sense that where the generating function is nonanalytic, so is the full correlation function. The nonanalytic structures seen in eq. (61) can therefore be related to some of the nonanalytic structures sketched for the retarded correlation function in the introductory Fig. 1. In particular, in the first line of eq. (61), the two terms ∝ ω +k G(ω +k, ρ) and ∝ ω −k G(ω −k, ρ) have a logarithmic branch cut for negative imaginary values ofω +k andω −k, respectively. This corresponds to the two nonhydrodynamic cuts depicted in Fig. 1. Moreover, the term in the second line of eq. (61) is nonanalytic in the entire strip Imω < 0 and −k < Reω <k due to the explicit appearance of Imω. This corresponds to the grey-shaded area of nonanalyticity in Fig. 1. We note that this nonanalytic contribution becomes nonperturbatively small for smallω due to the factor ∼ e 1/Imω in G(ω, ρ = 1). Therefore, the analytic region at Imω ≥ 0 is reached very smoothly, whereas the generating function is discontinuous when crossing the (Reω) 2 =k 2 lines. In contrast to poles and branch-cuts, the analyticity in this strip is also mild in the sense that a contour integral around a region of area A is proportional to A. The converse of the above statement is not true: the full retarded correlation functions can show additional nonanalytic features that are not visible in the generating function H. There are singular points, arising from the zeroes of the denominators of eqs. (40) and (41). These special points are embedded in the strip of mild nonanalyticity, but they give rise to pole-like structures in the sense that they give a finite contribution even when A goes to zero, provided that the special point lies within A (see Fig. 1). Some of these correspond to the hydrodynamical modes in the model. Indeed, the location of such a special point, in shear channel for example, is given for small k by as expected from hydrodynamical gradient expansion. We note that for this result, as for any expression derived in a gradient expansion, the nonanalytic parts of the generating function H cannot contribute because of the nonperturbative suppression factor.

B. Ambiguities in the analytical structure
To obtain correlation functions in the time domain, an inverse Fourier transformation needs to be taken This expression is typically evaluated by completing the contour of the ω-integration along a path at negative complex infinity, and writing the result as the sum of contour integrals around the nonanalytic structures in the negative complex half plane. In the present case, however, this standard strategy seems difficult to follow as instead of simple cuts and poles, the generating function H in eq.(62) and, a fortiori, the retarded correlation functions are nonanalytic in an entire two-dimensional region as sketched in Fig. 1.

The analytically continued generating function Ha
A better strategy for calculating (63), that is more practical and more physically revealing is to note that the correlation function is analytic in the upper complex half-plane and along the contour of integration in eq. (63). Therefore, for the purposes of calculating measurable quantities like (63), we may replace the correlation function in the lower complex half-plane with the analytic continuation of the function from the upper complex half-plane. The nonanalytic structure of the correlation functions G αβ,γδ (ω, k) and their generating function H in the lower complex half-plane are thus ambiguous to the extent to which the nonanalytic structures arising in H can be substituted by an analytic continuation from the upper half-plane.
As the nonanalytic part of the function has already been separated in eq. (61), an analytic continuation of H from the upper half plane is found simply by removing the nonanalytic part form eq. (61), Here, the subscript a stands for analytic continuation.
The function H a contains incomplete gamma functions with logarithmic branch cuts whose paths are arbitrary as long as their endpoints are fixed toω = ±k and to negative complex infinity. Here, we adopt the simplest, but ambiguous choice of continuing the complex gamma function to the full complex plane, resulting in branch cuts atω = ±k + iȳ, for realȳ ≤ 0. So, H a shows the nonhydrodynamic cuts depicted in Fig. 1, but unlike H, these cuts do not bracket a two-dimensional strip of mild nonanalyticity. To visualize how the analytic structure of H a shapes that of retarded correlation functions, we plot in Fig. 5 the real and imaginary part of the shear channel G 0x,0x R , calculated from H a . This correlation function clearly shares with H a the two branch cuts that run in the negative imaginary half plane along Re (ω) = ±k from zero to complex negative infinity. Closer inspection also reveals that the discontinuity across these branch cuts is expo-nentially small for small Im (ω), as expected from the factor exp [iρ/ω] in (58). In addition, there is a prominently visible structure of neigboring peak and trough close to Re (ω) = 0 at negative Im (ω), whose orientation is rotated by π/2 between the real and imaginary part of G 0x,0x R . This is the tell-tale signature of a simple pole ∝ 1/(ω + i const) in the complex plane. The precise location of this hydrodynamic pole will be discussed in the following. In the gradient expansion, it is given of course by (62).

Deforming the branch cuts
The purpose of this section is to show that in general, the presence or absence of hydrodynamic poles in the lower imaginary half plane of G αβ,γδ (ω, k) is not indicative of the onset or disappearance of fluid dynamic behavior.
To set the stage of this discussion, we note first that the same physical response G αβ,γδ R (t, k) in the time domain can be encoded in different analytical structures G αβ,γδ R (ω, k) in the complex frequency domain. This was illustrated already by showing that constructing G αβ,γδ (ω, k) from the generating function H in eq. (61) or from H a in eq. (64) yields physically identical responses G αβ,γδ R (t, k) while the analytic structure of G αβ,γδ (ω, k) is qualitatively different for both cases in the sense that it has a two-dimensional region of mild nonanalyticity if constructed from H, but not if constructed from H a . In the present section, we consider formulations of the latter kind, for which G αβ,γδ R (ω, k) is given in terms of branch cuts and poles only. In particular, the construction of G αβ,γδ R (t, k) from the generating function H a is technically advantageous, since the contour of the integration (63) can be closed by encircling the branch cuts going from ±k to ±k − i∞ and encircling any hydrodynamical poles ω i that may be found in the given channel, As we shall illustrate in the following with an explicit construction, only the sum of the pole and cut contributions on the right hand side of (65) is physical. The relative weight of both terms depends on the orientation of the branch cuts in the lower complex half plane, which is a purely technical choice without unambiguous physical interpretation.
If calculated from H a , all integral moments entering the correlators G αβ,γδ R (ω, k) can be expressed in terms of rational functions and rational functions times G + = arithmic branch cuts of Γ 0, ī ω±k . To be specific, we consider now a particular deformation of these branch cuts, sketched in Fig. 6 and defined by the replacement Here, R reg denotes the regular part of the Γ-function. In the replacement (66), the logarithm in the second line of (66) cancels part of the branch cut of the Γ-function, such that only the segment a) in Fig. 6 remains. The two logarithms in the third line of (66) combine to the segment b) in Fig. 6, and the logarithm in the last line corresponds to segment c). We deform the branch cut of G − symmetrically (see Fig. 6), so that both branch cuts meet atω = −iσ on the imaginary axis, and are then continued on top of each other up to complex imaginary infinity. This deformation leaves the generating function unchanged for Re (ω) ≥ 0 and it therefore encodes the same physics. In Fig. 7, we plot the real and imaginary parts of the retarded correlation function in the shear channel for this choice of branch cuts. 1 Depending on the depth −iσ in the complexω-plane at which the two branch cuts 1 We note that our construction of these branch cuts in (66) involves pairs of logarithmic cuts that cancel each other outside a finite segment. For instance, the two terms in the third line of (66) extend both toω = −iσ + ∞ but they cancel each other for Re (ω) >k. The numerical evaluation shown in Fig. 7 does not are joined, the shear pole is either clearly visible (left hand side of Fig. 7), or it disappears under the branch cut. We emphasize that while both choices of σ lead to qualitatively different features in the analytical structure of G αβ,γδ R (ω, k), they are physically equivalent in the sense that they give rise to identical physical responses G αβ,γδ R (t, k) in the time domain. In this sense, the appearance or disappearance of a hydrodynamic-like pole is related to purely technical and physically ambiguous choice of branch cut and it therefore cannot be related to the onset of fluid dynamic behavior. 3. Differences between the cases ξ = 0 and ξ > 0 As explained in Appendix A, eq. (A3), the integral moments (31) that define retarded correlation functions for the case of a scale-independent relaxation time, ξ = 0, can be written in terms of rational functions ofω andk, and in terms of rational functions times the difference of logarithms (67) [for the case ξ = 0] This is consistent with the qualitative argument leading to (14). As a consequence, for ξ = 0, the retarded correlation functions share the nonanalytic structure of (67).
According to the standard definition, the branch cuts of the logarithms in (67) start at ω = −i/t R ± k and they run parallel to the real axis to ω = −i/t R −∞. Therefore, they cancel each other outside the range −k ≤ Re ω ≤ k, and this gives rise to the nonanalytic segment sketched in Fig. 4. However, the two logarithmic branch cuts of (67) could also be deformed to run parallel to the imaginary axis from ω = ±k − i/t R to negative complex infinity, ω = ±k − i∞.
These two ways of orienting the branch cuts of (67) are reminiscent of the two choices of branch cuts for H a depicted in Fig. 7 and discussed for ξ = 1 in the previous subsections. However, there are marked physical differences between the cases ξ = 0 and ξ > 0: First, for ξ = 0, the branch cuts can be oriented such that for sufficiently small k, hydrodynamic poles are the unique nonanalytic structure closest to the real axis, thus determining the late-time behavior of retarded correlation functions, see eq.(65). In contrast, for ξ = 1, the branch cuts start always atω = ±k, and for a gradient expansion aroundk = 0, poles and the starting point of branch cuts are not separated. This observation is related to the finding that the gradient expansion for the position of the pole converges for the case ξ = 0 (for instance, ω shear (k)| ξ=0 = −i t R + ik tan(k) [31]), while it is an attribute values to these lines along which logarithm contributions cancel each other, even though the correlation function is regular there. Second, for ξ = 0, the branch cuts in (67) can cancel each other outside a finite segment. As illustrated in Fig. 8, this is not possible for the case ξ = 1. If one deforms the branch cuts of H a so that they lie on top of each other fromω = −iσ up toω = −∞, they will not cancel exactly. Rather, along the line of overlapping branch cuts, there will be a discontinuity where H Right In this section, we utilize our understanding of the nonanalytic structures of G αβ,γδ R (ω, k) in the frequency domain for a discussion of the physical response G αβ,γδ R (t, k) in the time domain. The connection between both is given by eq.(65).
In general, with small but increasing k, the pole contributions to G αβ,γδ R (t, k) in (65) move deeper into the complex plane and they start being cancelled more efficiently by the discontinuities from the branch cuts. While only the sum of these nonanalytic contributions has unambiguous physical meaning, the separate determination of both, the poles and their residues, and the discontinuities along the branch cuts is needed in practice for a discussion of the full physical response in the time domain. In the following, we discuss these nonanalytic contributions separately for the specific choice of the gener- ating function H a in (64) with branch cuts taken alonḡ ω = ±k + i y t R , y ∈ [0, −∞].
A. The location of the hydrodynamic poles in the shear and sound channel

The pole in the shear channel
The poleω shear (k) of the retarded correlation function G 0x,0x (ω,k) is defined implicitly in terms of the zero of the nontrivial denominator in eq. (40), This equation can be solved numerically without any recourse to the gradient expansion. Alternatively, it can be solved by determining the first N coefficients b i in a gradient expansion In Fig. 9, the exact solution is compared with this gradient expansion. With increasing orders ∝k 2N , the gradient expansion is seen to deviate from the exact result at smaller and smallerk. This illustrates that the gradient expansion is an asymptotic expansion that does not possess a finite radius of convergence.
For largek, the hydrodynamic pole moves deep into the complex planē Re(ω) k Im(ω)

The sound channel
In close analogy to the discussion of the pole in the shear channel, the poles in the sound channel can be determined in terms of the zeros of the denominator of (41). While the pole in the shear channel is purely imaginary, the pair of sound poles start at finite real values ω sound (k = 0) = ±c s = ± 1 √ 3 before diving into the negative imaginary half plane. The full numerical solution is shown in Fig. 10.
We note that branch cuts can be chosen such that hydrodynamic poles disappear below the cut in one channel while they do not disappear in another channel. Here, this is the case for the choice of branch cuts in H a along the imaginary axis. For this choice, the shear pole will remain visible for allk, while the sound pole disappears atk = 4, see Fig. 10. This is yet another illustration of the general statement that there is no unambiguous relation between the existence of hydrodynamic poles in the retarded correlator and the persistence of fluid dynamic behavior.
We further observe with curiosity that the positions of the sound poles move first away from the real axis, before they move closer to the real axis again, see Fig. 10. We note that other cases are known in the literature where a pole moves closer to the real axis with increasing k, see e.g. Ref. [44]. The asymptotic large-k behavior is given byω B. Contributions of the branch cuts to G αβ,γδ R (t, k) We now combine the information gathered about the nonanalytic structure of G αβ,γδ R (ω, k) to arrive via eq. (65) at a qualitative understanding of the timedependence of the physical response G αβ,γδ R (t, k). For the shear channel, this time dependence is illustrated with the numerical results in Fig. 11 that display the three characteristic stages of hydrodynamization, hydrodynamic evolution and dehydrodynamization. The following discussion aims at providing an analytic understanding for how these features arise.
For notational simplicity, we work in the following with (73)

The limit t → 0 of the retarded correlation functions
In the kinetic theories studied here, the physical response at time t = 0 starts always from This can be seen by expanding the exponent in the Fourier transform for small t, One checks explicitly for each channel of interest that G αβ,γδ R (ω, k) falls off like ∝ 1/ω 2 or faster for large ω. Therefore, the first term in the expansion (75) can be obtained by closing the integration contour in the positive imaginary half plane where integrals along closed contours vanish due to the analyticity of G αβ,γδ R (ω, k). This implies that the small-t expansion of G αβ,γδ R (t, k) starts with a positive power of t and that eq. (74) is satisfied.
According to eq. (65), the pole contribution to G αβ,γδ R (t, k) at time t = 0 is a sum of residues which is nonzero. To satisfy (74), this pole contribution must therefore cancel exactly the contribution from the cut at time t = 0. To see this cancellation explicitly at work, we consider the shear channel that has one single pole, and we focus for simplicity on large k. In this limit, the residue of the shear pole ofḠ 0x,0x Therefore, the pole contribution to the retarded correlation function (65) diverges for large k and small t. In the same limit, the cut contribution is sharply peaked around the location of the pole DiscḠ 0x,0x such that the contribution from the discontinuity for large k and small t reads So, indeed, cut and pole contribution cancel exactly for t = 0. Since both are continuous in t, they weill cancel partially for short times t > 0.

Hydrodynamization
In applications of hydrodynamics, it is often assumed that hydrodynamic behavior dominates the evolution of near-equilibrium perturbations on time scales t > τ π . In the kinetic model studied here, this hydrodynamic shear relaxation time (49) is τ π = 6t R .
According to eq. (65), the timescale over which the cut contribution dies out exponentially is inversely proportional to the depth y in the complex plane where the discontinuity becomes sizeable. The physics is particularly clear in the limit k → 0, where one is dealing with one single cut and avoids issues related to the partial cancellation between different cut contributions. In this limit, the shear viscous correlation function takes the form (We note as an aside that the first nontrivial order of the shear correlator is ∝ k 2 as a homogeneous shear perturbation corresponds to a boost and does not create shear flow.) In Fig. 12, we have plotted the suitably normalized imaginary part of the discontinuity DiscḠ 0x,0x R (k + iȳ,k)/k 2 | k=0 as a function of negative Im (ω). One finds that this function peaks indeed close to 1/τ π , thus indicating that the cut contribution to the retarded correlation function (65) will be governed initially by an exponential decay time close to τ π .
In summary, simple physics arguments, the numerical inspection of the imaginary part of the cut discontinuity, and the numerical calculation of the retarded correlation function shown in Fig. 11 all indicate that the physical response to perturbations starts being dominated by hydrodynamics on time scales t > τ π = 6t R . We emphasize, however, that it is difficult to make this numerical observation analytically precise. The kinetic theory studied here allows for physics on different momentum scales to relax on different time scales.

Late time limit of the correlation function
The late time behaviour of the correlation function is determined by the nonanalytic structures closest to the real axis which are the cuts running to the real axis at ω = ±k. In the physical responseḠ αβ,γδ R (t,k) in eq. (65), the cut discontinuity DiscḠ αβ,γδ R (k + iȳ,k) at distancē y = y t R from the real axis is weighted with an exponentially suppression eȳ t/t R . For the study of the late time behavior t 1/k and for sufficiently long wavelengths 1/k t R , i.e.k 1, it is therefore sufficient to expand this discontinuity around the "on-shell" pointω =k.
To be specific, let us consider the shear channel correlation functionḠ 0x,0x R where the expansion of the branch cut discontinuity around the on-shell point yields The corresponding contribution to the retarded correlation function (67) in the time domain reads This contribution to the retarded correlation function is clearly nonhydrodynamic. It is an oscillating function with subexponential decay, and it will therefore dominate at late times over any contribution from hydrodynamic poles. Eq. (81) confirms in an explicit calculation for ξ = 1 the parametric estimates obtained for arbitrary ξ in section II B, see eq. (18). To estimate the scale at which this dehydrodynamization takes place, we require that the negative exponent of the pole contribution in (65) is much larger than the nonexponential factor in (81), tIm (−ω i ) 2 t/t R . Since the imaginary part of the fluid dynamic poles ω i starts ∝ k 2 for small k, we therefore conclude that in the scale-dependent relaxation time approximation investigated here, the kinetic theory dehydrodynamizes for arbitrarily small k at sufficiently late times, This dehydrodynamization is visible in an oscillatory subexponential late-time decay of retarded correlation functions, as can be seen in Fig. 11. According to eq.(82), the timescale at which dehydrodynamization occurs varies strongly with the momentum k. While Fig. 11 shows a wide window of close-tohydrodynamic evolution fork = 0.4, this window closes ifk is increased to values larger than unity. As seen in Fig. 13, already fork = 2, the oscillatory late-time behavior is visible at all time-scales and a window of closeto-hydrodynamic behavior does not exist. FIG. 13: Same as Fig. 11, but now for a larger momentum scalek = 2 for which the hydrodynamic pole does not dominate the time evolution on any timescale.

VI. ASYMPTOTIC NATURE OF GRADIENT EXPANSION AND BOREL SUMMABILITY
We discuss now the use of Borel techniques to resum the divergent gradient series of the correlation functions. To simplify the discussion and to arrive at analytical expressions, we consider the shear channel correlation function (79) in the limit of vanishing k as an explicit example. Its hydrodynamical gradient expansion corresponds to a Taylor expansion in ω Comparing the full expression to different orders of this gradient expansion, one sees from Fig. 14 that the expansion in powers of ω is an asymptotic. For a given value of ω, inclusion of higher order terms does not improve the approximation but instead makes it worse. This poor convergence of the series is caused by a factorial growth in the Taylor coefficients and is a consequence of the cut of the Γ-function extending to the expansion pointω = 0. A standard trick for improving the convergence of the series near non-analytic structures is to replace the Tay-lor series by a Padé approximant which as a rational polynomial can account for nonanalytic structures. Indeed, as is evident from Fig. 14, the 25th order Padé approximant (that is, approximating the function with rational polynomial whose numerator and denominator are 25th order polynomials in ω, and whose Taylor expansion coincides with that of the original function up theω 50 ) performs vastly better numerically.
Whereas the non-analytic structure of the correlation function is a cut, the only non-analytic structures present in the Padé approximant are poles. The way the cut is mimicked by the Padé approximant is in term of an alternating string of poles and zeroes where the original cut lies, such that the poles become denser as the order of approximation is increased, see Fig. 15.
In order to gain further improvement, one may try to use Borel's trick of writing factorials in integral represen- The art is then to perform the Borel sum in the integrand of (85) which can be convergent since it has factorially suppressed coefficients. As typically one has information only of finite set of Taylor coefficients b j , the standard practice is to again approximate the Borel transform using a Padé approximant. In our case it turns out that the Borel transform is itself a rational function, and therefore the Padé approximation is exact once a required amount of terms are taken into account Of course, if we had access only to a finite number of Taylor coefficients, we could not know for sure that we have fully reconstructed the Borel transform. But in our case, we may simply compute the inverse transformation of eq. (85), and indeed we recover back the original expression (79). It is remarkable how using the Borel resummation we have been able to recover the non-analytic features of the correlation function with only perturbative information about the gradient series.
It has been suggested that the non-analytic features in the Borel transform arise from physics of nonhydrodynamical modes. In the example at hand, it is easy to see that the essential singularity at the origin arises from the residue of the only pole of the Borel transform We find it curious that the exponent in the previous equation, or the location of the nonanalyticity of the Borel transformation are not directly related to the location of the nonhydrodynamic mode with the smallest imaginary part. This is in contrast to the analogous problem in an expanding background, where the system is driven out of equilibrium because of longitudinal expansion instead of a external metric source. It has been suggested in [1,2,4] that in this case the location of the first nonanalytic structure in the Borel plane is given by the slowest decaying nonhydrodynamic mode. In our case, the nonhydrodynamic mode with the smallest imaginary part has always vanishing imaginary part, and indeed the nonanalytic behaviour arises from the combined effect of all nonhydrodynamic modes.
To contrast this picture with a case where the nonhydrodynamic modes are well separated from the expansion point ofω = 0, consider the correspoding shear channel correlation function at vanishingk in the case of ξ = 0 In this trivial case the gradient expansion is well behaved and the pole located atω = −i sets the radius of convergence. In this case the Borel transformation reads which is a complete function with only an essential singularity at large s. We also note that, from the point of view of Borel summation, the cases xi = 0 and ξ = 1 are both speacial. For ξ = 0, the gradient expansion is convergent series. For 0 < ξ ≤ 1, its Borel sum is convergent while the gradient expansion itself is asymptotic. As seen, e.g., from eq. (46), the coefficients of the gradient expansion for ξ > 1 grow faster than factorial, making also the Borel sum nonconvergent.

VII. CONCLUSIONS
Generically, the path to equilibration in relativistic systems described by Boltzmann transport is governed by an interplay of collective hydrodynamic and noncollective particle excitations. The present study allowed us to expose this interplay in detail. Generically, there is no sharp onset of hydrodynamic behavior. On all time and length scales, both hydrodynamic and nonhydrodynamic modes are present. To which extent the one dominates over the other can be at best a quantitative statement that changes gradually with scale. Also, the appearance of poles in the first (physical) Riemann sheet of retarded correlation functions is a matter of choosing a particular analytical continuation and thus cannot be related unambiguously to the onset of fluid dynamic behavior. Still, even if the pole can be made disappear from the physical Riemann sheet by utilizing the ambiguity in analytic continuation, its weight is translated unambiguously to other non-analytic structures in that sheet. In this sense, the relative closeness of hydrodynamic poles to the real axis carries quantitative information about the onset of hydrodynamic behavior irrespective of whether they are visible. The hydrodynamic behavior is fully characterized by the coefficients of a gradient expansion. As we showed for a generic kinetic theory, this expansion is asymptotic already for retarded correlation functions, since the starting point of the branch cut approaches the origin for small k. This is in marked different to results obtained for strong coupled field theories and in the standard scale-independent relaxation time approximation, where the gradient expansion for retarded correlation functions converges. Remarkably, however, the latter theories if pushed out of equilibrium by longitudinal expansion exhibit a time-dependent energy density whose gradient expanion (in powers of inverse time) is asymptotic. We note that non-linear transport coefficients appear in this expansion, while the above-mentioned gradient expansions of retarded correlation functions involve linear transport coefficients only. It would be interesting to understand the relation between the analytic structures of the higher n-point functions that give rise to non-linear transport coefficients, and the qualitatively different convergence properties of the above-mentioned gradient expansions.
Borel summation is employed in attempts to extract physically meaningful information from non-convergent asymptotic series. This technique is often advocated with the seemingly contradictory claim that it can reveal nonperturbative information from analysis of purely perturbative input. By explicitly resumming the Borel series of the gradient expansion of a retarded correlator, we demonstrated in section VI how this can function. To the best of our knowledge, this is the first time that an explicit Borel transformation of a hydrodynamizing nonequilibrium system has been fully performed.
Our study could be extended on several fronts. The present discussion remained limited to linear response and it could be extend within the present set-up to nonlinear response and, in line with the remarks above, to systems undergoing expansion. It may also be interesting to supplement the kinetic theories studied here with thermal fluctuations that via the fluctuation-dissipation theorem are known to give rise to characteristic long-time hydrodynamical tails. Furthermore, it would be interesting to observe, e.g., in numerical simulations, the features identified here in kinetic theories whose collision kernels are derived directly from quantum field theory. Finally, as mentioned in the introduction, a full quantum field theoretical treatment contains interference effects that go beyond simple kinetic theory and become relevant at higher orders in perturbation theory.
Appendix A: Calculation of the integral moments I abc for ξ = 1 In this appendix, we provide further information on how to evaluate the integral moments I abc of eq. (31), which can be written for ξ = 1 in the form I a,b,c = 1 24 dp p 5−a e −p dφ 2π Here,ω ≡ t R ω,k ≡ t R k and x denotes the cosine of the angle between v and k. The φ-integration leads to trivial prefactors. Only integral moments with even integer index b are non-vanishing. To bring the p-and x-integrations into a simpler form, we proceed as follows: We first observe that for b = c = 0, the elementary xintegral returns a logarithm For arbitrary positive integers b, c, the corresponding integral can be shown to be of the form For the components relevant for our calculation, we have tabulated the functions T b,c 1 (kp,ωp) and T b,c 2 (kp,ωp) in Table I. To obtain the moments I a,b,c in (A1), it then remains to perform the integral I a,b,c (k,ω) = 1 24 dp p 4−a e −p T b,c 1 (kp,ωp) + T b,c 2 (kp,ωp) For all moments that enter the retarded correlation functions (39), (40) and (41), the products p 4−a T b,c 1 (kp,ωp) and p 4−a T b,c 2 (kp,ωp) in the integrand of eq.(A4) are explicitly known polynomials in p that include only positive powers up to p 4 . The first term in (A4) is then easily integrated, using dp p n e −p = Γ[n + 1] . (A5) The second term in (A4) requires calculating for n = 0, 1, 2, 3, 4 the expression 1 24 dp p n e −p −i 2k L = i 24 dp p n e −p where H(ω,k, ρ) is the analytically known generating function defined in (56), or an analytically continued function H a that agrees with H along the real axis and the positive imaginaryω-half plane.
In this way, somewhat lengthy but explicit expressions for all relevant integral moments I a,b,c (ω,k) are obtained by inserting into (A4) the explicit terms given in table I, writing these terms in powers of p, and performing the p-integrals with the help of eqs.(A5) and (A6).