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

We will solve the interpolation problem between late and early times for an ODE describing the evolution of the effective temperature ⁴ 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. The energy momentum tensor of our system is given by,

$$\begin{equation} T^{\mu \nu} = \mathcal{E} \, u^\mu u^\nu + (g^{\mu \nu} + u^\mu u^\nu) p(\mathcal{E}) + \Pi^{\mu \nu}. \end{equation} \tag{ 1 }$$

where $g^{\mu \nu} = \text{diag}(-1,1,1,1)^{\mu \nu}$ is the flat Minkowski metric, $\mathcal{E}$ is the energy density in the rest frame of the fluid, $p(\mathcal{E}) = \mathcal{E}/3$ is the pressure of a perfect conformal fluid, the vector u^µ is the four-velocity of the fluid, and $\Pi^{\mu \nu}$ is the shear-viscosity tensor. Conservation of energy and conformal symmetry imply:

$$\begin{equation} \nabla_\mu T^{\mu \nu} = 0; \quad \quad {T^\mu}_\mu =0; \quad \quad \mathcal{E} \sim T^{\,4} \,. \end{equation} \tag{ 2 }$$

The third equation above states that the energy density is proportional to the fourth power of temperature $T = T(\tau)$. The most straightforward approach towards solving (2) for the temperature $T(\tau)$ is to expand the shear-stress tensor $\Pi^{\mu \nu}$ 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 $\Pi^{\mu \nu}$ 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 $(\tau,T)$, where T stands for the temperature, it is more convenient to work with the variables (w, f) defined by ⁵ :

$$\begin{equation} w = T \tau; \quad \quad \quad f = \frac{3 \tau}{2 w}\frac{\mathrm{d} w}{\mathrm{d} \tau}. \end{equation} \tag{ 3 }$$

The variable w measures proper time in units of inverse temperature, and the quantity f is closely related to the pressure anisotropy ⁶ . The differential equation describing the evolution of f(w) in MIS theory is:

$$\begin{equation} w f(w) f^{\prime}(w) +4 f(w)^2+ f(w) \left( -8+A w \right)+\left(4-\beta-A w \right) =0 . \end{equation} \tag{ 4 }$$

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 use the values for A and β as in [5] ⁷ :

$$\begin{equation} A = \frac{3 \pi }{2-\log 2}; \quad \quad \quad \beta = \frac{1}{2 (2-\log 2)}. \end{equation} \tag{ 5 }$$

Let us consider the solutions of equation (4). In the solutions plot shown in figure 1 the real solutions along the real axis are displayed. There are two distinct solutions, represented by the red and black curves in figure 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 = +\infty$, and $f_+$ (respectively $f_-$) is the attractor (respectively repellor) at $w = -\infty$. This means that all other solutions, represented by green curves in figure 1, become exponentially close to either $f_+(w)$ or $f_-(w)$ as $w \to +\infty$. An important feature of the solutions is the presence of square root branch points, whose locations we shall denote by w_s. It can be shown that equation (4) admits solutions of square root type, and admits the following analytic expansion in powers of $(w-w_\text{s})^{1/2}$:

$$\begin{align} f(w)=\sum_{n=1}^\infty h_n(w_\text{s})\left(w-w_\text{s}\right)^{n/2}, \quad\textrm{with}~~ h_1^2(w_s)\!=\!\frac{2Aw_s+2\beta-8}{w_s},~~ h_2(w_s)\!=\!\frac{16-2Aw_s}{3w_s},~~\ldots.\nonumber\\ \end{align} \tag{ 6 }$$

Zoom In Zoom Out Reset image size

The locations w_s of the branch points depend on the initial conditions we impose on f(w). The presence of these square root branch points implies that the natural domain of f(w) is a non-trivial Riemann surface. Note that the summation in (6) starts at n = 1, hence, f(w) is zero at the branch points. It is easy to check that the only possible regular zero for solutions of (4) is at $w = (4-\beta)/A\approx 0.5$, the point of intersection of red and blue curve in figure 1. Let us now analyse how the real solutions are related to each other by considering their expansions around the origin and around infinity.

Around w = 0 we have the following convergent expansions:

• (a)  
  Two finite solutions at the origin:

  $$\begin{equation} f_\pm(w)= \left(1\pm\sqrt{\beta}/2\right)+\mathcal{O}\left(w\right) , \quad w \to 0. \end{equation} \tag{ 7 }$$
• (b)  
  A one-parameter family of solutions diverging at the origin

  $$\begin{equation} f_C(w) =C w^{-4} +2+ \mathcal{O}\left(w\right), \quad w \to 0. \end{equation} \tag{ 8 }$$

We shall now analyse the solutions $f_C(w)$ and their relation to $f_+(w)$. In the solutions plot of figure 1 the green curves above the graph of $f_+(w)$ (in black) correspond to $f_C(w)$ for C > 0. As C becomes smaller, the green curves in figure 1 approach $f_+(w)$. In the limit $C\to 0^+$, $f_C(w)$ converges pointwise to $f_+(w)$ for w ≠ 0. Hence $f_+(w)$ can be understood as the $C \to 0^+$ limit of $f_C(w)$, in which the divergence at the origin disappears. For $0\gt C\gt C_{\text{split}}\approx -0.0874$, $f_C(w)$ has a square root branch point on the positive real axis. At $C = C_\text{split}$ this square root singularity splits into two singularities, one above and one below the real axis. The corresponding function $f_{C_\text{split}}(w)$ is represented by the blue curve in figure 1. For $C\lt C_\text{split}$, the function $f_C(w)$ 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 figure 1.

Around $w\,=\,\infty$ we also have two distinct 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_+(w)$ as $w\to +\infty$ can be described with the following transseries expansion [5]:

  $$\begin{equation}\mathcal{F}(w,\sigma) = \sum_{n=0}^{\infty} \sigma^n w^{\beta n} \mathrm{e}^{-n A w} \Phi^{(n)}(w). \end{equation} \tag{ 9 }$$

  The transseries $\mathcal{F}$ describes a one-parameter family of solutions converging to the finite value $\mathcal{F} \to 1$ as $w \to +\infty$. Hence all the green curves in figure 1 which approach the black curve $f_+(w)$ have a transseries parameter σ assigned to them. The parameter σ is not determined by the equation and has to be matched with the early-time behaviour of f(w) around w = 0, where we know all the solutions as convergent series expansions (7) and (8). We will show that the value of σ corresponding to the black and blue curves in figure 1 are approximately $\sigma_+ = -0.3493 + 0.0027 \mathrm{i}$ and $\sigma_\textrm{blue} = -14.4111 + 0.0027 \mathrm{i}$, 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 $\Phi^{(n)}(w)$ is the divergent, asymptotic series of the nth non-perturbative sector (i.e. the nth non-hydrodynamic mode):

  $$\begin{equation} \Phi^{(n)}(w) = \sum_{k=0}^{\infty} a^{(n)}_k w^{-k}. \end{equation} \tag{ 10 }$$

  The coefficients $a_k^{(n)}$ 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^{(1)}_0 = 3/2$ ⁸ . The perturbative sector $\Phi^{(0)}(w) = 1 + \frac{\beta}{A}w^{-1} + \mathcal{O}\left(w^{-2}\right)$ describes the hydrodynamic series and defines the attractor at $w\rightarrow +\infty$. Due to the factor $\mathrm{e}^{-n A w}$ multiplying the non-hydrodynamic series $\Phi^{(n)}(w), \, n\geqslant 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_-(w)$ as $w \to +\infty$ admit the following transseries expansion ⁹

  $$\begin{align} \Psi(p_4,w) & = \sum_{k=-1}^{3} p_k \,w^{-k} + \sum_{k=4}^{9} w^{-k}\left( p_k +q_k \log w\right)\nonumber\\ & \quad +\sum_{k=10}^{14} w^{-k} \left( p_k +q_k \log w+r_k \log^2 w\right) +\dots , \end{align} \tag{ 11 }$$

  with $p_{-1} = -A /5$. The first four coefficients, $n = \, -1,\,\, 0,\, 1,\, 2,\, 3,\, 4$, are uniquely determined by the MIS equation (4) alone. The coefficient p₄ is not determined by (4), and all other coefficients generically depend non-linearly on the coefficient p₄. Hence the transseries Ψ from (11) represents a one-parameter family of solutions. The red curve in figure 1, that is $f_-$, corresponds to $p_4 = -0.3474942558$. In figure 1 it can be seen that as $w\to-\infty$ the regular solutions have a transseries expansion of the form (11).

The current work concerns the solutions of finite limit, also known as attractor solutions, given by the transseries (9). This transseries can be regarded as an expansion with a two-scale structure, the perturbative variable w⁻¹, as well as an exponential variable:

$$\begin{equation} \tau \equiv \sigma\, w^{\beta} \mathrm{e}^{-Aw}. \end{equation} \tag{ 12 }$$

In general, the transseries (9) is formally written such that the outer sum is performed over powers of τ, with coefficients $\Phi^{(n)}(w)$ depending on the variable w⁻¹. In this work we will present different summation approaches of the asymptotic functions $\Phi^{(n)}(w)$. We will also explore an alternative way of summing the transseries (9) called transasymptotic summation, in which the order of summation is reversed [13], such that the coefficients in the w⁻¹-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) is an expansion around $w = +\infty$, the transasymptotic approach allows us to access different regimes where τ is no longer small.

In section 2 we perform the exponentially accurate ($\text{error} \sim \mathrm{e}^{-2|A w|}$ ) 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 $~ \mathrm{e}^{- |A w|}$ from the 1-parameter family of of solutions at initial time. We explain how the matching function $\sigma(C)$ can be used to illustrate the convergence of the initial solutions $f_C(w) \to f_+(w)$ as C → 0 (see figure 2). This matching is performed at a chosen matching point, and the analytic continuation $f_\text{ac}(w)$ 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 the of square-root branch points, and explain the differing asymptotic expansions in two different regions of our complex domain $(w\rightarrow \pm\infty$) as a direct consequence of the change in sign of the exponents $\log(\tau^n) \sim - n A \, w$ of our non-perturbative exponentials.

Zoom In Zoom Out Reset image size

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 order to compute an approximation for the solution $f(w_0)$ to the MIS ODE (4) at a finite 'matching time' w₀, we will keep terms up to linear order in the parameter σ from late time transseries solution (9). This corresponds to calculating level-one hyperasymptotics, for which we need to compute terms of the transseries sectors $\Phi^{(0)}(w)$ and $\Phi^{(1)}(w)$ from (9) (see appendix A for the computation). In level-one hyperasymptotics, the optimal number of terms N_Hyp at which the series expansions must be truncated is a function of the resummation point w at which we wish to resum the transseries [31]:

$$\begin{equation} N_\text{Hyp}(w) = 2 \left \lfloor \big| A w\big| \right \rfloor , \end{equation} \tag{ 13 }$$

where $\left \lfloor \cdots \right \rfloor$ is the floor function. Thus we need to compute the terms of the power series $\Phi^{(0)}(w)$ and $\Phi^{(1)}(w)$ to sufficiently high order (we used a maximum of 100 terms ¹⁰ for all our approximations).

The level-one hyperasymptotic summation is then given by [32] ¹¹

$$\begin{equation} f_\text{Hyp}(w_0) = f_{\text{Hyp},0}(w_0) + \sigma\,f_{\text{Hyp},1}(w_0). \end{equation} \tag{ 14 }$$

where the hyperasymptotic summations for the perturbative sector and the first non-perturbative sector are given by,

$$\begin{align} f_{\text{Hyp},0}(w_0) &= \sum_{m=0}^{N_\text{Hyp}(w_0)-1} a^{(0)}_m w_0^{-m} \nonumber\\ & \quad + w_0^{1-N_\text{Hyp}(w_0)} \frac{S_1}{2 \pi \mathrm{i}} \sum_{m=0}^{N_\text{Hyp}(w_0)/2-1} a_m^{(1)} F^{(1)}\left(w_0;{\genfrac{}{}{0pt}{}{N_\text{Hyp}(w_0)+\beta-m}{-A}}\right); \nonumber\\ f_{\text{Hyp},1}(w_0)& = \mathrm{e}^{-A w_0} w_0^\beta \sum_{m=0}^{N_\text{Hyp}(w_0)/2-1} a^{(1)}_m w_0^{-m}. \end{align} \tag{ 15 }$$

The function $F^{(1)}$ in (15) is called hyperterminant and defined in terms of incomplete gamma functions via [33]:

$$\begin{equation} F^{(1)}\left(w;{\genfrac{}{}{0pt}{}{M}{a}}\right) = \mathrm{e}^{a w+i \pi M} w^{M-1}\, \Gamma (M) \Gamma(1-M,a w). \end{equation} \tag{ 16 }$$

The quantity S₁ in (15) is the so-called Stokes constant, which may be defined as the change in the transseries parameter σ of (9) upon crossing the Stokes line, which in our case is the positive real axis. The constant S₁ has been calculated in previous work [23, 24] and is given by,

$$\begin{equation} S_1 \approx 5.4703\times 10^{-3}\,\mathrm{i}. \end{equation} \tag{ 17 }$$

This Stokes constant can also be determined using hyperasymptotics, see appendix B, where we also provide S₁ with more precision. Note that contributions of order $\mathcal{O}\left(\sigma^2\right)$ and above in the transseries (9) are not included in level-one hyperasymptotics. The error in (15) is therefore of order $\mathrm{e}^{-2 |A w_0|}$ [31].

In order to match the late time approximation with the early time solution, we need to bring our solution at early times (equations (7) and (8)) to the finite value w₀. 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(w_0)$ as ¹²

$$\begin{align} f_{\mathrm{ac}}(w_0) : = \mathrm{numerical}\, \mathrm{analytic}\, \mathrm{continuation}\, \mathrm{of} \,\,\, f(w)\, \mathrm{from}\,\, \mathrm{the}\, \mathrm{origin},\,\, \mathrm{evaluated}\, \mathrm{at}\, \mathrm{the}\, \mathrm{time}\, w_{0} .\nonumber\\ \end{align} \tag{ 18 }$$

By requiring $f_\text{ac}(w_0) = f_\text{Hyp}(w_0)$, we obtain the following approximation for σ:

$$\begin{equation} \sigma \approx \frac{f_\text{ac}(w_0)-f_{\text{Hyp},0}(w_0)}{f_{\text{Hyp},1}(w_0)}. \end{equation} \tag{ 19 }$$

By decreasing the step size and increasing the order of the Taylor expansions in the calculation of $f(w_0)$, 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 $\mathrm{e}^{-|A w_0|}$.

Do note that the approximation for $f(w_0)$ from the late-time transseries solution can easily be extended to higher orders in the transseries parameter σ, by computing more non-hydrodynamic sectors $\Phi^{(n)}(w)$ in (9). In figure 2 the results of the early-to-late-time matching $C \leftrightarrow \sigma$ are plotted for C > 0. Our results are consistent with the observation in section 1 that as $(C,\,\sigma) \to (0^+,\sigma_+ = -0.349261+0.00273515 \mathrm{i} )$, the solutions $f_C(w)$ converge pointwise to the solution $f_+(w)$, which is finite at the origin. The function $\sigma(C)$ in figure 2 is roughly linear (left plot) except for a tiny region around the origin C = 0, where the convergence towards $\sigma_+$ is very slow and is best visualised on a log-linear plot.

Another way of approximating $f(w_0)$ is through Borel resummation (see e.g. [8] for a review). The Borel transform of a series $\Phi(w) = \sum_{j\ge0}\,a_j\, w^{-j}$ is given by ¹³

$$\begin{equation} \mathcal{B}\left[ \Phi\right](\xi)= a_0\,\delta(\xi)+\sum_{j=0}^{+\infty} \frac{a_{j+1}}{j!} \xi^j. \end{equation} \tag{ 20 }$$

We truncate the series in (20) after N₀ terms ¹⁴ and calculate its Padé approximant $\text{BP}_{N_0}\left[ \Phi\right]$, i.e. we approximate the resulting truncated sum by a rational function $\text{BP}_{N_0}\left[ \Phi\right]$ with a numerator/denominator of order $\left \lfloor N_0/2 \right \rfloor$.

The Borel-Padé resummation method then consists of taking the inverse Borel transform of $\text{BP}_{N_0}\left[ \Phi\right]$, which is given by the Laplace transform:

$$\begin{equation} \mathcal{S}_{N_0,\theta} \Phi(w) = a_0 + \int_0^{\mathrm{e}^{\mathrm{i} \theta}\infty} \mathrm{d} \xi \, \mathrm{e}^{-w \xi}\text{BP}_{N_0}\left[ \Phi\right](\xi). \end{equation} \tag{ 21 }$$

The resurgence properties of the transseries (9) directly translate to the existence of pole singularities of the Padé approximant $\text{BP}_{N_0}\left[ \Phi\right](\xi)$ in equation (21) along the positive real axis—the Stokes line—and these singularities reflect the branch cuts of the Borel transform (20), starting at all $\xi = nA,\,\,n\in\mathbb{N}$ (one for each exponential in the transseries). Thus we cannot resum $\mathcal{S}_{N_0,\theta} \Phi(w)$ along the positive real axis, and we need to choose integration contours that have θ slightly away from this axis, either above or below. Although this ambiguity in the choice of integration contour gives rise to an imaginary contribution for each summed sector $\mathcal{S}_{N_0,\theta}\Phi^{(n)}(w)$, 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—this summation procedure it the the so-called median summation [34, 35]. To do so we pick a small negative angle $\theta\,=\,-\varepsilon \lt0$ for the integration (21), and require the imaginary value of σ in the following way (see appendix B for some more details):

$$\begin{equation} \mathrm{i} \, \mathbb{I}\text{m}(\sigma) = \frac{S_1}{2}. \end{equation} \tag{ 22 }$$

We then obtain an approximation for $f(w_0)$ to first order in the transseries parameter σ ¹⁵ , given by,

$$\begin{equation} f_\text{B}(w_0) \equiv \mathcal{S}_{N_0,-\varepsilon} \Phi^{(0)}(w_0) +\sigma w^{\beta}\,\mathrm{e}^{- A w_0} \mathcal{S}_{N_0,-\varepsilon} \Phi^{(1)}(w_0). \end{equation} \tag{ 23 }$$

In analogy to (19), using the Borel resummation method we arrive at the following expression for σ:

$$\begin{equation} \sigma \approx \frac{f_\text{ac}(w_0) - \mathcal{S}_{N_0,-\varepsilon} \Phi^{(0)}(w_0)} {w^{\beta}\,\mathrm{e}^{-A w_0}\, \mathcal{S}_{N_0,-\varepsilon} \Phi^{(1)}(w_0)} . \end{equation} \tag{ 24 }$$

Notice that for both equations (19) and (24) we only went up to linear order in σ in the approximation of $f(w_0)$. To obtain more accurate results, we could have included higher powers of σ, which amounts to including extra exponential orders ¹⁶ . For the Borel summation method we would only need to numerically compute the integrals (21) for the higher-order hydrodynamic sectors $\Phi^{(n)}(w)$ in (9), while the generalisation of the hyperasymptotic summation is a bit less straightforward. It can nonetheless be done, and we refer the reader to the literature [11, 30, 31, 36]. 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 w₀.

Once the parameter σ has been matched to a given initial condition ¹⁷ , the transseries (9) can be used to find an approximation of f(w) everywhere, via a 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 figure 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 w₀.

Zoom In Zoom Out Reset image size

In figure 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 $\mathrm{e}^{-2|A w|}$ (the order of the first exponential we have neglected). We can also clearly see that the approximations given by each summation method are still 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 figure 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 ¹⁸

$$\begin{equation} f_\text{opt}(w) = \sum_{n=0}^{N_\text{opt}(w)-1}\,a^{(0)}_n w^{-n},\quad \text{where}\qquad N_\text{opt}(w) = \left \lfloor | A w | \right \rfloor . \end{equation} \tag{ 25 }$$

The accuracy of this optimal truncation is approximately $\mathrm{e}^{-|A w|}$, in agreement with the plots in figure 3.

Now that we have discussed the interpolation between late and early times using Borel resummation and hyperasymptotics, the next section will be devoted to the transasymptotic summation method.

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—using numerical analytical continuation from the origin to the matching point $w\,=\,w_0$ to obtain the numerical approximation $f_\text{ac}(w_0)$ (see (18)). We next compute an approximation for $\tau(w_0)$, from which the transseries parameter σ can directly be calculated using the definition of $\tau(w_0)$, equation (12). The idea is the following: we want to solve for the function $\gamma(w)$ obeying:

$$\begin{equation} \mathcal{F}(\gamma(w), w) = \sum_{n \geqslant 0}\, F_r\left(\gamma(w)\right) w^{-r} =f_\text{ac}(w_0)= \text{constant},\,\textrm{for }\,w\gg1, \end{equation} \tag{ 31 }$$

which will be equal to $\tau(w_0)$ when evaluated at the point w₀, i.e. $\gamma(w_0) = \tau(w_0)$. The function $\gamma(w)$ satisfying (31) admits a perturbative, divergent asymptotic expansion in w⁻¹:

$$\begin{equation} \gamma(w) = \sum_{k=0}^{+\infty}\, \gamma_k\,w^{-k}\,, \end{equation} \tag{ 32 }$$

and determining $\gamma(w)$ will correspond to finding the coefficients γ_k . Truncating the above expansion at its first term $\gamma(w) = \gamma_0+\mathcal{O}\left(w^{-1}\right)$, we find from (31) that up to leading order:

$$\begin{equation} 1+W\left(\frac{3}{2}\gamma_0\right)=f_{\text{an}}(w_0). \end{equation} \tag{ 33 }$$

Then $\gamma_0 = \tau(w_0)$ (up to leading order in $w_0^{-1}$), and using the definition of $\tau(w)$ (12) we finally find:

$$\begin{equation} \sigma(w_0) = \big(f_\text{ac}(w_0)-1\big) \,w_0^{-\beta}\, \mathrm{e{\,}}^{f_\text{ac}(w_0)-1 +A w_0} \bigg(1+ \mathcal{O}\left(w_0^{-1}\right) \bigg). \end{equation} \tag{ 34 }$$

This result can be easily extended to higher orders in $w_0^{-1}$ by including higher orders in the ansatz (32) and matching powers of w⁻¹ in (31). The first four coefficients of the perturbative expansion of $\gamma(w)$ 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 $F_r(\tau(w))$ vary with w. The result is displayed in figure 3, where we can see that this approach slightly outperforms optimal truncation. Note that we only calculated the coefficient functions $F_r(w)$ up to $r\,=\,15$, and so the calculation is no longer optimal after the kink in the logarithmic error plot of figure 3. Furthermore, the kink happens at a higher value of w than we would expect from the resummation point w₀ corresponding to 15 terms with classical optimal truncation given by (25).

Transasymptotics can also be used to describe global properties of the function f(w) from (4), such as zeros, poles and 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 = +\infty$. Let us start by sketching out how the locations of the branch points may be obtained. Notice that in the solutions plot found in figure 1, the locations w_s of the square root branch points depend on the initial value problem for f(w) (given by f(0) for convergent solutions at the origin, or by the value of the continuous continuation of $w^4 f(w)$ at w = 0 for the divergent solutions). From the perspective of late-time asymptotics, this means that the locations w_s are functions of the transseries parameter σ. As already mentioned, all the coefficient functions $F_r(\tau)$ in the transasymptotic summation (26) can be expressed as rational functions of the Lambert-W function $W(\frac{3}{2}\tau)$, which has a square-root branch point at $\tau = -\frac{2}{3}\mathrm{e}^{-1}$. This branch point in the τ-plane translates to an infinite number of branch points in the w-plane, once we substitute $\tau = \tau(w)$ as in (12). Since the Lambert-W function appears in all the coefficient functions in the transasymptotic summation (26), we expect the function f(w) 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 $F_r(\tau)$ can be used to provide an approximation for the locations w_s.

All zeros of f(w) are square root branch point singularities (see the expansions in equation (6)) with the exception of a potential regular zero at $w = (4-\beta)/A \approx 0.502$. Hence we can solve for the branch points w_s by solving the equation,

$$\begin{equation} \mathcal{F}(w=w_\text{s},\sigma) = 0, \end{equation} \tag{ 35 }$$

for $w_\text{s}(\sigma)$, where $\mathcal{F}$ is the transasymptotic summation in equation (26). We find,

$$\begin{equation} w_\text{s}(t) \simeq w^\text{(approx)}_\text{s}(t)= \frac{t}{A}+\frac{\beta}{A}\, \log t+\frac{1}{A\, t}\left(\beta^2\,\log t+\beta^2-5\beta-\frac{3}{A}\right)\, , \quad \text{ as } t \to \infty, \end{equation} \tag{ 36 }$$

where we have introduced the variable:

$$\begin{equation} t\equiv t(n,\sigma) = \log\left(\frac{3\sigma \mathrm{e}}{2A^\beta}\right) + \pi \mathrm{i} (1+2n ), \quad \quad n \in \mathbb{Z}. \end{equation} \tag{ 37 }$$

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 equation (36) is only a good approximation for the branch points/zeros of f(w) when $|t|$ is large enough. In particular, $w^\text{(approx)}(t)$ becomes more accurate for large values of the discrete parameter n, because the auxiliary variable t grows as an affine function of n. Since the leading order approximation $w^\text{(approx)}(t)$ from equation (36) grows linearly in t, the branch points which lie far from the origin are the ones best approximated by $w^\text{(approx)}(t)$.

Numerically, we can compute the zeros of f(w) by initially guessing the position of the branch point using (36) ²⁰ and then using a contour integral to find an approximation for the exact location. We start by choosing a value for the transseries parameter σ and use the hyperasymptotic approximation equation (15) to find $f(w_0)$ (at e.g. $w_0 = 10$). We then analytically continue f(w) from w₀ to a point in the vicinity of our prediction (36) using the Taylor series method, $w_1 = w^\text{(approx)}_\text{s}(t) + \varepsilon$ (e.g. ε = 0.3). Next we analytically continue again to compute the data on the circle $|w- w^\text{(approx)}_\text{s}(t)| = \varepsilon$. Using the trapezoidal rule [40] we evaluate the contour integral of $\frac{w f^{\prime}(w)}{f(w)}$ to obtain the zeros of f(w) ²¹ . The approximate locations obtained with equation (36) as well as the numerical results are listed in table 1, and plotted in figure 4.

Zoom In Zoom Out Reset image size

Table 1. Approximations for the locations of the square-root branch points of (36) versus their numerically computed values for $\sigma=\frac23$.

	Approx. (36)	Numerical
$n\,=\,0$	$0.0975\,+\,0.6040\mathrm{i}$	$0.1147\,+\,0.5076\mathrm{i}$
$n\,=\,1$	$0.1555\,+\,1.416\mathrm{i}$	$0.1580\,+\,1.384\mathrm{i}$
$n\,=\,2$	$0.1817\,+\,2.276\mathrm{i}$	$0.1827\,+\,2.257\mathrm{i}$
$n\,=\,3$	$0.1991\,+\,3.143\mathrm{i}$	$0.1997\,+\,3.129\mathrm{i}$
$n\,=\,4$	$0.2122\,+\,4.012\mathrm{i}$	$0.2125\,+\,4.001\mathrm{i}$

There is yet another powerful application of transasymptotics, which is that it can be used to correctly predict the different asymptotic behaviour of our solutions in separate regions of the w-plane. Consider the attractor $f_+$ in the solutions plot of figure 1 (the black curve). At large, positive w, $f_+(w)$ 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 $\mathrm{e}^{-A w}$, see [41]). 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 $\tau\sim \mathrm{e}^{-A w}$, the large, negative w limit is associated with exponentially large values of τ. Thus when one flips the sign $w \to -w$, the powers of w⁻¹ in the transasymptotic summation (26) do not change size, while the exponential variable $\tau\sim \mathrm{e}^{-A w}$ becomes exponentially large instead of exponentially small. Since the coefficient functions in the transasymptotic summation are just rational functions of $W(\tau)$, and the large τ expansion of $W(\tau)$ is known (see [42] and appendix E), we were able to use transasymptotics to correctly derive the first four terms of the linearly growing expansion (11). The details of the calculations in this section are beyond the scope of this publication and will be explored in an upcoming paper [43].