Transseries for causal diffusive systems

The large proper-time behaviour of expanding boost-invariant fluids has provided many crucial insights into quark-gluon plasma dynamics. Here we formulate and explore the late-time behaviour of nonequilibrium dynamics at the level of linearized perturbations of equilibrium, but without any special symmetry assumptions. We introduce a useful quantitative approximation scheme in which hydrodynamic modes appear as perturbative contributions while transients are nonperturbative. In this way, solutions are naturally organized into transseries as they are in the case of boost-invariant flows. We focus our attention on the ubiquitous telegrapher's equation, the simplest example of a causal theory with a hydrodynamic sector. In position space we uncover novel transient contributions as well as Stokes phenomena which change the structure of the transseries based on the spacetime region or the choice of initial data.


Introduction
Understanding nonequilibrium phenomena with hydrodynamic tails has been a very active research direction of the past two decades. One motivation for this quest has been ultrarelativistic heavy-ion collisions at RHIC and LHC and the success of hydrodynamic modelling there [1][2][3][4]. Another set of motivations came from condensed matter and quantummany body physics [5]. There have also been more foundational questions driving this field, such as existence of bounds on transport [6], or the character of the hydrodynamic gradient expansion [7].
One of the driving forces in this endeavour was holography [8][9][10], in which the nonequilibrium dynamics of a certain class of quantum field theories is represented by time-dependent geometries involving black holes [11]. Results in this context are often based on numerical solutions of the equations of motion, with analytic or semi-analytic insights limited to basically two cases: linear response theory around equilibrium or highly symmetric dynamics. The former cases are based on Fourier transform techniques while the latter concern situations which can be effectively captured as comoving flows known from cosmology, or from heavy-ion collisions.
Our paper is directly motivated by such a description of ultrarelativistic heavy-ion collisions in terms of a Bjorken flow -one-dimensional expansion of matter, which looks the same in any coordinate frame boosted in the direction of expansion [12]. In the case of conformal models, such dynamics can be expressed as a single function A(w) which measures deviations away from local equilibrium. w is a dimensionless clock variable and the key motivating point for our paper is that A has an asymptotic late-time form of a transseries, i.e. a double expansion in powers of 1/w and an exponential suppression factor: µ n w −n + e −Ω w w γ ∞ n=0 ν n w −n + . . . . (1.1) The ellipsis in the above equation stands for other exponentially suppressed contributions to the sum. The first sum in (1.1) is the onshell hydrodynamic gradient expansion of the energymomentum tensor, Π ab = −η σ ab + η τ Dσ ab + λ 1 σ a c σ b c + . . . , (1.2b) where we focus on conformal theories and the ellipsis denotes three additional nonlinear terms with two derivatives of velocity U a which are irrelevant for our discussion, as well as higher order contributions, see, for example, [1] for details. In (1.1), the µ 1 /w term represents the shear viscosity η contribution and the µ 2 /w 2 term represents the combined effect of τ and λ 1 [13]. More generally, each term µ n /w n comes from the n th order of the hydrodynamic gradient expansion. The key aspect of [7,[14][15][16][17][17][18][19][20][21][22] was the ability to explicitly calculate higher order contributions to the sum (1.1) with the conclusion that the gradient expansion evaluated on shell for the Bjorken flow diverges, i.e. µ n ∼ n! at sufficiently large n (see, however, [23]).
The second contribution in (1.1) is associated with transient, exponentially decaying phenomena known from linear response theory and dressed in the hydrodynamic variables [24,25], hence the other, also divergent gradient expansion with ν n /w n terms. As is the case in the mainstream application area of resurgence in theoretical high-energy physics [26,27], i.e. coupling constant expansions in interacting quantum mechanical systems, the hydrodynamic sector carries information about transient phenomena through the phenomenon of resurgence.
It is an important open problem to what extent the picture uncovered for Bjorken flow, i.e. divergent hydrodynamic expansion and spatiotemporal dependence as transseries, survives when symmetry assumptions are lifted. New light on this question was shed by our recent article [28], in which we combined linear response theory with the Fourier transform to investigate convergence of hydrodynamic gradient expansion in linearized conformal hydrodynamics in complete generality.
In linear response theory, a component of conserved currents ρ(t, x) acquires spatiotemporal dependence given by where ω q (k) are dispersion relations for different modes in the system. 1 Hydrodynamic modes are those for which ω q (k) = κ k z + . . . , κ ∈ C, (1.4) with Lifshitz exponent z > 0, which guarantees that dissipation accounted for by the imaginary part of κ can be made arbitrary small by giving initial data support at arbitrarily small k. The ellipsis in (1.4) denote terms with higher powers of k, which are a counterpart of the hydrodynamic gradient expansion (1.2). The aforementioned transients are simply contributions to the sum in (1.3) which do not share the property (1.4). The aim of the present article is to build on [28] to construct a transseries solution describing nonequilibrium processes going beyond the class of comoving flows represented by Bjorken dynamics. Our guiding principle will be to have hydrodynamic phenomena captured by the perturbative part of the transseries with nonperturbative transient phenomena captured by higher transseries sectors -in analogy with Bjorken flow, where the perturbative part is given by the 1/w expansion while the transient effects are expressed by terms exponentially suppressed in w. In any given model, there may be many different ways to construct such an expansion, e.g. by treating some microscopic parameter as small. The expansion we use here can be applied to any model.
At the technical level, the key idea is to introduce a formal expansion parameter for the transseries by rescaling space and time coordinates in Minkowski space where a nonequilibrium phenomenon of interest takes place. Since we expect hydrodynamics to be a late time phenomenon, we introduce a formal parameter through the rescaling with α > 0, and treat as small. The resultant effect on (1.4), through the corresponding scaling ω → α ω, k → k, is given by, For a given hydrodynamic sector parameterized by some Lifshitz exponent z, the natural choice is therefore the marginal one, α = z, preserving the hydrodynamic scaling (1.4). This choice focuses attention on the sector of interest, whilst not scaling as far as to render it trivial. It also has the desired effect of ensuring that all nonhydrodynamic modes appear nonperturbatively in a small expansion, since they scale as ω = O( ) −z . The outcome is that the spectral decomposition (1.3) becomes a transseries in the parameter with perturbative sectors corresponding to hydrodynamic mode contributions, and nonperturbative sectors corresponding to nonhydrodynamic ones. We emphasise that is only a formal parameter, and it takes the value = 1 at the end of the calculation. For definiteness, in the present work we focus on nonequilibrium phenomena described by the telegrapher's equation which at large distances and long times describes diffusion, or z = 2 hydrodynamic scaling (1.4). Later, unless we keep explicitly τ and D, we use their following numerical values The linear partial differential equation (1.7) is well-known in the literature and, as we review in appendix A, it arises as the description of shear channel perturbations in the Müller-Israel-Stewart (MIS) formulation of relativistic hydrodynamics [30]. 2 From a broader perspective, the telegrapher's equation features prominently in the context of quasihydrodynamics [33], where it provides the simplest example of a diffusion-to-sound crossover. Quasihydrodynamics is the natural generalization of standard hydrodynamics in the presence of weakly broken symmetries. Apart from MIS itself, examples of theories featuring a diffusion-to-sound crossover described 3 by the telegrapher's equation include quantum fluctuating superconductors [34], systems breaking spatial translations spontaneously in the presence of phase relaxation [35] and, in the AdS/CFT context, probe branes at finite temperature and large baryon density [36][37][38], models of momentum relaxation [39], higherderivative gravity [40,41], and constructions based on generalized global symmetries that describe dynamical electromagnetism in the boundary QFT [42,43] or viscoelastic media [44]. 4 Furthermore, with a straightforward modification, our methods also apply to the chiral magnetic waves in the presence of axial charge relaxation discussed in [46,47]. These observations suggest that the results we will derive in this work are potentially relevant for a wide range of distinct physical systems. 2 In the AdS/CFT context, the natural counterpart of this problem would be a shear channel fluctuation in all-order hydrodynamics as discussed in [31,32]. 3 It is worth remarking that the telegrapher's equation might only emerge in a suitable parametric limit. 4 See also [45] for an embedding of the telegrapher's equation into a field-theoretic context.
After these considerations, let us comment briefly on the mode structure of the telegrapher's equation. We have two modes in the sense of (1.3): Among these two modes, ω H is a hydrodynamic diffusion mode 10) while ω N H remains gapped in the same limit. As figure 1 illustrates, both modes are nonpropagating and purely decaying below a critical momentum , (1.11) which corresponds to a branch point where the hydrodynamic and the nonhydrodynamic modes collide in the sense introduced in [48] and developed in [49][50][51][52][53][54]. For τ and D given by (1.8), k 2 c = 1. Past the critical momentum, the modes acquire a propagating component. For asymptotically large |k|, both modes become purely propagating, with a linear dispersion relation ω → ± D/τ k. This should not come as a surprise, since in the end the telegrapher's equation is nothing but a dissipative wave equation. If one wants to impose relativistic causality with the speed of light set to unity, this requires D ≤ τ .  Another illustrative perspective is to view the mode collision represented by the telegrapher's equation as the simplest incarnation of the so-called k-gap phenomenon, which features widely across physics (see [55] for a review). Apart from the examples mentioned before, the existence of a k-gap in the dispersion relation of the transverse collective excitations is a crucial feature distinguishing liquids and solids [56][57][58]. Other instances of the k-gap phenomenon in the AdS/CFT context include p-wave superfluids [59] or plasmons [60][61][62].
In the special case of the telegrapher's equation, the expansion in powers of can equally well be regarded as an expansion in powers of the relaxation time τ . In consequence of our choice to seek a perturbative scheme in which hydrodynamic modes appear at leading order while nonhydrodynamic degrees of freedom enter as nonperturbative corrections, we observe that we are expanding around the acausal Navier-Stokes limit, and nonperturbative contributions are necessary to ensure causality. Likewise, physical effects of the propagating modes mentioned above enter only at the nonperturbative level. These observations are illustrated further in appendix D, where we apply our perturbative scheme to a particular example.
From a more mathematical perspective, we will be studying solutions of a partial differential equation (1.7) in the form of a transserieŝ where the hats indicate momentum space quantities. One of the crucial points of our work will be to make sense of this transseries expression in position space. Similar analyses of other partial differential equations have been carried out before in various contexts, for example in [63][64][65]. In holography, the large-w expansion of Bjorken flow [7,20,21] is also governed by partial differential equations, the fully nonlinear Einstein equations with negative cosmological constant. 6 One of the interesting conclusions in [63][64][65] is that the structure of the transseries may be different in different parts of space and time. As we will see, this is also the case in our studies. Furthermore, in our study we want to stress the dependence of the transseries on initial conditions. Transseries have also been studied in the context of attractors in Bjorken flow [14]. The non-perturbative contributions quickly decay, leaving a universal perturbative piece identified as the attractor. Similarly, in the present context there are many different solutions that differ only by the behaviour of the transients. For that class of solutions, the perturbative part acts as an attractor in the same sense. 7 However, the Bjorken flow attractor arises from far from equilibrium behaviour involving the fast expansion at early times [66] while in this paper, because of linearity, all dynamics can be understood as small perturbations away from equilibrium. On the other hand, we are able to study more general flows, but leave open the question of what happens when non-linearities are introduced. Note on this front that less symmetric, non-linear flows have been studied numerically in [66][67][68][69].
Finally, let us come back to our initial motivation, i.e. the (conformal) Bjorken flow, and discuss the expansion and (1.12) in relation to w-transseries in (1.1). While the Bjorken flow is a genuinely nonlinear one-dimensional expansion, nothing stops us from applying an -rescaling (1.5) also to this case. Choosing x to be the expansion direction, the boost-invariance forces the dynamics to depend on proper time τ = √ t 2 − x 2 only (not to be confused with the relaxation time τ ). As a result, the natural rescaling preserving 6 When the microscopics is given by kinetic theory models, see [18,19], the situation is even richer, as the collisional kernel in the Boltzmann equation generically involves integration over momenta. 7 The series is generically divergent, so it must be properly resummed or optimally truncated.
the character of the proper time is the homogeneous one, i.e. α = 1, which simply takes τ → τ/ . The clock variable is defined as (1.13) where T (τ) is the local effective temperature associated with the local energy density. The homogeneous expansion forces large proper time expansion of w in powers of 2/3 /τ 2/3 starting with the term 2/3 /τ 2/3 −1 . This reorganizes the transseries (1.1) from a transseries in w, whose each perturbative contribution corresponds to a given order of hydrodynamic gradient expansion evaluated on-shell, to a transseries in 2/3 /τ 2/3 . In the present work we will be working with the analogue of the latter expansion, whereas the former was discussed in full generality in linearized hydrodynamics in our recent paper [28]. In appendix C we discuss the link between the hydrodynamic gradient expansion and the perturbative part of the transseries.
2 The momentum-space transseries and its Fourier transform

Constructing the transseries
In this section we construct transseries solutions to the telegrapher's equation (1.7) in the small formal parameter as defined in (1.5). Our starting point is thus the following -rescaled telegrapher's equation, with (1.7) and its actual solution recovered by setting = 1. The point of introducing is that this is the approach that one can adopt generally, for example in holography or kinetic theory, while the τ -expansion is specific to the telegrapher's equation and related systems. Due to the spatial translational invariance of (2.1), it is convenient to work in momentum space, where (2.1) reduces to a linear ODE, 8 One then immediately recognises an appropriate transseries ansatz forρ(t, k) as → 0 as given by (1.12). Plugging it into (2.2) and considering terms order-by-order in , reveals the following pair of nested ODE systems where we takeâ −1 =b −1 = 0. We see thatâ n+1 (t, k) obeys a heat equation sourced by the previous order, whileb n+1 (t, k) is governed by the time-reversed equation. At this stage, the two equations are decoupled from one another. In order to solve (2.3), we supplement the equations with initial data at t = 0, and denote the corresponding Fourier-transformed functions asû(k),v(k). For the sake of simplicity of presentation in what follows, we are going to focus on the u(x) = 0 case. At the level of the expansion coefficients, this initial condition reduces toâ where the last equality applies only for n > 0. The initial conditions couple the coefficients arising in the perturbative series to those in the nonperturbative series. Finally, it is possible to find closed-form expressions forâ n (t, k) andb n (t, k) that solve (2.3) and obey (2.5). They are given bŷ Let us now turn our attention to position space. Our guiding principle here will be to define the expansion coefficients of a given sector of the position space transseries as the inverse Fourier transform of the corresponding momentum space coefficients, as long as this inverse Fourier transform exists for positive real t. This is always the case for the perturbative sector of (1.12) (i.e. theâ n coefficients). In fact, the perturbative sector proceeds straightforwardly, and we can immediately obtain closed-form results which we present in the remainder of this section. Nonperturbative contributions in position space are both subtle and interesting, and later sections of this paper are devoted to this topic.
The position space coefficient, can be computed in closed form, with the result where * represents a convolution in the spatial coordinate, and the (n+1)-th kernel K n+1 (x) is given by An alternative representation of the same result is given by where G 0 (t, x) is the propagator for the heat equation, Withâ n+1 (t, k) and a n+1 (t, x) now computed, we can arrange them to produce a piece of the full solution forρ and ρ respectively. These we label with a superscript 'H', and since we obtained them by inverse Fourier transform of a series in we have added a subscript to denote the fact that they are not exact in . In the light of (2.10), each term in ρ (H) (t, x) corresponds to a gradient series in ∂ 2 x acting upon a solution of the heat equation that depends on the initial data. In this sense, as advocated in the Introduction and expanded upon in Appendix C, ρ (H) (t, x) provides a particular reorganization of a gradient expansion construction of the contribution of the hydrodynamic mode to ρ(t, x), hence the label 'H'. The convergence of ρ (H) (t, x) will be the focus of the next section, and we refer the reader to Appendix D for an analysis of how well ρ (H) (t, x) -when we set = 1 -reproduces the exact microscopic ρ(t, x) in a particular example.

Large-order behavior
In this section we analyse the large-order behavior of ρ (H) (t, x) =1 . We first provide a fully general, model-independent condition for the convergence of this object, which relies only on the support of the initial data in momentum space. Then, we discuss the large-order behavior of the series for initial data where this condition fails. We begin by splitting the full microscopic ρ(t, x) into hydrodynamic and nonhydrodynamic mode contributions, ρ(t, x) = ρ (H) (t, x) + ρ (N H) (t, x), as defined by individual contributions to the Fourier integral along a path γ, 9 where (2.13) with analogous expressions for the nonhydrodynamic mode, N H. As written, these expressions are exact in , and notationally this is indicated by the lack of an subscript. For our choice of initial data ((2.4) with u = 0), we have (2.14) By series expanding the exactρ (H) (t, k) around = 0 we recover the perturbative sector of the momentum space transseries,ρ (H) (t, k), computed earlier (2.12). The existence of 9 Note that the splitting is only unequivocally defined once a particular integration path γ is specified.
While, for entire initial dataû(k) andv(k), ρ(t, x) is the Fourier transform of an entire function, this is not the case for the individual hydrodynamic and nonhydrodynamic contributions as defined in the text, due to the branch points of the mode frequencies at k 2 = k 2 c . Each individual contribution is well-defined only after a particular γ to go around the corresponding branch cuts has been provided. two branch points in ω H ( k) and f H (k) at k = ±|k c | implies that, for |k| > |k c |,ρ (H) is a divergent series. On the other hand, in defining ρ (H) we assumed that the momentum space integral commutes with the infinite sum inρ (H) , in such a way that the individual series expansion coefficients,â n and a n , were directly related by the Fourier transform. Hence, From the expression above it is immediate to see that, unless the momentum space support of the initial datav is restricted to 10 the integral contains contributions from the momentum space region whereρ (H) diverges.
Modulo fine-tuned cancellations, the natural expectation to draw from this analysis is that the series ρ (H) does not converge if the initial data does not satisfy the condition (2.16).
It is important to note that, even if we have obtained (2.16) for a particular class of initial data for the telegrapher's equation, its applicability is not restricted to this example. As long as the hydrodynamic mode frequency ω H (k) has a complex singularity at |k| = |k c | and f H (k) is analytic for |k| ≤ |k c |, equation (2.16) will hold. It can be argued that the existence of this complex singularity in ω H (k) is a necessary condition for a microscopic theory to behave causally [28]. From this viewpoint, ρ (H) would be a divergent series in any theory that respects relativistic causality for generic initial data with unrestricted momentum space support.
In order to illustrate how the convergence condition (2.16) applies to our case, we consider compactly supported initial data in momentum space of the form where Θ denotes the Heaviside step function. These initial data correspond to a regularized δ-function in position space. For x = 0, a n (t, 0) can be explicitly computed To check whether the series expansion defined by the coefficients above is convergent, we use the ratio test and compute numerically q = lim n→∞ q n , q n = |a n+1 /a n |. (

2.19)
A sufficient condition for convergence is that q < 1. In figure 2 (left), we plot q as a function of A at t = 1. We always find that q = A 2 to exceedingly good accuracy. This implies that ρ (H) (t, x) =1 only converges for A < 1, i.e., for initial data with no support 10 Recall that one should set the formal parameter = 1 at the end of the analysis. for |k| > |k c | = 1 in k-space. If the support exceeds this bound, we obtain a divergent asymptotic series. This statement is time-independent, as figure 2 (right) shows. It is worth mentioning that, as long as |k c | < A < ∞, the series divergence is just a geometric one. A factorial growth of the a n (t, 0) coefficients only appears in the 20) and it follows that the ratio between successive coefficients increases linearly with n, implying a factorial divergence. We have found that the correlation between the infinite support of the initial data in k-space and the factorially divergent character of ρ (H) (t, x) =1 always extends beyond the particular case of δ-function initial data discussed above. This is the behavior to expect for generic initial data: initial data with a sharp cutoff in momentum space represent finetuned states from a position space perspective, in the sense that any modification of the initial data which is localized in position space -for instance, by any Gaussian -would make the momentum space support unbounded. Let us illustrate further the generic factorial divergence of the perturbative series by focusing on the example provided by Gaussian initial data of the form 11 As it happened with our compactly-supported initial data, this case is also analytically solvable at x = 0. We obtain These coefficients, which diverge factorially as n → ∞, show a more intricate behavior than their δ-function initial data counterparts. As n → ∞, where the function µ depends only the scaling variable and is given by This implies that, for Gaussian initial data, q n only shows the behavior of δ-function initial data for times greater than a critical time. In figure 3, we plot how r evolves with t for s = 1/4 (left) and s = 16 (center). These two curves, when re-scaled by 8Dτ /s 2 and expressed as functions of 2Dt/s 2 , collapse to a universal function, which agrees with µ(η) (right, in blue).    Let us close this section by mentioning that the r parameter can be equivalently extracted by means of Borel transforms. Given a factorially divergent asymptotic series its Borel transform f B (z) given by and, by definition, has a finite convergence radius in z. This convergence radius is set by the singularity which is closest to the origin in the complex z plane. In general, the analytical continuation of f B (z) to complex z cannot be computed in closed form and some approximation technique needs to be invoked. A common choice is to use Padé approximants. The (m, n)-Padé approximant of f B (z) is the unique rational function A well-known property of Padé approximants is that branch cuts in f B (z) manifest themselves as lines of pole condensation in P(z). For the Gaussian initial data (2.22), including their δ-function limit s → 0, we always find a line of pole condensation along the positive z-axis starting at z c = r, where r is given by (2.24). Therefore, we conclude that the exact analytically-continued Borel transform has a branch point at this location, as expected on the basis of the large order behaviour of the series. Note, however, that for Gaussian initial data the location of the branch point depends on t/t c , and is given by t/τ only for t > t c . For t < t c the branch point location also depends on the initial data through its dependence on s. Note that a condensation of poles can hide more than a single branch cut, 12 as happened, for example, in [14]. We will come back to this point in the next section, where we discuss the physical origin of these singularities and their dependence on t.

Saddle point analysis
In this section we take a different perspective and address the → 0 limit of the hydrodynamic mode contribution ρ (H) (t, x) by means of a saddle point analysis. As we will show, and as expected on general grounds, this procedure allows us to understand the physical origin of the branch points in the Borel plane.
To start, let us consider a schematic integral and focus on its behavior for |λ| → ∞. For us, the parameter λ will be simply 1 2 and we introduce it purely for notational convenience. This behavior can be obtained from a 12 For example, for fB(z) = √ 1 − z + √ 2 − z, the Padé approximation (2.29) is going to display a condensation of poles emanating from 1 towards larger values of z on the real axis. The poles with 1 ≤ z < 2 will correspond to the branch cut associated with √ 1 − z and the poles with z ≥ 2 will include also contributions from the branch cut associated with √ 2 − z.
saddle point analysis. The relevant analysis can decomposed intro three subsequent steps. First, one finds the stationary points u s of the 'action' S(u), Second, the steepest descent contours emanating from these saddle points are determined. The steepest descent contour γ s associated to the u s saddle point is the path emanating from u s along which Re S(u) decreases the fastest. It obeys Im λS(u(ξ)) = Im λS(u s ).
Finally, one decomposes the original integration path γ into steepest descent contours γ s and calculates the corresponding integrals in a |λ| → ∞ expansion. The steepest descent path associated to a saddle u 0 can collide with another saddle u 1 when arg λ is such that Im λS(u 0 ) = Im λS(u 1 ). These saddles are known as adjacent saddles [70] and play a prominent role in controlling the large-order behavior of the |λ| → ∞ expansion around the original u 0 saddle. In particular, when considering the Borel transform of this |λ| → ∞ expansion, adjacent saddles manifest themselves as branch points in the Borel plane [71], located at Let us apply this line of reasoning to understand the large-order behavior of ρ (H) (t, 0) we reported in the previous section for the case of Gaussian initial data (2.22). In this case, we have that with action The analogous expressions for ρ (N H) , G N H and S N H are obtained by flipping the sign of the square root in (3.6). The most important point to draw from the expression above is that the initial data contribute nontrivially to the relevant action and, as a result, also to its saddle points. Up to this point, our expressions are exact in . When x = 0, the saddle points of S H , S N H are as follows. For S H , we have a single saddle point located at k = 0 with vanishing action. On the other hand, S N H has three saddle points. The first one is located at k = k 0 = 0 at all times and has action The positions of the two remaining ones are time dependent Re S NH (t,k,ϵ=e ⅈθ ) with t c given by (2.25) and lead to actions At t = 0, the k ± saddles start at the critical momenta ±k c . As t grows, they approach each other along the real k-axis, until colliding with the k 0 saddle at t = t c . Past this time, they recede from the origin in opposite directions along the imaginary k-axis.
The perturbative sector of our transseries, ρ (H) (t, x), corresponds to the saddle point expansion around the hydrodynamic saddle. 13 Crucially, for complex the steepest descent path emanating from this saddle crosses the branch cuts in the hydrodynamic dispersion relation and continues on the nonhydrodynamic sheet. Once there, this path can collide with adjacent nonhydrodynamic saddle points at specific values of θ. An example of this behavior is provided in figure 4, where we also show the steepest descent contours associated to the three nonhydrodynamical saddles. lines, from lighter to darker, represent the values θ = 2.5 × 10 −3 , 2.5 × 10 −4 and 2.5 × 10 −5 , while dashed red lines represent their negatives. As θ approaches zero, the steepest descent path gets progressively closer to the k ± saddles before veering off to infinity. Furthermore, θ = 0 marks a discontinuous change in the path behavior, which is the reason why we considered θ as a parameter to vary; for instance, the quadrant in which the red curve extends to infinity undergoes a sudden change. These observations are compatible with the hypothesis that k ± are the adjacent saddles controlling the divergence of the perturbative series expansion. On the other hand, as illustrated in figure 5 (right), the collision of nonperturbative saddle points at t = t c causes the nature of the adjacent saddles to change: past t c , k 0 becomes a new adjacent saddle for the hydrodynamic steepest descent contour.
We can now relate this saddle point analysis to the singularities we observed in the Borel plane in section 2.2. Since for t > t c we have that |S 0 | < |S ± |, the above observations entail that the branch point of the Borel transform which is closest to the origin should go from being located at z c = −S ± for t < t c to being located at z c = −S 0 for t > t c . This prediction matches precisely the behavior of r reported at the end of section 2.2.
Let us emphasize that for t > t c the branch point associated to the S ± adjacent saddles, z ± = −S ± , does not disappear. What actually happens is that, since both z 0 = −S 0 and z ± are real positive quantities with z ± > z 0 , the branch cut associated to z 0 superimposes with the branch cuts associated to z ± , causing the latter to be superficially hidden in the Borel plane (as we illustrated using a simple example in footnote 12).
To expose the hidden z ± branch points, we proceed along the lines of [72]. We take our original Padé approximant P(z) and introduce the variable change z = z(w), with z(w) analytic at w = 0. Then, we series expand P(z = z(w)) around w = 0. Finally, we compute the Padé approximant of the resulting series, P(ω). By a suitable choice of z(w), the images of our original branch points lead to non-superimposing branch cuts in the w Borel plane. A convenient choice of variable change is as follows. Define with z c being, as before, the branch point closest to the origin in the z Borel plane. The variable change (3.11) maps a point z ∈ (z c , ∞) to two complex conjugated images on the right half of the |w| = 1 unit circle, with z c being mapped to w = 1. For t < t c , the pole structure of P(w) is as shown in figure 6 (left). We see three lines of pole condensation: two of them are complex conjugated and emanate from the singular points of the map (3.11), located at w = ±i, and are therefore unphysical; the third one starts at w = 1 and runs along the positive real axis. This latter line corresponds to the original branch cut starting at z c = z ± . For t < t c we see no trace of additional branch points associated to z 0 in the w plane, confirming that for t < t c k 0 is not an adjacent saddle.
On the other hand, for t > t c , this state of affairs changes. As figure 6 (right) shows, besides the branch point at w = 1 corresponding now to z c = z 0 , we also find two additional, complex conjugated lines of pole condensation emanating from the unit circle. It is immediate to check that these points are nothing but the images of z ± under the map (3.11). These images are represented by the upper and lower blue stars in the figure.

Superexponential decay of the nonhydrodynamic saddle points
From a physical standpoint, the existence of the k + , k − saddles bears crucial consequences for the late-time behavior of the nonhydrodynamic mode contribution ρ (N H) for the Gaussian initial data, (2.22). To expose these consequences, for the rest of this section we focus on the original = 1 telegrapher's equation. We take the integration path γ, used to compute ρ (N H) , to lie strictly below the real axis in the complex k-plane. With this choice of path, ρ (N H) is real. Linearity of the telegrapher's equation means that ρ (N H) is now a fully-fledged solution to the telegrapher's equation on its own, associated to the initial data,û in such a way that f H (k) = 0. 14 In the remaining part of this section, we will denote this solution simply as ρ.
The chosen integration path γ can be deformed into the steepest descent contour associated with the k − saddle-point (see figure 7). Therefore, the natural expectation is that the late-time behavior of the initial data (4.1) at x = 0 is controlled by the S − action, which predicts a late-time decay faster than e − t τ . Two questions arise naturally at this point. The first one is whether the faster-thanexponential decay we have uncovered extends to finite x. The second one is whether the existence of these additional saddle points extends to generic initial data with Gaussian asymptotic behavior as k → ∞.
Both questions can be answered affirmatively by the following argument. For initial data with Gaussian asymptotic behavior, the relevant action to consider in the nonhydrodynamic sheet at finite x is (we take = 1 from now on) By solving the saddle point equation in a late-time expansion, we find a saddle k * located in the lower half complex k-plane given by 15 with action As the expression above shows, the finiteness of x does not change the leading order latetime behavior of S * we encountered before for x = 0. Moreover, only the asymptotic behavior of the initial data at large k matters in reaching this conclusion.
To test whether the late-time behavior predicted by (4.4) is actually realized, we select the logarithmic derivative ∂ t log ρ(t, x) as our probe. According to (4.4), we must have that In order to check whether equation (4.5) holds, we restrict ourselves to initial data in which we multiply the Gaussian appearing in (4.1) by a polynomial P (k) (taken to be a function of k 2 with real coefficients), fix s, and compute numerically ∂ t log ρ(t, x) for a range of spatial positions. In figure 8, we show results for P (k) = 1 (left) and a fourth-order P (k) with random coefficients drawn from the interval [−1, 1] (right), for s = 1 and at x = 0, 1 and 2. These results are in agreement with the hypothesis that the nonhydrodynamic contribution displays a faster-than-exponential decay. The physical origin of the faster-than-exponential decay is propagation of the data with a Gaussian tail rather than an effect of dissipation governed by the imaginary part of the mode frequency. To see this, note that (4.4) can be suggestively rewritten as  seen in figure 1. Figure 9 shows the utility of this line of reasoning. In this figure we plot a dense set of snapshots of ρ -determined by direct numerical integration -as a function of x. We observe that the initial data (4.1) indeed evolve as a propagating wave that recedes from x = 0 as time grows. Furthermore, as the argument above suggested this wave gets progressively damped, but just exponentially so. -20 -For Gaussian initial data, figures 5 and 6 show that both the saddle point structure and the Borel plane structure take a different character depending on whether t is larger or smaller than the critical time t c = s 2 /2D. The transseries also undergoes an abrupt change, known as a Stokes phenomenon [27]. It occurs when crossing a point where the imaginary part of the action is the same, i.e. Im (S 0 − S ± ) = 0. (5.1) If we consider the action as a function of complex t, with all other parameters real 16 and x = 0, this condition is realized when As the Stokes line at s 2 2D is crossed along the real axis, the contribution of the S 0 saddle is turned off and the contributions of S − and S + are turned on.
In terms of the coefficients b n (t, x), this is seen as follows. Recall that for Gaussian initial data of the form (2.22),b n (t, k) = −â n (−t, k), whereâ n (t, k) is given by (2.6a). For k ∈ R and as |k| → ∞, we have thatb n+1 (t, k) behaves aŝ Therefore, the Fourier integral that computes b n (t, x) exists for positive real t as long as t < t c . The result is b n (t, 0) = −a n (−t, 0) with the latter given by (2.23). This changes dramatically for t > t c , where the Fourier integral diverges along the standard Fourier contour. The contour must be deformed to the k − saddle point. For x = 0 this series expansion can be computed in closed form, with the end result that where η = t/t c > 1. This series is factorially divergent. We find that the ratio test behaves as where r is a function given by the difference between the k 0 and k − saddle point actions, 16 There is no obstruction against considering complex values of s, which give rise to oscillating solutions.
We provide an example of this behavior in figure 10 (left). In line with this result, the Padé approximant to the Borel transform of the asymptotic expansion (5.4a) displays a line of pole condensation along the negative real axis, starting at z c = S − − S 0 (see the right plot in figure 10). This location of the branch cut follows from the fact that k = 0 is an adjacent saddle for the steepest descent contour emanating from k − when arg 2 = π.  To summarize, for these initial data, we take our nonperturbative transseries sector to be defined by We see that the form of the nonperturbative transseries sector depends on the spacetime location. This is the reason why we started with the transseries in momentum space, see (1.12), which is uniquely defined in terms of the modes in the system. This kind of Stokes phenomena in spacetime is not unique to this model. It has been investigated previously in, for example, references [63][64][65]. Those studies show that nonlinearity can also be handled and that a more intricate higher-order Stokes phenomenon can occur.
To close this section, let us emphasize the role that the initial data plays in the spacetime picture. Rather than the Stokes phenomenon occurring as t is varied, one can equivalently regard it as occurring as the initial data is varied. It arises here because the transseries coefficients depend on the initial data, which we expect is generic and holds also in nonlinear models. If a richer set of initial conditions were considered, a richer set of transseries sectors could be the result. The takeaway lesson here is that both the spacetime location as well as the initial data must be taken into account when formulating the transseries.

Discussion
The main motivation for our work is understanding nonequilibrium phenomena with a hydrodynamic tail by expressing them as transseries with resurgent relations connecting their various sectors. This point of view originates from the studies of expanding matter in ultrarelativistic heavy-ion collisions described by the paradigmatic example of Bjorken flow. Our work proposes and explores an approach that allows one to frame more general examples of dynamics in the same kind of language. We focus on the linearized regime but, reminiscent of [73], transseries techniques can in principle be applied when deviations from equilibrium are large. It would very interesting to study this question in detail, and to make contact with far from equilibrium attractors, which we leave for future research.
To describe a nonequilibrium process with a transseries one needs to define a small parameter. Guided by results from Bjorken flow, we introduce a formal parameter based on rescalings of spacetime coordinates (1.5), that organizes the hydrodynamic and nonhydrodynamic contributions into different transseries sectors. While the momentum space picture is straightforward, when passing to coordinate space we see new features which are as yet unseen in Bjorken flow and other expanding plasma systems.
In particular, we find that the initial conditions affect the form of nonperturbative contributions in in the spacetime picture. In our work we focused on a particular simple yet rich and widely encountered equation of motion -the telegrapher's equation (1.7)and a few classes of initial conditions. The nonhydrodynamic sector of the transseries took on two different forms. One is as a nonpropagating transient mode evaluated at zero momentum, similar to what was found in the Bjorken flow. However, when the initial conditions produce propagating wave packets, the receding tails of these wave packets gave rise to new transseries sectors. We have seen then that the decay of the nonhydrodynamic data is not only governed by the transient mode at zero momentum, but also by the form of the initial data and the dispersion relations at finite momentum.
While we observed this phenomenon for the telegrapher's equation, its ingredients seem to originate from the underlying causality of the system. This is necessarily shared by all the models of relativistic matter, in particular holography, which suggests that it is ubiquitous. 17 In studies of the transition to hydrodynamics in relativistic heavy-ion collisions using holography, the dominant theme has been the decay of transient modes as the mechanism governing it. Here we see that contributions from the nonhydrodynamic sectors can propagate away from a given spatial location, which effectively may render them zero. It would be very interesting to understand possible phenomenological implications of this observation.
Finally, let us comment on the utility of transseries solutions. The transseries allows one to organize hydrodynamic and nonhydrodynamic phenomena using a unified mathe- 17 Perhaps similar phenomena can even occur in the holographic Bjorken flow. One can view the gravity dual to the Bjorken flow as a set of nonlinear wave equations with constraints in two variables. This is not too dissimilar from what we considered here if the transseries analysis is extended into the bulk. matical language. This is crucial when the hydrodynamic gradient expansion diverges and requires resummation. The transseries provides a framework to resum it yielding a unique answer for a nonequilibrium solution. Furthermore, transseries provides a way to encode different asymptotic behavior in different spacetime regions, or for different initial data. Transitions between these behaviors are described by the Stokes phenomena. In our case, these considerations allowed us to uncover a new interesting physical effect in the context of hydrodynamization.
A Müller-Israel-Stewart in the shear channel MIS theory is the simplest phenomenological model of stress-energy tensor equilibration that agrees with relativistic Navier-Stokes hydrodynamics at leading order in the gradient expansion and, at the same time, respects causality.
For a conformal fluid, the construction of MIS theory proceeds as follows. We start from the Landau-frame constitutive relations of first-order viscous relativistic hydrodynamics in d-dimensional Minkowski space (1.2) with the equations of motion of the theory being nothing but the conservation of the full energy-momentum tensor, ∂ a T ab = 0. The algebraic relation between Π ab and σ ab implied by (1.2) entails that first-order relativistic viscous hydrodynamics violates causality. In MIS theory, this problem is overcome by promoting Π ab to a set of new independent degrees of freedom that obey a relaxation equation, in such a way that the original constitutive relation (1.2) is recovered at times sufficiently larger than a new time-scale set by a relaxation time τ , The operator D a is a Weyl-covariant derivative originally introduced in [74].
In this work, we consider infinitesimal fluctuations of this theory away from thermal equilibrium. We thus write where E 0 is the equilibrium energy density, and treat δE/E 0 and u 2 as infinitesimally small. Moreover, we also make the symmetry assumption that the hydrodynamic variables are independent of x 1 , ..., x d−2 . Defining x d−1 ≡ x, this corresponds to δE = δE(t, x), u i = u i (t, x). This ansatz can be viewed as the minimal generalization of a boost-invariant flow, for which the hydrodynamic variables would also be functions of t and x, but only through the combination √ t 2 − x 2 . As mentioned in the main text, the telegrapher's equation emerges when considering a shear channel fluctuation, for which δE(t, x) = 0 and u i (t, x) = u 1 (t, x)δ i,1 with no loss of generality due to rotational invariance. Defining and linearizing in the velocity fluctuation amplitude, the equation for energy-momentum conservation and the dynamical constitutive relation (A.1) can be expressed as where the diffusion constant is D = η/(sT ). As mentioned in the Introduction, the linear PDE system (A.4)-(A.5) is well-known in the literature [30,33]. Combining both equations, we recover (1.7).

B The causality of the telegrapher's equation
In this appendix we show that the telegrapher's equation respects causality (see also [30]). The basis of our proof is the following theorem [75].
Theorem (Paley-Wiener). Let f (x) ∈ L 2 (R) be supported in x ∈ [−A, A]. Then its Fourier transformf (k) belongs to L 2 (R) and is an entire function of exponential type A.
We remind the reader than an entire function is a function which is analytic everywhere in the finite complex plane; an entire function of exponential type A is an entire function that obeys the bound |f (z)| ≤ Ce A|z| , ∀z ∈ C.
Consider the most general solution to the telegrapher's equation, and suppose that the initial data are supported only between −R and R. Causality demands that, at time t, ρ(t, x) is supported at most in the interval |x| ≤ R + t. Therefore, we have to show that ρ(t, k) is an entire function of exponential type at most R + t.
The most general square-integrable solution of the telegrapher's equation can be written as We start by noting that an entire function of the square root of a complex number is itself entire if the original function is even. Therefore, both cosh ∆t 2τ and 1 ∆ sinh ∆t 2τ are entire functions of k. It then follows thanρ(t, k) is entire in k, sinceû,v are entire by assumption, and the product of two entire functions, as well as their sum, are themselves entire.
To show that (B.3) is of exponential type at most R + t, we proceed as follows. First, we note that both cosh ∆t 2τ and 1 ∆ sinh ∆t 2τ are of exponential type D/τ t. This also applies to their sum. Next, we recall that the product of two functions of exponential types σ 1 and σ 2 is at most exponential type σ 1 + σ 2 . Therefore, the type of each term in the sum (B.3) is at most R + D/τ t. Finally, since the type of the sum of two functions of exponential types σ 1 and σ 2 is smaller or equal than max(σ 1 , σ 2 ), it follows that the type of ρ(t, k) is at most R + D/τ t. As long as D/τ ≤ 1, we see that the system respects relativistic causality.
It is worth pointing out that this result conforms with the expectation that, in any local quantum system, causality bounds the diffusion constant in terms of the effective light-cone speed and the local equilibration time [76] (see also [77][78][79]). In the case at hand, the effective light-cone speed is to be identified with the speed of light, and the local equilibration time with τ .

C The large-scale expansion and hydrodynamics
In our previous work [28], we described the most general construction of linearized hydrodynamics for a neutral conformal fluid in arbitrary number of spacetime dimensions. Specializing to a shear channel fluctuation, the hydrodynamic description of the microscopic field ρ is provided by the conservation equation (A.4) in combination with the purely-spatial gradient expansion of the constitutive relations. For the MIS case, and in the notation of appendix A, the hydrodynamic gradient expansion reads The transport coefficients c n are extracted from the microscopic shear hydrodynamic mode by a matching computation. In MIS, they are given by where C n is the n-th Catalan number.
As we have illustrated extensively in the main text, the perturbative series ρ (H) (t, x) corresponds to the hydrodynamic shear mode contribution to the exact ρ(t, x). It is then natural to ask whether the effective description of MIS theory in terms of classical hydrodynamics contains exactly the same physical information as the perturbative sector of our transseries. Let us view (C.1) as a formal series, making no assumptions about the relative size of subsequent terms, and perform the rescaling (1.5) both at the level of the gradient expansion (C.1) and the conservation equation. Then, it turns out that the perturbative sector of our transseries, ρ (H) (t, x), solves the resulting system order-by-order in an expansion around = 0. Equivalently, if we ignored the specific values of the transport coefficients c n but somebody handed to us ρ (H) (t, x), we could fix the former by demanding that the latter is a solution order-by-order in an expansion around = 0. This procedure could then be viewed as the position space counterpart of the matching computation performed in [28].
It should be noted that while ρ (H) (t, x) is a solution, the initial conditions we imposed on the a n (t, x) coefficients in the main text are completely unnatural from the perspective of the hydrodynamic description of the system. These boundary conditions relied on the existence of the nonperturbative sector of the momentum space transseries, which allowed us to enforce simultaneously thatρ(0, k) =û(k), ∂ tρ (0, k) =v(k). This nonperturbative sector is absent now, reflecting the fact that the algebraic relation between J(t, x) and ρ(t, x) in the hydrodynamic description implies the loss of the nonhydrodynamic degree of freedom.
We now proceed to explain the natural choice of initial conditions from the perspective of the hydrodynamic description. After the spacetime rescaling, our equations of motion are given by The ansatz results in the following nested ODE system with an equivalent expression for the odd coefficients. Again, the n-th term in the series expansion (C.4) is a solution of the heat equation sourced by the n − 1 previous orders. Let us assume that, at a time slice t = t 0 , ρ(t 0 , x) = ρ 0 (x) is known. In this situation, the natural boundary conditions to impose on (C.4) are that, at t = t 0 , the leading-order term a 0 (t 0 , x) agrees with ρ 0 (x), with the remaining higher-order terms vanishing. Due to the structure of (C.5), these boundary conditions result in the vanishing of the odd order terms in (C.4) for all times.
With this choice of boundary conditions, (C.5) can be explicitly solved in closed-form. In Fourier space, we find that 18 a 2n (t, k) = δ n,0ρ0 (k) − (−1) n c n t k 2(n+1) 1 F 1 (1 − n, 2 + n, Dk 2 t) e −Dk 2 tρ 0 (k). (C.6) The expansion coefficients can be explicitly computed in position space. The zeroth-order one is the solution of the heat equation given by a 0 (t, x) = (G 0 * ρ 0 )(t, x), where G 0 is given in equation (2.11). The remaining ones are a 2n (t, x) = c n n−1 q=0 Γ(n)Γ(n + 2) Γ(n − q)Γ(n + 2 + q)Γ(q + 1) x a 0 (t, x). (C.7) As it happened with our previous choice of initial conditions, each term is a gradient series in ∂ 2 x acting on a particular solution of the heat equation.

D Applicability of the truncated perturbative series
The practical usefulness of the perturbative piece of the asymptotic expansion developed in the previous section is that, when truncated to low order, it provides a good description of the exact velocity field ρ in some spacetime regions. Owing to our general discussion in the Introduction, where we introduced our large-scale expansion, we expect that these regions correspond to those in which the nonhydrodynamic mode contribution has significantly decayed. Let us illustrate this by considering the time evolution of δ-function initial data of the form u(x) = 0, v(x) = δ(x). In this case, ρ(t, x) corresponds to the propagator of the telegrapher's equation [30], The expression above clearly demonstrate that the telegrapher's equation is causal: for any t > 0, ρ(t, x) vanishes if |x| > D τ t. We compare the exact ρ(t, x) above with the with the perturbative expansion ρ obtained by means of equation (2.10) truncated to first and third nontrivial order. The results can be found in figure 11. We clearly see that, for any x, ρ =1 never provides an accurate description of ρ at early times. This is due to the fact that, in this regime, the nonhydrodynamic mode contribution, which is necessary to enforce the initial condition ρ(0, x) = 0 we chose, is still sizable. This state of affairs changes at later times. In particular, focusing on a fixed x as t grows, we eventually observe a very good agreement between ρ and the low-order truncated ρ Another important point to be drawn from figure 11 is that ρ (H) =1 is never a good approximation of the exact microscopic ρ outside the causal cone of the system. While the exact ρ vanishes there, ρ (H) =1 does not. The reason behind this difference is that ρ (H) is solely built out of the hydrodynamic shear mode, and is therefore blind to the nonhydrodynamic contribution needed to enforce the causal response of the system. This fact provides a nice illustration of the general lesson that the nonhydrodynamic sector is essential to guarantee causality [1].

E Large-order behavior for Gaussian initial data
In the main text, we showed that the large-order behavior of the perturbative series is controlled by the adjacent nonperturbative saddle points. In this appendix, we elaborate further on this connection for the case of Gaussian initial data. Following Berry and Howls [70], we can express the a n+1 (t, 0) coefficient of the perturbative series expansion as the following contour integral, a n+1 (t, 0) = Γ n + 1 2 with Γ a positively-oriented path enclosing k = 0 in the hydrodynamic sheet. It can be demonstrated that (E.1) can be alternatively expressed as a contour integral in the nonhydrodynamic sheet, a n+1 (t, 0) = Γ n + 1 2 where Γ is a path starting at ∞ − i0 below the right branch cut, going around the right branch point, and ending at ∞ + i0 above the right branch cut. 19 In order to analyze the behavior of (E.3) when n → ∞, it is natural to decompose Γ into the steepest descent contours associated to the adjacent saddles discussed in the main text, and employ a saddle point approximation afterward. 20 The relevant steepest descent contours involved in computing (E.3) in a large n asymptotic expansion are depicted in figure 12. For t < t c , only the k + saddle contributes to a n (t, 0). On the other hand, for t > t c , we have contributions both from k + , k − and k 0 . Let us focus on the former case first. At leading order, it is immediate to show that from which it follows that r = lim n→∞ n a n (t, 0)/a n+1 (t, 0) = −S + , as reported in the main text. A plot of the deviation of the ratio , (E.6) 19 In writing (E.3), we have taken into account that left and right branch cuts contribute equally to the total result. 20 Since at k * = k+, k− or k0 we have that S N H (k * ) = 0, but SNH (k * ) = 0, these saddles are also stationary points of log(−SNH (k)).
from one, where a exact n+1 is given by (2.23) and a s.p. n+1 by (E.5), can be found in figure 13 (left). We have considered several different times before t c . It is readily seen that, as n → ∞, r n → 1 plus a O(1/n) correction, confirming the validity of equation (E.5). On the other hand, for t > t c , we have three separate contributions to consider. We find that at leading order the k 0 saddle contributes as while the first nontrivial term of the combined contribution of the k + , k − saddles is given by Γ n + Since |S 0 | < |S + |, the k ± contribution is exponentially suppressed with respect to the k 0 one, 21 which now governs the divergence of the perturbative series expansion. The behavior of r n for t > t c with a s.p. n (t, 0) given by (E.7) is illustrated by figure 13 (right).

F Large-order behavior for Lorentzian initial data
Another family of initial data that allows to compute in closed-form the coefficients of the perturbative expansion at x = 0 is that of Lorentzian initial data, v(x) = α π(x 2 + α 2 ) ,v(k) = 1 2π e −α|k| , (F.1) where we find a n+1 (t, 0) = α 2n+1 τ n+1 4π 2 D n+1 t 2n+1 Γ n + As it happened for Gaussian initial data, we find that, for n → ∞, q n = lim n→∞ | a n+1 (t,0) an(t,0) | is linear in n, implying that the series is factorially divergent. The large n behavior determined numerically is compatible with irrespectively of the value of α. See figure 14 for an example of this asymptotic behavior in the case α = 1. G Reconstructing ρ from the transseries: Gaussian initial data case In this appendix we illustrate how Borel resummation can be employed to reconstruct the exact ρ(t, x) from the different transseries at our disposal. The agreement between the Borel resummation and the exact solution determined by direct numerical integration provides a nontrivial consistency check of the results presented in the main text. For the sake of brevity, we focus on Gaussian initial data of the form (2.22). Before starting, let us recall that the splitting of ρ(t, x) into hydrodynamic and nonhydrodynamic contributions is only unequivocally defined once a particular integration path γ is provided, due to the branch points of ∆(k). 22 To keep track of which integration path we are employing, we will denote by γ (+,−) the equivalence class of integration paths starting at −∞ + i0 + above the left branch cut and ending at ∞ − i0 + below the right one, with an analogous interpretation for γ (−,+) , γ (−,−) and γ (+,+) .

G.1 The perturbative sector
As mentioned before, the Padé approximant of the Borel-transformed asymptotic series displays a line of pole condensation along the positive real z-axis starting at z c = −S − and, as a consequence, we have to resort to lateral Borel resummations. In turns out that the negative (positive) lateral Borel resummation corresponds to the choice of γ (+,−) (γ (−,+) ) integration contour in the original integral, with the average of both corresponding to the γ (−,−) , γ (+,+) integration paths. We plot an example of this agreement in figure 15.

G.2 The nonperturbative sector
For the nonperturbative sector, we have a discontinuous change in the nonperturbative series when t = t c ; correspondingly, we discuss the t < t c and t > t c cases separately: • t < t c . In this regime, the Fourier integral defining b n (t, x) in terms of b n (t, k) converges. Therefore, b n (t, 0) = −a n (−t, 0), with a n (t, 0) given by (2.23). The relation between the lateral Borel resummations and the γ integration contours is the same as for the perturbative series, with the real result being given by the average of both lateral resummations. We plot this average in figure 16 (left), where we compare it with the ρ (N H) (t, 0) determined by a direct numerical integration, finding very good agreement between both quantities.
• t > t c . In this regime, the nonperturbative sector of our transseries is given by (5.4a)-(5.4b). The absence of singularities on the positive real axis in the Borel plane implies that the Borel transform of ρ (N H) (t, 0) is Borel resummable; hence the integration along the real axis agrees trivially with the average of the lateral resummations. In figure 16 (right), we compare the ρ (N H) (t, 0) determined by a direct numerical integration with the resummation result, finding again excellent agreement.