The large proper-time expansion of Yang-Mills plasma as a resurgent transseries

We show that the late-time expansion of the energy density of $\mathcal{N}=4$ supersymmetric Yang-Mills plasma at infinite coupling undergoing Bjorken flow takes the form of a multi-parameter transseries. Using the AdS/CFT correspondence we find a gravity solution which supplements the well known large proper-time expansion by exponentially-suppressed sectors corresponding to quasinormal modes of the AdS black-brane. The full solution also requires the presence of further sectors which have a natural interpretation as couplings between these modes. The exponentially-suppressed sectors represent nonhydrodynamic contributions to the energy density of the plasma. We use resurgence techniques on the resulting transseries to show that all the information encoded in the nonhydrodynamic sectors can be recovered from the original hydrodynamic gradient expansion.


Introduction
The discovery of quark-gluon plasma (QGP) at RHIC and the ongoing studies of its properties both there and at the LHC have lead to a burst of activity related to the theoretical description of this new state of matter. Tiny drops of QGP created in these heavy-ion collision experiments appear initially in highly non-equilibrium states. Nevertheless, after a short time this system reaches a state amenable to an effective description formulated in the language of hydrodynamics, which continues up until the medium hadronizes as the local effective temperature drops below the confinement scale. Recent activity aimed at understanding the equilibration of QGP has lead to significant progress concerning the foundational aspects of relativistic hydrodynamics (for reviews see e.g. Refs. [1,2]), such as the regime of applicability of hydrodynamics [3][4][5][6][7], the role and meaning of the hydrodynamic gradient expansion [8] and attractor behaviour far from equilibrium [9][10][11]. The difficulties of treating real-time evolution far from equilibrium in QCD have provided strong motivation to look for other systems where such an analysis can be more tractable. Since many of the fundamental questions concern relativistic hydrodynamics itself rather than its specific application to QCD, the study of related model systems has proved both fruitful and practical.
An important framework where many of these ideas were developed is N = 4 supersymmetric Yang-Mills theory (SYM), where the AdS/CFT correspondence provides an effective method of carrying out ab-initio calculations of highly non-equilibrium phenomena in a strongly-coupled quantum system by relating non-trivial observables at infinite coupling to solvable problems in classical Einstein gravity [12,13]. While supersymmetric Yang-Mills theory is very different from QCD, these differences are less pronounced at finite temperature and some of the results obtained using AdS/CFT appear to give a qualitatively useful picture even when interpreted in the context of QGP. Most importantly, since this approach has allowed for reliable calculations of equilibration, it has provided a priceless theoretical laboratory where the emergence of a hydrodynamic behaviour could be explored.
A kinematic situation of phenomenological interest is that of boost-invariant longitudinal expansion -Bjorken flow, which mimics the behaviour of matter during the early stages of an ultra-relativistic heavy ion collision [14]. Here, under the assumption of conformality, symmetry constraints impose that the expectation value of the 3 + 1 dimensional energy-momentum tensor can be expressed in terms of the energy density E as a function of a single parameter, the proper time τ elapsed since the collision event. This leads to an amazingly simple model, yet one which preserves much of the essential complexity of the original problem, making it possible to build a detailed picture of the transition from a highly non-equilibrium initial state to hydrodynamics.
A key approximation method which makes analytic calculations possible on the gravity side of the duality is the large proper-time expansion [15]. The bulk geometry obtained in this way implies that the energy density of N = 4 SYM at late times takes the form where u = τ 2/3 , with τ the proper time. Following Ref. [8], here and in the following we have suppressed a dimensionful parameter, effectively choosing to measure the proper-time in units set by the scale of the energy density. This expansion is directly related to the hydrodynamic gradient expansion of Bjorken flow, which is the expansion in powers of τ E 1/4 . In fact, given the coefficients of the late proper-time expansion one can calculate those of the gradient expansion and vice-versa -in other words, these two expansions contain the same information. With this understanding we will refer to Eq. (1) as the hydrodynamic expansion.
The dimensionless expansion coefficients ε (0) k are calculable from the gravity solution. The leading 240 coefficients were computed in Ref. [8], and further 140 were obtained in Ref. [16]. The series appearing in Eq. (1) is asymptotic, and it is natural to ask how to best interpret it.
Already in Ref. [8] it was observed that the analytic continuation of the Borel transform of this series has branch-point singularities at locations corresponding to frequencies of the least-damped nonhydrodynamic quasinormal modes of the AdS black-brane. This connection was subsequently studied at the level of hydrodynamics, where series expansions such as Eq. (1) can be generated as asymptotic solutions to MIS-type differential equations such as those introduced in Refs. [17,18]. These studies [9,19] lead to the realisation that the proper setting for the hydrodynamic gradient expansion is in fact a transseries [20,21], which systematically incorporates the contributions from nonhydrodynamic modes. These contributions take the form of terms which are exponentially suppressed at large times, but become significant at earlier times, where they encode the initial state data. 1 In this paper we present the first calculations which directly reveal the transseries structure at the level of the microscopic theory. This is done by finding a generalisation of the late proper-time expansion guided by the form of asymptotic solutions of the hydrodynamic model studied in Ref. [19]. However, in contrast to this hydrodynamic model which features just two conjugate nonhydrodynamic modes, in the case of N = 4 SYM there is an infinite number of nonhydrodynamic modes, naturally ordered by the rate of exponential suppression following from the complex values of the black-brane quasinormal mode (QNM) spectrum. Each such mode introduces a separate sector -a separate asymptotic series weighted by an independent exponential involving the QNM frequencies. The bulk solution will be described in detail in Section 2. Here we will present the ensuing form of the energy density of N = 4 SYM.
In order to describe this structure explicitly let us introduce some notation which we will employ throughout this paper. With each pair of nonhydrodynamic QNMs (characterised by a pair of complex frequencies A k , A k with k > 0) we associate a pair of unit vectors e k , e k , whose components vanish apart from a one in slots 2k − 1, 2k respectively. These unit vectors form a basis of an infinite-dimensional semi-lattice space N ∞ 0 . Each QNM is labelled by a vector n which is equal to one of these basis vectors and introduces an asymptotic expansion of the form, with some characteristic exponent β n , and expansion coefficients ε (n) m , with the understanding that the leading coefficient ε (n) 0 = 0. One can view this series as a hydrodynamic "dressing" of the individual QNMs. 1 It is worth noting that different asymptotic analysis of QNM frequencies of black holes can be found in the context of stability of black holes and their infinitely damped QNMs, see e.g. [22][23][24].
The contribution of each pair of QNMs comes with a pair of complex normalisation constants, the transseries parameters, which will be denoted by σ k , σk. These numbers can be interpreted as integration constants -they contain information about the initial state. We find that the energy density at large proper time τ (or, equivalently, for large u) has an expansion as a transseries of the form: where σ n ≡ σ n 1 1 σ The exponential weights appearing in Eq. (3) are expressed in terms of the vector of QNM and vectors n = (n 1 , n 1 , n 2 , n 2 , · · · ) with non-negative integer entries. Note that the vector of QNM frequencies in Eq. (5) determines the form of Eq. (4), which encodes the initial data of the boundary theory through the values of the transseries parameters. The collections of contributions to Eq. (3) with n equal to one of the basis vectors e k , e k will be referred to as fundamental sectors -each such sector corresponds to a specific QNM. Note however that the sum appearing in Eq. (3) also includes terms with vectors n corresponding to linear combinations of the basis vectors with non-negative coefficients. These contributions will be referred to as mixed sectors. They can be thought of as a reflection of QNM coupling. We will be referring to all of these as nonhydrodynamic sectors. The contribution to Eq. (3) corresponding to n = 0 (with β 0 = 2) is the original hydrodynamic expansion seen in Eq. (1).
This type of (resurgent) transseries structure appears often in studies of asymptotic expansions, and the expansion coefficients appearing in different sectors are known to be related by intricate consistency conditions, the so-called large-order relations, thoroughly described bý Ecalle's theory of resurgence [25] (see also e.g. the review [21] and references therein). These remarkable relations in principle allow us to extract the full content of the nonhydrodynamic sectors from the original hydrodynamic series. In practice, one can use them as a consistency check of the nonhydrodynamic sectors directly determined from the bulk solution. Such relations provide conclusive evidence that the transseries solution, supplemented by the relevant initial conditions, captured the full non-perturbative behaviour of the fluid. These large-order relations have successfully been used to predict novel non-perturbative physics [26][27][28][29][30][31][32] as well as in consistency checks of perturbative and non-perturbative sectors appearing in transseries solutions [19,[33][34][35][36][37][38][39][40][41][42][43].
Our findings are reminiscent of what is known about such series in the context of couplingconstant perturbative expansions in quantum mechanics, where these non-perturbative, exponentially suppressed sectors correspond to specific instanton solutions and are usually referred to as instanton sectors. It is well known that one needs more than just the instanton sectors as non-perturbative corrections to perturbative calculations of energy levels in quantum-mechanical models -one also needs multi-instanton sectors, which correspond to nonlinear effects arising from the quantisation condition (see e.g. Refs. [44,45] and Refs. within [21]). In the case at hand, the analogue of multi-instanton sectors are the mixed nonhydrodynamic sectors described above. Their content corresponds to going beyond the linearised approximation which gives rise to the QNM spectrum. From the point of view of resurgence, these additional sectors are unavoidable and are constrained by the aforementioned large-order relations. Verifying that they are satisfied gives us a very high level of confidence in the consistency of the transseries ansatz in the present context.
In Section 2 we construct a transseries solution for the gravitational dual to the Bjorken expansion of N = 4 SYM plasma by supplementing the large proper-time expansion with exponentially suppressed contributions, and in Subsection 2.3 we extract the coefficients of the dual gauge theory for several transseries sectors of Eq. (3) to high orders. The asymptotic expansions of these sectors show similarities with those arising in the context of the hydrodynamic model of Ref. [18], whose asymptotic behaviour was investigated in Ref. [19]. In particular, the large-order behaviour of the coefficients in these series involve various factorial and exponential growths with complex parameters (related to the exponential weights and characteristic exponents), which make it impossible to apply standard tools used in many recent studies of asymptotic series. The study of the large-order behaviour of these coefficients can be found in Section 3, where we introduce a novel method for reliably extracting resurgent information from an asymptotic series whose leading large-order relations depend on multiple complex factorial growths. The numerical checks showing the necessity of the full transseries solution Eq. (3) to describe the energy density of the N = 4 SYM plasma can be then found in Secs. 4 (for evidence of inclusion of fundamental sectors) and 5 (for the inclusion of mixed sectors).
This work offers a full non-perturbative picture of the late time expansion of the N = 4 SYM plasma's energy density, paving the way to use these new techniques in more general situations. In Section 6 we summarise our main results and discuss future applications of our work.
For ideal Bjorken flow, the local effective temperature of the plasma at the 4-dimensional Minkowski space boundary should behave as T ∼ τ − 1 3 at late times. We therefore fix s = 1 r τ − 1 3 3 Einstein's equations with cosmological constant Λ = −6 are given by R µν + 4g µν = 0.
so that in the late time limit the naive location of the horizon in the s co-ordinate will remain finite. Further, we note that gradient corrections should go like 1 τ T ∼ τ − 2 3 , so we fix u = τ 2 3 and expand the metric functions in inverse powers of u. Note also that 0 < s < 1, u > 0, and that late times correspond to large u.
Having in mind the hydrodynamic expansion as well as presence of the transient modes we consider a transseries ansatz for the metric functions of the form where Here (as already described in the introduction) n = (n 1 , n 1 , n 2 , n 2 , · · · ) ∈ N ∞ 0 is an infinite dimensional vector of non-negative, integer components which will be used to define the space of solutions. We will use the field equations to determine the quantities A = A 1 , A 1 , A 2 , A 2 , · · · and α which correspond to particular nonhydrodynamic sectors. 4 We will denote the hydrodynamic sector by 0 and define special unit vectors e k which have all zero entries except the 2k − 1 entry which will be set equal to 1. Similarly, e k is defined to have only the 2k entry non-zero (and equal to 1). These are defined so that the scalar products with A select the corresponding elements, in the sense that e k · A = A k , and e k · A = A k .
After substituting our ansatz into Eqs. (8) to (12), the linear independence of Ω n (u) will imply an infinite hierarchy of equations. These equations can again be expanded in inverse powers of u to find linear ODEs for the functions h We use residual gauge invariance to set the radial co-ordinate r to keep the horizon fixed at s = 1 at every order. This will amount to choosing h

Hydrodynamic and nonhydrodynamic sectors
The case of n = 0 will be referred to as the hydrodynamic sector. Any power of Ω n =0 (u) that enters into the equations of motion will remain linearly independent from terms proportional to Ω 0 (u). Therefore h can be found independent of the solutions to any other sector. 4 Note that we will assume that each component of A is not a rational multiple of another component.
The zero-th order solution (in our chosen gauge) that preserves flatness and bulk regularity will be given by a boosted black-brane written as [52], The cases of n = e k and n = e k will be called the fundamental nonhydrodynamic sectors. 5 Equations linear in Ω e k (u) can depend on the hydrodynamic sector and solutions within its own nonhydrodynamic sector. The zero-th order solution in each case is given by, where Z e k (s) satisfies and where we have defined ω = − 2 i 3 A to put Eq. (24) into the standard form of the Quasinormal Mode (QNM) equation given in infalling EF co-ordinates [53,54]. Imposing flatness and regularity in Z e k (s) turns Eq. (24) into an eigenvalue problem with an infinite number of solutions ω. We take the eigenvalues ω with negative real part to correspond to the unit vector e k , and those with positive real part to correspond to e k , with index k ordering them naturally by the negative imaginary part of each eigenvalue. The infinite dimensional vectors e k and e k can now be understood as spanning the space of QNM frequencies.
Each eigenfunction Z e k (s) is determined up to an arbitrary integration constant associated with the choice of Z e k (s = 1). This overall normalisation we will later identify as the transseries parameter. This freedom comes from the fact that (along with imposing flatness at the boundary) we only require Z e k (s) to be regular in the bulk, allowing a family of solutions which satisfy this condition. We can interpret each associated sector as an independent nonhydrodynamic excitation. The vector α is fixed through a similar eigenvalue problem at order i = 1. Our numerically computed values of α and A satisfy α = A 6 with high precision. 6 The cases of n = e k and n = e k will be called the mixed sectors. Once we have fixed ω in Eq. (24), if we replaced e k by a general vector n, then this equation would have no non-trivial solutions obeying the chosen boundary conditions for n = e k . All functions h for n = e k , e k and i ≥ 0 will be fully determined by solutions in the hydrodynamic and the fundamental 5 In the discussion that follows one can exchange e k for e k where appropriate. 6 A heuristic argument for this value is suggested by Refs. [8,52]. At late times Ω e k (u) is expected take the form of a decaying mode with frequency ω k in a slowly evolving plasma of effective temperature T ∼ E 1/4 given by, Equating A k = 3 i 2 ω k for every k, we can identify α = A 6 .
nonhydrodynamic sectors, and so will contain no free integration constants. In this way these sectors can be viewed as a cascade of interactions of the fundamental nonhydrodynamic modes rather than independent solutions.

Series solution and numerical implementation
Subsequent equations of motion for i ≥ 1 where n = 0, e k or e k , (and for i ≥ 0 where n = e k and n = e k ) can be found by directly substituting Eqs. (18) to (20) into Eqs. (8) to (10). They can be written respectively as where j d,n i , j h,n i and j b,n i are source terms, which are sequentially given in terms of solutions at lower orders in i, and solutions at the same order which are determined when solving the equations in the order (26) to (28). The linear operators above have relatively simple forms, In our chosen gauge, Eq. (12) can be written as a constraint at s = 1, with J (n) i a number determined by lower order solutions. Eq. (11) is redundant, being implied by the other four equations.
The implementation of the numerical methods used in this paper are a continuation of the techniques used in Ref. [8]. To integrate Equations (26) to (28) we use Chebyshev spectral methods [55] with 400 grid points over the s ∈ [0, 1] radial co-ordinate with a minimum of 240 digits of numerical precision. 7 As mentioned before, in our gauge we have chosen the warp factor h(r, t) to vanish at the naive location of the horizon (s = 1) which results in the boundary condition that h vanish at s = 0 and are regular at s = 1) fully determines the system for i ≥ 1 for the hydrodynamic sector, i ≥ 2 for the fundamental nonhydrodynamic sectors, and i ≥ 0 for the mixed sectors. 7 In each case the computation could be done on a typical laptop to high orders within several hours.
A non-obvious implementation is that of i = 1 for the fundamental nonhydrodynamic sectors, where the parameter α must be tuned so as to allow a solution for b (e k ) 1 which vanishes at s = 0. For α = A 6 this condition is satisfied for any finite value of b at the horizon (s = 1) must be chosen so as to fix b (e k ) i+1 (s = 0) = 0. This was employed in our code via a shooting method.

Energy density of the dual theory
It was shown in [15] that imposing the symmetries of conformal Bjorken flow and the conservation of energy and momentum leads to the following form of the energy-momentum tensor: The particular energy-momentum tensor expectation value in the state of the boundary theory which is dual to our bulk geometry can be evaluated via holographic renormalization [56], which determines the function f (τ ) appearing in Eq. (33) in terms of the expansion Re-expressing our solution in terms of u and s variables we can write, where the coefficients f (n) i can be read from our gravity solution as The factor of σ n (see Eq. (4)) is introduced into Eq. (35) so that the dependence on the choice of initial condition is explicit in the energy density. For generic n, some number of the first coefficients f in the sum (36) will be zero. For convenience we define ε (n) 0 to be (up to normalisation) the first non-zero coefficient of f i , i > 0 given by the subsequent coefficients), and absorb the shift of the index into the definition of the characteristic exponent, which will be given by β n . One then arrives at the final formula for the energy density in the form of a transseries, given in Eq. (3), which we repeat here for convenience, These coefficients ε (n) k have been included with this submission for the sectors Φ n with n = 0, e 1 , e 2 , 2e 1 , (e 1 +e 1 ), along with their corresponding exponential weights A i and characteristic exponents β n . For the hydrodynamic sector Φ 0 coefficients for the hydrodynamic expansion were taken from [16].
The normalisations for the hydrodynamic series Φ 0 and each of the fundamental sectors Φ e k , Φ e k ( associated to each k th QNM frequency) are not fixed, and have been chosen such that With this choice of normalisation, the mixed sectors have no freedom in their respective coefficients. The list of sectors we have included are as follows: • The first 250 coefficients of the fundamental sector Φ e 1 , and • The first 200 coefficients of the fundamental sector Φ e 2 , and β e 2 = − A 2 6 +3 with A 2 = i 3 2 ω 2 where ω 2 ≈ 5.1695 − i 4.7636 is the second lowest nonhydro QNM frequency.
• The first 100 coefficients of the mixed sector Φ e 1 +e 1 , with β e 1 +e Note that the fundamental sectors Φ e 1 , Φ e 2 are just complex conjugate of the sectors Φ e 1 and Φ e 2 , respectively. Furthermore, direct calculation shows that Φ 2e 1 is complex conjugate of Φ 2e 1 . Their respective characteristic exponents β n are also related by complex conjugation (as are their exponential weights). 8

The resurgent transseries of the energy density
The coefficients appearing in the different sectors of the transseries for the energy density (3) calculated in the previous section can be seen to grow factorially at large-order. Thus, each sector is formally defined as an asymptotic series. It is well known that one can remove the factorial growth by performing a Borel transform in each of the asymptotic sectors Φ n , which is the first step to the summation of these series. To actually perform such summation one needs to know the singularity structure of the analytic continuation of the Borel transform. In the case at hand this large-order behaviour has an oscillatory component which prevents the application of standard techniques used to extract information from it. Some methods of approaching this problem were presented in Ref. [19], but we have developed a new, powerful approach explained in detail in Appendix B, which we will apply below.

The singularities in the Borel plane
In general terms, the Borel transform corresponds to mapping a divergent series in the variable u 1 to a series in the Borel variable ξ via the map u −α → ξ α−1 Γ(α) . The series Φ n (see Eq. (2)) is mapped to This series has a finite radius of convergence and we can analyse the singularity structure of its analytic continuation 9 , which is necessary if we eventually wish to carry out the Borel summation procedure. Fig. 1 shows the poles of the analytic continuation of the Borel transform for the leading asymptotic sectors of our transseries: the hydrodynamic sector Φ 0 (top plot), the fundamental sectors Φ e 1 , Φ e 2 (middle plots) and the mixed sectors Φ 2e 1 and Φ e 1 +e 1 (bottom plots). The condensation of poles is taken to be indicative of a branch point singularity. We can directly check there that the positions of the inferred branch points are determined by the different exponential weights appearing in our transseries solution (3). We can see the appearance of the different fundamental sectors (shown in the figure as filled circles) as well as mixed sectors (shown as filled purple diamonds).
However, not all the mixed sectors appearing in our transseries (3) will contribute to the singularity structure of a given asymptotic sector. Applying resurgence techniques to our transseries, one can predict all the branch point singularities for the Borel transform of each sector. The procedure is thoroughly explained in section 5 of [21] (for multi-parameter transseries). The full list of the branch point singularities for each sector can be found in Appendix A, and these were accordingly marked in Fig. 1 in the form of circles (fundamental) or diamonds (mixed). Analysing the figure further, we note that there are several predicted branch points without any singular behaviour (grey dots). This is due to the fact that we are only taking a finite and limited number of terms from the original series and approximating the analytic continuation of the Borel transform by means of Padé approximants. Depending on the number of terms available the effectiveness with which we can detect branch point singularities in the Borel plane diminishes with distance as one moves away from the origin (or the empty circle in the shifted Figure 1: The singularity structure (grey dots) in the complex ξ-plane (Borel plane) associated to the Borel transforms of sectors Φ n , with n = 0, e 1 , e 2 , 2e 1 , e 1 + e 1 . This singularity structure is obtained via the poles of the respective Borel-Padé approximant (of order N ) BP N [Φ n ]. The Borel-Padé approximants used were: The Borel planes have been shifted to be represented on a common axis, with the expansion points shown by a hollow circle and the hydrodynamic expansion located at the origin (in dark yellow). The predicted branch points for each sector (listed in Appendix A) have been superimposed on each Borel plane and are divided into fundamental (solid circles) and mixed sectors (purple diamonds). The fundamental sectors include: A 1 and A 1 (blue), A 2 and A 2 (red), A 3 and A 3 (green). Any branch points associated to sectors of the transseries which do not give origin to branch points are shown as orange crosses. plots of Fig. 1).
The Borel plane singularities shown in Fig. 1 have an unexpected feature, particularly evident in the middle left plot (the Borel plane singularities of sector Φ e 1 ): the hydrodynamic series has a cut associated to it, at the origin of the plot. This is very surprising: typically, in a resurgent transseries, the perturbative sector has no transseries parameter associated with it (such as the σ i appearing in (3)). The appearance of the cut associated to the hydrodynamic series suggests otherwise. The resolution of this puzzle is very simple: the energy density as a function of proper time contains in fact a dimensionful parameter, which we have set to a convenient value (following [8]). This parameter is a vestige of the initial conditions and as such it leads to a branch-point singularity and a cut in the Borel plane. We have checked that there are no such cuts if instead of the energy density we consider the pressure anisotropy (defined by A = p T −p L /3 , see also Eq. (33)) as a function of τ E 1/4 whose hydrodynamic gradient expansion is known to be universal at late times, and so does not involve any integration constants. 10 Each of the branch cuts appearing in the analytically continued Borel transforms contains all the information pertaining to the sector associated with it. Naturally, one can expect this information to be encoded in the asymptotic behaviour of the the Borel transform of the hydrodynamic sector, B [Φ 0 ], as shown on the top plot of Fig. 1. In other words, the large-order (factorially divergent) behaviour of the coefficients of this series will encode information about all of the fundamental sectors and mixed sectors. This can be most efficiently and systematically studied through the tools of resurgence, which when applied to our transseries solution (3) give rise to the so-called large-order relations, which express precise relations between the coefficients of the different transseries sectors. These large-order relations allow us to very accurately check the consistency of the the transseries ansatz (3), including the presence of mixed sectors reflecting the non-linear QNM coupling. A comprehensive explanation of how to determine these relations can be found in [21] (sections 2 and 5).
Taking as an example the hydrodynamic series, the large-order behaviour of its coefficients, 10 From the independence of the pressure anisotropy on initial conditions we can apply resurgence techniques to determine the exact behaviour of the branch cut appearing at the origin of the middle left plot of Fig. 1: the asymptotic series associated to it is related but not exactly equal to the hydrodynamic series. To do a full resurgent analysis of the transseries describing the energy density, including the cut at the origin, one would have to recover the dependence on this dimensionful parameter. For the purposes of this work we will keep this parameter fixed, thus obtaining the solution in Eq. (3).
implied by the assumption that our transseries is resurgent, is schematically given by: Each line above corresponds to different exponentially suppressed contributions (top line: leading; middle line: first exponentially suppressed; third line: second exponentially suppressed). The proportionality constants S 0→n are the so-called Borel residues, which are well known functions of the Stokes constants associated with the problem (see [21] for the exact relations). 11 Naturally, if a particular sector Φ n does not appear as a singular branch point in the Borel plane of B [Φ 0 ], then it won't contribute to the large-order behaviour of ε (0) k , because the corresponding Borel residue vanishes, S 0→n = 0. The expansions χ 0→n (k), which depend solely on the coefficients of the sector Φ n and the leading power of the hydrodynamic expansion β 0 , are asymptotic expansions in k 1 and are generally defined as Expanding the above expression at large-order we thus obtain an asymptotic series, e.g.
Recalling that the coefficients of the hydrodynamic series are all real, from these large-order relations we can predict certain conjugation properties of the Borel residues. Consider, for instance, the leading contributions to the large-order behaviour (41): these come from the complex conjugate sectors Φ e 1 and Φ e 1 , whose contributions are of the same order. To sum these contributions into a real result, we need to have S 0→e 1 = −S 0→e 1 .
A question now arises: is it possible to retrieve information concerning the sub-leading, exponentially suppressed contributions, such as (42) and (43)? To reach these exponentially smaller terms, we need to subtract the first line (i.e. Eq. (42)) from the asymptotic behaviour of the original hydro series coefficients ε (0) k . However, as mentioned before, χ 0→n (k) is itself an asymptotic series and we need to perform a resummation of it for each value of k needed to carry out the subtraction. We do so via the so-calledÉcalle-Borel-Padé procedure, wherein we determine its Borel transform, analytically continue it via a Padé approximant BP N [χ 0→n ], and then perform the inverse transform 12 to obtain the summed result S 0 + χ 0→n (k): where we assumed > 0 small, and k ≥ 1 (we cannot use this procedure to determine the summation at k = 0). We can then define the subtracted coefficients as 13 The subtracted series defined by these coefficients k u −k is again asymptotic, and the singularity structure in the respective Borel plane can be seen in Fig. 2. As before, in blue are shown the singular points associated to the actions A 1 , A 1 , while in red we can see the singular branch points associated to the actions A 2 , A 2 . Comparing to the top plot of Fig. 1, it is evident that indeed the singularities associated to s = A 1 , A 1 have been subtracted, and now the leading singularities are at s = A 2 , A 2 .  Fig. 1 we can see that the singularity structure of the leading sector can be consistently removed from the coefficients of Φ 0 , and the information of the higher sectors can be recovered. The slight asymmetry in the figure above is due to the particular choice of lateral resummation where we have integrated along the real line from above.
One can now numerically check that the large-order behaviour of the subtracted coefficients δ 1 ε (0) k is just given by the second (leading) and third lines (sub-leading) (42) and (43). From explicitly examining the numerical coefficients, one can notice that Re δ 1 ε |2A 1 | k , which means we only need the real part of this summation procedure to analyse the contributions from the fundamental sectors Φ e 2 , Φ e 2 , appearing in (42). Given that the sectors Φ e 2 and Φ e 2 are complex conjugates, with β e 2 = β e 2 , then we again expect S 0→e 2 = −S 0→e 2 . This conjugacy relation between Borel residues should extend to every fundamental sector: The resurgent properties of our transseries also impose several constraints on Stokes constants and consequently on the Borel residues. Of particular interest to our analysis is the expected constraint S 0→n e k = (−1) n+1 (S 0→e k ) n (see [21] for more details). Given the conjugacy relations just mentioned between Borel residues of conjugate fundamental sectors, we obtain From direct gravity calculations we know that Φ 2e 1 and Φ 2e 1 are complex conjugates of each other. We could then worry that the large-order contribution coming from the third line (43) is not real, as one could have expected. However, the resummations from the previous lines (41) and (42) will give a non trivial imaginary contribution of the same order. 14 The above discussion was framed in the context of the hydrodynamic sector, but an analogous large-order analysis can be done for any sector of the transseries. The upshot of this analysis is that despite the appearance of an infinite number of fundamental QNM sectors one can really analyse each contribution individually due to the hierarchy of exponential damping. We have thus verified that our original transseries ansatz is consistent at this level. In the following section we will develop the techniques necessary for a quantitative and systematic testing of the resurgence relations between coefficients appearing in different sectors.

A Borel analysis of large-order relations
In many cases one can directly use the large-order relations to determine the coefficients of exponentially suppressed sectors from the values of the hydrodynamic gradient series using numerical acceleration methods such as Richardson transforms [19,37]. However, in the case of interest here the oscillatory nature of our problem, due to the presence of complex conjugate sectors with complex characteristic exponents β n , complicates the task exceedingly. This type of issue has already been encountered in the analysis of large-order behaviour of the hydrodynamic model of [18] which was studied in Ref. [19]. The equations considered there were designed specifically to model the leading QNMs of N = 4 SYM and they served as a testing ground for the analysis undertaken here. Due to the lack of suitable techniques Ref. [19] tested the large-order relations using values of the expansion coefficients appearing in the various sectors. The methods used there started from the large-order relations such as Eq. (41), and performed the resummation of the expansions χ 0→e 1 , χ 0→e 1 (using Eq. (46)). This way one obtains as the leading contributions to the large-order behaviour. In the above equation everything is known except the Borel residues S 0→e 1 , S 0→e 1 , which we know to be related by Eq. (48) due to the reality condition which our solution satisfies. We can therefore solve for the modulus and argument of the Borel residue S 0→e 1 . While the result formally depends on k, the large-order relations imply that it should saturate at large values of k, which can easily be verified, leading to an accurate estimate of the Borel residue. The result of such a procedure -implementing precisely the techniques of Ref. [19] but now using the series calculated for N = 4 SYM leads to the results shown in Fig. 3. This provides further support for our approach. While the arguments presented so far provide strong evidence that we are dealing with a resurgent transseries, we would still like to extract values of the coefficients appearing in the exponentially suppressed transseries sectors from the hydrodynamic gradient expansion alone. The difficulties discussed in the previous paragraph preclude the use of standard techniques, but nevertheless, an equivalent analysis can be carried out at the level of the respective Borel transforms. This new method is presented in some detail in Appendix B, but we briefly summarise it below for use in the following Sections.
We consider the hydrodynamic series, Φ 0 (u), whose coefficients ε (0) k have the asymptotic large-order behaviour shown in Eqs. (41) - (43). As we saw previously, this large-order behaviour is intimately linked to the existence of branch cut singularities on the Borel plane, and resurgence analysis explicitly reveals the analytic properties of this link. It is well known 15 that if we multiply our asymptotic series Φ 0 (u) by an overall factor u −β for a specific value of β such that the Borel transform removes the exact leading factorial growth of its coefficients, then this Borel transform will present a logarithmic cut at its leading singularity. For the particular case of the hydrodynamic sector and the leading singularity at ξ = A 1 (corresponding to the first large-order contribution of (41)) we can write This means that if we analyse each branch cut of the Borel plane separately, we can recover the coefficients of the sector associated to that branch cut (in this case Φ e 1 ). It is vital that we know the factorial growth corresponding to each of the branch points (as given in the large-order relations) with as high an accuracy as possible, as this determines the exact type of branch cut one finds in the Borel plane. This is necessary to be able to transform it into the logarithmic cut shown in Eq. (51) above. From the knowledge of Eq. (51), we can transform the logarithmic behaviour into a leading square root branch cut by further multiplying the original sector by u −1/2 This result can be demonstrated through a simple extension of an analogous proof contained in Section 4 of [21] (for the case of the quartic integral), see also Appendix B.
Various changes of variables have been used in the literature to more effectively analyse the Borel plane singularity structure, in particular through the use of conformal maps, see e.g. [58][59][60], as well as recent applications [61]. 16 To further simplify our results, we can now perform a change of variables effectively transforming the branch cut into a simple pole. To this end we define ξ = A 1 − (ζ − A 1 ) 2 . In terms of this new variable we can write the above equation From this expression we can see that by calculating the Borel transform B u βe 1 −1/2 Φ 0 (u) (ζ) and analysing its residue at ζ = A 1 we have direct access to the Borel residue S 0→e 1 . 17 Moreover, by subtracting the simple pole from the Borel transform: and multiplying everything by (ζ − A 1 ) −2 , we can once again use the residue theorem to predict the next coefficient ε of the nonhydro sector Φ e 1 . This procedure can be systematised to predict an arbitrary number of coefficients of the nonhydro sectors, as is shown in Appendix B.
In the following we will consider in some detail the singular behaviour of sectors Φ 0 , Φ e 1 , and how this behaviour is governed by the other fundamental sectors, as well as mixed sectors. We will observe the remarkable accuracy of the predictions obtained via resurgence from the analysis of the Borel planes. In Section 4 we look at leading contributions coming from fundamental sectors (associated directly with the AdS black-brane QNMs): we calculate the coefficients of fundamental sectors Φ e 1 and Φ e 2 from the Borel analysis of the hydrodynamic sector Φ 0 and compare them with the values obtained directly from the exponentially-suppressed contributions to the bulk gravity solution of Section 2. In Section 5 we will show that the mixed sectors, coming from non-linear couplings between QNMs, contribute non-trivially to the singular behaviour of the transseries sectors and cannot be ignored in a full description of the energy density. Our analysis will focus on the Borel analysis of sector Φ e 1 and the corresponding prediction of coefficients from sectors Φ 2e 1 and Φ e 1 +e 1 .
The results presented in Sections 4 and 5 provide significant further evidence that the energy density of N = 4 SYM can at late times be described by the resurgent transseries solution given in Eq. (3) and discussed in this section.

Resurgence of QNMs in the hydrodynamic expansion
We will now use the techniques explained in the previous section to analyse the resurgence relations between different sectors in the transseries shown in Eq. (3). We will focus on the hydrodynamic gradient expansion and the appearance of the QNMs in the large-order behaviour of its coefficients. The appearance of the these modes was already evident from the visual Borel 17 In fact we have access to the combination S 0→e1 ε (e1) 0 , but we normalise the fundamental sectors such that their first non-zero coefficient is 1. Figure 4: Comparison of the predicted results from resurgence techniques and the numerical ones obtained from solving Einstein equations. In the above plots we show this comparison for the coefficients of the fundamental sectors Φ e 1 and Φ e 2 , as predicted by the large-order behaviour of the coefficients of the perturbative sector Φ 0 (see Eq. (55)).
analysis shown in the top plot of Fig. 1. By applying the methodology introduced in the last section (and thoroughly explained in Appendix B) to the perturbative sector Φ 0 , we can predict the values of the coefficients of the sectors associated the leading (least damped) QNMs, associated to the transseries fundamental sectors Φ e 1 and Φ e 2 . In Fig. 4 we can find the normalised error of the resurgence prediction and numerical gravity calculation of the first coefficients associated to these two sectors. This error is given by where we defined ε (m) k | numerical as the coefficients associated to the sector Φ m numerically determined from the gravity calculation of Section 2, and ε (m) k | n−predicted being the same coefficients as predicted by the resurgence analysis of sector Φ n , which is outlined by Eqs. (53) and (54). The comparison of predicted and numerical coefficients for k ≥ 1 follows the calculation of the respective Borel residues S 0→e 1 and S 0→e 2 , as detailed below.

Behaviour of the hydrodynamic series at the leading Borel singularity
As mentioned earlier, the singularity structure of the Borel transform of the hydrodynamic (perturbative) series Φ 0 (u) can be found in the top plot of Fig. 1, where we have plotted the zeros of the denominator of the (diagonal) Borel-Padé approximant BP N [Φ 0 ], of order N = 189. We clearly observe pole condensation indicative of branch point singularities at ξ = A i , A i with i = 1, 2, 3 shown in blue, red and green (filled) circles, respectively. These are the branch points associated to the first fundamental sectors Φ e i , Φ e i with i = 1, 2, 3. The purple diamonds correspond to mixed sectors, which will be the subject of Section 5. 18 The two leading nonhydrodynamic, fundamental sectors appearing in this plot are Φ e 1 , Φ e 1 , responsible for the singular points ξ = A 1 , A 1 . We can then perform the Borel plane analysis introduced in the previous section and thoroughly described in Appendix B for any of these singular points. Consider for instance the sector associated to ξ = A 1 , given in Eq. (2) with the characteristic exponent β e 1 given on page 12. The leading growth of the coefficients of the hydrodynamic series ε (0) k associated to this sector is given by β = β 0 − β e 1 , as shown in the first line of the large-order relations Eq. (41). Then from a direct application of Eq. (53) we can determine the Borel residue associated with this singular point: The convergence plots of Fig. 3 show exactly this result. But we can now go further and analyse the sub-leading behaviour of the perturbative series at the singular point ξ = A 1 , as predicted by the procedure shown in Eq. (54). 19 The comparison between values predicted using resurgence and what was numerically calculated directly from the bulk gravity solution in sector Φ e 1 can be inspected on the left plot of Fig. 4 for the first 11 coefficients ε (e 1 ) m . We can see that the error remains extremely small, thus confirming the resurgent properties of our transseries solution (3).
As we have seen in the previous section, the Borel residue associated to S 0→e 1 corresponding with the singularity ξ = A 1 is related to S 0→e 1 by Eq. (48): This was directly checked by repeating the above calculation for the singularity ξ = A 1 .

Sub-leading exponentially suppressed behaviour of the hydrodynamic series
We have now verified the resurgent behaviour of the perturbative series in the nonhydrodynamic sector corresponding to the least damped QNM. We can ask if we can go further and analyse the next pair of complex-conjugate singularities of the hydrodynamic series at ξ = A 2 , A 2 , which correspond to the second QNM. The sectors associated with these singularities are Φ e 2 and Φ e 2 , respectively, given in Eq. (2) with the characteristic exponents given on page 12. As explained before, in order to analyse the resurgent behaviour at these singularities we need to first subtract the leading contributions coming from the first QNM. This can done by subtracting the resummed leading large-order behaviour (through a Borel-Padé-Écalle lateral summation), shown in (41), from the coefficients of the hydrodynamic series. We obtain the subtracted coefficients given in Eq. (47), which now obey the large-order relation given in Eqs. (42), (44): These large-order relations allow us to predict the coefficients ε (e 2 ) k , through the same Borel analysis applied before to the leading singularities, but now applied to the series The singularity structure on the Borel plane for this asymptotic series was already shown in Fig.  2, where we can clearly notice the absence of the leading singularities at ξ = A 1 , A 1 .
Note that we have an extra shift of the characteristic exponent β 0 : this was because we cannot define a resummation procedure for k = 0 to determine δ 1 ε (0) 0 . Nevertheless, the asymptotic properties are encoded in the large-orders, and taking this shift into account, we can easily determine the Borel residue associated to the second exponential weight (i.e. the second QNM): Also, from the direct analysis of the branch cut at ξ = A 2 we have checked that the relation between Borel residues S 0→e 2 = −S 0→e 2 is satisfied as expected.
We can now apply the method proposed in the previous section to predict the subleading resurgent behaviour of the hydrodynamic series' coefficients at the singular point ξ = A 2 . This corresponds to the resurgent prediction of the coefficients ε (e 2 ) m . 20 The right plot of Fig. 4 shows the comparison between prediction based on resurgence and the direct numerical calculation starting from the bulk gravity solution for the first 8 of these coefficients. We can see that the error remains small (although we have lost accuracy due to the resummation process), confirming once more the resurgent properties of the transseries solution (3).

Nonhydrodynamic sectors and effects of QNM coupling
In the previous section we have analysed the resurgence relations revealed by the asymptotic hydrodynamic expansion of the energy density and the exponentially suppressed fundamental sectors, directly associated with QNMs of the AdS black brane. This analysis already showed 20 These coefficients are explicitly predicted by Eq. (84). very strong evidence that the energy density has a representation as the resurgent transseries (3). However, a much more striking part of this transseries is the existence of mixed sectors, corresponding to non-trivial non-linear coupling between (different) QNMs, such as Φ 2e 1 or Φ e 1 +e 1 . The occurrence of these mixed sectors has already been seen in the gravity calculation of Section 2; here we show that they are indeed essential to complete the resurgent picture of our transseries for the energy density.
The easiest way to show the non-linear coupling between quasinormal modes is by studying the large-order behaviour of the coefficients of the fundamental sector associated to the first QNM Φ e 1 . The singularity structure of its Borel transform can be found on the middle left plot of Fig. 1. Here we already have direct evidence of the contribution of fundamental sectors such as Φ e 2 , Φ e 1 (via the branch cut starting at ξ + A 1 = A 2 and ξ + A 1 = A 1 , respectively), as well as the mixed sectors Φ 2e 1 and Φ e 1 +e 1 (via the branch cuts starting at ξ + A 1 = 2A 1 and ξ + A 1 = A 1 + A 1 ). The resurgent properties of these contributions can be explicitly written in the form of the large-order behaviour of coefficients ε (e 1 ) k (analogous to Eq. (44), see [21]): The values of the characteristic exponents appearing above can be found on page 12. The first line of Eq. (61) gives the leading growth of these coefficients, governed solely by the sector Φ e 2 (closest singularity on the Borel plane), while the second line shows the sub-leading exponentially suppressed growth governed by the mixed sectors Φ 2e 1 and Φ e 1 +e 1 . 21 Applying the methodology introduced in Section 3 (and explained in detail in Appendix B) to the coefficients of sector Φ e 1 we have predicted the values of the coefficients of the the fundamental sector Φ e 2 , as well as of the mixed sectors Φ 2e 1 , Φ e 1 +e 1 . The comparison between the resurgence predictions and the numerical gravity calculations can be found on the plots of Fig. 5, where we have used the notation for the comparison of coefficients introduced in Eq. (55). The following subsections discuss these developments in more detail.

Resurgence of different fundamental sectors
The first evidence of non-trivial relations between different QNMs comes from the fact that the leading large-order behaviour of the fundamental sector Φ e 1 associated to the least damped QNM is governed by the coefficients of the fundamental sector Φ e 2 , associated to the second Figure 5: Comparison of the predicted results from resurgence techniques and the numerical ones obtained from solving Einstein equations. In the above plots we show this comparison for the coefficients of the fundamental sector Φ e 2 and mixed sectors Φ 2e 1 and Φ e 1 +e 1 , as predicted by their relation to sector Φ e 1 (following Eq. (55)).
QNM. The precise relation can be analysed by focusing on the leading singularity ξ = A 2 − A 1 of the Borel transform of Φ e 1 . Applying an approach almost identical to Eq. (53) we can determine the Borel residue associated with this singular point: This Borel residue appears as the proportionality constant for the leading large-order growth shown in Eq. (61). Using the procedure shown in Eq. (54) we can further analyse the sub-leading large-order behaviour of the series Φ e 1 at the singular Borel plane point ξ = A 2 − A 1 and obtain a prediction of the coefficients of sector Φ e 2 . The comparison between the values of the first 11 coefficients ε (e 2 ) m predicted using resurgence and their values calculated from the bulk solution is found in the top plot of Fig. 5. The error observed there is very small, confirming the presence of nontrivial resurgence relations between different fundamental sectors, and consequently different QNMs. This provides further evidence for the correctness of the picture based on resurgence theory.

Mixed sectors
Until this moment we have dealt solely with the resurgent relations between fundamental sectors appearing in our transseries for the energy density, Eq. (3). However, the full transseries (3) also includes mixed sectors which do not correspond to single QNMs: they are in fact expansions in sectors whose exponential weight is given by linear combinations of the QNM frequencies as a result of the coupling between different fundamental sectors. The presence of such mixed sectors was already shown by the direct gravity calculation of Section 2, and it can be further justified via their resurgent relations to the fundamental sectors.
To determine the sub-leading, exponentially suppressed, behaviour of the Φ e 1 sector shown on the second line of Eq. (61), we want to analyse the singularities on the Borel plane at ξ + A 1 = 2A 1 and ξ + A 1 = A 1 + A 1 . As explained in Section 3 and applied in Section 4, we can do so by performing a resummation of the contributions from leading singularity ξ + A 1 = A 2 , shown on the first line of Eq. (61), and determine a subtracted series: where the subtracted coefficients are given by the laterally resummed series whose large-order behaviour is now given by the second line of Eq. (61). Through the analysis of the leading singularities of B [δ 1 Φ e 1 ] (ξ) at ξ + A 1 = 2A 1 , A 1 + A 1 , we can now determine the Borel residues associated with contributions of the mixed sectors Φ 2e 1 , Φ e 1 +e 1 to the large-order behaviour of Φ e 1 : S e 1 →2e 1 = 2S 0→e 1 ; S e 1 →e 1 +e 1 = S 0→e 1 = −S 0→e 1 .
These relations between Borel residues were expected as they are some of the many such relations which can be found from the direct resurgent analysis of the transseries in Eq. (3) (see Secs. 2 and 5 of the review [21] for more details). We can continue our procedure further and predict the coefficients ε (2e 1 ) k and ε (e 1 +e 1 ) k of mixed sectors Φ 2e 1 , Φ e 1 +e 1 , respectively. The two bottom plots of Fig. 5 show the comparison between resurgence prediction and numerical calculation for the first 8/9 of these coefficients. While the errors are somewhat larger here (with some loss of accuracy due to the resummation process), they are still small enough to be confident that the presented picture, including the presence of mixed sectors in the full solution for the energy density (3), is fully consistent with expectations based on the theory of resurgence.

Outlook
Through the AdS/CFT correspondence the problem of plasma equilibration can be viewed equivalently as the process by which a dynamical black object evolves toward its static state. The dissipative character of the black-brane horizon translates into the dissipative properties of the late-time asymptotic solution whose leading term was found in Ref. [15]. Some of the subleading terms were determined analytically in Refs. [46][47][48] and many more numerically in Refs. [8,16] through a power series in the variable u = τ 2/3 . In this paper we presented a novel transseries solution for the bulk geometry which supplements the asymptotic power-series by including contributions which are non-perturbative (exponentially damped) in the late time expansion. These additional transseries sectors can be interpreted either as linearised perturbations of the late-time bulk solution or as non-linear couplings of these perturbations. They become most relevant at early times and naturally encode information about the initial data, which is "lost" as the system evolves toward equilibrium. These perturbations can be identified with quasinormal modes of the static AdS black-brane solution, modified by gradient corrections. We were able to explicitly compute the bulk solution expanded around several transseries sectors to high orders in the late-time expansion. Since the AdS/CFT correspondence maps QNMs of the black brane to nonhydrodynamic modes of the dual gauge theory plasma, from the bulk solution we derived a transseries for the energy density of N = 4 SYM plasma undergoing Bjorken flow.
On the field theory side, our bulk solution provides an asymptotic form of the energy density as a function of proper time together with corrections which depend on the initial state. The physical meaning of this object is the same as at the level of hydrodynamics [9,19]. Even though the energy density depends on a scale parameter (which we have fixed as in Ref. [8]), it is independent of other features of the initial state, which enter only through the exponentially suppressed contributions. Thus, the transseries for the energy density has a similar interpretation as the transseries for quantities which are universal in the sense of Ref. [63] -it describes the dissipation of initial state information in the process of hydrodynamization [9,19,63,64].
From a mathematical perspective our results provide a novel example where resurgence is at work: the late-time expansion of a strongly coupled QFT computed in a microscopic theory. The intricate large-order relations which express the resurgent properties of our transseries solution provide strong links between the original hydrodynamic series and all the non-hydrodynamic sectors. The predictive power of these relations is such that one can in principle recover all nonhydro sectors directly from the large-order behaviour of the coefficients of the hydro sector (the limiting factor being the actual number of coefficients determined). One can ask how is this possible, given that the Borel plane of the hydro sector (top plot of Fig. 1) does not directly show all the mixed sectors? The answer is that this is an iterative process: the hydro series returns full information about all fundamental sectors, which in turn provide all the information about the mixed sectors which was previously not visible.
The asymptotic series which arise in this way show similarities with those arising in the context of the hydrodynamic model of Ref. [18], whose asymptotic behaviour was investigated in Ref. [19]. In particular, the large-order behaviour of these series involve multiple, competing factorial growths dependent of different complex parameters (exponential weights and characteristic exponents), which make it impossible to apply standard tools used in many recent studies of asymptotic series. We have however succeeded in developing an equivalent method for extracting information from asymptotic series directly at the level of the Borel plane singularities, which has made it possible to effectively calculate Borel residues (Stokes constants) and verify large-order relations. In Secs. 4 and 5 presented extremely accurate numerical evidence that both nonhydro fundamental and mixed sectors need to be included in the full solution in order to have a complete picture of the fluid's energy density. Although natural from the point of view of resurgent transseries, the existence of these mixed sectors, interpreted as the result of coupling between different QNMs, is completely novel.
Finally, we comment on the applicability of our approach to early time dynamics. Our transseries solution, once supplemented by initial data, provides a full non-perturbative description of the energy density of the fluid. Thus it naturally extends the regime of validity of the original hydrodynamic expansion. Although the transseries was formally defined for large-proper times, one can use a summation prescription on every (asymptotic) sector to obtain the energy density for smaller values of proper time. Naturally the smaller the proper time, the more nonhydro sectors one will need to consider to obtain an accurate result. Of possible summation prescriptions the most widely used in the so-called Borel-Padé-Écalle summation procedure (see [21]) as discussed, for example, in [16,18,29,[65][66][67][68]. Such results may help to shed light on the question which initial conditions correspond to the early-time attractor in N = 4 SYM [10,11]. Alternative summation procedures can also be adopted to analyse the early-time regime, such as trans-asymptotics and analytic transseries summation [69][70][71].
Finally, let us note that the explicit identification of nonhydrodynamic modes with the quasinormal mode spectrum suggests that a generalisation of the gradient expansion through the fluid-gravity duality [72] should exist. Such a solution would systematically map exponentially suppressed corrections in the bulk to nonhydrodynamic mode contributions in the dual Yang-Mills theory even without the imposition of the symmetries of Bjorken flow.

B.1 Behaviour of an asymptotic series at the first Borel singularity
We assume an asymptotic series expansion (for u 1) of the form whose coefficients F (0) n show large-order factorial growth. This growth will contain possibly multiple contributions, each of which growing as Γ(n + β) A −n−β for n 1 for some possibly complex parameters β 0 , β, A. The parameters β are the so-called characteristic exponents, defining the type of branch point singularities that we will find in the Borel plane, while A determines location of these singularities. In the context of the present paper we could just write ε n instead of F n here, but the actual method is clearly more general than the present context.
Given the factorial growth of the coefficients in the hydrodynamic (in general: perturbative) series in Eq. (66), it is convenient to define a one-parameter family of series 22 whose Borel transform (defined via the rule u −δ → ξ δ−1 /Γ(δ)) will have a singularity structure dependent on the choice of α (see Section 4 of [21] for more on this subject). In the vicinity of the singular point ξ = A the Borel transform for the particular case of α = 0 has the standard behaviour of simple resurgent functions: 23 In this equation S 0→1 is a a Borel residue, given by the first Stokes constant of the problem, and Φ (1) 0 is defined via Eq. (67) with This asymptotic series will be associated with the first nonhydrodynamic (in general: nonperturbative) sector appearing in our transseries (with |A 2 | > |A|) The key point of this procedure is that given the behaviour in Eq. (68), we can transform the logarithmic behaviour into a square root behaviour. To see this, it will be very convenient to define which, as we now proceed to show, has a square root branch point singularity at ξ = A.
We begin by observing that in the vicinity of that point one has The validity of this relation can be verified using a simple generalisation of the methods described in Section 4 of the review [21] for the case of the quartic integral. 24 On the other hand, by definition, for general values of β i this will have a complicated leading behaviour at ξ = 0. However, when Φ (0) and Φ (1)  where we also used that As anticipated, Eq. (75) has a square root branch point singularity. The next step involves changing the Borel plane variable in order to more easily analyse the behaviour around each singular branch point. Such changes of variables have been used in the 24 In general, one can obtain Eq. (68) from Eq. (72) by evaluating the semi derivative associated with B Φ (0) 1 2 (ξ) term by term, expanding the result about (ξ − A) and isolating the contributions that are linear in log (ξ − A). The reverse will also be true, so Eq. (72) will always hold. literature, in particular through the use of conformal maps, see e.g. [58][59][60], as well as recent applications [61]. In our case, it will be useful to change variables to transform the square root behaviour into a pole. This in turn will make it possible to extract information from Eq. (66) by calculating residues. The required change of variable is Note that the the point ξ = A corresponds to ζ = A. With this change of variables it follows from (75) that where The C It is now straightforward to see that evaluating residues of Eq. (71) will give us direct access to these coefficients, as well as the Borel residue (or, equivalently, the Stokes constant). To effectively perform this analysis, we need to use the hydrodynamic expansion coefficients to calculate the series Ψ (0) (ξ) and perform the change of variable in Eq. (77). Expanding the result around ζ 0 = A − √ A (which is the regular point of the Borel transform corresponding to ξ = 0) one finds Assuming that we have N terms for this series, we define an analytic continuation of this truncation using the Padé approximant BP N 2 Ψ (0) (ζ). By virtue of the arguments presented earlier, this quantity should have an isolated pole at ζ = A. This can be verified directly, and we can then calculate the residue at that point. The result can be compared with the value obtained from the asymptotic calculation (78), which leads to the conclusion that This gives explicit access to the Borel residue S 0→1 , since all the remaining quantities appearing in Eq. (82) are known (recall that we normalize the leading coefficients in the fundamental nonhydrodynamic sectors to 1). Once the Borel residue is known, we can gain access to higher coefficients F (1) m of the non-perturbative series Φ (1) (which first appear in C (0→1) m of (78)), by calculating the residue of Indeed this residue will give Res ζ=A B m Φ (0) = C (0→1) m .
The case m = 0, corresponds to the calculation of the Borel residue found in Eq. (82). With the value of this Borel residue one can predict the value of the coefficients F B.2 Behaviour of an asymptotic series at its second Borel singularity We have analysed the properties of the leading branch cut (closest to the origin) appearing in the Borel plane of the perturbative series. The position of the branch point was ξ = A. In order to differentiate between different branch point in the Borel plane, let us redefine A = A 1 . We would now like to analyse the behaviour around the the singularities in the Borel plane further away from the origin, such as the one at ξ = A 2 (where |A 2 | > |A 1 |). As it was explained in Section 3, in order to analyse the branch cut at this singularity, we need to first subtract the contributions from the closer one at ξ = A 1 . We saw that this can done by removing the leading large-order behaviour from the coefficients of the perturbative series. In much the same way as in Eqs. (41)- (43), this leading behaviour is given by (see [21] for further details): where h .
The sum appearing in χ 0→1 (n) is to be taken as a series in n −1 for large-order n, and is asymptotic. We can then perform a Borel-Padé-Écalle resummation for each value of n. 25 However, given that the singularities of the corresponding Borel transform also fall on the path of integration (the real axis), we will need to perform a lateral summation S 0 + χ 0→1 (n), as defined in Eq. (46) (for a Padé approximant of order N ). The subtracted coefficients obey new large-order relations governed by the coefficients F (2) k of the next asymptotic sector Φ (2) (u) in the transseries solution to our problem (70).
To check the large-order behaviour of these subtracted coefficients we perform the analysis described above, but now applied to Borel transform of the series and compare the residues Res ζ=A 2 B m δ 1 Φ (0) to the corresponding coefficients C (0→2) m . Note that the the resummation process can only be done for positive n, and we cannot determine the coefficient δ 1 F (0) 0 . Thus, the leading term of the subtracted series δ 1 Φ (0) (u) is now u −β 0 −1 .