The late to early time behaviour of an expanding plasma: hydrodynamisation from exponential asymptotics

We use exponential asymptotics to match the late time temperature evolution of an expanding, conformally invariant fluid to its early time behaviour. We show that the rich divergent transseries asymptotics at late times can be used to interpolate between the two regimes with exponential accuracy using the well-established methods of hyperasymptotics, Borel resummation and transasymptotics. This approach is generic and can be applied to any interpolation problem involving a local asymptotic transseries expansion as well as knowledge of the solution in a second region away from the expansion point. Moreover, we present global analytical properties of the solutions such as analytic approximations to the locations of the square-root branch points, exemplifying how the summed transseries contains within itself information about the observable in regions with different asymptotics.


Introduction
Viscous relativistic hydrodynamics is a long-wavelength effective theory which has been traditionally thought to be valid only near local thermal equilibrium. Surprisingly, hydrodynamic models can be successfully applied to certain physical systems which are far from equilibrium, such as an expanding quark-gluon plasma created from heavy-ion collisions at relativistic energies [1]- [4]. In those cases the hydrodynamic model contains within itself emergent, non-hydrodynamic degrees of freedom which are non-perturbative in nature and decay exponentially in time toward a hydrodynamic attractor [5]. This process is known as hydrodynamisation (see e.g. [6]). These non-hydrodynamic modes play a major role during the early times of the expanding plasma, and are quite sensitive to the different initial conditions. At the hydrodynamisation time, the system is still far from equilibrium and its pressure is quite anisotropic, but nevertheless the different initial solutions all become exponentially close to each other, and the evolution of the system towards equilibrium is effectively described by viscous hydrodynamics, via the same power series expansion in small gradients valid at late times.
From the point of view of asymptotics, however, such behaviour is expected. Mathematically, the latetime attractor is described by a divergent, asymptotic perturbative series, whose resurgent properties encode all the information about the non-perturbative modes. The information about the initial conditions is instead uniquely encoded in a set of parameters determining the strength of the (exponentially small) nonhydrodynamic modes. 1 The full description of the system can then be achieved via a so-called resurgent transseries.
Thus, when solving for the time-evolution of an observable (generally determined by some ODE/PDE within the hydrodynamic model), the initial-conditions (constrained by the physics) can dramatically change the behaviour at early-times by fixing the strength of the non-hydrodynamic modes dominating this earlytime regime, while at late times the non-hydrodynamic modes are exponentially suppressed and thus negligible in terms of their numerical magnitude, washing away the information on the initial conditions, with only the hydrodynamic power-law decay towards the attractor remaining.
Having access to the behaviour of our system at an initial time, as well as a description of its latetime asymptotic behaviour, one is naturally left with a few questions: how can we match the late-time behaviour to any given initial condition? Beyond a purely numerical analysis of the observable, how can we use this matching to describe the system at all times? Can we hope to describe the analytic behaviour of our observable?
The factorial divergent nature of the late-time expansions and their resurgent properties provides a path to answer these questions. Unlike many of their convergent counterparts, it is well known that these asymptotic expansions converge to their expected results quite quickly -in fact keeping just a few terms provides a very precise approximation, which can be extended far beyond the original expansion point, in our case large time (see e.g. [7]). Moreover, there are well established asymptotic summation methods, based on the underlying resurgent properties (see e.g. [8] and references therein), which provide such an approximation with exponential accuracy, thus effectively distinguishing between the exponentially close solutions at late-times [9]- [11]. As we will see, some of these methods even allow us to study global analytic properties of the asymptotic observable in its domain of interest, such as existence of poles or branch points [12]- [15] Naturally one should start with how to interpolate our late-time solution with a given initial condition. Unlike previous work discussing this matching in the context of relativistic hydrodynamics [16], [17], where the interpolation was done using a numerical least-square fit, our approach will involve various resummation methods based on the resurgent properties of the late-time solution. We will show that these resummation methods can be used to calculate the residual parameter labelling the exponentially close late-time solutions, being also highly effective at computing the solutions with exponential accuracy. Hence they are excellent approximations for most times outside of some region at very early times, where all orders of the nonhydrodynamic exponential modes are of notable size and drive the behaviour of the system.

A simple model of hydrodynamics
We will solve the interpolation problem between late and early times for an ODE describing the evolution of the effective temperature 2 of a conformal fluid in d " 4 dimensions undergoing a boost-invariant expansion. The model can be regarded as a toy model for the expansion of a strongly-coupled Quark-Gluon-Plasma created after a collision of two heavy ions beams. We assume rotational and translational invariance transverse to the collision axis. Further, we assume boost invariance with respect to boosts along the collision axis (Bjorken flow, see [18]), which is a reasonable approximation at high energies in the central rapidity region. Hence all observables in our system only depend on the proper time τ of some inertial observer and we may write T " T pτ q for the temperature. The energy momentum tensor of our system is given by where g µν " diagp´1, 1, 1, 1q µν is the flat Minkowski metric, E is the energy density in the rest frame of the fluid, ppEq " E{3 is the pressure of a perfect conformal fluid, the vector u µ is the four-velocity of the fluid, and Π µν is the shear-viscosity tensor. Conservation of energy and conformal symmetry imply (the symbol T in the third equation stands for temperature): The most straightforward approach towards solving (2) for the temperature T pτ q is to expand the shearstress tensor Π µν from (1) by summing all the allowed derivative terms up to a given order. However, the equations one obtains are not hyperbolic, and hence the model is acausal. An alternative way of dealing with (2) is to upgrade the shear-stress tensor Π µν to an independent field satisfying a relaxation-type differential equation. This approach, called Müller-Israel-Stewart (MIS) theory, [19]- [22] results in a causal model and is the one we will use in this work (see [5], [23], [24]). Instead of using the variables pτ, T q, where T stands for the temperature, it is more convenient to work with the variables pw, f q defined by 3 : The variable w measures proper time in units of inverse temperature, and the quantity f is closely related to the pressure anisotropy 4 . The differential equation describing the evolution of f pwq in MIS theory is The parameters A and β depend on phenomenological constants. If the microscopic theory behind a physical system is known, it can be used to derive these parameters. Our analysis can be performed in exactly the same way for any values of A and β. However, in this work we will only work with the following values for A and β 5 : Let us consider the solutions of Eq. (4). In the solutions plot Fig. 1 the real solution along the real axis are displayed. There are two distinct solutions, represented by the red and black curves in Fig. 1, which are finite at the origin. We will call the solution represented by the black curve f`, and the one represented by the red curve f´. The functions f`and f´are special solutions because they are attractors at w "`8, and f`{f´is the attractor/repellor at w "´8, respectively. This means that all other solutions, represented by green curves in Fig. 1, get exponentially close to either f`pwq or f´pwq as w Ñ`8. An important feature of the solutions is the presence of square root branch points, whose locations we shall denote by ws. It can be shown that Eq. (4) admits solutions of square root type, and admits the following analytic expansion in the variable pw´wsq 1{2 : 2Aws`2β´8 ws , h2pwsq " 16´2Aws 3ws ,¨¨¨. (6) The locations ws of the branch points depend on the initial conditions we impose on f pwq. The presence of these square root branch points implies that the natural domain of f pwq is a non-trivial Riemann surface. Note that the summation in (6) starts at n " 1, hence, those solutions are zero at the branch point. It is easy to check that the only possible regular zero for solutions of (4) is at w " p4´βq{A « 0.5, the point of intersection of red and blue curve in Fig. 1. Let us now analyse how the real solutions are related to each other by considering their expansions around the origin and around infinity.

Solutions around the origin
Around w " 0 we have the following convergent expansions: There is a relationship between f`pwq and fC pwq. In the solutions plot Fig. 1 the green curves above the graph of f`pwq (in black) correspond to fC pwq for C ą 0. We can see that as C becomes smaller, the use the term 'temperature' to refer to T . 3 Our definition of f pwq in (3) differs from the convention in [5] by fours " 3 2 f theirs . We nonetheless chose this normalisation because it leads to simpler equations. 4 The pressure anisotropy A is related to f by A " 8 pf´1q [25]. 5 Regarding our choice of parameters: The three phenomenological parameters defining the second-order transport coefficients which are relevant for the MIS dynamics are C τ Π , C λ 1 , and Cη. Assuming the microscopic theory is N " 4 SYM these parameters have been derived using holography and are given by [26], [27]: The ODE (4) is obtained by setting C λ 1 " 0, and identifying A " 3 2C τ Π and β " Cη C τ Π . We followed [5] and chose C λ 1 " 0 because this special case leads to a very interesting mathematical structure. The other two phenomenological parameters are chosen as , leading to (5). Note that Eq. (4) is only correct in the case C λ 1 " 0. green curves in Fig. 1 get closer to f`pwq. In the limit C Ñ 0`, fC pwq converges pointwise to f`pwq for w ‰ 0. Hence f`pwq can be understood as the C Ñ 0`limit of fC pwq, in which the divergence at the origin disappears. For 0 ą C ą C split «´0.0874, fC pwq has a square root branch point on the positive real axis. At C " C split this square root singularity splits into two singularities, one above and one below the real axis. The corresponding function fC split pwq is represented by the blue curve in Fig. 1. For C ă C split , the function fC pwq has no square root branch points along the real axis and admits a real solution for all w ą 0. These solutions are represented by the green curves in the bottom right corner of Fig. 1.

Solutions around infinity
Around w " 8 we also have two distinct transseries expansions, depending on whether the solutions converge or diverge at infinity.
(a) Solutions of finite limit: the solutions which converge to the hydrodynamic attractor f`pwq as w Ñ`8 can be described with the following transseries expansion [5]: The transseries F describes a one-parameter family of solutions converging to the finite value F Ñ 1 as w Ñ`8. Hence all the green curves in Fig. 1 which approach the black curve f`pwq have a transseries parameter σ assigned to them. The parameter σ is undetermined by the equation and has to be matched with the early-time behaviour of f pwq around w " 0, around which point we know all the solutions as convergent series expansions (7) and (8). We will determine later in the paper that the σ corresponding to the black/blue curve in Fig. 1 are approximately σ`"´0.3493`0.0027i and σ blue "´14.4111`0.0027i, respectively. The particular form of (9) implies that we know the amplitudes of all the non-perturbative modes once we know the transseries parameter σ. The expression Φ pnq pwq stands for the divergent, asymptotic series of the n-th non-perturbative sector or non-hydrodynamic mode: The coefficients a pnq k above can be determined from recurrence equations obtained by using the ansatz (9) into the MIS ODE (4), and matching equal powers of σ (see Appendix A). We use the convention a p1q 0 " 3{2. 6 The hydrodynamic series Φ p0q pwq " 1`β A w´1`O`w´2˘describes the perturbative sector and defines the attractor. Due to the factor e´n Aw multiplying the non-hydrodynamic series Φ pnq pwq, n ě 1, the convergence of the solutions to the attractor is exponentially fast. (b) Growing solutions: the solutions which are linearly growing to leading order and asymptotically approximate f´pwq as w Ñ`8 admit the following transseries expansion 7 Ψpp4, wq " w´k pp k`qk log wq`1 4 ÿ k"10 w´k`p k`qk log w`r k log 2 w˘`. . . , (11) with p´1 "´A{5. The first four coefficients, pn,´1 ď n ď 3, are uniquely determined by the MIS equation (4) alone. The coefficient p4 is undetermined by (4), and all other coefficients generically depend non-linearly on the coefficient p4. Hence the transseries Ψ from (11) represents a one-parameter family of solutions. The red curve in Fig. 1, that is f´, corresponds to p4 "´0.3474942558. It should be obvious from Fig. 1 that as w Ñ´8 the regular solutions have a transseries expansion of the form (11).
The exponential transseries (9) can be regarded as an expansion with a two-scale structure, the perturbative variable w´1, as well as an exponential variable Notice the form of the transseries (9): the outer sum is performed over powers of τ , with coefficients Φ pnq pwq depending on the variable w´1. In this work we will present different summation approaches of the asymptotic functions Φ pnq pwq. We will also explore an alternative way of summing (9) called transasymptotic summation, in which the order of summation is reversed [13]: the coefficients in the w´1-expansion are then functions of τ (defined by a convergent Taylor series in τ ). Thus the divergence of the transseries is not caused by the large-order behaviour of the exponential scales, but instead by the divergent asymptotic expansions at each order of the non-perturbative exponential. Although the transseries (9) was an expansion around w "`8, the transasymptotic approach allows us to access different regimes where τ is no longer small. 6 With this convention our Stokes constant and transseries-parameter normalisation is the same as in [5], [23], [24], and choosing a different value for a p1q 0 corresponds to a rescaling of σ. 7 Note that the transseries Ψ from Eq. (11) is constructed from the basis monomials w´1 and log w, whereas F is constructed from the basis monomials w´1 and e´A w .

Outline
In Section 2 we perform the exponentially accurate (error " e´2 |Aw| ) interpolation between late and early times with two different asymptotic methods: hyperasymptotics and Borel resummation. In particular, we show how to compute the transseries parameter σ from (9) with accuracy e´| Aw| from the 1-parameter family of of solutions at initial time. We explain how the matching function σpCq can be used to illustrate the convergence of the initial solutions fC pwq Ñ f`pwq as C Ñ 0 (see Fig. 2). This matching is performed at a chosen matching point, and the analytic continuation facpwq from the origin to the matching point is performed numerically using the Taylor series method. In Section 3 we introduce the transasymptotic summation and derive an asymptotic expansion for σ in closed-form at the matching point. Although the accuracy is worse with respect to hyperasymptotics and Borel resummation, it allows us to obtain analytic results which are useful far beyond the interpolation problem between late and early times and which can be employed to deduce global properties of the solutions. For example, one can derive an asymptotic formula for the location of square-root branch points, and explain the differing asymptotic expansions in two different regions of our complex domain pw Ñ˘8) as a direct consequence of the change in sign of the exponents logpτ n q "´nA w of our non-perturbative exponentials.

Interpolation between late-times and early-times
In the previous section we described the behaviour of the solutions to the MIS equation (4) both for earlyand late-times. We found that there exists a one-parameter family of solutions with a finite limit at infinity. These solutions converge exponentially fast to a hydrodynamic attractor described by a perturbative series. We saw that this series could be upgraded to the transseries (9) by including decaying exponential terms at large times. The transseries parameter σ from (9) was identified as a proxy for the amplitudes of the nonperturbative exponential modes. In the early-time regime near the origin, we found another representation of said family of solutions (8), labelled by the leading-order coefficient C of their Laurent expansion around the origin. Linking the magnitude of the non-perturbative modes of the late-time asymptotic transseries to the early-time behaviour can be very useful, and has previously been done by numerical fitting [17], [28]. However, the fitting approach does not exploit the vast possibilities arising from the rich late time asymptotics of the solutions. In particular, the difficulty with the fitting method lies in the exponential proximity of any two distinct solutions at late times, and hence a significant deviation from the desired solution is weakly penalised at late times, while at early times the function is not accurately captured by the fit model due to the finite truncation of the transseries (9). Fortunately, given that our solution at late times is divergent asymptotic, we have a range of tools at our disposal to do the matching, whose exponential accuracy provides a way to differentiate the behaviour of the different solutions. The matching between late and early times will be achieved in three main steps: (i) we will sum the factorial divergent expansion at late times, using exponentially accurate methods, keeping not only the perturbative series but also a non-perturbative, exponentially small part (effectively keeping the exponential accuracy). We will then evaluate this sum at a finite but large enough time w0.
(ii) we will analytically continue the solution at the origin to the same value w0.
(iii) the two approximations we will find depend on their respective parameters (C representing early times and σ late times) and their relation can be obtained via direct comparison.
After having found the transseries parameter for a given solution at early times, we can use the asymptotic summation methods to find exponentially accurate interpolations in the regime between early-times and infinity.

Hyperasymptotic summation
Hyperasymptotics is a resummation method which exploits the asymptotic properties of the transseries to approximate the value of a function by truncated sums [10], [11], [29], [30]. In computing our approximation for the solution f pw0q to the MIS ODE (4) at a finite 'matching time' w0 from the late time solution (transseries), we will keep terms up to linear order in the transseries parameter σ from (9). This corresponds to calculating level-one hyperasymptotics, for which we need to compute terms of the transseries sectors Φ p0q pwq and Φ p1q pwq from the transasymptotic summation (9) (see Appendix A for the computation). The optimal number of terms NHyp at which the series expansions arising in level-one hyperasymptotics must be truncated is a function of the resummation point w at which we wish to resum the transseries [31]: where t¨¨¨u is the usual floor function. Thus we need to compute the terms of the power series Φ p0q pwq and Φ p1q pwq to sufficiently high order (we used a maximum of 100 terms 8 for all our approximations). 8 Using 200 terms allows us to use the hyperasymptotic approximation with optimal precision up to w " 7.
The level-one hyperasymptotic summation is then given by [32] 9 fHyppw0q " fHyp,0pw0q`σ fHyp,1pw0q, where the hyperasymptotic summations for the perturbative sector and the first non-perturbative sector are given by the function F p1q in (15) is called hyperterminant and defined in terms of incomplete gamma functions via [33] The quantity S1 in (15) is called Stokes constant, which may be defined as the change in the transseries parameter σ from (9) upon crossing the Stokes line, which in our case is the positive real axis. The constant S1 has been calculated in previous work [23], [24] and is given by This Stokes constant can also be determined using hyperasymptotics, see Appendix B, where we give many more digits. Note that contributions of order O`σ 2˘a nd above in the transseries (9) are not included in level-one hyperasymptotics. The error in (15) is therefore of order e´2 |Aw 0 | [31].
In order to match the late time approximation with the early time solution, we need to bring our solution at early times (Eqs. (7) and (8)) to the finite value w0. This is done by analytical continuation with the numerical Taylor series method (see Appendix F). Let us denote the numerical approximation we obtain for f pw0q as 10 facpw0q :" numerical analytic continuation of f pwq from the origin, evaluated at the time w0 (18) By requiring facpw0q " fHyppw0q, we obtain the following approximation for σ: By decreasing the step size and increasing the order of the Taylor expansions in the calculation of f pw0q, we can achieve arbitrary accuracy, such that the error in the approximation (19) is determined by limitations of the hyperasymptotic approximation. Hence the parameter σ in (19) is accurate up to an error of order e´| Aw 0 | . Do note that the approximation for f pw0q from the late-time transseries solution can easily be extended to higher orders in the transseries parameter σ, by computing more non-hydrodynamic sectors Φ pnq pwq from (9). In Fig. 2 the results of the early-to-late-time matching C Ø σ are plotted for C ą 0. Our results are consistent with the observation in Section 1 that as pC, σq Ñ p0`, σ`"´0.349261`0.00273515iq, the solutions fC pwq converge pointwise to the solution f`pwq, which is finite at the origin. The function σpCq in Fig. 2 is roughly linear (left plot) except for a tiny region around the origin C " 0, where the converge toward σ`is very slow and is best visualised on a log-linear plot.

The Borel resummation
Another way of approximating f pw0q is through Borel resummation (see e.g. [8] for a review). For a series Φpwq " ř jě0 aj w´j, the Borel transform of Φpwq is given by 11 B rΦs pξq " a0 δpξq``8 9 f Hyp,0 pw 0 q is not the same as the level-0 hyperasymptotic approximation or optimal truncation, since the number of terms at which the series is truncated must be increased as more non-perturbative sectors are included in the calculation. 10 Note that facpw 0 q depends on which solution we pick around the origin from the set tf`, f´, f C |C P Cu. 11 As usual with Borel transforms, any finite number of powers w j , j ě 0 need to be addressed separately, see e.g. [8].  (7); For the imaginary part of σ we always have Impσq " Imp S1 2 q. The convergence σpCq Ñ σ`as C Ñ 0 shows the pointwise convergence f C pwq Ñ f`pwq.
We truncate the series in (20) after N0 terms 12 terms and calculate its Padé approximant BPN 0 rΦs, i.e. we approximate the resulting truncated sum by a rational function BPN 0 rΦs with a numerator/denominator of order tN0{2u.
The Borel-Padé resummation method then consists of taking the inverse Borel transform of BPN 0 rΦs, which is given by the Laplace transform The resurgence properties of the transseries (9) directly translates to the existence of singularities of the integrand BPN 0 rΦs pξq in Eq. (21) along the positive real axis -the Stokes line -and the singularities reflect the branch cuts of the Borel transform (20), starting at all ξ " nA, n P N, one for each exponential in our transseries. Thus the value of the resummation S N 0 ,θ Φpwq depends on the choice of the angle θ from the positive real axis. Although this ambiguity in the choice of integration contour gives rise to an imaginary contribution for each summed sector S N 0 ,θ Φ pnq pwq, there is a natural way of summing the resurgent transseries (9) such that the final result is unambiguous and real for real positive values of w: the so-called median summation [34]. To do so we pick a small negative angle θ "´ε ă 0 for the integration (21), and require the imaginary value of σ in the following way (see Appendix B for some more details): i Impσq " S1 2 .
We can now give an approximation for f pw0q to first order in the transseries parameter σ 13 : In analogy to (19), we arrive at the following expression for σ for the Borel resummation method: Notice that for both Eq. (19) and (24) we only went up to linear order in σ in the approximation of f pw0q. If we wanted to obtain more accurate results, we could have included higher powers of σ, which amounts to including extra exponential orders 14 . For the Borel summation method we would only need to numerically compute the integrals (21) for the higher-order hydrodynamic sectors Φ pnq pwq in (9), while the generalisation of the hyperasymptotic summation is a bit less straightforward. It can nonetheless be done, 12 We used N 0 " 100, which allows us to perform the Borel resummation with optimal accuracy up to w " 7. 13 We are using the transseries (9) and throwing away all the terms of order O`σ 2˘a nd above. 14 A similar matching was already done in [5] for the solution f`using Borel resummation with two exponentials. and we refer the reader to the literature [11], [30], [31], [35]. However, one can obtain the same accuracy if instead of increasing the number of exponentials/powers of sigma, we would just increase the value of the matching time w0.
Once the parameter σ has been matched to a given initial condition, 15 the transseries (9) can be used to find an approximation of f pwq everywhere, via some summation technique such as hyerasymptotics and Borel summation described above. The hyperasymptotic method does not require computing numerical integrals, but has the disadvantage of yielding discontinuous approximations to the summed transseries: it provides a piecewise analytic approximation (which is clear from the left plot of Fig. 3). On the other hand, the Borel summation integrals (21) must be computed as numerical approximations at each evaluation point, but the method has the advantage of giving a continuous function of w0.
In Fig. 3, we can see how different resummation methods compare with each other: in terms of accuracy the hyperasymptotic summation and the Borel resummation method are equivalent outside of a very small region near the origin, both giving an exponentially small error of approximately " e´2 |Aw| (the order of the first exponential we have neglected). We can also clearly see that the approximations given by each summation method are quite accurate at very early times even though we have only included a single exponential mode -to obtain accurate results at earlier times one would need to include further exponentials and their respective asymptotic expansions from (9). Also in Fig. 3 one can find results corresponding to a transasymptotic resummation, which will be discussed in the next Section 3. Let us also briefly mention the optimal truncation method, which consists of truncating the power series of the perturbative sector before the least term 16 foptpwq " The accuracy of the optimal truncation method is approximately " e´| Aw| , which agrees with our plots in Fig.3. Now that we have discussed how to perform the interpolation between late and early times using Borel resummation and hyperasymptotics, the next section will be devoted to the transasymptotic summation method.

Transasymptotic summation
We have seen in Section 2 that approximating the transseries (9) by keeping only the perturbative and the first non-perturbative sector gives excellent approximations of exponential accuracy for the function f pwq outside a small region near the origin. However, truncating the transseries in this way only works if the exponentials are small. Along the negative axis, the exponential monomial τ " e´A w defined in (12) grows arbitrarily large, and it is clear that truncating the transseries (9) can no longer work since all orders of τ contribute significantly towards the sum in that regime. This raises the question whether the transseries is of any use at all in regions where the exponential monomial is large enough. The answer is yes: we can exploit the fact that the divergent behaviour in the transseries comes only from large orders of the perturbative variable w´1, whereas the large order behaviour of the exponential variable τ is convergent. All we need to do is change the order of summation in (9): Frpτ qw´r; The coefficient functions Frpτ q are analytic at τ " 0 , and we will see that it is possible to systematically calculate them in closed form. This approach is called the transasymptotic summation [12], [13], and has been shown to be a powerful tool in the study of non-linear problems [14], [17], [36]. The special form of (26) allows us to compute the functions Frpτ q by treating τ and w as independent variables. Let us start with the lowest order approximation Then F0 obeys the ODE´1`F 0pτ q`1´τ F 1 0 pτ q˘" 0, which is solved by F0pτ q " 1`W p 3 2 τ q 17 , where W stands for the branch W0 of the Lambert-W function (see Appendix E). We can go further and calculate all Frpτ q. For r ě 1 we find the following differential 15 Value of the function at w " 0 for the solution f`p0q, or for the solutions which diverge at the origin f C pwq " Cw´4, the value of C is used as an initial value. 16 The formula for Nopt in (25) is a good approximation for the least term. 17 The general solution to (28) is F 0 pτ q " 1`W pcτ q. The integration constant c is found by matching the transasymptotic expansion to the transseries (9), and depends on the choice for a A`τ F0pτ qF 1 r pτ q``τ F 1 0 pτ q´1˘Frpτ q˘"p4´βqδr,1´8Fr´1pτ q`9´r 2 Note that in (29) all the derivative terms come multiplied by the variable τ , and that the variable τ does not appear other than as a multiplier of the derivatives. This motivates the convenient variable transformation τ Ñ W " W p 3 2 τ q. The derivatives transform as With the transformation (30) it is possible to rewrite the original recursive set of ODEs (29) and integrate them exactly. The details of this calculation as well as the method of fixing the integration constants are given in Appendix (C). It turns out that all the Fr are rational functions in W and can be computed exactly [37], [38]. Let us now see how the functions Frpτ q can be used to solve the interpolation problem between early and late times.

Interpolation with transasymptotics
We want to find an approximation for the transseries parameter σ corresponding to a given solution around the origin ((8) or (7)) using the transasymptotic summation (26). The first step of our approach is the same as in Section 2: we use numerical analytical continuation from the origin to the matching point w " w0, obtaining the numerical approximation facpw0q (see (18)). In a second step, we compute an approximation for τ pw0q, from which the transseries parameter σ can directly be calculated using our definition of τ pw0q, Eq. (12). The idea is the following: we want to solve for the function γpwq obeying Fr pγpwqq w´r " facpw0q " constant, for all w, which will be equal to τ pw0q when evaluated at the point w0, i.e. γpw0q " τ pw0q. The function γpwq satisfying (31) admits a perturbative, divergent asymptotic expansion in w´1: and determining γpwq will correspond to finding the coefficients γ k . Truncating the above expansion at its first term γpwq " γ0`O`w´1˘, we find from (31) that up to leading order Then also up to leading order in w´1 0 , we have γ0 " τ pw0q, which together with the definition of τ pwq (12) returns: Eq. (34) can be easily extended to higher orders in w´1 0 by including higher orders in the ansatz (32) and matching powers of w´1 in (31). The first four coefficients of the perturbative expansion of γpwq are given in Appendix D. The transasymptotic summation (26) can also be used to re-sum the transseries by truncating the series at the term of least magnitude. The difference with respect to the classical optimal truncation is that coefficients Frpτ pwqq vary with w. The result is displayed in Fig. 3. We can see that this approach outperforms optimal truncation. Note that we only calculated the coefficient functions Frpwq up to r " 15, and so the calculation is no longer optimal after the kink in the logarithmic error plot of Fig. 3. Furthermore, the kink happens at a higher value of w than we would expect from the resummation point w0 corresponding to 15 terms with classical optimal truncation given by (25).

Analytic results: branch points and global behaviour
Transasymptotics can be used to describe global properties of the function f pwq from (4), such as zeros, poles, branch points or to link distinct expansions in different asymptotic regimes. This is quite remarkable given that the transasymptotic summation was derived as a local expansion around the point w "`8. Let us start by sketching out how the locations of the branch points may be obtained. Notice that in the solutions plot Fig. 1 the locations ws of the square root branch points depend on the initial value problem that f pwq solves. From the perspective of late-time asymptotics, this means that the locations ws are a function of the transseries parameter σ. As already mentioned, all the coefficient functions Frpτ q in the transasymptotic summation (26) can be expressed as rational functions of the Lambert-W function W p 3 2 τ q, which has a square-root branch point at τ "´2 3 e´1. This branch point in the τ -plane translates to an infinite number of branch points in the w-plane if we substitute τ " τ pwq from (12). Since the Lambert-W function appears in all the coefficient functions in the transasymptotic summation (26), we expect the function f pwq to have an infinite number of square root branch points as well. The analytic information about the non-perturbative exponentials encoded in the coefficient functions Frpτ q can be used to provide an approximation for the locations ws.
Do note that all zeros of f pwq are square root branch point singularities (see the expansions in Eq. (6)) with the exception of a potential regular zero at w " p4´βq{A « 0.502. Hence we can solve for the branch points ws by solving the equation Fpw " ws, σq " 0 approximately for wspσq, where F is the transasymptotic summation from Eq. (26). We find where we have introduced the following variable: t " tpn, σq " logˆ3 σe 2A β˙`π ip1`2nq, n P Z; The integer n in (37) parameterises the sequence of branch points. Note that (36) is the partial sum of a divergent asymptotic expansion in t, and thus Eq. (36) is only a good approximation for the branch points/zeros of f pwq when |t| is large. In particular, w (approx) ptq becomes more accurate for large values of the discrete parameter n from (37), since the auxiliary variable t grows as an affine function with n. Since the leading order approximation w (approx) ptq from Eq. (36) grows linearly in t, the branch points which lie far from the origin are are best approximated by w (approx) ptq. Numerically, we can compute the zeros of f pzq by initially guessing the position of the branch point using (36) 18 and then using a contour integral to find a good approximation for the exact location. We start by  Figure 1: Active Stokes curves (the straight lines), higher order Stokes curve (the unit circle) for the asymptotic terms of (??) and the curve where |f 1 (z)| = |f 2 (z))| (the cardioid and teardrop). 1 Figure 4: Branch cuts of the f pwq for the case σ " 2 3 . Green dots: numerically computed branch-points; blue dots: approximations (36) to the locations of the branch points obtained from the transasymptotic summation of the late time solution Fpwq (9) (compare Table 1). Red dots: poles for the Padé approximant (about w " 1 2`5 2 i, shown as ‹) representing the branch cuts of the solution.
choosing a value for the transseries parameter σ and use the hyperasymptotic approximation Eq. (15) to find f pw0q (e.g. w0 " 10). We then analytically continue f pwq from w0 to a point in the vicinity of our prediction (36) using the Taylor series method, w1 " w (approx) s ptq`ε (e.g. ε " 0.3). Next we analytically continue again to compute the data on the circle |w´w (approx) s ptq| " ε. Using the trapezoidal rule [39] we evaluate the contour integral of wf 1 pwq f pwq to obtain the zeros of f pwq 19 . The approximate locations obtained with Eq. (36) as well as the numerical results are listed in Table 1, and plotted in Fig. 4. approx. (36) numerical n " 0 0.0975`0.6040i 0.1147`0.5076i n " 1 0.1555`1.416i 0.1580`1.384i n " 2 0.1817`2.276i 0.1827`2.257i n " 3 0.1991`3.143i 0.1997`3.129i n " 4 0.2122`4.012i 0.2125`4.001i Table 1: Approximations for the locations of the square-root branch points of (36) versus their numerically computed values for σ " 2 3 .
Let us now turn to another powerful application of transasymptotics: it can be used to correctly predict the different asymptotic behaviour of our solutions in separate regions. Consider the attractor f`in the solutions plot Fig. 1 (the black curve). At large, positive w, f`pwq converges to a finite value, while at large, negative w the same solution grows linearly with w. Therefore we have have two different asymptotic expansions, the transseries (9) at large positive w and the linearly growing expansion (11) at large negative w (which is also a transseries, but with log w-monomials instead of exponentials e´A w , see [40]). This is not surprising given the presence of square root branch points in the domain of our solutions. But it also raises an interesting question: can we somehow relate the two expansions to one another? The answer is yes: the great power of the transasymptotic approach lies in the possibility of analytically accessing regions in which the non-perturbative exponentials are no longer small. While the large, positive w limit corresponds to exponentially small values of τ " e´A w , the large, negative w limit is associated with exponentially large values of τ . Since the coefficient functions in the transasymptotic summation (26) are just rational functions of W pτ q, and the large τ expansion of W pτ q is known (see [41] and Appendix E), we were able to use transasymptotics to correctly derive the first four terms of the other, linearly growing expansion (11). The reason transasymptotics is so powerful in this case is that in flipping the sign w Ñ´w, the powers of w´1 in the transasymptotic summation do not change size, while the exponential variable τ " e´A w changes its regime, and becomes exponentially large instead of exponentially small. The details of the calculations in this section are beyond the scope of this publication and will be explored in an upcoming paper [42].

Summary/Discussion
The main focus of this work was to solve the problem of late time to early-time matching for arbitrary solutions of the ODE (4). We have a one-parameter family of solutions in two different regions of our domain: at late times, the ODE (4) admits formal transseries solutions consisting of the hydrodynamic perturbative sector as well as non-hydrodynamic sectors incorporating positive integer powers of the nonperturbative exponentials e´A w in the variable w representing time. In the early time regime near w " 0 there is a one-parameter family of divergent solutions which behave asymptotically as " w´4, as well as two finite solutions which are special limits of the one-parameter family (see solution plots Fig. 1). The different exponentially small contributions appearing at late times can be expected to be the leading contributions at early times. Beyond the MIS case, one expects to find similar transseries solutions in other hydrodynamic systems which observe a factorially divergent late time behaviour (see e.g. [43]) The resummation methods we used are hyperasymptotics, Borel-summation, and transasymptotics, and are all well-established. However, previous work did not exploit their strengths to do the parameter-matching and relied instead on less accurate procedures such as numerical least square fits [17], [28]. We carried out an analysis of said methods, and have shown that they are very effective tools for the parameter-matching. In terms of accuracy, the hyperasymptotic approach and the Borel resummation perform best. Both give an exponentially small error " e´2 |Aw| in the variable w.
The hyperasymptotic approximation has discontinuities since the number of terms which are included in the series varies with w, but requires no intricate numerical computations other than determining the series coefficients of the perturbative and first non-perturbative sectors.
On the other hand, the Borel resummation is a continuous function of w. Both resummation methods can be extended to include an arbitrary number of exponentials. Hence the method can be made arbitrarily accurate by increasing the number of non-perturbative modes we include in the approximation. However, the calculation of the Laplace transform (21) in going from the Borel-plane to the complex plane of our original variable w requires the numerical computation of an integral. As a consequence, Borel resummation is more computationally expensive than the hyperasymptotic summation, especially since said integral must be computed to exponential accuracy e´A w for the method to perform as well as the hyperasymptotic summation.
The transasymptotic summation has originally been used very effectively in the analysis of solutions of non-linear problems [13]. While the transasymptotic summation is less accurate in performing the interpolation, giving an error of " e´| Aw| (as opposed to e´2 |Aw| for the other methods), it is an extremely useful tool in the study of the global analytic properties of the solutions. The power of the transasymptotic approach lies in encoding the behaviour of the non-perturbative exponentials in analytic closed-form expressions, the transasymptotic coefficient functions. We have provided a systematic way of calculating these functions and used them to derive intricate analytic results such as analytic approximations to the locations of the square-root branch points as well as a way of linking distinct asymptotic expansions in two different regions of the domain to each other. These results have only been sketched out in this work, and a larger analysis will be presented in an upcoming paper [42].
The matching procedure we used is quite general and can be used beyond relativistic hydrodynamics. In fact, we can apply it to any interpolation problem between two different regions (e.g. late-time to earlytime, strong/weak coupling, large charge to small charge), where the solutions in one region are described by resurgent, asymptotic perturbative expansions, and where the behaviour in the other regime is known analytically (e.g. [17], [44]- [49]).
With the series ansatz we obtain the recurrence relations for Φ p0q and Φ p1q : ff , for j ě 2; The coefficient a p1q 0 is undetermined by (39), and any redefinition of a p1q 0 can be absorbed into the transseries parameter σ.

B The Stokes constant S 1 and median summation
An approximation for S1 relying on hyperasymptotics is given by [50]: We did compute S1 with (42) with an accuracy of O`10´6 5˘u sing N0 " 200. Eq. (42) requires knowledge of the coefficients of both the perturbative and the first non-perturbative sector. Note that it is possible to compute S1 without knowing the coefficients of the first non-perturbative sector using the so-called largeorder relations a p0q n " Γpn`βq A n`β S1´a p1q 0`O`n´1˘¯, as n Ñ 8 The leading order behaviour in (43) provides a sequence which converges to S1 as O`n´1˘and involves only the free coefficient a " 3{2, which we did so value of the Stokes constant is the same as in [5], [23], [24]. There is a connection between the value of the Stokes constant and the ambiguity in the value of the parameter σ. The positive real axis is a Stokes line, meaning that the Borel transform B rΦs pξq has branch-cut singularities at the locations ξ " A, 2A, 3A . . . . Therefore the definition (21) is ambiguous in the choice of angle θ.
When we move the integration path across the Stokes line from below and thus increase the angle θ in (21) from θ´"´ε to θ`"`ε we get a discontinuity in the result of the Borel resummation (21). Crossing the Stokes line in (21) while keeping the value of the transseries parameter σ from (9) constant corresponds to moving from one Riemann sheet to the other. Alternatively, we can alter the value of the transseries parameter as σ Ñ σ´S1 in order to cancel the discontinuity. We require the result of the resummation for the whole transseries (9) to be real-valued on the positive real axis, which is known as Median-resummation. The reality constraint fixes the imaginary part of the transseries parameter σ. Median resummation requires 20 i Impσq "˘S 1
If we choose a convention on the path along which we carry out the integration in (21) (below/above the real axis in the Borel plane), the only degree of freedom that is left is the real part of the parameter σ, which makes sense given that we have a one-parameter family of real solutions.

E Lambert-W function
The Lambert-W function (see [51, §4.13]) is defined as the solution to the equation The function W pzq has infinitely many branches, which are known as W k pzq, where k is an integer. Only two of those branches, W´1pzq and W0pzq, return real values on subsets of the real line. In the case of our problem, the MIS equation (4), the Lambert-W function appears in the context of the transasymptotic summation (26), where the leading-order contribution in w´1 is given by F0pτ pwqq " 1`Wˆ3 2 σw β e´A w˙.
As w Ñ`8 we require f pwq Ñ 1. This means that W p. . . q Ñ 0 in (57). For k ‰ 0 the branches W k pzq diverge as z Ñ 0. Therefore, we need to choose the branch W0 at w "`8, which admits the Taylor expansion W0pzq " z`. . . around z " 0 and is hence consistent with the behaviour of f pwq near w "`8.
For large arguments, the branch W0 admits the following expansion [41]: where L1 " log w ; L2 " logplog wq ; The expression Stir denotes Stirling circle numbers of the first kind. The presence of logarithmic terms in the expansion 58 explain how logarithmic terms arise in the transseries Ψ in (11) from the transseries F in (9) when going from w "`8 to w "´8. Note that the magnitude of the exponential scale τ " e´A w changes from small to large when the sign of w is flipped from p`q to p´q, which makes it necessary to use the expansion (58). Let us also note that the Lambert-W function has a square root branch point at z "´e´1.

F Taylor-series method
In the Taylor-series method (see [51, §3.7(ii)]) we combine, at a regular point w " w0, the Taylor series f pwq " ř 8 n"0 bn pw´w0q n with our differential equation (4) and obtain the recurrence relation w0pn`1qb0bn`1 "Aδn,1`pAw0`β´4qδn,0´1 2 w0pn`1q n ÿ m"1 bmbn`1´ḿ With this method it is very easy to 'walk' in the complex w-plane. Once we know b0 " f pw0q (either from a local expansion at the origin, or a branch-point, or from the asymptotic expansion) we can compute many coefficients in the Taylor-series expansion, and use this Taylor series to make a small step in the complex w-plane, that is, compute f pw0`hq and use this as the new b0.