Relativistic Hydrodynamics: A Singulant Perspective

There is growing evidence that the hydrodynamic gradient expansion is factorially divergent. We advocate for using Dingle's singulants as a way to gain analytic control over its large-order behaviour for nonlinear flows. Within our approach, singulants can be viewed as new emergent degrees of freedom which reorganise the large-order gradient expansion. We work out the physics of singulants for longitudinal flows, where they obey simple evolution equations which we compute in M\"uller-Israel-Stewart-like models, holography and kinetic theory. These equations determine the dynamics of the large-order behaviour of the hydrodynamic expansion, which we confirm with explicit numerical calculations. One of our key findings is a duality between singulant dynamics and a certain linear response theory problem. Finally, we discuss the role of singulants in optimal truncation of the hydrodynamic gradient expansion. A by-product of our analysis is a new M\"uller-Israel-Stewart-like model, where the qualitative behaviour of singulants shares more similarities with holography than models considered hitherto.

There is growing evidence that the hydrodynamic gradient expansion is factorially divergent.We advocate for using Dingle's singulants as a way to gain analytic control over its large-order behaviour for nonlinear flows.Within our approach, singulants can be viewed as new emergent degrees of freedom which reorganise the large-order gradient expansion.We work out the physics of singulants for longitudinal flows, where they obey simple evolution equations which we compute in Müller-Israel-Stewart-like models, holography and kinetic theory.These equations determine the dynamics of the large-order behaviour of the hydrodynamic expansion, which we confirm with explicit numerical calculations.One of our key findings is a duality between singulant dynamics and a certain linear response theory problem.Finally, we discuss the role of singulants in optimal truncation of the hydrodynamic gradient expansion.A by-product of our analysis is a new Müller-Israel-Stewart-like model, where the qualitative behaviour of singulants shares more similarities with holography than models considered hitherto.

I. INTRODUCTION
The past two decades have been a true golden age for our understanding of dissipative relativistic fluids.This surge of interest has been primarily driven by the interplay between growing sophistication in hydrodynamic modelling of the newly discovered quark-gluon plasma at RHIC (and later also at LHC) [1,2] and unprecedented progress in studying nonequilibrium phenomena with hydrodynamic tails at strong and weak interaction strength using, respectively, holography and relativistic kinetic theory [3,4].Today the domain of relativistic hydrodynamics encompasses also time-dependent black hole phenomena via the holographic fluid-gravity duality [5], neutron star modelling in the context of gravitational wave physics [6][7][8], as well as condensed matter phenomena [9], including electron flow in graphene [10].Furthermore, recent developments indicate that relativistic hydrodynamics can apply also to far-from-equilibrium situations, which is a subject of the hydrodynamics attractors research program [11][12][13].
The key notion underlying relativistic hydrodynamics is that of constitutive relations.As an effective field theory of transport of conserved currents, the fundamental object of interest is the expectation value of the energy-momentum tensor T µν .Hydrodynamic constitutive relations provide an ansatz for this quantity as an infinite series in gradients where E and U µ (U µ U µ = −1) are the local energy density and the local velocity encapsulating slow degrees of freedom and Π µν is a gradient-expanded dissipative part of the energy-momentum tensor.Here and in the following, is a formal bookkeeping parameter counting the number of gradients, and equalities involving infinite power series in are to be understood in the sense of a formal ansatz.To remove ambiguities in the gradient expansion (1) associated to field redefinitions, we work in the Landau frame and hence demand that Π µν is transverse to the fluid velocity, Π µν U ν = 0. We refer the reader to Refs.[3,4] for contemporary reviews on relativistic hydrodynamics.
The gradient expansion (1) can be constructed as long as the energy-momentum tensor at the spacetime point being considered has a single real timelike eigenvector [14].Provided that this condition is met the gradient expansion (1) can be written down, and thus it is natural to ask whether it can capture processes arbitrarily far away from global thermal equilibrium.The answer to this question depends crucially on whether the gradient expansion is a convergent series.If convergent, the gradient expansion (1) can in principle provide a full description of the underlying microscopic theory as it pertains to the expectation value of the energy-momentum tensor.If this were the case, there would be no obstruction in obtaining an approximation to T µν with an error as small as desired: one just needs to truncate the gradient expansion at a sufficiently high order.On the other hand, if the gradient expansion is a divergent series, this is not possible, even in principle.In this situation, the first step becomes that of finding the optimal truncation order.
In this paper we propose to systematically apply the idea of the singulants to the gradient expansion (1b) regardless of the presence of any flow symmetries.This leads to the large-order ansatz (2) Π (n)  µν (t, x) ∼ A µν (t, x) Γ(n + α(t, x)) χ(t, x) n+α(t, x) , with the singulant itself becoming a scalar field in spacetime (see also Refs.[40,41]).Again, A µν (t, x) and α(t, x) play no role at leading order at large n.Our analysis will focus on longitudinal flows [34,42], where we explicitly confirm the validity of the ansatz (3) in a number of models.We will see that, in general, there will be multiple singulants which contribute additively to Eq. (3).From the point of view of the hydrodynamic gradient expansion alone, singulants can be thought of as a novel emergent phenomenon in relativistic hydrodynamics, not necessarily related to a single contribution to the expansion at a given order, but rather to a reorganisation of the whole series.
For systems near equilibrium (or other special solutions such as attractors) quasinormal modes describe the decay of non-hydrodynamic degrees of freedom and the approach to equilibrium.Far from equilibrium however, quasinormal modes can no longer be used and it is not obvious what concept or object replaces them, if anything.We advocate that singulants can fill this role.The reason for this is that singulants are fundamentally nonhydrodynamic collective fields whose contribution to the energy-momentum tensor decays over time.To see that H < l a t e x i t s h a 1 _ b a s e 6 4 = " L w S X b 7 s L L T K i p 8 I D E V K K Z 2 o f H i 8 = " > A A A B 8 X i c b V B N S w M x F H y p X 7 V + V T 1 6 C R b B U 9 m V g n o r e O m x g m 3 F d q n Z N N u G J t k l y Q p l 6 b / w 4 k E R r / 4 b b / 4 b s 2 0 P 2 j o Q G G Z e e P M m T A Q 3 1 v O + U W F t f W N z q 7 h d 2 t n d 2 z 8 o H x 6 1 T Z x q y l o 0 F r G + D 4 l h g i v W s t w K d p 9 o R m Q o W C c c 3 + R + 5 4 l p w 2 N 1 Z y c J C y Q Z K h 5 x S q y T H n q S 2 J G W W W P a L 1 e 8 q j c D X i X + g l R g g W a / / N U b x D S V T F k q i D F d 3 0 t s k B F t O R V s W u q l h i W E j s m Q d R 1 V R D I T Z L P E U 3 z m l A G O Y u 2 e s n i m / v 6 R E W n M R I Z u M k 9 o l r 1 c / M / r p j a 6 C j K u k t Q y R e e L o l R g G + P 8 f D z g m l E r J o 4 Q q r n L i u m I a E K t K 6 n k S v C X T 1 4 l 7 Y u q X 6 t e 3 9 Y q 9 c f 6 v I 4 i n M A p n I M P l 1 C H B j S h B R Q U P M M r v C G D X t A 7 + p i P F t C i w m P 4 A / T 5 A / d g k W 0 = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " Z j t 5 g p 8 T Y L Y D 1 + h 8 A Y g r 5 F N 7 h P g = " > A A A B 6 3 i c b V B N S w M x E J 3 U r 1 q / q h 6 9 B I v g q e y K o N 4 K X j x W s B / Q L j W b Z r u h S X Z J s k J Z + h e 8 e F D E q 3 / I m / / G b N u D t j 4 Y e L w 3 w 8 y 8 M B X c W M / 7 R q W 1 9 Y 3 N r f x n B 1 k k P X Y H t U I v s e t I r V / y q P w V a J s G c V G C O e q / 8 1 e 0 r k g o q L e H Y m E 7 g J z b M s L a M c D o p d V N D E 0 x G e E A 7 j k o s q A m z a e I J O n F K H 8 V K u y c t m q q / f 2 R Y G D M W k Z v M E 5 p F L x f / 8 z q p j S / D j M k k t V S S 2 a I 4 5 c g q l J + P + k x T Y v n Y E U w 0 c 1 k R G W K N i X U l l V w J w e L J y 6 R 5    they are non-hydrodynamic we note that our large-order ansatz is intimately related to the appearance of contributions to the transseries generalisation of the hydrodynamic gradient expansion which are non-perturbative in .In all examples considered in this paper, the value of |χ| grows after enough time has passed.This leads to a growth in the optimal truncation order and a decrease in the truncation error, thus gradually extending the applicability of the hydrodynamic expansion over time.
Our results indicate that singulants (χ) provide a new perspective on collective states of matter that complements and connects in a novel way the existing paradigms of the hydrodynamic gradient expansion (H) and the amplitude expansion exemplified by linear response theory (A).This is represented by Fig. I.In this work, we primarily consider the region i, where the large-order behaviour of the gradient expansion (1) for nonlinear fluids is governed by singulants.One of our most relevant findings is a duality between the singulant dynamics, as computed in region (i), and a particular linear response theory problem defined in region (ii).Region (iii) corresponds to the realm of linearised hydrodynamics, where one considers the gradient expansion (1) for infinitesimal perturbations away from global thermal equilibrium, but does not focus on its large-order behaviour.Finally, in region (iv), when utilising the three approaches-singulants, the gradient expansion and linear response theory-we gain complete analytic control of the singulants.Appendix D in the Supplemental Material provides an example of such a study.
This paper is organised as follows.In Sec.II, we introduce longitudinal flows and discuss several results pertaining to the large-order behaviour of the gradient expansion (1) and singulants in this context.Then, in Sec.III, we test the general results for longitudinal flows put forward in Sec.II in a series of phenomenological models of the Müller-Israel-Stewart (MIS) class.These include the Baier-Romatschke-Son-Starinets-Stephanov (BRSSS) model [43], the Heller-Janik-Spaliński-Witaszczyk (HJSW) model [44], and a new model we introduce for the first time in this work.These and related models originate from the need to have equations of motion for relativistic fluids that give rise to causal and stable evolutions [45,46].In particular, the original MIS model is the workhorse for the vast majority of hydrodynamic simulations of quark-gluon plasma in nuclear collisions.In the three MIS-like models we consider, we demonstrate that the singulant dynamics can be mapped to a linear response theory problem that consists in computing the poles of a momentum-dependent sound attenuation length γ s .In Sec.IV, we explore the largeorder behaviour of the gradient expansion for longitudinal flows in the context of holography and show that, if a factorial divergence is present, then the singulant equation of motion is also determined by the poles of γ s .Sec. V discusses singulants in kinetic theory, focusing on the gradient expansion of the distribution function.We test our analytic predictions by explicit numerical computations for Bjorken flow in the relaxation time approximation (RTA).Sec.VI explores the interplay between singulants and optimal truncation in the context of BRSSS theory.The paper closes with a discussion of the physical interpretation of singulants and key open problems in Sec.VII.
Several computations supporting the results presented in the main body of the paper have been relegated to Supplemental Material as appendices.Appendix A provides the demonstration of the general results introduced in Sec.II, while Appendix B discusses the large-order behaviour of the gradient expansion (1) beyond longitudinal flows in the context of MIS-like models.Appendix C contains the detailed computation of γ s in holography.Appendix D-which has been mentioned before-explores the large-order behaviour of the gradient expansion (1) in linearised hydrodynamics.Finally, in Appendix E we discuss the causality and stability of the new MIS-like model we introduced in this work, at the level of linearised perturbations of global thermal equilibrium.

II. SINGULANTS IN LONGITUDINAL FLOWS
In this paper we develop the paradigm put forward in the Introduction for a particular class of fluid flows which we refer to as longitudinal.Longitudinal flows simplify drastically the technical aspects of our analysis and, at the same time, keep intact the essential features of far-from-equilibrium fluid dynamics.In particular, as originally shown in Ref. [34], for longitudinal flows in MIS-like theories it is feasible to compute the gradient expansion (1b) numerically up to an order sufficiently large to assess its asymptotic behaviour.
A longitudinal flow in d-dimensional Minkowski spacetime is defined by singling out one spatial direction, x, and demanding translational invariance and isotropy in the hyperplane spanned by the remaining spatial coordinates .This requirement implies that that the nontrivial dynamics is confined to the plane spanned by t and x, which we refer to as the longitudinal plane.With these symmetry restrictions, both the fluid velocity U µ and any two-tensor A µν which is symmetric, transverse to U µ , and traceless can be parameterised in terms of a single degree of freedom.Specifically, where P µν T is the projector into the transverse hyperplane and u, A depend solely on t and x.In this work, we consider only longitudinal flows in conformal theories, in such a way that T µν is traceless.This entails that the equation of state of the fluid is given by P(E) = E/(d − 1) and that Π µν is also traceless.Using the energy density we define T (t, x) = (E(t, x)/E 0 ) 1 d and refer to this field as the effective temperature.In equilibrium T (t, x) becomes the temperature of the system, as determined by the equation of state.
In a longitudinal flow, the gradient expansion (1b) becomes a gradient expansion for Π , and the asymptotic ansatz (3) reads The singulant field, χ(t, x), controls the subleading geometric correction to the leading-order factorial growth of the gradient expansion.
In this work, we explore the gradient expansion (5) in the context of MIS-like models (Sec.III), holography (Sec.IV) and RTA kinetic theory (Sec.V).We note the following.
• In MIS-like theories, we work directly at the level of Π and the gradient expansion (5).
• In holography, we work with the gradient expansion of the bulk metric, from which the gradient expansion (5) descends by holographic renormalisation.In Eq. ( 7), X = (r, x) and r parameterises the radial direction of the higher-dimensional geometry.The singulant field governing the large-order behaviour of the gradient expansion (5) follows from the singulant field governing the large-order behaviour of the gradient expansion (7): Note that we have assumed that the singulant field is r independent.We will show that this assumption is self-consistent in Sec.IV.
• In kinetic theory, we work with the gradient expansion of the distribution function, from which the gradient expansion (5) descends by computing the second-order moments.In Eq. ( 9), p is the momentum.The singulant field governing the large-order behaviour of the gradient expansion (5) follows from the singulant field governing the large-order behavior of the gradient expansion (9): In every theory being considered in this work, the exact Π can be computed in terms of the energy density E and fluid velocity U ≡ U µ ∂ µ associated to a given out-of-equilibrium state by solving a nonlinear system of integro-PDEs (partial differential equations).This nonlinear system of integro-PDEs is specific to the theory in question.One has the following.
• In holography, it corresponds to a subset of the Einstein equations governing the higher-dimensional bulk geometry [cf.Eq. (62) in Sec.IV].
Given the hydrodynamic fields E and U , we determine the coefficients of the gradient expansions ( 5), (7), and (9) through the following systematic procedure.
(i) We introduce the bookkeeping parameter into the nonlinear system of integro-PDEs determining Π mentioned above by means of a homogeneous rescaling of the longitudinal plane coordinates: (ii) We plug in the ansatz between Eqs. ( 5), ( 7), (9) appropriate for the theory in question.
(iii) We solve in a small-expansion.
The end result of this procedure is a set of recursion relations that fix the nth order coefficient of the relevant gradient expansion in terms of the previous lower-order ones.Explicit examples are provided in the rest of the paper whenever necessary.We emphasise that, once the recursion relations have been written down, any ambiguity in the gradient expansion associated with the usage of the conservation equations is removed.The singulant equation of motion is determined by the recursion relations, and it follows straightforwardly from the requirement that the asymptotic ansatz (6) [( 8) and (10) resp.] for the coefficients of the gradient expansion (5) [ ( 7) and ( 9) resp.]solves the associated recursion relations at leading order at large n.A fact that will be crucial for our analysis is that, for a factorially divergent gradient expansion, the recursion relations simplify drastically in this large n regime.The two major simplifications taking place at large n are the following.
(a) The recursion relations become linear and, (b) terms in the recursion relations associated with gradients of the hydrodynamic fields, E and U , drop out.
Hence, rather than being sensitive to every possible term in the recursion relations, the singulant dynamics only depends on a subset of dominant terms, which are characterised by points (a) and (b) above.
To illustrate the general discussion of the previous paragraph, we consider the case of MIS-like models.The cases of holography and kinetic theory are similar and are discussed in Secs.IV and V, respectively.In MIS-like models, when one introduces the large-order ansatz (6) into the recursion relations and takes the n → ∞ limit, one finds that the recursion relations simplify to In accordance with property (a), Eq. ( 12) is linear in Π ; in accordance with property (b), the tensors f only depend on the local values of the hydrodynamic fields at the spacetime point being considered.Terms involving gradients of these fields and/or nonlinear in Π result in subleading contributions at large n and are therefore not included in Eq. (12).The singulant equation of motion follows immediately from the fact that the factorial-over-power ansatz (6) has to solve the simplified form of the recursion relations (12) and the observation that, at leading order at large n, the following identity holds: One important consequence of the linearisation of the recursion relations is that the most general large-order behaviour of a factorially divergent gradient expansion of the form (5) (resp.( 7) and ( 9)) is described by a linear combination of contributions of the form (6) (resp.( 8) and ( 10)), where each χ q satisfies the same equation of motion but with different initial conditions.Note that, in every theory being considered, the reality of the gradient expansion coefficients appearing in the left hand side of Eq. ( 14) implies that every singulant contribution appearing in the right hand side is either real or accompanied by a complex-conjugated partner.In the specific models where we computed the gradient expansion numerically, we always find that these complex-conjugated singulant pairs are present.When multiple singulant contributions are present, we refer to the one with the smallest norm as the dominant singulant.
A complementary viewpoint on the singulant dynamics is provided by the observation that, when upgrading the perturbative series (5) to a transseries, the singulants weight the nonperturbative transseries sectors, where the ellipsis represents possible nonperturbative transseries sectors associated to nonlinear interactions between different singulants χ q .In the holography and kinetic theory cases, there exist transseries analogous to Eq. ( 15) for g AB and f , from which Eq. ( 15) descends naturally.
The transseries approach allows for an alternative derivation of the singulant equation of motion.In this alternative derivation, one starts with the nonlinear system of integro-PDEs appropriate for the theory in question after the bookkeeping parameter has been introduced but, rather than introducing the relevant gradient expansion among (5), (7), and (9), one plugs in its associated transseries ansatz and expands around = 0. Demanding that the coefficient multiplying e − χq vanishes at leading order in the → 0 limit gives directly the singulant equation of motion.From this perspective, the singulant equation of motion can be interpreted as the eikonal equation coming from a WKB analysis of the nonlinear integro-PDE system determining Π .
As we discuss in detail in Appendix A in the Supplemental Material, the procedure described in the previous paragraph is equivalent to (i) take the original nonlinear system of integro-PDEs and linearise it around the zeroth-order term of the corresponding gradient expansion, (ii) neglect terms associated to gradients of the hydrodynamic fields, and (iii) perform the replacement, The end result of (i)-(iii) is again the singulant equation of motion.One of the main benefits of the transseries approach, as embodied in these three steps, is that, in the cases of MIS-like models and holography, it allows one to identify in a straightforward way a linear response theory computation dual to the singulant equation of motion.Indeed, steps (i)-(iii) are respectively equivalent to first taking the nonlinear system of integro-PDEs that determines Π , then setting the hydrodynamic fields, T and U , to spacetime-independent constants, T 0 and U 0 , and finally finding the dispersion relation for linearised perturbations of Π : There exists a map involving the identifications, that transforms the equations determining the dispersion relation of the linearised perturbation ( 17) into the equations of motion determining the singulant dynamics.This map is discussed in further detail in Appendix A in the Supplemental Material and we will see explicit incarnations of it in Secs.III and IV.
A linear response theory problem that features prominently in studies of nonequilibrium physics in MIS-like theories, holography and kinetic theory is the computation of the modes of the system [43,[47][48][49].A mode is a singularity of a retarded thermal two-point correlator in Fourier space, and can be thought of as being located at a frequency ω ∈ C that depends nontrivially on the the spatial momentum k ∈ C d−1 as determined by the dispersion relation ω = ω(k).In linear response theory, hydrodynamic modes are associated with long-lived and slowly varying excitations such that ω(k) → 0 as |k| → 0, while nonhydrodynamic modes are associated with excitations that do not obey the latter property.
A fact that the reader should keep in mind is that the linear response theory problem dual to the singulant equation of motion may not correspond to a mode computation.The reason is that, in a mode computation, the hydrodynamic fields T and U are treated dynamically.These two problems are therefore equivalent if and only if the dynamics of the hydrodynamic fields decouple from the dynamics of δΠ .For general longitudinal flows, this will be the case in the MIS-like models discussed in Secs.III B and III C, but not in the in the new model introduced in Sec.III D, nor in holography.
To conclude our overview of singulants in longitudinal flows, we comment briefly on how singulants can be extracted from the numerical values of the gradient expansion coefficients.For illustrative purposes, we focus on a gradient expansion of the form (5) with the largeorder behaviour (14), and note that identical strategies exist for the gradient expansions of the metric (7) and the distribution function (9).
First, in the restricted case in which one is solely interested in the absolute value of the dominant singulant, that, asymptotically, a root-test plot of gradient expansion coefficients results in a straight line of slope 1/(e|χ d |), in such a way that |χ d | is fixed once this slope is known. 1n the other hand, if one wishes to determine the values of the subdominant singulants as well, the best strategy is to consider the analytical continuation of the Borel transform of the gradient expansion.In this case, each singulant contribution χ q appears as a singularity in the Borel plane located at ζ = χ q , where ζ is a complex-valued variable used to parameterise the Borel plane throughout this paper.

III. MIS-LIKE MODELS A. Introduction
In this section, we compute the singulants contributing to the large-order behaviour of the gradient expansion for longitudinal flows in a class of phenomenological MIS-like models.As mentioned in the Introduction, these are the BRSSS model [43], the HJSW model [44], and the new model introduced for the first time in this work.The rationale behind these models is embedding hydrodynamics in a framework compatible with relativistic causality.To achieve this, the dissipative tensor Π µν is promoted to a set of independent dynamical degrees of freedom obeying their own equation of motion.While such models are inspired by hydrodynamics, it is crucial to note that unlike hydrodynamics they are exact, in the sense that they are not defined using a perturbative expansion in the number of gradients.
Different models in this class are distinguished by the different equation of motion obeyed by Π µν which we detail case by case in the sections that follow; however, each case obeys the same current conservation equations, ∇ µ T µν = 0, which for longitudinal flows are where As a final comment, we note that the MIS-like models we consider in this work are causal and stable at the level of linear response around global thermal equilibrium, however the reader should keep in mind that these are not sufficient conditions for their causality and well-posedness at the fully nonlinear level [46,50].

B. BRSSS model
In the BRSSS model the exact dynamical equation governing Π µν , when specialised to longitudinal flows, is given by from which the series ( 5) can be obtained by solving the following recursion relation: This recursion relation is solved at large orders by the factorial-over-power ansatz ( 14), provided the singulant equation is obeyed, In earlier work [34] we established the factorial growth of Π (n) for general nonlinear longitudinal flows in BRSSS, by numerically solving Eqs. ( 20)-( 22) and evaluating Eq. ( 23) on them.An example of this is given in the top panel of Fig. 2 where the n! behaviour is illustrated. 2The presence [Eq.( 14)] and motion [Eq.( 24)] of the singulants governing the large-order behaviour of such flows can be readily confirmed for such solutions, and we now present two such ways of doing so.
First, according to Eq. ( 14) the singulant field χ q with the smallest |χ q | at any given spacetime point dominates the large-order behaviour yielding the prediction ∼ n e|χq| where χ q obeys Eq. ( 24).Solving Eq. ( 24) on a given background is unique up to a choice of complex integration constant per flow line.Two such solutions are presented as the blue and red curves in the middle panel of Fig. 2 alongside the black dots which correspond to fitting a straight line to the numerical series data, |Π , showing excellent agreement.Here we have chosen to determine the singulant integration constants by adjusting for the best fit on the entire flow line. 2 The initial data corresponds to a periodic overdensity in T that locally resembles a Gaussian, see Ref. [34] for more details on this function and the numerical method.We have set E/T 4 = 1, η/s = 1/(4π) and τ Π T = 1/4.Second, the presence of χ q can be seen as singularities appearing in the Borel transform of Eq. ( 5).This is demonstrated in the bottom panel of Fig. 2.This is a snapshot at a time labelled by the dashed vertical line in the middle panel, where the blue singulant dominates.Given the integration constants as determined in the previous paragraph, there is a precise match between χ q and the location of a singularity in the Borel plane inferred by a Padé approximant of the Π (n) series.Along a flow line these singulants move from left to right in the Borel plane according to Eq. ( 24) as indicated by the arrows in the lower panel.There is a time for which they are at their point of closest approach to the origin, corresponding to the maxima in the middle panel of Fig. 2. Similarly there is a time for which |χ red | = |χ blue | corresponding to an exchange of dominance as seen by the crossing of the red and blue singulant trajectories in the middle panel.

C. HJSW model
The HJSW model [44] is a generalisation of conformal BRSSS theory such that the equation of motion for Π µν includes second-order derivatives along the fluid velocity U , where Ω = Ω R + iΩ I ∈ C, c σ ∈ R, and we are neglecting possible nonlinear terms in Π µν .As the temperature T is the only dimensionful scale of the theory, the shear viscosity takes the same functional form as it did in BRSSS theory.
The physical motivation behind the construction of the HJSW theory was upgrading the original conformal BRSSS theory to a model with a nonhydrodynamic sector closer to the AdS/CFT one, in the sense of having two nonhydrodynamic sound modes, ω (±) NH (k), with opposite real parts at zero momentum.Indeed, in the HJSW model, ω For linearised perturbations around thermal equilibrium (and in the local rest frame of the fluid) the HJSW model is known to be causal and stable provided that the parameters Ω, η/s and c σ are chosen appropriately.Linear stability, in particular, always requires a finite c σ [44].
Demonstrating the causality of this model at the fully nonlinear level is a challenging open problem.It should be kept in mind that even in the case of MIS-type theories such a proof has only become available recently [46].That being said, exploratory studies of numerical evolution in the HJSW model, with appropriate parameters, have not exposed any pathologies.
It is immediate to compute the recursion relation obeyed by Π (n) once Eq. ( 25) is known.The only structural difference with respect to the BRSSS case is that, since Eq. ( 25) is second order, Π (n) now depends both on Π (n−1) and Π (n−2) for n > 2. In the following, we work in d = 4.We have that where θ is the expansion of the flow, θ = ∇ α U α .Empirically, we always find that the gradient expansion defined by the recursion relations (27a)-( 27c) is factorially divergent when evaluated on any longitudinal flow.In Fig. 3, we provide an example for initial data of the form with σ = 1.We have chosen E/T 4 = 1, η/s = 1/(4π), Ω R = 2, Ω I = 4, and c σ = π.The root test applied to the gradient expansion associated to these initial data at x = 0 and t = 0.5 is represented in the upper panel.The asymptotic behaviour at large n is clearly linear, implying that the gradient expansion is factorially divergent.In the lower panel, we represent the singularities of the Padé approximant to the Borel transform of the gradient expansion as black dots.There are several well-defined lines of pole condensation, grouped in complex-conjugated pairs.The poles from which these lines emanate are candidate singulants, and have been highlighted as stars.
The singulants with the smallest norm correspond to the cyan stars.
According to the general results presented in Sec.II, the singulant equation of motion is determined by the first three terms in Eq. (27c).One finds that Hence, just as it happened in the conformal BRSSS model, the singulant trajectory is determined by ω NH (k = 0) evaluated at the effective temperature.Equation (29) implies that, when one moves along a particular flow line, the trajectories described by the singulants are inclined lines in the Borel plane, with slopes of magnitude |Ω R /Ω I |.This is in stark contrast with the BRSSS case, where the singulant trajectories were given by horizontal lines.Let us select the flow line passing through x = 0 at t = 0, which corresponds to the t axis in the longitudinal plane.In Fig. 4, we plot the real and imaginary parts of the three singulants of smallest norm as a function of time along this particular flow line.We have selected the singulants with positive imaginary part.These singulants, which have been computed by means of a Padé approximant to the Borel transform of the gradient expansion, have to be compared with the prediction of Eq. ( 29): Since we cannot determine independently the initial values of the singulants, we have chosen to match the prediction (30) to the numerical results at t = 0.5.The final outcome of our analysis is that, for every singulant under consideration, the prediction of Eq. ( 30), which is represented as a solid line, does an excellent job in describing the time evolution of the singulants we have computed numerically.

Physical motivation
The analysis performed so far in the BRSSS and HJSW models revealed that the large-order behaviour of the gradient expansion is related to the sound channel nonhydrodynamic mode frequencies evaluated at k = 0 through the singulant equation of motion.Schematically, This is not a coincidence.As we mentioned in Sec.II, both the singulant equation of motion (31) as well as the equation that determines the mode frequencies are obtained through a linearisation procedure and, in the BRSSS and the HJSW models, the former equation can always be mapped onto the latter one evaluated at zero momentum.This rests crucially on the fact that the equation of motion for Π µν only involves comoving derivatives along a flow line.In light of this, one should not expect the relation between the singulant equation of motion and the sound channel nonhydrodynamic modes evaluated at zero momentum to hold when the equation of motion for Π µν features derivatives along directions orthogonal to U .As we will see later in Sec.IV, there is strong evidence that holography does not display the structure indicated in Eq. ( 31): derivatives in directions orthogonal to U appear in the singulant equation of motion.Motivated by this, we introduce a new MIS-like model that also displays this feature.We achieve this by extending Eq. ( 25) to include derivatives along spacelike directions orthogonal to U .

Description of the model
We work in d = 4.Our new model is defined by the following the equation dictating the spacetime evolution of Π µν , where the new term LΠ µν corresponds to the symmetric, transverse and traceless part of ( Note that, in principle, adding this term is only allowed by hyperbolicity once we go beyond the first-order equation of motion obeyed by Π µν in BRSSS theory.In this work, we will not establish that the initial value problem associated to the conservation equation and Eq. ( 32) is well-posed at the fully nonlinear level.Despite this, we emphasise that (i) as we discuss in detail in Appendix E in the Supplemental Material there exists a parametric regime for η/s, c L , c σ , and Ω such that linearised perturbations around thermal equilibrium are causal and stable in a general reference frame3 and (ii) when performing numerical simulations of the nonlinear problem within such parametric regime, we have seen no issues arise.
We conclude our presentation by clarifying that Eq. ( 32) is not the most general equation of motion for Π µν featuring second-order derivatives that one can write down.For instance, the term D µ D α Π αν can certainly be incorporated into the left-hand side.Since Eq. ( 32) already serves the purpose we outlined, the question of finding out what the most general version of this model is will not concern us here, although it might be relevant for the phenomenological modelling of scenarios like Ref. [51], where MIS or BRSSS approaches are known to fail.

Gradient expansion
In model (32), the recursion relations obeyed by Π (n) are related to the HJSW ones in a straightforward manner.To wit, Π (1) and Π (2) obey the HJSW recursion relations while, for n > 2, one has that where [HJSW] represents the right-hand side of Eq. ( 27c) and As in the BRSSS and HJSW cases, in this model we also find that the gradient expansion is factorially divergent for all our choices of parameters, initial data and spacetime point.To characterise the large-order behaviour of the gradient expansion, we turn now to the question of the singulant trajectory.

Singulant trajectory
According to the factorial-over-power ansatz and Eq. ( 34), the singulant equation of motion in the new model is where This equation features derivatives along Z and, as a consequence, U (χ) stops being related to the nonhydrodynamic sound modes at zero spatial momentum.We now turn to verifying that Eq. ( 35) governs singulant motion in our numerical examples.In order to solve Eq. ( 35) it is no longer enough to specify integration constants per flow line, due to the appearance of derivatives transverse to the flow line Z(χ).This term did not appear in the BRSSS and HJSW models.Now singulant initial data must be specified for some portion of a Cauchy surface, preventing us from testing Eq. ( 35) by fitting one or two complex numbers.To sidestep this issue, we attempt to extract Z(χ) 2 , U (χ), U (χ) 2 for the dominant singulant and test that Eq. ( 35) holds as an algebraic relation at each spacetime point.To extract Z(χ) 2 , U (χ), U (χ) 2 we utilise the following relations, , (36) and similarly for Z(χ) 2 as n → ∞, and confirm Eq. ( 35) in Fig. 5, up to some scatter associated with this numerical procedure.Note that in this figure the new term Z(χ) 2 makes a contribution of similar magnitude as those of U (χ), U (χ) 2 .are extracted from the large-order behaviour of Π (n) and its derivatives, as described in the text.Here the black disks are rendered partially transparent to convey their density.

Singulant trajectory and linear response theory
Despite the fact that the singulant trajectory is not controlled by the nonhydrodynamic sound modes at zero momentum, the lesson that the singulant equation of motion can be understood in terms of a linear response theory problem stands.Let us set the hydrodynamic fields to their values at thermal equilibrium in the rest frame of the fluid, and consider linearised perturbations of Π , δΠ , around this static state.At leading order, one finds that δΠ satisfies For a plane-wave perturbation of the form Eq. ( 38) reduces to a second-order polynomial equation for ω, and the map transforms Eq. ( 40) into the singulant equation of motion (35).The above analysis is nothing but a particular illustration of the general results from Sec. II.A crucial observation that has to be emphasised is that, just as we foresaw at the end of that section, the roots of Eq. ( 40), which determine the singulant trajectory, do not correspond to the sound channel nonhydrodynamic modes when the momentum is finite.The latter (including the hydrodynamic modes) are given by the roots of the following fourth-order polynomial in ω, and only include the roots of Eq. ( 40) when k = 0. Despite the fact that they do not correspond to nonhydrodynamic modes, the roots of Eq. ( 40) still allow for a physical interpretation.This physical interpretation follows from the observation that, due to rotational invariance, the sound channel dispersion relation is constrained to take the form [52] where γ s is a momentum-dependent sound attenuation length.Indeed, the sound channel dispersion relation (42) can be put in the form (43) with the identifications H(ω, k 2 ) = − 1 3 and It is immediate to see that the roots of Eq. ( 40) correspond precisely to the poles of γ s .It is important to keep in mind that, for a general CFT in the linear response regime, γ s can be defined without resorting to Eq. ( 43) by the relation where E 0 is the equilibrium energy density and δ û represents an infinitesimal perturbation of the fluid velocity in the sound channel: With the definition (45), the equation that the sound channel modes have to obey always takes the form (43) with H = −1/3.We conclude our analysis with two comments.The first and most important one is that the connection between the singulant equation of motion and the poles of γ s implied by the map (41) is not restricted to the case at hand: it applies to the BRSSS and HJSW models as well.However, in the BRSSS and HJSW cases, the poles of γ s are independent of the spatial momentum k.In this situation, these poles have to coincide with the sound channel nonhydrodynamic modes evaluated at zero momentum.This observation is in agreement with our analysis of the singulant equation of motion in the BRSSS and HJSW theories.
The second comment is that, when one restricts to Bjorken flow [53], the singulant equation of motion is always insensitive to the difference between the nonhydrodynamic sound modes and the poles of γ s .The reason is that for Bjorken flow χ is a function of the proper time τ = √ t 2 − x 2 alone.This implies that Z(χ) = 0 and, according to the map (41), that the relevant poles of γ s are the ones with k = 0, which agree with the zero-momentum sound channel nonhydrodynamic modes.

A. Introduction
In this section, we examine the large-order behaviour of the gradient expansion for a longitudinal flow in the context of holography. 4 We do not attempt to perform the numerical computation of the gradient expansion evaluated on a particular state.Rather, we take advantage of the asymptotic ansatz (14) to prove that a large-order factorial growth is allowed.We remind the reader that finding out that the asymptotic ansatz ( 14) is consistent in holography does not imply that the gradient expansion is necessarily factorially divergent; rather, this observation just demonstrates that such large-order behaviour is, in principle, possible.
Our study builds on previous developments in the AdS/CFT context.The first one is the original construction of fluid/gravity duality [70], which we employ to build the gradient-expanded constitutive relations at the fully nonlinear level.The second one is the analysis of the exact constitutive relations of the microscopic CFT in the linear response regime [71,72].The bridge between both approaches is the linearisation enacted by the asymptotic ansatz (14).
The key result that we show in this section is as follows.If a large-order factorial growth is present, then the singulant equation of motion in holography is given by the poles of the momentum-dependent sound attenuation length, γ s , under a map analogous to that which we saw in the MIS-like models (41).

B. Longitudinal flows in holography
To construct the geometry dual to a longitudinal flow, we follow Ref. [73] and put forward the following metric ansatz, where x µ represents the boundary coordinates, U µ as elsewhere in the text is the unit-normalised fluid velocity and G µν is transverse, It is convenient to employ a curvilinear coordinate system to describe the longitudinal flow.We focus on longitudinal flows in four-dimensional Minkowski spacetime.We take our boundary coordinates to be where τ and σ are, respectively, a timelike and a spacelike coordinate that parameterise the longitudinal plane.We choose this curvilinear coordinate system in such a way that the boundary metric is diagonal, and the fluid velocity reads The orthonormal longitudinal vector field Z is thus given by Imposing flatness of the boundary metric (49) leads to an equation linking a and b.Finally, to comply with the symmetry restrictions of the flow, we take where V τ , V σ , Σ, and B are functions of τ , σ, and r.
The transversality condition G µν U ν = 0 is automatically satisfied since G τ τ = 0.The energy-momentum tensor of the dual CFT is dictated by the holographic dictionary [74] in terms of the coefficients of the series expansion of the metric around the asymptotic boundary located at r → ∞.In our coordinate system for a given quantity f , we denote the coefficient of r −n in the near-boundary (large-r) expansion as f n .A straightforward computation shows that 5 where the only nonzero components of t µν are given by and we have defined To conclude our analysis, we impose the Landau frame condition by demanding that the fluid velocity U = e −a ∂ τ is the only timelike eigenvector of T µ ν .This can only be achieved provided that Once the Landau frame condition ( 56) is imposed, the eigenvalues of T µ ν are in one-to-one correspondence with its diagonal components.One has that the energy density E, the longitudinal pressure P and the transverse pressure P ⊥ are given by Hence, Π = B 4 − 1 2 Φ.This relation can be simplified further by recalling that the metric ansatz (47) is invariant under and that, in the r → ∞ limit, implying that Σ 0 (x) can be tuned to any desired value by a suitable choice of f (x).In our case, the gauge choice 5 We work with the normalisation L = 4πG = 1.
sets Φ = 0 and implies that Π can be directly identified with B 4 : The existence of this relation is one of the main advantages of our metric ansatz and coordinate choice, as it implies that the gradient expansion of Π can be obtained from the near-boundary behaviour of the gradient expansion of B in a straightforward manner.

C. Gradient expansion
In order to construct the gradient expansion of Π we follow the strategy we described in Sec.II.As it is standard in fluid/gravity duality, the first step is decomposing the Einstein equations into dynamical equations and constraint equations.The dynamical equations can be employed to express Π µν as a functional of E and U ; they can be interpreted as the AdS/CFT counterpart of the equation of motion for Π µν in MIS-like theories.In our case, the dynamical equations are On the other hand, the constraint equations enforce the conservation of the energy-momentum tensor of the dual CFT.
Once the dynamical equations have been identified, we take the metric (47) and introduce the bookkeeping parameter by means of a homogeneous rescaling of the boundary coordinates: Then, we express V τ , V σ , Σ, and B as asymptotic series in , insert them into the dynamical Einstein equations, and demand that they are a solution order by order in an expansion around = 0.This procedure transforms the dynamical equations into recursion relations for σ , Σ (n) , and B (n) .The zeroth-order solution reads where we taken into account that the bulk spacetime is asymptotically AdS with boundary metric (49).Note that a, b and r h are functions of τ and σ only.
The gradient expansion of Π can be obtained from Eq. ( 64) through the relation (61), provided that appropriate boundary conditions are satisfied by the nth order solution.These boundary conditions are divided into two sets, infrared and ultraviolet, depending on whether they are imposed at the event horizon of the zeroth-order solution (65), located at r = r h , or at the asymptotic boundary, located at r = ∞.Our infrared boundary conditions are regularity of σ , Σ (n) , and B (n) at r = r h .Since we are working in Eddington-Finkelstein coordinates, demanding regularity of our metric functions at r = r h corresponds to imposing infalling boundary conditions.As for the ultraviolet boundary conditions: 1. We demand that the metric is asymptotically AdS 5 .
Since this boundary condition has already been fully taken into account by the zeroth-order solution (65), it follows that for n ≥ 1.We note that, when imposing Eq. ( 66), the dynamical Einstein equations imply that, in the r → ∞ limit, V = O(r −2 ) for n > 2 and B (n) = O(r −4 ) for n > 0. This observation will be important later.
2. We require that the Landau frame condition (56) is obeyed at every order in the gradient expansion.This requires that for n ≥ 0. When combined with the last observation performed in point (1) above, Eq. ( 67) implies that 3. We demand that the energy density as computed from the zeroth-order solution, 3 4 r 4 h , agrees with the actual energy density of the system E.This requires that for n ≥ 1. r h is thus related to the effective temperature T as πT = r h .
4. We enforce the condition (60).Since this relation is first order in gradients, we must impose it at the level of Σ (1) .This implies that for n = 1.

D. Large-order behaviour of the gradient expansion
In order to assess the large-order behaviour of the gradient expansion (64), we assume that the singulant fields corresponding to V τ , V σ , Σ, and B are equal and, subsequently, find that this assumption is self-consistent.Hence, we take Vτ,q (r, τ, σ) Γ(n + α q (r, τ, σ)) χ q (τ, σ) n+αq(r,τ,σ) , (70) with analogous expressions holding for V σ , Σ, and B. In the end, this restriction is equivalent to the assumption that the gradient expansion of the full spacetime metric (7) takes the asymptotic form (8).
To simplify the subsequent argument we introduce Ψ, which stands for any of the functions V τ , V σ , Σ and B. According to the general analysis presented in Appendix A in the Supplemental Material-and whose most important take-home points were mentioned in Sec.II-the terms in the dynamical equations that set the singulant equation of motion are of the form where p, m are non-negative integers.Since the dynamical equations are second order in spacetime derivatives, one has that p + m ≤ 2. Furthermore, it is very important to stress that, due to the functional form of the dynamical equations, the large-order ansatz (8) implies that χ is independent of the radial coordinate r.
Let us discuss E rr first.One finds that, at leading order in n as n → ∞, This relation, when combined with the boundary conditions spelled out before, entails that Taking into account this fact, the remaining dynamical equations imply that, at leading order in n as n → ∞, where f (0) = r 2 − r 4 h r −2 .Note that Vτ decouples from Vσ and B.
Let us assume that the singulant field χ is known at the τ = 0 slice.In order to integrate the singulant motion, we need to find ∂ τ χ.This can be achieved by solving the eigenvalue problem for ∂ τ χ posed by Eqs.(74a) and (74b) with boundary conditions given by Eqs. ( 66) and (67).We note that the singulant dynamics is ultralocal in the boundary coordinates: to compute ∂ τ χ at a given spacetime point, the inputs that we need are the values of ∂ σ χ and a, b, and r h -or, equivalently, the energy density and velocity of the fluid -at the point.This ultralocality is in line with the general results discussed in Sec.II.Furthermore, in the next section, we also confirm that the eigenvalue problem spelled out above has a natural counterpart in linear response.

E. Singulant equation of motion and linear response theory
In Refs.[71,72] it was shown that, in the linear response regime, it is possible to find the exact constitutive relations that express the dissipative tensor as a functional of the hydrodynamic fields in the Landau frame.Consider infinitesimal perturbations of the hydrodynamic fields and Π µν around a reference thermal state of energy density E 0 and fluid velocity U 0 = ∂ t .These infinitesimal perturbations, which we denote as δE, δ u, and δΠ µν , are defined as where and we have enforced the Landau frame condition.Latin indices range from one to three and refer to spatial directions of the boundary.The main result of Refs.[71,72] is obtaining the constitutive relations that express δΠ ij as a functional of δE and δu i in closed form.In momentum space, where δ Πij is the Fourier transform of and η(ω, k 2 ), ξ(ω, k 2 ) are momentum-dependent transport coefficients.As explained in Refs.[71,72], these momentum-dependent transport coefficients are computed by solving a system of four coupled radial ODEs (ordinary differential equations) in a black brane background.
For a sound wave δ u is parallel to k and one can employ rotational invariance to set δ ûi = δ û δ i,1 , k i = k δ i,1 with no loss of generality.This implies that πij = − 1 2 k 2 σij and, hence, Recalling the general definition of γ s given in Eq. ( 45), we have that Equation ( 80) shows that γ s is a linear combination of the dynamical transport coefficients η and ξ originally defined in Refs.[71,72], and indicates that the method put forward there can be straightforwardly modified to compute γ s directly.Both this computation and the results it leads to are described in detail in Appendix C in the Supplemental Material.For our purposes here, it suffices to mention the following.
1.For fixed k ∈ R, γ s is a meromorphic function of ω ∈ C, with an infinite number of simple poles, Ω (±) q (k), symmetric around the imaginary ω axis, Ω q (k) is always negative.
2. These simple poles can be computed by solving the following eigenvalue problem: where f = r 2 − µ 4 r −2 .The boundary conditions to be imposed on P and Q are regularity at the black brane horizon located at r = µ and the nearboundary behaviour The two observations above are sufficient to state the main result of this section: if the gradient expansion grows factorially and the singulant field is independent of the metric component under consideration, then the singulant equation of motion is determined by the poles of γ s , just as it was the case in the model discussed in Sec.III D. This fact follows immediately from the realisation that the map, transforms Eqs.(81a) and (81b) into Eqs.(74a) and (74b) and that the infrared and ultraviolet boundary conditions obeyed by P, Q and V (n) σ , B (n) are the same for n > 2. The existence of this map conforms to the general expectations put forward in Appendix A in the Supplemental Material and summarised in Sec.II.

V. KINETIC THEORY
In kinetic theory, we take a factorial ansatz for the distribution function itself, and derive constraints on the singulant from the Boltzmann equation.We cross-check our results in RTA kinetic theory for Bjorken flow, where we can compute precise numerical solutions.We also reanalyse the 1/w expansion in this new perspective, extending results from previous work on the large-order behaviour of moments [18,20,29].

A. Singulants of the distribution function
The Boltzmann equation is which together with the ansatz (9) implies that f (n) is of order n in gradients of f (0) .With the factorial-over-power ansatz (10) the resulting linearisation implies that up to subleading corrections in 1/n, where C is the linearised collision operator: Differently from the other theories, the factor A/χ n does not cancel out, so the large-order equation seems to depend explicitly on n.
Let us sidestep the issue of defining what A/χ n converges to.In the limit, we see that the action of C must be identical to the action of −p µ ∂ µ χ(x, p).The latter does not couple different values of p. Thus, the equation of motion for the singulant can be determined by finding those distribution functions where C effectively acts diagonally in p.This is automatically fulfilled in the simplest RTA model, 6 where leading to We can write the solution as where χ FS is any solution to the free-streaming Boltzmann equation.

B. Moments
Consider some moment, whose nth term arises from the nth term of the distribution function as With the singulant ansatz for f (n) , this is To evaluate this in the limit when n → ∞, we can try a saddle point integration.Then the leading contributions come from those points p s where ∇ p ln χ(x, p)| p=ps = 0.This heuristic argument suggests that the large-order behaviour of the moment is governed by the saddle point singulant χ s (x) ≡ χ(x, p s ).This implies that where α s = α(x, p s ) − j/2.A s gets contributions from A, the Hessian of χ, and the measure [76].
We also note with this argument, the tensor structure is not sensitive to n, which is consistent with the general ansatz in Eq. ( 3).

C. Gradient expansion in Bjorken flow
We now consider Bjorken flow in RTA kinetic theory, where the singulant equation of motion can be verified using the results of Ref. [29].With a relaxation time τ rel = γ/T the Boltzmann equation is where p T is the magnitude of the transverse momentum, v = tP L − zE, and There is a simple recursion relation for f (n) : ) and the singulant satisfies which is just Eq. ( 88) specialised to this model.Equation ( 92) can be integrated to give an equation for the temperature T (τ ) [77,78], which can be used to get precise numerical solutions.Together with the recursion relation this allows access to large orders of f (n) .
In Fig. 6, we verify that the singulant equation of motion is satisfied, and that there are momentum-dependent singulants.for three values of the momentum (v, pT ).One singulant seems to be approximately momentum independent, and is the dominant one.Another depends on their ratio.Bottom panel: test of the singulant equation of motion.χ can be measured numerically at each point in time by calculating the gradient expansion.Alternatively, we can take the dominant singulant from the top panel and evolve it by Eq. ( 94).The figure shows that these procedures agree.

D. 1/w expansion revisited
Another expansion, used in previous works [20,29], is to expand the pressure anisotropy in inverse powers of the variable w ≡ τ T .The essential difference between these expansions is that the latter one uses the conservation equations to remove all appearances of ∂ τ T ; see Ref. [34] for details.This expansion leads to different singulants from those in previous section.
It is convenient to change variables to and we now expand as with f0 = e − √ p1+p2 .Here tildes denote series coefficients in the 1/w expansion.The pressure anisotropy is also expanded as The Boltzmann equation leads to the following recursion relation for the expansion coefficients: This recursion can be used to calculate several hundred orders of fn and ãn analytically.The ãn coefficients were observed to diverge factorially in Ref. [20].Numerical studies of these series show that the coefficients fn (p 1 , p 2 ) diverge factorially, with a slope that depends on the values of p 1 and p 2 .We are therefore led to the ansatz with a momentum-dependent singulant.The recursion relation implies that χ satisfies where the terms involving the coefficients of the pressure anisotropy have dropped out.This equation has the solutions where G is an arbitrary function.G is in principle determined by matching a sum of singulants to the zeroth order of the perturbative series, but we do not know how to do this in practice.However, even without knowing the function G, the values of χ(p 1 , p 2 ) on a curve where p 2 √ p 1 is constant can be calculated from any other.We verified this using the explicit large-order computation; see Fig. 7.
For the pressure anisotropy, let us apply the saddle point argument from above.Equation (101) immediately implies that χ * = 3 2γ at a saddle point.Thus, even though G is unknown, and therefore the saddle points p * cannot be determined, the value of the singulant at the saddle point is uniquely fixed.However, it is also known that there are contributions with a nonzero imaginary part [20], which do not correspond to physical excitations.We do not know why they do not show up in this analysis, but a possibility is that they come from end point contributions of the saddle point integral.
Finally, note that introducing A(w) simplifies a lot of analysis in the Bjorken flow and makes some features  102), the value of the singulant χ(p1, p2) can be predicted along a certain trajectory in the space of (p1, p2).We verify this by calculating the singulants for two points on such a curve.These singulants are shown in the Borel plane and the two sets of parameters correspond to the blue and green points in the figure .From the singulant revealed by one set, we predict the location for the other parameters and confirm this prediction using the explicit large-order computation.
manifest.For example, it led to the notion of hydrodynamic attractors.Therefore, it has been an interesting and potentially important question if A(w) acquires a natural generalisation when symmetries of the flow are relaxed.While generalisation of A can be trivially obtained as long as the local rest frame exists, the w clock variable is less obvious.Earlier attempts in this direction include [79,80] and our work adds singulants as a particularly natural possibility.The reason for it is twofold: the leading singulant governs the asymptotics of the gradient expansion and, through resurgence, also the most significant exponentially decaying term supplementing optimally truncated (see the next section) hydrodynamic constitutive relations, as w also does.

VI. OPTIMAL TRUNCATION OF THE GRADIENT EXPANSION
Having established that the gradient expansion ( 5) is factorially divergent in MIS-like theories, the natural question that arises is what its practical usefulness is.There are two different ways in which one can employ a factorially divergent gradient expansion: a fixed-order truncation or an optimal truncation.In this section, we explore the second option.We remove the bookkeeping parameter by setting it to one.
The optimal truncation of the gradient expansion at a given spacetime point is the partial sum closest to the actual value of Π (t, x).Note that, due to this definition, the order of optimal truncation n opt is expected to be spacetime dependent.
Our main objective in this section is to put forward a criterion for estimating n opt that relies exclusively on the gradient expansion itself.Our choice is the following.Let |χ d (t, x)| correspond to the absolute value of the dominant singulant at a given spacetime point.We propose to estimate the order of optimal truncation at that spacetime point by the relation where the brackets instruct us to take the integer part of the quantity they enclose.
To explore the consequences of Eq. ( 104), let us consider the case in which the large-order behaviour of the gradient expansion is of the form (6)  is the smallest coefficient of the gradient expansion, and, furthermore, that this smallest coefficient is exponentially suppressed in |χ d |: The standard expectation is that S (nopt,est) provides a representation of the actual Π which is accurate up to such exponentially small term: In the case that the dominant singulant is not real but rather corresponds to a complex-conjugated pair, the norm of Π (n) displays additional oscillations at large n with frequency arg χ, where we assumed that α is real.In this situation, the definition (104) singles out the gradient expansion coefficient for which the envelope of these oscillations is the smallest.
To test our estimate, we consider the following initial data in BRSSS theory: which correspond to an initial overdensity in the effective temperature such that |T (0,0)−T (0,∞)| T (0,∞) = 1.We have set E/T 4 = 1, η/s = 1/(4π) and τ Π T = (2 − log 2)/(2π).We evaluate the gradient expansion along the flow line passing through x = 0 at t = 0. Our results will be parameterised in terms of the dimensionless variable At late times, a complex-conjugated singulant pair, whose trajectory corresponds to the solid red line, takes over.The discontinuous solid black line corresponds to nopt,est, as given by Eq. (104).Middle panel: coefficients of the gradient expansion at three selected times.The ones associated to the corresponding nopt,est have been highlighted in red.Lower panel: absolute error of the optimal truncation (black crosses) and the truncation at order nopt,est (red dots).Crucially, the latter error provides an upper bound for the former that displays the same time dependence.The dashed (dotted) black line represents the absolute error of the first-order (secondorder) order truncation.
In the upper panel of Fig. 8, we represent the time evolution of the dominant singulant and the corresponding n opt,est as we move along this flow line.This dominant singulant has been extracted with the help of the method put forward in our analysis of the BRSSS model in Sec.III.The reader can find where the gradient expansion coefficient Π (nopt,est) picked by our procedure sits at three different times in the middle panel of Fig. 8.
Finally, in the bottom panel of Fig. 8, we plot the time evolution of the error of our estimate for optimal truncation and compare it with the error of the actual optimal truncation as originally defined.The latter is always upper bounded by the former and, furthermore, both errors decrease exponentially as ξ grows.For reference, we also include the absolute errors incurred by the first and second-order truncations of the gradient expansion.

VII. OUTLOOK
We have established that there exists a deep connection between far-from-equilibrium nonlinear relativistic hydrodynamics and linear response around global thermal equilibrium.This connection is manifest in the large-order behaviour of the hydrodynamic gradient expansion and is naturally encoded in singulants.When viewing this duality from the gradient expansion side, one can regard singulants as emergent excitations coming from a dramatic reorganisation of the series at large orders.We have outlined the principles governing singulant dynamics in longitudinal flows for MIS-like models, holography and kinetic theory, and checked their validity by explicit numerical computations in a number of cases.
In the course of this work we have also investigated the large-order behaviour of the gradient expansion and singulants in non-relativistic theories along the lines of [81,82], finding results analogous to the ones presented here.This is beyond the scope of the current paper and we hope to report on it elsewhere [83].This suggests that singulants are a useful language for hydrodynamics in general.
Our work opens new research directions which come with new technical and conceptual challenges to be addressed.On the technical front, it would be interesting to solve for nonlinear longitudinal flows as an initial value problem in holography and RTA kinetic theory and confirm the predicted behaviour of singulants by analysing the gradient expansion at large order (just as we have done here for MIS-like models).On the conceptual front, we singled out four questions that we believe are of particular importance and discuss them in the sections below.

A. Singulants and initial conditions
While we have succeeded in determining the time evolution of individual singulants, we have not discussed how the number of singulants and their initial values are related to the initial out-of-equilibrium state under consideration.The complete answer to this problem remains unknown to us; we have, however, made partial progress on this front by working in the linear response regime.In this case, the problem of finding initial conditions for the singulant fields can be solved, and we refer the reader to Appendix D in the Supplemental Material for a worked-out example in BRSSS theory.
Revisiting this question for nonlinear flows is of fundamental importance in order to understand how to promote the gradient expansion (1) to a full-fledged transseries representation of the energy-momentum tensor.If there were a method to extract the values of the singulant fields at t = 0 from the initial data, then the singulant equation of motion would allow us to determine the dominant singulant at any given point without having to compute the gradient expansion itself.In light of the results presented in Sec.VI, the natural expectation that arises is that full knowledge of the spacetime profile of the dominant singulant would allow us to single out the spacetime region where relativistic hydrodynamics, when truncated at low order, will not be applicable-for instance, because |χ d | ∼ 1 and no optimal truncation even exists.Finally, it would be interesting to find out whether making progress on this front could shed light on the necessary conditions that the initial data have to obey for a factorial growth to be present, and the singulant ansatz (3) to be applicable.Further discussion on necessary conditions for factorial growth can be found in [34].

B. Singulants beyond longitudinal flows
Another important question we have not discussed up to this point is whether the results presented in this work generalise beyond longitudinal flows.The main simplification brought by restricting ourselves to this class of fluid flows was that Π µν could be described in terms of a single scalar field, Π .We have seen that a longitudinal flow is nothing but a nonlinear sound wave.Because of this, and the linearisation entailed by the asymptotic ansatz (14), the singulant equation of motion was determined by a linear response theory problem formulated in the sound channel.In the absence of symmetry constraints, a d-dimensional conformal energy-momentum tensor is associated with a Π µν that encompasses d(d−1)/2−1 degrees of freedom.It is natural to wonder whether, in this completely general situation, the leading large-order behaviour of Π (n) µν can be expressed as a linear superposition of contributions defined in independent channels7 and, furthermore, if these conjectured channels can be put into one-to-one correspondence with the tensor, vector, and scalar channels that arise when decomposing an infinitesimal perturbation of the energy-momentum tensor around thermal equilibrium.
In Appendix B in the Supplemental Material we provide the first steps to address this question in the context of MIS-like theories.The analysis presented there shows that, for a decomposition along the lines of Eq. ( 112) to be possible, one would need to work with a basis for the tangent space at a given point that depends explicitly on the singulant fields themselves.A fully fledged numerical analysis beyond the realm of longitudinal flows in these models along the lines of Ref. [34] is required to find out whether this conclusion is correct.This computation will be numerically more costly than the longitudinal flow one, but we expect it to be feasible.
We would like to point out that, when it comes to this issue, kinetic theory stands aside the other models we have considered.In the kinetic theory case, we know that the large-order behaviour of the gradient expansion of the distribution function does not decompose into channels.It would be interesting to understand whether, in spite of this, computing the second-order moments to get the energy-momentum tensor leads naturally to a channel decomposition.
While this discussion might sound a bit technical, it does have interesting physical implications.The key reason for this is the idea of rheology, which aims to improve hydrodynamic predictions by making transport coefficients depend on gradients of fluid variables, see in particular [84] for a pioneering effort in this direction.Recently it was revisited in [12,79,[85][86][87] by including contributions from transients in redefined transport coefficients.The true potential of the latter approach depends on correctly accounting for transient effects, which in our language are simply the singulant fields.In particular, the extent to which this is feasible in general depends crucially on understanding the detailed tensor structure of the nonhydrodynamic sectors of the transseries representing the energy-momentum tensor.Our work opens this possibility.

C. Singulants in other gradient expansions
Finally, we would like to point out that Eq. (1b), while defining classical relativistic hydrodynamics as an effective field theory, is not the only gradient expansion one can work with.Our analysis of the 1/w expansion in RTA kinetic theory already illustrates that singulants provide novel insights beyond the realm of the gradient expansion (1b), and we expect this observation to be completely general.
To further illustrate this point, let us stay within the realm of longitudinal flows and consider a scenario in which, besides Π , both E and U are gradient-expanded, with the conservation equation ∇ µ T µν = 0 being solved in a small-expansion.In this alternative gradient expansion, E(0) and Ů (0) represent the solution of ideal hydrodynamics associated with the initial spatial profiles of the energy density and the velocity field, while higherorder terms encode dissipative corrections to the ideal solution which vanish at t = 0. We refer to the gradient expansion (113) as the ideal expansion.In the context of Bjorken flow, the ideal expansion is precisely the expansion in inverse powers of the proper time τ , and the 1/w expansion arises as a (partial) resummation of it.
Let us assume that the large-order behaviour of E(n) , Ů (n) , and Π(n) are governed by a common singulant field.The general rules we put forward in Appendix A in the Supplemental Material entail that the singulant equation of motion is connected to a linear response theory problem through the map (A22).In this regard, the fundamental differences between the gradient expansions (1) and (113) are twofold.First, the role that the microscopic E and U had in the map associated to the gradient expansion ( 1) is taken over by E(0) and Ů (0) when the ideal expansion is considered.Second, the particular linear response theory problem associated to the ideal expansion is different from the one relevant for the gradient expansion (1): rather than being determined by the poles of γ s , the singulant equation of motion in the ideal expansion is set by the transient sound channel modes.The reason is that, to compute the ideal expansion, the conservation equation is also gradient expanded.
Another expansion to which one can apply some of the techniques developed in this work is the Cauchy data expansion.In this case one employs the conservation equations to systematically replace longitudinal derivatives in the constitutive relations by transverse ones.The Cauchy data expansion is known up to third order in gradients at the fully nonlinear level [88,89], and its large-order behaviour was analysed in Ref. [35] for conformal fluids in the context of linearised hydrodynamics, where it only involves spatial derivatives-hence our nomenclature.The aforementioned replacement, while inconsequential from the effective field theory perspective, makes a difference when the constitutive relations are evaluated on a fluid flow.Given the results of Ref. [35], we expect the Cauchy data expansion to be factorially divergent at the nonlinear level, and it would be interesting to explore in detail the link between singulant dynamics and linear response theory in this case.
Finally, we want to mention that our general expectation is that different frame choices, i.e., different choices of collective fields, would result in singulants with different dynamics.This should not come as a surprise, given the previous discussion on the ideal expansion.In fact, as the causal completions of relativistic hydrodynamics recently put forward by Bemfica, Disconzi, Noronha [90] and Kovtun [91] (BDNK) illustrate, there exist phenomenological models for which even the very existence of the gradient expansion ( 1) is contingent on the choice of collective fields.To make meaningful comparisons between gradient expansions across different models, one has to keep the choice of collective fields invariant.In this regard, it might be worthwhile to perform a field redefinition in the BDNK theories (or the theory put forward in Ref. [92]) to put their gradient-expanded constitutive relations in the Landau frame and check whether the singulant dynamics in longitudinal flows is still controlled by the poles of γ s .

D. Physical interpretation of singulants
For systems near thermal equilibrium a particularly useful theoretical language is provided by linear response theory.To the leading order in the amplitude of distortions away from equilibrium, the dynamics of the system is determined by the properties of retarded two-point functions of local operators.Singularities of these correlators in Fourier space determine the dynamics of the system.In particular, black hole quasinormal modes in holography arise as single pole singularities of boundary correlators.MIS-like theories also have such single pole singularities, whereas in kinetic theory it is known that this singularity structure is richer and includes also branch points.
In holography (with a straightforward generalisation to MIS-like theories), short-lived quasinormal modes describe the decay of non-hydrodynamic degrees of freedom with the subsequent approach to equilibrium being governed by a few long-lived quasinormal modes realising linearised hydrodynamics.It is fair to say that quasinormal mode picture underlies much of progress achieved in studying strongly coupled systems in the course of the past two decades.However, outside the realm of linear response theory and in the absence of high degree of symmetries, it has never been clear what object, if any, arises as a natural generalisation of a quasinormal mode.We advocate that singulants can fill this role.
They are quasinormal-mode-like because they are excitations which describe non-hydrodynamic decay, reflected in the fact that in all examples considered in this paper, the value of |χ| grows after enough time has passed.We have not observed any examples where χ orbits the origin of the Borel plane, which is fundamentally tied to this interpretation.At the same time, they do not call for a proximity of global equilibrium.They differ from standard quasinormal modes in the sense that they apply arbitrarily far from equilibrium.
Finally, for singulants to fill the role of quasinormal modes for far-from-equilibrium systems, then what is their imprint on a given flow?We will now outline a class of measurements that can be performed, at least in principle.The data we assume is a measurement of the energy-momentum tensor T µν in a given spacetime region.Given this data one can solve the eigenvalue problem T µν U ν = −EU µ at each spacetime point to obtain the fields E and U .Taking derivatives of E and U one can assemble the structures we have given in detail and directly observe singulant dynamics, just as we have done through numerical simulation in the preceding sections.
In the special case of Bjorken flow, see Sec.V D, the hydrodynamic gradient expansion is the same for all nonequilibrium states.Therefore, another way of accessing the singulant dynamics in this case is to construct A(w) for two different flows and subtract them from one other.In this case the hydrodynamic sector of the transseries cancels completely, leaving singulant physics visible to the naked eye [93].
As we discussed in section VII C one has the freedom to vary the type of gradient expansion.This leads to different singulant dynamics which can be observed by constructing combinations of derivatives of E and U appropriate for the gradient expansion in question.In this way, the freedom to choose the gradient expansion leads to an infinite family of related self-consistent predictions.
Our goal is constructing the five-dimensional geometry dual to a sound wave of infinitesimal amplitude propagating on a background thermal state.Since this sound wave is a longitudinal flow, we take the metric ansatz (47) and write where f = r 2 − µ 4 r −2 is the blackening factor of the background black brane.We parameterise the longitudinal plane with Minkowski coordinates t and x.The perturbations δV t , δV x , δΣ, and δB are functions of r, t, and x, while λ is a bookkeeping parameter that reminds us that these perturbations are infinitesimally small.Due to the translational invariance of the background black brane, it is convenient to work in momentum space.Hence, we take and similarly for δV x , δΣ, and δB.With our choices, the Einstein equations read at leading order in λ.As dynamical equations, we take the following linear combinations of δ Êµν , The requirement that the boundary metric is the Minkowski metric entails that A 1 = 0. Furthermore, it is always possible to shift A 2 by the coordinate change r → r + λΨe −iωt+ikx .Therefore, A 2 can always be set to zero and, as a consequence, we take δ Σ = 0 with no loss of generality.Once the latter condition is imposed, the remaining dynamical equations read Note that δ Vt decouples from δ Vx and δ B. According to the dynamical equations, the asymptotic series expansion of the metric is where we have chosen the boundary conditions at r = ∞ in such a way that the bulk spacetime is asymptotically AdS.From these asymptotic series expansions, holographic renormalisation gives a dual energy-momentum tensor of the form T µν = T µν + λδ Tµν e −iωt+ikx , with and with the remaining components vanishing.The Landau frame condition demands that, up to O(λ), the fluid velocity U µ = (1, λδ ûe −iωt+ikx , 0, 0) is a timelike eigenvector of T µ ν .This condition can only be satisfied provided that Once that the relation above has been imposed, we can finally identify Recalling Eq. ( 45), we observe that γ s can be obtained directly from the near-boundary series expansion of δ B. The poles of γ s play a central role in the analysis presented in Sec.IV.To compute them, we plug the ansatz into Eqs.(C7a) and (C7b) and expand around ω = Ω p (k).The leading order result corresponds precisely to Eqs. (81a) and (81b) in the main text.These equations, on the other hand, are nothing but Eqs.(C7a) and (C7b) with δ û = 0 after the identification Ω p → ω, Q → δ B and P → δ Vx .This is physically consistent: at a pole of γ s , Eq. ( 45) predicts that Π can be finite even if δ û = 0. Taking into account the asymptotic series expansions (C8b) and (C8c) and the Landau frame condition (C11), one recovers the near-boundary behaviour quoted in the main text: P = O(r −3 ) and Q = O(r −4 ) as r → ∞.

Numerical results
To compute γ s , we need to solve the coupled ODEs (C7a) and (C7b).We proceed as follows.First, we take into account the r → ∞ series expansions (C8b) and (C8c) and perform the field redefinitions, where the Landau frame condition (C11) has been imposed.Note that, according to Eq. (C12) and the definition (45), The field redefinition (C14) transforms equations (C7a) and (C7b) into where we have set µ = 1 and traded the radial coordinate r for In Eqs.(C16), the primes denote derivatives with respect to u.The reader should keep in mind that, due to the choice µ = 1, ω and k are measured in units of πT .
To solve (C16), we resort to pseudospectral methods.We introduce a Chebyshev-Gauss-Lobatto grid in the variable u, transforming the two coupled ODEs (C16) into a single algebraic equation.There are no boundary conditions to be imposed on a(u) and b(u) on the discretisation grid, neither at the asymptotic boundary, located at u = 0, nor at the black hole horizon, located at u = 1.The reason is that the correct u → 0 behaviour has been already taken into account by the field redefinitions (C14), while the requirement of a(u) and b(u) to be regular at u = 1 is automatically incorporated by our choice of spectral decomposition.
Once the numerical solution corresponding to a given pair (ω, k) is known, we can read the corresponding γ s by employing the relation (C15).In Fig. 9, we plot γ s for ω, k ∈ R. Due to the symmetry of equations (C16) under k → −k, we restrict ourselves to k ≥ 0. We see that, for a given k, the real part of γ s is even in ω, while the imaginary part is odd.We refer the reader to Refs.[71,72] for a discussion of earlier numerical results on η and ξ-linked to γ s by Eq. ( 80)-for real ω and k.
The relation between the singulant trajectory and the poles of γ s we uncovered in Sec.IV entails that our main interest does not lie in the behaviour of γ s for real wavevectors, but on its singularities in the complex ω-plane.
We start by demonstrating that γ s is not an entire function of ω at fixed k.To this aim, in the top panel of Fig. 10 we represent the (logarithm of) the norm of γ s as a function of ω ∈ C for k = 0. Our results show that γ s , while being analytic in the upper-half complex ω-plane, has poles in the lower one.At zero spatial momentum, these poles coincide with the nonhydrodynamic sound channel quasinormal modes.
The poles Ω (±) p describe trajectories in the complex ω-plane as k varies.While, for a given k, Ω (±) p could be found from Eqs. (C16) by means of a binary search algorithm that minimises − log(|γ s |), it is more efficient to compute them through Eqs.(81a) and (81b).To solve these equations numerically, we perform the field redefinitions which ultimately transform them into Eqs.(C16) with vanishing source terms on the right-hand-side.The pseudospectral discretisation proceeds in the same way as described before, resulting in a generalised algebraic eigenvalue problem whose solution returns Ω (±) p .In the bottom panel of Fig. 10, we plot the trajectories followed by the three lowest-lying poles in the complex ω-plane.At large k, we always find that Re Ω (±) p → ±k, while Im Ω (±) p approaches zero, in such a way that Ω (±) p approach the real ω axis from below.In this Appendix, we illustrate how the singulant field can be completely determined in the linear response regime.For definiteness, we focus on conformal BRSSS theory, which has been discussed at the fully nonlinear level in Sec.III.
To access the linear response regime, we place ourselves in the fluid rest frame, write Re Ω (±) (orange), Ω (±) 2 (red) and Ω (±) 3 (brown) in the complex ω-plane as k is increased from 0 to 50 in units of 0.5.The poles at k = 0 are represented by the black stars.We have employed a collocation grid of 80 points in order to get these results.
and take the λ → 0 limit.At leading order in λ, the recursion relations simplify to The equations of motion for BRSSS theory in the linear response regime can be solved in closed-form by working in momentum space.One finds that δu(t, x) = R dk e ikx q=+,−,NH δu q (k)e −iωq(k)t (D3) where ± and NH refer respectively to the contributions of the two hydrodynamic modes and the nonhydrodynamic one.δu + (k), δu − (k) and δu NH (k) can be expressed in terms of the Fourier transforms of the initial data, δT (0, x), δu(0, x) and δΠ (0, x).Expressions analogous to (D3) hold for δT and δΠ as well.
Given the form of the linearised recursion relations, it is immediate to see that contributions of individual modes decouple.Focusing on the qth contribution one finds that, in Fourier space, δΠ (1) (k) = 2 3 η(T 0 )ikδu q (k), (D4a) The recursion relations above are immediately solved by Let us assume that the analytical continuation of the Borel transform of the gradient expansion, has a well-defined Fourier transform, δΠ (B) (t, k; ).Then, we can employ (D5) to compute the contribution of the qth mode to δΠ (B) (t, k; ), δΠ ,q (t, k; ), with the final result that δΠ (B) ,q (t, k; ) = 2η(T 0 )(e i τΠ(T0)ωq(k) −1) 3τ Π (T 0 )ω q (k) kδu q (k).(D7) Our main working hypothesis is that, when the Fourier integral ,q (t, k; )e −iωq(k)t (D8) ceases to exist for a particular , the original analytically continued Borel transform δΠ (B) (t, x; ) becomes singular.Hence, it is natural to identify the 's at which this phenomenon takes place with the singulants controlling the large-order behaviour of the gradient expansion.
In the light of expressions (D7) and (D8), the convergence of the Fourier integral (D8) depends on two inputs: • The large-k behaviour of the initial data, which determines the large-k behaviour of δu q (k).
• The large-k behaviour of the mode frequencies.
Regarding the latter, it is straightforward to show that, as k → ∞, in a way compatible with relativistic causality, we note that in the |k| → ∞ limit with k x = α|k| → ∞, |α| ≤ 1, one has that ω ± = ±|k| − iT 0 Ω I (cosh(u 0 ) ∓ α sinh(u 0 )) + . . ., (E8) in such a way that the group velocity is exactly one.The limit in which |k| → ∞ with fixed k x is obtained from the expression above by setting α = 0 and leads to identical conclusions.
We performed a numerical scan of Im ω ± in the threedimensional parameter space spanned by |k|, α and u 0with the remaining parameters as in the numerical simulation discussed in the main text-and did not encounter any point where this quantity became positive.In particular, note that in the local rest frame one always has that Im ω ± = −T 0 Ω I ≤ 0. Therefore, our stability results conform to the findings of Refs.[50,94,95] showing that stability is a Lorentz-invariant notion in a causal theory.

Shear channel
Shear channel modes correspond to the roots of the following cubic polynomial in k • U 0 , and are divided into one hydrodynamic mode and two nonhydrodynamic ones.In the |k| → ∞ limit with k x = α|k|, |α| ≤ 1, these modes behave as the three modes are asymptotically stable.These asymptotic stability conditions are the identical to the ones obtained in the local rest frame of the fluid, since cosh(u 0 ) ∓ α sinh(u 0 ) ≥ 0 for |α| ≤ 1.Finally, we note that in the regime in which |k| → ∞ with fixed k x , the behavior of the shear channels modes is described by (E10) with the replacements α → 0, α|k| → k x .
We performed a numerical scan of Im ω q in the threedimensional parameter space spanned by |k|, α and u 0with the remaining parameters as in the numerical simulation discussed in the main text-and did not encounter any point where this quantity became positive.The shear modes for our parameter choice, real |k|, and u 0 = 0 are plotted in Fig. 12.We have set cL = 1, cσ = 0, η s = 1 4π , ΩI = 2 and ΩR = 1, and worked in T0 = 1 units.The two nonhydrodynamic modes are associated to red curves, while the hydrodynamic mode is associated to black ones.
Regarding asymptotic stability, the modes ω 3,4 impose no new condition since cosh(u 0 )β ± α sinh(u 0 ) ≥ 0 (E16) for |α| ≤ 1 and u 0 ∈ R while the modes ω 1,2 demand that Since Ω I is a positive real number, the condition above is more stringent than the second entry in (E11) and overtakes it.We performed a numerical scan analogous to the shear channel one and did not encounter any instances of Im ω q becoming positive.The sound modes for our parameter choice, real |k|, and u 0 = 0 are plotted in Fig. 13.
J 2 Z W d 3 b / + g e n j U N k m m K W v R R C S 6 G x L D B F e s Z b k V r J t q R m Q o W C c c 3 x Z + 5 4 l p w x P 1 Y C c p C y Q Z K R 5 x S m w h 9 W n M B 9 W a V / d m w K v E X 5 A a L N A c V L / 6 w 4 R m k i l L B T G m 5 3 u p D X K i L a e C T S v 9 z L C U 0 D E Z s Z 6 j i k h m g n x 2 6 x S f O W W I o 0 S 7 U h b P 1 N 8 T O Z H G T G T o O i W x s V n 2 C v E / r 5 f Z 6 D r I u U o z y x S d L 4 o y g W 2 C i 8 f x k G t G r Z g 4 Q q j m 7 l Z M Y 6 I J t S 6 e i g v B X 3 5 5 l b Q v 6 v 5 l / e b + s t Z 4 b M z j K M M J n M I 5 + H A F D b i D J r S A Q g z P 8 A p v S K I X 9 I 4 + 5 q 0 l t I j w G P 4 A f f 4 A H 6 C O n A = = < / l a t e x i t > A < l a t e x i t s h a 1 _ b a s e 6 4 = " J H G Y j A f T m I Q E d 3 B i Y E A E h p T 1 F o o = " > A A A B 8 X i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s y I o O 4 q b l x W s A 9 s h 5 p J M 2 1 o H k O S E c r Q v 3 D j Q h G 3 / o 0 7 / 8 Z M 2 4 W 2 H g g c z r n h n n u i h D N j f f / b K 6 y s r q 1 v F D d L W 9 s 7 u 3 v l / Y O m U a k m t E E U V 7 o d Y U M 5 k 7 R h m e W 0 n W i K R c R p K x r d 5 H 7 r i W r D l L y 3 4 4 S G A g 8 k i y u 8 e c Z 7 8 d 6 9 j 9 l o w Z t X e A h / 4 H 3 + A O y 9 k W Y = < / l a t e x i t > i < l a t e x i t s h a 1 _ b a s e 6 4 = " F e e Q o o i m 9 2 9 / r 5 P u V d 1 r 1 G 8 f G r V m u 6 i j D G d w D p f g w T U 0 4 R 5 a 0 A E C E p 7 h F d 4 c 7 b w 4 7 8 7 H c r T k F D u n 8 A f O 5 w 8 S + J F C < / l a t e x i t > ii < l a t e x i t s h a 1 _ b a s e 6 4 = " N X b z x 6 9 e 4 8 2 / M t L P Q 1 g O B w z n 3 k n N P m H C m j e t + O 6 W t 7 Z 3 d v f J + 5 e D w 6 P i k e n r W 1 T FIG. 1. Hydrodynamics (H) and linear response theory (A) are well-known techniques to study collective states of matter.Our paper introduces in this context a new perspective based on singulants (χ).The main text discusses the meaning of different overlaps between these domains represented as disks.

FIG. 2 .
FIG. 2. Singulants in the BRSSS model for a nonlinear longitudinal flow.Top panel: root-test for a numerical evaluation of the series Π (n) (black), and the slope predicted by the dominant singulant at this point χ blue (blue).Middle panel: fitted slope of the root-test plot evolving along a fluid flow line (black), with singulant predictions (red and blue).The crossing corresponds to an exchange of singulant dominance.The dashed line marks t used in the other panels.Bottom panel: singularities in the Borel plane ζ, indicated by Padé poles (black), illustrating agreement with singulant values (red and blue circles).

FIG. 3 .
FIG.3.Upper panel: large-order behaviour of the gradient expansion at x = 0, t = 0.5 for initial data of the form(28) as quantified by a root-test plot.The factorial growth is manifest.Bottom panel: singularities of the Padé approximant of the Borel transform of the gradient expansion at x = 0, t = 0.5.The points at which a line of pole accumulation starts are singulants and have been highlighted as stars.Cyan, orange, red, green, purple and brown stars correspond to candidate singulant pairs χ1, . . ., χ6 of progressively increasing norm.

FIG. 4 .
FIG. 4. Comparison between the time evolution of the three singulants of smallest norm as determined numerically (open circles) and the prediction (30) (solid lines).The matching has been performed at t = 0.5.The colour coding of the lines coincides with the colour coding of the stars in the bottom panel of Fig. 3. Upper panel: real parts.Lower panel: imaginary parts.

c 2 |Ω| 2 T 2 FIG. 5 .
FIG. 5. Numerical confirmation of the singulant equation of motion Eq. (35) along a flow line in the model Eq.(32).This model is distinguished from BRSSS and HJSW by the appearance of transverse derivatives in the singulant equation (controlled by cL), a feature it shares with holography and our main motivation for studying it.The values of Z(χ) 2 , U (χ), U (χ) 2 are extracted from the large-order behaviour of Π(n) and its derivatives, as described in the text.Here the black disks are rendered partially transparent to convey their density.

FIG. 6 .
FIG.6.Singulants in the gradient expansion in RTA kinetic theory in Bjorken flow.Top panel: Borel plane at τ ≈ 0.6 for three values of the momentum (v, pT ).One singulant seems to be approximately momentum independent, and is the dominant one.Another depends on their ratio.Bottom panel: test of the singulant equation of motion.χ can be measured numerically at each point in time by calculating the gradient expansion.Alternatively, we can take the dominant singulant from the top panel and evolve it by Eq. (94).The figure shows that these procedures agree.

FIG. 7 .
FIG.7.Borel plane for the distribution function in RTA kinetic theory in Bjorken flow, for the 1/w expansion.According to Eq. (102), the value of the singulant χ(p1, p2) can be predicted along a certain trajectory in the space of (p1, p2).We verify this by calculating the singulants for two points on such a curve.These singulants are shown in the Borel plane and the two sets of parameters correspond to the blue and green points in the figure.From the singulant revealed by one set, we predict the location for the other parameters and confirm this prediction using the explicit large-order computation.

FIG. 8 .
FIG.8.Upper panel: evolution of the norm of the dominant singulant |χ d | as a function of ξ along the flow line located at x = 0 for initial data (109).A real singulant, whose trajectory is denoted by the solid blue line, dominates at early times.At late times, a complex-conjugated singulant pair, whose trajectory corresponds to the solid red line, takes over.The discontinuous solid black line corresponds to nopt,est, as given by Eq. (104).Middle panel: coefficients of the gradient expansion at three selected times.The ones associated to the corresponding nopt,est have been highlighted in red.Lower panel: absolute error of the optimal truncation (black crosses) and the truncation at order nopt,est (red dots).Crucially, the latter error provides an upper bound for the former that displays the same time dependence.The dashed (dotted) black line represents the absolute error of the first-order (secondorder) order truncation.

kFIG. 9 .
FIG. 9. Real (upper panel) and imaginary (bottom panel) parts of γs, for ω ∈ [10, 10] and k ∈ [0, 10].ω and k are measured in units of πT = 1.These numerical results have been obtained from a collocation grid of 60 points, with ω and k forming a square lattice of spacing 0.125.

3 γs poles at k = 0 FIG. 10 .
FIG. 10.Upper panel: log(|γs|) as a function of ω for k = 0.The γs has been obtained using a collocation grid of 60 points, with the real and imaginary parts of the frequency forming a square lattice of spacing 0.125.Bottom panel: trajectories described by Ω (±) 1

FIG. 12 .
FIG.12.Real (upper panel) and imaginary (bottom panel) parts of the shear channel modes in the new MIS-like model in the local rest frame of the fluid as a function of |k|.We have set cL = 1, cσ = 0, η s = 1 4π , ΩI = 2 and ΩR = 1, and worked in T0 = 1 units.The two nonhydrodynamic modes are associated to red curves, while the hydrodynamic mode is associated to black ones.

FIG. 13 .ss
FIG.13.Real (upper panel) and imaginary (bottom panel) parts of the sound channel modes in the new MIS-like model in the local rest frame of the fluid as a function of |k|.We have set cL = 1, cσ = 0, η s = 1 4π , ΩI = 2 and ΩR = 1, and worked in T0 = 1 units.The two nonhydrodynamic modes are associated to red curves, while the two hydrodynamic modes are associated to black ones.
and that |χ d | 1.In this situation, it is immediate to demonstrate that Π