Resurgence and Hydrodynamic Attractors in Gauss-Bonnet Holography

We study the convergence of the hydrodynamic series in the gravity dual of Gauss-Bonnet gravity in five dimensions with negative cosmological constant via holography. By imposing boost invariance symmetry, we find a solution to the Gauss-Bonnet equation of motion in inverse powers of the proper time, from which we can extract high order corrections to Bjorken flow for different values of the Gauss-Bonnet parameter $\lambda_{GB}$. As in all other known examples the gradient expansion is, at most, an asymptotic series which can be understood through applying the techniques of Borel-Pad\'e summation. As expected from the behaviour of the quasi-normal modes in the theory, we observe that the singularities in the Borel plane of this series show qualitative features that interpolate between the infinitely strong coupling limit of $\mathcal{N}=4$ Super Yang Mills theory and the expectation from kinetic theory. We further perform the Borel resummation to constrain the behaviour of hydrodynamic attractors beyond leading order in the hydrodynamic expansion. We find that for all values of $\lambda_{GB}$ considered, the convergence of different initial conditions to the resummation and its hydrodynamization occur at large and comparable values of the pressure anisotropy.


Introduction
In recent years there has been increasing interest in understanding the emergence of hydrodynamic behaviour in relativistic theories. In addition to many new theoretical advances, the strong multi-particle correlations observed in heavy ion collisions both at RHIC [1][2][3] and the LHC [4][5][6] and its successful description via hydrodynamical modelling [7][8][9][10][11][12] has provided a testing ground to explore how this collective behaviour arises from a microscopic theory. One of the most surprising empirical insights that this type of modelling of subnuclear dynamics has revealed is the fact that hydrodynamics can describe the bulk properties of the system even for extreme pressure gradients, which at face value, question the applicability of this long distance effective theory to those collisions. The recent observation of collective phenomena also in even smaller systems, such as p-p collisions [13][14][15] together with the success of the same hydrodynamic modelling in describing them [16][17][18] poses new challenges to our understanding of the applicability of this theory.
In conjunction with the above phenomenological observations, in recent years several numerical experiments have been performed to test the validity of hydrodynamics in different ultra-reltavisitic scenarios. Both in the infinitely strongly coupled limit of N = 4 SYM, described via holography, [19][20][21][22] and in the weakly coupled (perturbative) limit of gauge theories, described via kinetic theory [23], the direct comparison of the full stress tensor in different out of equilibrium processes with the hydrodynamic expectation showed that hydrodynamics can provide an accurate description of the evolution of the system even where the gradient terms are as large as the leading order terms. These experiments demonstrate that hydrodynamics can be applied even if the system under consideration is very far away from local thermal equilibrium and with strong deviations from the equation of state, as explicitly demonstrated in the analysis of non-conformal theories [24] .
Complementary to these numerical studies, the convergence of the hydrodynamic series has been recently analysed. In a seminal paper [25], the late time behaviour of boost invariant expansion of N = 4 SYM was analysed in a power series of the inverse proper time up to very high order. This series may be viewed as an expansion in the Knudsen number, and the coefficients of the series are controlled by increasing orders in the hydrodynamic expansion. The analysis of these large order perturbations showed that the hydrodynamic expansion is an asymptotic series, exhibiting a factorial growth of the series coefficients. Similar behaviour was found in different ultraviolet completions of second order hydrodynamics [26][27][28] and kinetic theory in the relaxation time approximation (RTA) [29] (see also [30] for a complementary analysis of the convergence of the hydrodynamic series). Interestingly, the analysis of these series via Borel-Padé techniques showed that these large order gradient expansions are sensitive to non-hydrodynamic modes, which play an equivalent role to non-perturbative corrections in perturbation theory.
Numerical analyses of these same boost invariant flows led Heller and Spaliński to suggest the existence of a hydrodynamic attractor [26] which may be thought as an extension of hydrodynamics beyond local thermal equilibrium [31]. These are special time dependent configurations to which all other boost invariant evolutions of the system converge at different times . These solutions have been found in different theories, such as Israel-Stewart (IS) [32] and Baier-Romatschke-Son-Starinets-Stephanov (BRSSS) [33] hydrodynamics [26,34], N = 4 SYM, kinetic theory [31] or anisotropic hydrodynamics [35]. It has been further argued that these special solutions may be found, or at least well approximated, via a trans-series solution of the system evolution [26,36], that non-perturbatively completes the gradient expansion. These types of solutions have also been recently analysed in less symmetric situations, including non-conformal theories and less symmetric flows [37,38]. Quite remarkably in all those cases the hydrodynamic attractor is very well approximated up to unexpectedly large values of the gradient by first order hydrodynamics, providing a dynamical understanding to the unexpected success of this theory.
The emergence of these solutions has only been studied in the two extreme cases of infinitely strong and perturbatively weak coupling. However, to better connect these theoretical advances with phenomenological applications it is important to understand how these special types of solutions behave at intermediate coupling. Starting from the infinite coupling limit, finite coupling corrections are studied via the gauge/gravity duality by introducing higher curvature terms in the dual gravitational theory [39,40]. For the gravity dual of N = 4 SYM, the relevant correction that affects the dynamics of the field theory stress tensor are expressed in terms of the Weyl tensor and are quartic in the curvature. This implies that those higher derivative terms can only be studied perturbatively, to avoid the emergence of instabilities, ghosts and other pathologies. Within this limit, these types of corrections have been vigorously studied in the past, exploring the correction of many different quantities (see [41] for a recent compilation of results). Recently, the relaxation of small stress tensor fluctuations in the thermal ensemble of N = 4 SYM has been studied in detail [42,43], (see also [44,45] for previous studies). In these studies, a new set of purely dissipative modes have been found, which are an intrinsic consequence of finite coupling.
Another higher-derivative theory which has also received significant attention in the holographic context is Gauss-Bonnet gravity. The action of this theory includes both quadratic and quartic curvature terms, governed by a single parameter λ GB ; nevertheless, as a Lovelock theory, its classical equations of motion contain only up to quadratic derivatives of the metric, which in principle allows the non-perturbative treatment of the high-derivative terms. Unfortunately, the dual field theory to Gauss-Bonnet gravity is unknown. Nevertheless, in this holographic construction it is possible to extract static (thermodynamic) and transport properties of the putative dual field theory, as well as the relaxation of small non-thermal perturbations [42,43], its off-equilibrium dynamics [46,47] and corrections to its hydrodynamic expansion [48]. The comparison of those analysis and the results obtained from finite coupling corrections of N = 4 SYM show that for negative λ GB these two theories share many common qualitative aspects, which makes Gauss-Bonnet holography an interesting laboratory with which to explore finite coupling effects to holographic theories. Note however that causality, positivity of energy and hyperbolicity considerations constrain the range of values of λ GB [49][50][51][52][53][54][55][56] (see [43] for a detailed discussion on these limitations).
As a first step towards intermediate coupling, in this paper we analyse the convergence of the hydrodynamic series in the gradient expansion of matter dual to Gauss-Bonnet holography. As in the previously mentioned examples, we find that the series is asymptotic as a consequence of non-perturbative (in gradient strength) contributions, given by the quasinormal modes of the dual black-hole. Remarkably, these characteristic modes dictate that the structure of singularities of the Borel transform associated to the hydrodynamic expansion interpolate between the known examples of infinitely strongly and weakly coupled theories. After characterising the analytic properties of the Borel transform, we resum the hydrodynamic series and constrain the dynamics of the hydrodynamic attractor. For all values of λ GB considered, the resummation approaches first order hydrodynamics even when gradient corrections to the stress tensor are large and at comparable values of the pressure anisotropy. Therefore, even though as the viscosity increases the approach to hydrodynamics occurs at a decreasing value of the gradient, this hydrodynamization processes takes place at comparable values of a viscosity-rescaled gradient, as suggested in [57]. Analysing the leading contributions in the trans-series, we estimate how close the resummation is to the attractor solution. By varying the amplitude of those non-perturbative corrections we observe that the convergence of different initial configurations towards the attractor occurs at values of the viscosity-rescaled gradient for which the resummation has hydrodynamized. We also perform an identical analysis of the hydrodynamic series of RTA kinetic theory [29] and find that in this model both hydrodynamization and convergence to the attractor solution occurs at smaller viscosity-rescaled gradients than in the strong coupling computations. This paper is organised as follows: in Section (2) we review the main properties of boost invariant flow and holographic Gauss-Bonnet to set up the notation, define the main quantities and outline the strategy to compute large order gradient corrections in this higher-curvature theory. In Section (3) we perform a Borel transformation of the hydro expansion and analyse the phenomenon of resurgence, i.e. the emergence of non-perturbative information in the perturbative series, in this context. After analytically continuing the Borel series, in Section (4) we resum the hydrodynamic expansion and constrain the dynamics of the attractor by analysing the leading order terms in the trans-series. Finally, in Section (5) we put our results into context and conclude.
2 Boost Invariant Flow in Holographic Gauss-Bonnet Gravity

Boost Invariant Hydrodynamics
Boost invariant flow, also known as Bjorken flow, is a particular class of solutions of ultrarelativistic hydrodynamics obtained by assuming that the dynamics of the fluid in 3+1 dimensions is independent of boost transformation in one of the space directions, z. This solution was first described by Bjorken [58] as a model for particle production in high-energy hadronic collisions. This solution was motivated by the observation that particle productions in this type of collision exhibits relatively flat momentum-rapidity distributions; while deviations of this independence are observed, boost invariant solutions of hydrodynamics are routinely used in the phenomenology of heavy ion collisions (see [59] for recent review).
Hydrodynamics may be viewed as an approximation of the stress tensor in a gradient expansion around the local thermal equilibrium at each point, characterised by the local energy density , and the fluid velocity U µ , where p( ) is the equation of state. For conformal theories, such as the ones we will describe in this paper, this is given by p( ) = 1 3 . The tensor Π µν encodes all the deviations of the stress tensor in a given state from local thermal equilibrium. In hydrodynamics, this tensor is expanded in space gradients 1 of the hydrodynamic fields as where η and ζ are the first order transport coefficients, known as shear and bulk viscosities, and σ µν is the shear tensor, constructed from the symmetrised and traceless space gradient of the velocity field. In this expression the ellipses indicate additional gradient corrections which, as at first order, can be expressed as a set of (independent) derivatives of the hydrodynamic fields multiplied by a set of transport coefficients which depend only on the local energy density. The complete set of second order coefficients and their operators can be found in [33,60] and [61] for conformal and non-conformal theories respectively. Following Bjorken, a solution to leading order hydrodynamics equation (Π µν = 0) can be easily found assuming that the hydrodynamic fields depend only on one spatial direction, z, and that the fluid is boost invariant. Under this assumption, independently of the microscopic symmetries of the theory, entropy conservation of ideal hydrodynamics imposes that the entropy density behaves as s 0 (τ ) ∝ 1/τ , with τ = √ t 2 − z 2 the proper time. Focusing on conformal theories further simplifications can be made; in particular, dimensional analysis dictates that the local temperature and energy density are given respectively by T (τ ) ∝ τ −1/3 , ∝ τ −4/3 . Similarly the size of the gradients can be expressed in terms the temperature T (τ ), the only scale that characterises the state at a given proper time τ . Therefore, the hydro expansion can be expressed as series of negative fractional powers of the proper time as with u = (Λτ ) −2/3 , Λ is a characteristic scale of mass dimension one, which in the hydrodynamic limit encodes all the information on the initial conditions of the flow and i are dimensionless constants which are solely dependent on the degrees of freedom and transport coefficients, and are therefore a property of the theory, as opposed to the state. In particular, the constant we may choose g * = eq /T 4 with eq the equilibrium energy density; the leading gradient correction is controlled by the shear viscosity and similarly all additional coefficients are a combination of the transport coefficients up to the corresponding order.
In the absence of transverse dynamics, boost invariance further imposes that the velocity field of the fluid U µ is given by U τ = 1. Another remarkable aspect of this effectively 1-Dimensional flow is that once the functional form of the energy density is determined, stress tensor conservation and conformal symmetry completely fix the stress tensor of the theory. Defining the longitudinal (P L ) and transverse (P T ) pressures as the diagonal components of the stress tensor in the fluid rest frame in the direction of expansion and perpendicular to it respectively, these are given by Note that these expressions are valid independently of whether the system is far from local thermal equilibrium, and that in thermal equilibrium both pressures must be the same. Therefore, we may gauge how far the system is from local thermal equilibrium by monitoring how anisotropic the system becomes with respect to the average pressure, which we can determine by computing the anisotropy parameter where the equilibrium pressure has been used to normalise the pressure anisotropy. The anisotropy parameter R can be used to monitor the approach of a system to hydrodynamics [19]. Following the hydrodynamic expansion Eq. (2.3), for a conformal theory the anisotropic parameter R becomes a function of the dimensionless gradient, where out of equilibrium we can identify the effective temperature as where we have defined g * below Eq. (2.3). In the hydrodynamic limit, R admits an expansion in inverse powers of w as R = n=1 r n w n , (2.8) where, once again, the leading order term is controlled by the viscosity to entropy ratio r 1 = 8 η s and, to first order in hydrodynamics, the anisotropy function is given by Finally, for later reference, another dimensionless quantity which has been used in the literature for identifying attractor solution is the logarithmic proper time derivative of w f ≡ τ w (2.10) The two function f and R are not independent of each other, since R(w) = −12 + 18f (w).

Gauss-Bonnet Holography and Boost Invariant flow
The Gauss-Bonnet gravity action is given by where κ 2 5 is proportional to the five dimensional Newton constant, L is the AdS radius of the λ GB = 0 theory, and λ GB is a dimensionless coupling which controls the magnitude of the higher derivative corrections. Without loss of generality, in the rest of this work we will set L = 1. An important feature of this theory is that in spite of incorporating higher derivatives terms, black-brane solutions, dual to thermal ensembles of the associated gauge theories, can be found for non-perturbarive values of the Gauss-Bonnet parameter. From these the equation of state of the dual field theory can be extracted [39,40] = 3 8 (2.12) As expected the λ GB → 0 limit of Eq. (2.12) agrees with the equation of state of N = 4 SYM. This expression indicates that not all values of λ GB are physical, since only λ GB ∈ (−∞, 1/4] yield real energy densities. As already mentioned in the introduction, for arbitrary values of λ GB this theory posseses causality problems associated with the superluminal propagation of high-momentum modes as well as negativity of the energy flux. [52,62]. These considerations impose further constraints in the allowed values 2 of −7/36 < λ GB ≤ 9/100. Nevertheless, since these constraints concern the ultraviolet behaviour of the theory, we may still consider values of λ GB beyond this region to explore its infrared dynamics, such as the approach towards hydrodynamics of the theory, as already done in [42,43]. Note also that as in a strongly expanding system as boost invariant flows, the early time dynamics, w 1 are sensitive to these high-momentum pathological modes. For this reason in this work we will not explore the Bjorken flow dynamics at arbitrariy early times.
In addition to the equation of state, the transport coefficients of the holographic dual have also been analysed. In particular, the ratio of shear viscosity to entropy density is given by [62] η Second order transport coefficients of this theory have also been analysed in [43]. From this expression we can observe that positive values of λ GB yield smaller values of η s than N = 4 SYM [49]. However, negative values of λ GB yield larger viscosity to entropy density ratios, as expected from finite t'Hooft coupling corrections of the infinite coupling limit in N = 4 SYM [63]. The analysis of the relaxation of small fluctuations of the thermal state via the computations of the quasi-normal mode spectrum of the dual black-branes indicates that many qualitative features of finite coupling corrections to N = 4 SYM are captured by Gauss-Bonnet holography with negative λ GB [42]. For this reason, in this paper we will only consider negative values of this parameter.
Numerical solutions of the Einstein equations (with no higher derivative corrections, λ GB = 0) with boost invariant symmetry from initial data at τ = 0 have been found in [19]. These solutions first showed the success of viscous hydrodynamics to describe the evolution of strongly coupled N = 4 SYM even when gradient corrections are large, later confirmed in less symmetric solutions. Boost invariant solutions of the Gauss-Bonnet equations of motion can also be found starting with this same ansatz. Imposing AdS asymptotics (with radius L c ) leads to the following near boundary (r → ∞) expansion of the different metric functions where the boundary value of the functions b and c are chosen such that the metric has AdS asymptotics with radius L c . The functions A (4) (τ ), b (4) (τ ) c (4) (τ ) cannot be determined from the near boundary expansion and additional infrared conditions, such as regularity, must be imposed. However, these are not all independent, since the power series expansion imposes Energy-momentum conservation, which emerges from the boundary expansion as well, relates these functions to A (4) (τ ) which must be extracted from the numerical computation. From these solutions the dual field theory stress tensor can be extracted after holographic renormalisation, which has been performed for Gauss-Bonnet gravity in [68,69]. In terms of those functions, the stress tensor is diagonal and with the conventions (2.19) This same ansatz has been used to obtain the holographic equivalent of a gradient expansion [25]. Motivated by the expansion of the energy densities in powers of u = τ −2/3 , Eq. (2.3), the different metric functions can be expanded in a power series in this variable. After introducing the new holographic coordinate s = 1/(rτ 1/3 ), a series solution of the Einstein equation can be found by expanding With this assumption, the solutions of Einstein's equations become a set of ordinary differential equations (ODE's) in s that can be solved by imposing AdS asymptotics at s → 0 and regularity at the horizon, which using reparametrisation invariance, can be set 3 at s = 1.
The equations of motions for Gauss-Bonnet gravity can also be solved by using the same expansion, which has been recently used to determined the first few (three) orders of the expansion in [48]. Imposing that AdS asymptotics and the the metric has a horizon at s = 1 the leading order solution is given by which coincides with the black-brane metric in Gauss-Bonnet gravity expressed in the Eddington-Finkelstein gauge 4 . Recalling the definition of s, this leading order solution may be interpreted as a black-brane falling in the holographic direction at a rate given by the temporal change of the temperature scale, as defined in Eq. (2.7). Starting from this solution, all higher orders can be computed by solving a set of subsequent linear ordinary differential equations; in Appendix A we describe a formal solution to all orders which can be used to organise the perturbative solution. From this set of differential equations, the energy density can be computed by analysing the s → 0 limit of the functions A i . Given the boundary expansion Eq. (2.15) and using that A 0 (s = 0) = 1/L 2 c , all metric coefficients A i for i ≥ 1 must vanish at the origin. Comparing the expression of the holographic stress tensor, Eq. (2.19) with the gradient expansion Eq. (2.20), the energy density of the dual field theory is given by Comparing this with the expression for the equation of state Eq. (2.12), we can determine the late proper time expansion of the effective temperature. Combining this expansion with the definition of the anisotropy function R, Eq. (2.5), we can use the series to determine the coefficients r i , as defined in Eq. (2.8).
Following the procedure outlined above, 240 orders in the gradient expansion of the energy density Eq. (2.3) for strongly coupled N = 4 SYM, were determined in [25]. In this work we have extended this computation up to 380 orders and extended it to fixed negative values of λ GB . Since the Gauss-Bonnet equations of motion contain many more terms, the computation of higher order expansion coefficients becomes much more numerically demanding than for N = 4 SYM. Similarly, the presence of a singularity, increasingly close to the horizon as λ GB becomes more negative also makes the numerical computation more challenging. The analysis in this paper is based on the determination of N coefficients = 94, 86, 80, 66 coefficients for λ GB = −0.1, −0.2, −0.5, −1 respectively. These coefficients, as defined in Eq. (2.3) and Eq. (2.26) can be found in the arXiv submission on this paper.

Resurgence
One of the main conclusions of the analysis of boost invariant flows is that the hydrodynamic expansion does not converge and behaves instead as an asymptotic series [25]. This conclusion is based on the factorial growth of the coefficients in the large order hydrodynamic series. As a consequence, for any fixed gradient strength, increasing orders in the gradient expansion lead to larger contributions to the hydrodynamic functions. This behaviour has been observed both in N = 4 strongly coupled SYM [25] and in kinetic theory in the RTA approximation [29], as well as in phenomenological completions of hydrodynamics [26]. As expected, this same behaviour is also observed in Gauss-Bonnet gravity. In Fig. (1) we show the growth of the magnitude of the series coefficients for the anisotropy function Eq. (2.8).
The factorial growth of the series coefficients indicates that the hydrodynamic series may be Borel summable. As is standard (see [70,71] for recent reviews in resurgence), we may define the series expansion of the Borel transform of the anisotropy function R, as Unlike the hydrodynamic series, since the leading n! growth has been removed, the Borel transform defined above has a finite radius of convergence, which is controlled by the asymptotic large n slope of the growth of the coefficients shown in Fig. (1). The slopes of these curves grow as λ GB becomes more negative, which in turns means that the radius of convergence of each series decreases with decreasing λ GB . From the point of view of finite coupling corrections this may be a natural expectation, since at finite coupling we expect the magnitude of the gradient corrections to grow. In the rest of this section we will explore the precise dynamics behind this observation. The finite radius of convergence indicates that the Borel transform of the anisotropic function possesses singularities in the complex-ξ plane. The Borel transform and the original series are related via a Laplace transform. In order to be able to perform this integral, we will need to first analytically continue the Borel transform beyond its radius of convergence. A standard method to do this is to approximate the Borel transform via a Padé approximant as where all N + M + 1 coefficients are fixed by demanding that the power series of the Padé agrees with Eq. For comparison, in the lower right panel we show the Borel plane for the same analysis 5 in kinetic theory within the RTA, using the coefficients tabulated in [29]. For all these cases additional poles exist for very large values of |ξ|; however these have a strong dependence on the Padé order, which indicates that they are numerical artefacts. The singularity structure of the Borel Plane is particularly interesting. As a first observation, the locations of the poles control the convergence of the Borel series, since the distance of the closest pole to the origin is proportional to the inverse of the slope of coefficients shown in Fig. (1). Furthermore, in all cases, the Padé approximant exibits an accumulation of alternating poles and zeroes, starting at a well defined point in the borel plane. This concentrated sum of simple poles indicates the emergence of a branch cut [73]. Nevertheless, the structure of poles at finite λ GB is qualitatively different to that of N = 4 SYM at infinite coupling. While in the latter case all poles are complex, for all finite λ GB new branch cuts emerge along the real axis. For small negative λ GB these new branch cuts are far from the origin, but as λ GB becomes more negative these poles move closer to ξ = 0, and eventually dominate the radius of convergence for the Borel Transform. This behaviour qualitatively interpolates the structure of the infinitely coupled limit of N = 4 SYM with the expectation from perturbation theory as obtained via kinetic theory in the RTA approximation.
We can note from Fig. (1) that for λ GB = 0 and λ GB = −0.1 the leading behaviour of the coefficients at large n follows the form of an oscillating factorial function (r n ∼ n! cos(an)), in a similar fashion to scenarios noted in [28,36]. As we vary λ GB to decreasing values we find that the oscillating behaviour becomes suppressed and the coefficients tend to follow r n ∼ n! as in [29]. This is consistent with the dominant contribution to the large n coefficients for the hydrodynamic series transitioning from two dispersive non-hydrodynamic modes, to a single dissipative non-hydrodynamic mode.
Having analytically continued the Borel transform beyond the power expansion, we can determine the anisotropy function beyond the power series via the inverse Borel transform where C is a contour in complex plane which connects ξ = 0 and ξ = ∞. The presence of singularities in R B shows that different choices of contour C yield different answers. Since we require this to be an analytic continuation for every complex value of ω or ξ, this implies that all choices of C must yield identical results. The theory of resurgence indicates that the anisotropy function cannot be simply approximated by a gradient expansion, but must  also incorporate non-perturbative contributions in the gradient strength. Denoting by ξ (α) 0 the origin of each independent branch cut in R B (each of which lead to an independent trans monomial) this trans-series is given by [25,28,72], 4) where N denotes the number of independent non-perturbative modes, and the functions Φ (n 1 n 2 ...,n N ) (w) admit power series in inverse powers of w at large w. The non-perturbative behaviour in the gradient strength is encoded in the functions Ω α , given by where γ (α) is constant for each branch cut which may be determined from the analysis of residues along the branch cut and C α are Stokes parameters, which must be chosen such that the non-perturbative ambiguity obtained in the Borel-summation of Φ 0 (w), is exactly cancelled by the next terms in the trans-series, yielding a real final result [28,70]. These parameters will jump discontinuously every time the contour C crosses a singularity in the Borel plane. However, this reality condition does not completely fix these complex numbers [28]. While in the hydrodynamic limit all the information about initial conditions reduces to an overall scale, the different values of these constants correspond to additional information on the initial state of the evolution, which controls the magnitude of the non-perturbative modes.
As noted in [25] the form of this trans-series coincides with the expected contribution of the evolution from non-hydrodynamic perturbations of the system away from local equilibrium in a boost invariant expanding medium of the equilibrium state. As it is well known, at strong coupling these non-hydrodynamic excitations are characterised by a set of characteristic complex frequencies, which in the dual theory coincided with the quasi-normal modes of the associated black-branes. In the adiabatic approximation, each of these these excitations relax according to the local relaxation where ω (α) is the characteristic frequency of each mode. Since the system under consideration is conformal, the τ -dependence of those frequencies is controlled by the effective local temperature. From the late time T-dependence, the position of the branch cut can be related the frequencies of the quasi-normal modes as ξ We have checked that even at finite λ GB , at leading order in gradients, the resulting ODEs become independent of γ (α) and coincide with the QNM equations of the static black-brane in [42], after the appropriate relation of ξ (α) 0 with the quasi-normal mode frequency 6 . From the above result, the observed qualitative differences between the Borel planes of N = 4 SYM and Gauss-Bonnet gravity can be traced back to the structure of quasi-normal modes. As noted in [42], this higher-derivative theory possesses a new set of dissipative (imaginary) quasi-normal modes in addition to the characteristic discrete complex modes of N = 4 SYM. These purely imaginary poles are not an artifact of this particular higherderivative theory. As explicitly shown in [42], the higher-derivative term responsible for finite coupling corrections to N = 4 SYM also lead to this new type of relaxation mode; and the t'Hooft coupling dependence of these poles is qualitatively similar to the λ GB dependence as long as λ GB is negative. We can therefore infer that the structure of the Borel plane singularities for N = 4 will also be qualititatively similar to the one observed in our analysis. It is rewarding to realise that these corrections seem to interpolate between the perturbatively weak and infinitely strong coupling limits.
To explicitly show the relation between the quasi-normal mode spectrum and the Borel plane singularities we show the positions of these characteristic frequencies, after an appropriate rescaling, in Fig. (2). In this figure the positions of the singularities associated with the first purely imaginary and complex QNM's (with smallest imaginary part, i. e. the smallest damping rate) are shown by the orange and red solid circles respectively. Note that from the relation above between the QNM frequency and the parameter ξ (α) 0 that these correspond to poles in the Borel plane with the smallest real part. In all panels, such poles coincide with the start of an accumulation of singularities in the Padé approximant, which may be interpreted as the origin of the branch cut. The singularities associated with higher QNM's must also be present in the Borel plane; in Fig. (2) we have shown the positions of these singularities associated with the second and third complex QNM's (with the next two smallest imaginary parts) by blue and green squares. Integer mutliples of all QNM frequencies above are given in squares of the same associated colour.
In all our finite coupling computations we do not identify poles coinciding with higher order modes. This however is likely an artefact of the limited number of coefficients we have been able to extract from our computations. For N = 4, where we are able to determine many more coefficients, these singularities indeed emerge, as already pointed out in [72]. Note that resonant singularities, associated to the product of trans-monomials in Eq. (3.4) are also visible; in Fig. (2), yellow triangles are used to identify the sum of the lowest two complex QNM frequencies (in the case of N = 4). We conclude this description by noting that these resonant singularities are also visible in the Borel plane of the RTA kinetic theory; this observation strengthens the significance of these structures, which have only been observed in the non-linear analysis of [29], but do not correspond to poles of the retarded correlator of stress tensor in the linear response analysis of RTA kinetic theory performed in [74].

Resummations and the Hydrodynamic Attractor
We now study the extension of the anisotropy function for small values of w by analysing the result of the inverse Borel transform, Eq. (3.3). As we stressed in the previous section the presence of poles in the Padé approximant, R B (ξ), will lead to an ambiguity of the inverse Borel transform, since depending on the choice of contour we will obtain different answers. This ambiguity can be fixed by demanding that the coefficients C α (Stokes parameters) will be discontinuous across each branch cut, as is known as Stokes phenomenon. However, as already mentioned, this procedure does not completely fix the value of these coefficients on a given contour, only its change across Stokes lines. The remaining ambiguity can only be determined from additional knowledge of the far-from-equilibrium early time dynamics of the system's evolution, since specifying different values is equivalent to selecting different choices of initial configurations for the evolution.
Among all the different possible time evolutions of the system, it has been recently proposed that there is a particular configurations which behave as an attractor in the space of initial conditions [26]. While currently there is not a precise definition of the attractor (see [38] for a recent attempt to provide such a definition based on the theory of non-autonomous dynamical systems), numerical analysis of different theories have identified well defined attractor solutions at all w, at which all time evolutions of the system converge. Given the previously mentioned difficulties in studying the very early time dynamics of this high-derivative gravity, in this section we will constraint the properties of the attractor solution in holographic Gauss-Bonnet by resumming the hydrodynamic series.
An obvious physical requirement for the choice of contour is that R must be real. In the case of N = 4 SYM this requirement is easily fulfilled by choosing, for example, the real axis as an integration contour, while setting all the coefficients C α to zero. This choice was recently analysed in [36], motivated by the results obtained in a hydrodynamic theory with a similar singularity structure in the Borel plane as that of N = 4 [28]. In that theory, the direct integration of the inverse Borel transform in the real axis coincided with the numerically computed attractor for of w > 0.3. We have tested that the additional coefficients computed in this work do not change the result of this resummation in this case. Recent analysis of exact solutions of IS hydrodynamics has shown [34] that in that model the attractor coincides with the direct resummation of the hydrodynamic series explicitly setting to zero all non-perturbative tails.
For finite λ GB , the presence of poles in the real axis complicates the extraction of the anisotropy function, since any contour that avoids those poles generically leads to an complex R. This feature is a clear manifestation of the need to include non-perturbative corrections in the form of a trans-series to determine the anisotropy function beyond the power expansion. trans-series corrections were indeed studied in [26] in the context of BRSSS hydrodynamics, which also exhibit real poles in the Borel plane. While the expansion Eq. (3.7) provides a clear starting point to complete this program in holographic theories, computing these corrections are numerically more challenging and goes beyond the scope of this work. Nevertheless, since these corrections are exponentially suppressed at large w we will use the leading term in the trans-series to constraint the dynamics of the resummation.
In this work we focus on the real part of the anisotropy function computed by performing the inverse Borel transform, Eq. (3.3) of the Padé approximant determined in the previous section. We choose a contour of integration given by a straight line in the complex plane z = ξe θ , with θ > 0 such that all the poles with positive imaginary part lie above the contour. For latter reference, we will refer to this contour as C + and it is shown in the upper right panel of Fig. (2). As already mentioned this choice leads to a complex R-value; however, as expected, at large w this imaginary part becomes increasingly small. Note that by choosing the real part we are making the computed R value independent of whether the integration is performed along C + or along an analogous contour C − obtained by reflexion along the real axis (see Fig. (2)), since the discontinuity across the real axis is purely imaginary. In essence, by this prescription we are effectively incorporating part of the first trans-series corrections. In fact, a procedure to determine these coefficients is precisely to choose the C α to cancel all imaginary contributions. However, this is not the full answer, since the functions Φ (n 1 ...,n N ) in Eq. (3.4) may contain additional imaginary parts which are cancelled only by higher order terms in the trans-series either associated to independent real poles further away from the origin (absent in Gauss-Bonnet) or by integer multiples of the leading real pole. Nevertheless, these contributions possess stronger exponential suppression factors, which make them only relevant at sufficiently small time. Note that this prescription implies, in particular, that all the trans-series coefficients C α associated to complex singularies are set to zero along this integration contour.
The results of this integration for all the different values of λ GB = 0, −0.1. −0.2, −0.5, −1 are given by the grey, green, blue, red and orange curves displayed in the left panel of Fig. (3). For all non-zero values of λ GB we have supplemented each curve with a band generated by adding and subtracting to the real part the imaginary part of the integral Eq. (3.3). When the band is narrow, this is a conservative estimate of the deviation of the trans-series from our prescription, since, as already argued, additional contributions are exponentially suppressed. As the width of the band increases, the sensitivity to the Padé order and the number of coefficients also increases. In the right panel of Fig. (3) we compare the results from our holographic computation to the resummation of the RTA kinetic theory coefficients from [29]. Following [57], to better compare the different theories we have rescaled the values of w by the viscosity to entropy density ratio, such that the first order hydrodynamic prediction R 1st hyd , shown by the dashed line, agrees in all theories by construction. Even though the RTA computation is performed with 200 coefficients of the hydrodynamics series, we find that the inverse Borel transform is much more sensitive to the Padé order, which prevents us from studying the very small ws/η regime.  The inspection of this figure shows that, after w is properly rescaled, the evolution of the anisotropy function is very similar, but not identical, for all cases considered, at least at sufficiently large values of w. All the resummations at fixed λ GB exhibit small imaginary parts for w < ∼ 4πη/s, a region where viscous corrections are large. Note also that the magnitude of the non-perturbaive corrections does not scale with the ratio of η/s, since the width of the different bands at similar values of the re-scaled variable is different. From the point of view of the gravitational dual, this is a consequence of the fact that the imaginary part of the non-hydrodynamic QNM's does not scale with the transport coefficients, at least for the values of λ GB considered. Similarly, the effect of non-perturbative modes is bigger for the RTA calculation, and significant deviations from the leading order term in the trans-series persist at w > ∼ 4πη/s. When the width of the bands is small, we can use the resummation to explore the process of hydrodynamization of the system. As already observed in [36] for N = 4, the result of these resummations quickly approaches the first order hydrodynamic predictions for all the values of λ GB . To better quantify this process, we will assume that the system has hydrodynamized at w hyd if for any larger value of w the anisotropy function satisfies where R hyd is the first order hydrodynamic expression Eq. (2.8). The values of w hyd and the corresponding anisotropy for the different theories are tabulated in Table ( Table 1. Inverse gradient size and anisotropy function at hydrodynamization for different theories. Note that for RTA the quoted range reflects the sensitivity of the resummation to Padé order and does not include the uncertainty associated with the imaginary part of the inverse Borel transform.
becomes more negative, the value of the temperature-normalised gradient w at which hydrodyninamization occurs increases, as expected by the fact that the dual fluid is more viscous. Nevertheless, as in other theories where the resummation has been performed [26,36], R hyd approximates the resummed result even when the value of this normalised gradient is comparable to the microscopic scale. At these small values of the inverse gradient, the anisotropy function is larger than 1, which means that the viscous contribution to the pressures is as large as the equilibrium pressure, demonstrating that the contribution of higher order terms is potentially large. This is once again a manifestation of hydrodynamization without isotropization as discussed in [19]. In fact, since the series is only asymptotic, it is easy to test that the corrections given by the truncated hydrodynamic series at orders greater than 10 give divergent and sign alternating contributions at those values of w. These corrections are, nevertheless, tamed by the resummation.
In this table we have also quoted the values obtained for RTA. As already mentioned, these results are much more sensitive to the Padé order and the extracted values reflect this sensitivity. This sensitivity hints towards a larger contribution of the non-perturbative corrections, which we will explore below in detail, making the hydrodynamization interpretation harder. Nevertheless, it is worth noting that all the computations performed via the gauge/gravity duality hydrodynamize at comparable values of the viscosity re-scaled gradient ws/η, and significantly earlier than in RTA kinetic theory. Since both RTA and λ GB may be viewed as oversimplified treatments of finite coupling effects in gauge theories, it would be interesting to investigate more realistic higher derivative corrections and collision kernels to explore whether the size of the re-scaled gradient at hydrodynamization shows consistent trends in these complementary approaches towards gauge theories at intermediate coupling.
We now turn to the relation between the resummation of the hydrodynamic series and the attractor. As we have stressed, in performing our resummations we have implicitly selected some particular values of the initial conditions, which tantamount to a specific selection of the constants C α in Eq. (3.4). It is therefore unclear whether this choice leads to the hydrodynamic attractor. In the simpler example of [26], where the trans-series program has been performed, non-trivial values for this constant, beyond the cancellation of imaginary parts, must be introduced (fitted) to describe the numerically computed attractor. Therefore, to fully determine the attractor numerical computations from an early initial time are needed. For N = 4, an attractor has been identified by studying the behaviour of different initial conditions at a very early initial proper time in [31]. As already stressed, performing these types of computations in Gauss-Bonnet holography introduced practical and conceptual difficulties which make them a challenge beyond the scope of this paper. For this reason, in this paper we will use information extracted from the resummation to constrain the position of the attractor.
To estimate the relaxation of different sets of initial conditions, we will focus on the dynamics of the leading non-perturbative corrections. From the point of view of holography, these may be interpreted as the effect of the least damped quasi-normal modes. As we have already mentioned, the w-dependence of these contributions could be obtained via the computation of a series expansion in gradients, analogous to Eq. (2.20), but supplemented with the non-perturbative prefactor Ω α for each mode. However, we can also determine the late time behaviour of this contribution by examining the discontinuities of the inverse Borel transform for different choices of the contour integration. By inspection of the Borel planes at finite λ GB , Fig. (2), we identify three representative contours of integration, which yield different answers for the inverse Borel transform. We have already used one of those contours, C + , to define the inverse Borel transform above. The second contour C − , is the reflection of the previous to the lower half plane. Finally, the third contour is a straight line in the upper half plane at angle above the argument of the start of the complex branch cut, C c . All these contours are shown in Fig. (2). Denoting by R + , R − and R c the results of integrating Eq. (3.3) over each of these contours, we define the discontinuities D ± is real and coincides with the imaginary part of R + while D c is complex. Note that we could have defined an equivalent discontinuity by reflecting both C c and C + to the lower half; however, this discontinuity is simply the complex conjugate of D c (w). Independence of the inverse Borel transform on the integration contour imposes that the constants C α must differ for each integration. Therefore, if the trans-series would contain only the contribution of the independent singularities ξ α with smallest damping rate, this contribution would be proportional to the discontinuity, so that the ambiguity could be cancelled. In more general cases the presence of additional non-perturbative modes as well as resonant contributions among poles imply that the cancellation is more subtle and the whole trans-series is necessary. For this reason, the discontinuities we have computed are not solely dependent on the leading singularity. But since those additional contributions occur at larger values of the conjugate variable ξ, the exponential suppression of these contributions at sufficiently large values of w, w > 1/ξ, is larger than those of the leading singularity, as inferred from the exponential contribution Ω α in Eq. (3.4).
Within the above approximation, the late time behaviour of the different initial conditions is given by where R is the result of the resummation described above and the coefficients a i contain information from the early proper time evolution of the system beyond the hydrodynamic approximation. The functional form of each of these discontinuities is shown in Appendix B.
Note that the edges of the band displayed in Fig. (3) corresponds to setting a 2 = a 3 = 0 and a 1 = ±1. Varying the values of this constant, we can estimate how different initial configurations deviate from our resummed result. In Fig. (4) we explore the effect of different initial conditions on the time evolution of boost invariant expansion of different theories. The discontinuities discussed before yield a characteristic magnitude of the size of the non-perturbative corrections needed to appropriately define the trans-series. We will vary the coefficients a i ∈ (−1, 1) to gauge the spread of typical initial conditions of such time evolution 7 . Following different extractions of attractor solutions in the literature, [26,31,37], the attractor may be identified by the so called "small roll" condition, which demands that the the time derivative f , with f defined in Eq. (2.10), is small all along the evolution of the system. For this reason, in Fig. (4) we show the logarithmic derivative of the energy density, τ ∂ τ ln , from which such a derivative may be inferred. In these plots, the solid thick line corresponds to the resummation, while the colourful thin lines correspond to different time evolutions obtained by varying a i . In all panels, the dashed line corresponds to the first order hydrodynamic prediction for this quantity. Finally, the vertical dotted line marks the hydrodynamization time extracted in Table (1) and the horizontal dotted lines indicate the values of τ ∂ τ ln which correspond to varying R by 10% around the resummation.
The inspection of this figure shows that in all holographic calculations the result of the resumation provides a good proxy to the attractor at hydrodyanamization time. For all values of λ GB , the variation within typical initial conditions of the evolution of the energy density is approximately within the hydrodynamization criterium used to determine the w hyd . This shows that our extraction of w hyd is trustable; in additions, assuming that the attractor is captured by this set of typical initial conditions, the attractor should also be approximated to better than 10% after this time. Furthermore, the spread of the different initial conditions also show that, while specific initial conditions may converge faster, for typical configurations, we only expect convergence towards the attractor when this special solution is well described by first order hydrodynamics. Note also that increasing the range of a i only make this conclusion stronger; for generic initial conditions the properties of the attractor prior to hydrodynamization do not strongly affect the time evolution of the system. We close this section by noticing that the RTA computation exhibits a much stronger dependence on initial conditions that the holographic computations 8 . For RTA not only our estimated hydrodynamization time occurs much later than for the holographic computations, but also at this larger value of the re-scale gradient the spread of initial conditions is large and many of the individual initial conditions are not well approximated by first order hydrodynamics. This implies that the relevance of the attractor for individual initial 7 We have checked that this procedure leads to a spread in anisotropy paramater comparable to that induced by the different initial conditions in N = 4 SYM reported in [36]. 8 We thank M. Heller, M. Spaliński and V. Svensson for private communication on recent analysis of RTA kinetic theory with a non-conformal relaxation time [75] which exhibits a trans-series structure with multiple independent contributions with identical exponential suppressions. This indicates that for conformal RTA the trans-series may be more complicated than what we have assumed in this paper.    Table (1) (including its uncertainties for RTA). Every curve corresponding to evolution in a Holographic theory is insensitive up to 2% to the choice of the Padé order N used for w > 0.3. For the case of N = 4 there is no visible change for the entire region plotted. All RTA curves plotted were insensitive to N at the level of 1%.
conditions becomes important at values of the gradient when the attractor is better approximated by hydrodynamics. As explicitly shown in Appendix B, the origin of this large spread is the ambiguity associated to the complex poles in the Borel plane, since these posses much larger residues than the real poles. Remarkably, the origin of these poles is unclear, since they do not appear in a linear analysis [29].

Discussion
Understanding the unexpected success of hydrodynamics to describe the off-equilibrium dynamics of interacting systems is an important challenge, not only theoretically but also with important practical applications to heavy ion physics and beyond. To address the success of this low energy effective theory much beyond its expected regime of validity, the emergence of special time-dependent configurations of the interacting theory, which act as attractors for all possible system evolutions and which generalise the hydrodynamic expansion beyond the limit of small gradients, has been suggested as a possible explanation. Motivated by this suggestion, in this paper we have applied the extension of hydrodynamics beyond the gradient expansion to the boost invariant flow of the field theory dual of Gauss-Bonnet gravity in 5D, which may be viewed as a laboratory to study finite coupling corrections to infinitely strongly coupled theories.
As we have already stressed, we have chosen to analyse Gauss-Bonnet holography since, at least in principle, it allows us to explore non-perturbative values of λ GB the parameter that controls higher-derivative corrections to Einstein gravity. We would like to remark once again that the holographic dual to this theory is unknown; and it is not even clear whether in that putative dual theory the terms included via non-vanishing λ GB correspond solely to t'Hooft coupling corrections or if they also include finite corrections in the rank of the gauge group N c . Nevertheless, the relaxation dynamics of non-hydrodynamic modes at finite (and negative) λ GB exhibits qualitative similarities to the effect of finite t'Hooft coupling corrections for those dynamics in N = 4 SYM. In particular, both theories exhibit purely dissipative relaxation channels which, from the point of view of holography, are due to higher curvature terms. Note also that our analysis has been performed for values of λ GB beyond the causality bounds of [52,62]. For this reason, we are not able to explore the boost invariant expansion in the τ → 0 limit. Nevertheless, since the unphysical behaviour of Gauss-Bonnet gravity occurs in the ultra-violate, we have concentrated our analysis around how different field configurations approach the hydrodynamic regime.
One of our main results is the analysis of singularities of the Borel transform of the hydrodynamic series in this theory. In accordance with the general theory of resurgence, and as already observed in N = 4 SYM, these singularities reflect the characteristic QNM frequencies that control the relaxation of small non-hydrodynamical excitations. A direct consequence of the new purely dissipative modes is the presence of singularities on the real ξ axes, the variable conjugate to the inverse gradient. These, together with the complex singularities associated to other QNM's of the dual theory make the analytical structure of the fixed coupling calculation richer than in the infinite coupling limit. But even more importantly, the structure of singularities qualitatively interpolates between the infinite coupling limit obtained via holography and the weakly coupled limit, obtained by kinetic theory in the RTA. This may be viewed as an additional motivation to study the large order gradient expansion in this higher-derivative theory.
To explore the effect of this analytic structure on the early time dynamics of the system, we have resummed the hydrodynamic series via Borel-Padé techniques. This allows us to extend the information in the large order gradient expansion to the large gradient region, for values of w such that the contribution of increasing orders in the gradient expansion lead to large, alternating contributions. Remarkably, as in all other examples studied in the literature, the resummation of the gradient expansion of the field theory dual to this high-derivative gravity is approximated by first order hydrodynamics at an unexpectedly early time, in a region where viscous effects are large. At this hydrodynamization time, the pressure anisotropy in the expansion is comparable in all strongly coupled computations, independent of λ GB , which implies that the hydrodynamization occurs at comparable viscosity-scaled gradients, ws/η. By comparison, our analysis of the RTA kinetic theory gradient expansion computed in [29] indicates that hydrodynamization occurs later, even in the viscosity-scaled gradient, at smaller values of the anisotropy parameter, although these are still large. Our results are consistent with the numerical solutions of RTA described in [29].
This resummation allows us to explore the dynamics of the hydrodynamic attractor in this holographic model. Certainly, resummation techniques cannot solely determine the behaviour of the attractor. To fully determine this configuration, analysis able to explore the w → 0 limit must be performed. However, at sufficiently late times, when all non-perturbative contributions have relaxed, the resummation must coincide with the attractor. To gauge how close the resummation is from the attractor, we have estimated the relaxation of transient behaviour by studying the discontinuities of the inverse Borel transform over different contours of integration. Since those discontinuities must be cancelled by non-perturbative contributions, these provide a natural scale for the magnitude of these corrections. By varying the magnitude of these modes we can gauge the deviation from the resumation of generic initial conditions. This procedure may be also understood as varying the contribution of the leading QNM over the evolving system. From this analysis we conclude that in all holographic computations, the expected deviation of generic initial conditions from the resummation at hydrodynamization time is comparable to the difference between the resummation and first order hydrodynamics. As a consequence, our resummation will be a good approximation to the attractor at hydrodynamization time; however, at earlier times this not may be the case 9 . We may therefore conclude that while individual configurations may converge to the attractor earlier, in all these strongly coupled computations the relaxation of generic initial conditions occurs whenever the system has hydrodynamized. Our analysis also suggests that the sensitivity of kinetic theory to initial conditions persists up to significantly smaller viscosity-rescaled gradients.
Finally, we would like to conclude with an intriguing observation. By analysing the magnitude of the discontinuities in different directions in the complex plane we can estimate the dominant source of initial data dependence in the late time transient behaviour. Surprisingly, for all the values of λ GB studied the dominant contribution is always associated with the complex QNM, which leads to complex singularities in the Borel Plane. The pure dissipative mode is always subleading, even for large negative values of λ GB , when the associated singularity is close to the origin and therefore does not possess an obvious suppression (see Eq. (3.4)). The numerically extracted discontinuities are shown in Appendix B. What is even more remakable is that an identical behaviour is observed in kinetic theory, where the dissipative poles are much closer to the origin that the complex ones. This is even more surprising after realising that in RTA it is only the pure dissipative mode that can be obtained from linear response, while the origin of the complex singularities is not yet understood. It would be interesting to explore the effects of this behaviour in other observables.

A Power Series Solution to all Orders
While constructing solutions for our bulk geometry we found that specific redefinitions of the metric coefficients allowed us to express our equations of motion as λ GB -independent linear operators sourced by λ GB -dependent functions. It is easy to see that these linear operators are in fact those of the λ GB = 0 case which have exact solutions in terms of Greens functions. We could not however find a closed form expression for each source for arbitrary order, and each solution generically has explicit dependence on solutions at all orders below it.
Starting from the ansatz given in Eq. (2.14) ds 2 = −r 2 A(τ, r)dτ 2 + 2drdτ + (rτ + 1) 2 e b(τ,r) dy 2 + r 2 e c(τ,r) dx 2 we make a further redefinition d(τ, r) = c(τ, r)+ 1 2 b(τ, r) and expand the unknown functions as power series' in u = τ −2/3 , where s = 1/(rτ 1/3 ). We can solve the Field equations of Gauss-Bonnet perturbatively at each order in u, for which we will find 3 e.o.m for A i (s), d i (s) and b i (s), and 2 constraint equations that we will evaluate at s = 1. The i = 0 solutions are given by a standard blackbrane metric solution stated in Eq.'s (2.23) to (2.25). For each order i ≥ 1 the functions A i (s), d i (s) and b i (s) must satisfy linear second order ODE's of the form where L f λ is a linear operator depending on λ GB which will act on function f i to give the source j f i . For all i ≥ 1 we impose that A i (s), b i (s) and d i (s) all vanish at the boundary (s = 0) and are regular at the horizon which we fix (through co-ordinate reparameterization invariance) to be at s = 1. A consequence of this choice of co-ordinates is that A i (1) = 0 for i ≥ 1 so that the constrain equations take the simple forms the equations of motion become where the linear operators L C 0 no longer have a dependence on λ GB and the source terms j C i are scaled asj The λ GB -independent linear operators are Here integration constant associated withb i has necessarily been fixed to 0 to ensure regularity of b i at the horizon.
One can retrieve all perturbative solutions iteratively by findingd i ,Ã i thenb i and substituting the results to find d i , A i and b i . 11 The first few sources are given bỹ 11 One can find bi through a an integral setting the integration constant such that bi(s = 0) = 0.

B Non-perturbative Modes
The discontinuity Re[D c ], Im[D c ] and D ± defined in Equation (4.2) are shown in Fig. (5) in blue, yellow and red respectively for all the models considered. For every case, we use several values of the Padé order to estimate the uncertainty of this integration. To avoid ambiguities associated with the oscillatory character of these functions, we estimate the error of each curve by computing the maximum value for different choices of the Padé order N = N − ∆ with N = N coefficients /2 and ∆ = 1, 2, 3 of the function (B.1) Using this procedure, we find that dispersive modes, Re[D c ] (blue) and Im[D c ] (yellow), deviate for no more than 5% for w λ GB = −0.1 and for RTA, convergence at 5% level is only achieved at w 4π η s > 0.85 and 2 respectively. Even though these uncertainties are large, since the dissipative mode is much smaller than the dispersive contribution, these errors do not alter the spread of initial conditions displayed in Fig. (4).
Remarkably, we find that in all cases the dissipative mode is suppressed by at least an order of magnitude relative to the dispersive modes, even when the poles along the real axis in the Borel plane are closer to the origin than the leading complex mode. This is surprising since inspection of the trans-series Eq. (3.4) suggests that the contribution of each of the non-perturbative modes to this discontinuity is controlled by the exponential suppression associated to the position of the corresponding singularity in the complex Borel Plane Eq. (3.5). This observation is even more striking for RTA, since the dissipative branch cut in this case is much closer to ξ = 0 than the complex one, and hence one would expect a larger suppression of those non-linear modes. This suggests that it will be necessity to understand the role of the complex modes in RTA kinetic theory to properly describe the evolution of the system at intermediate to late times.