Hydrodynamic Gradient Expansion Diverges beyond Bjorken Flow

The gradient expansion is the fundamental organising principle underlying relativistic hydrodynamics, yet understanding its convergence properties for general nonlinear flows has posed a major challenge. We introduce a simple method to address this question in a class of fluids modelled by Israel-Stewart--type relaxation equations. We apply it to (1+1)-dimensional flows and provide numerical evidence for factorially divergent gradient expansions. This generalises results previously only obtained for (0+1)-dimensional comoving flows, notably Bjorken flow. We also demonstrate that the only known nontrivial case of a convergent hydrodynamic gradient expansion at the nonlinear level relies on Bjorken flow symmetries and becomes factorially divergent as soon as these are relaxed. Finally, we show that factorial divergence can be removed using a momentum space cutoff, which generalises a result obtained earlier in the context of linear response.

Introduction-Hydrodynamics plays a pivotal role in the description of nonequilibrium phenomena, with applications ranging from condensed matter systems [1] to scenarios in astrophysics [2][3][4] or nuclear physics [5,6]. The reason is that hydrodynamics captures the infrared behavior of any medium endowed with conserved quantities. For a given set of conserved currents, the expression of hydrodynamic behaviour rests on the derivative expansion in the spirit of an effective field theory [7][8][9][10][11]. For a neutral relativistic fluid, the natural choice of dynamical variables are the energy density E(x) and the unit-normalized fluid velocity U (x) = U µ (x)∂ µ , with the conserved currents, T µν , given by the constitutive relation T µν = E U µ U ν + P (E)(g µν + U µ U ν ) + Π µν . (1) Here, the first two terms describe ideal flow with g being the Minkowski metric and Π µν captures dissipative effects organised as where Π (n) µν contains n spacetime derivatives of E, U and we have introduced as a formal derivative-counting parameter. The gradient expansion in Eq. (2) is defined up to redundancies associated with frame choice and current conservation ∇ µ T µν = 0.
Understanding the character of the expansion (2) constitutes a fundamental open problem. Is it convergent, in * michal.p.heller@aei.mpg.de † alexandre.serantes@ub.edu ‡ michal.spalinski@ncbj.gov.pl § viktor.svensson@aei.mpg.de ¶ b.s.withers@soton.ac.uk such a way that subsequent truncations are progressively more accurate? If not, how does its divergent nature relate to the empirical success of low-order truncations? Studies of comoving flows in Refs. [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31], in which all fluid flow lines can be mapped to each other under symmetry transformations, rendering the problem effectively (0+1)-dimensional, have been instrumental in advancing our understanding of the hydrodynamic expansion (2). Among these, Bjorken flow [32] in conformally invariant theories is the most thoroughly explored example due to its role in studies of quark-gluon plasma. In these cases a particular strategy to solve the dynamical equations is an expansion in the Knudsen number, 1/w [33]. It is possible to compute a sufficient number of terms to assess that these expansions are factorially divergent. The expansion in defined in Eq. (2) encapsulates the expansion in 1/w, as we review in the Supplemental Material. Another well-studied example of a comoving flow is the Gubser flow which is reached by Weyl transformation from a (0 + 1)-flow on dS 3 × R [34,35].
Outside the realm of comoving flows, the only generic result on (2) was restricted to the linear response regime [36]. It showed that depending on how the momentum space support k max of E and U µ compares to an intrinsic scale of the underlying microscopic theory k * [37], the gradientexpanded constitutive relations could either be convergent (k max < k * ), geometrically divergent (k * < k max < ∞) or factorially divergent (k max → ∞).
In this Letter, we break free both from the symmetry constraints of comoving flows and the conveniences of linearisation to address, for the first time, the largeorder behavior of the hydrodynamic gradient expansion beyond comoving flows at the fully nonlinear level. Specifically, we ask the following question: Given a generic nonequilibrium, nonlinear configuration of E, U µ arising on-shell, what is nature of the expansion in (2) when evaluated on this solution? We answer this question by introducing a simple method that allows one to calculate Eq. (2) up to high order on a desktop computer. In this Letter we illustrate it in two examples. The first is the model put forward by Baier, Romatchske, Son, Starinets and Stephanov (BRSSS) in Ref. [38], while the second is the model originally introduced by Denicol and Noronha (DN) in Ref. [27]. Both theories are representative examples within a broad class of models employing the Israel-Stewart approach to embed hydrodynamics in a framework compatible with relativistic causality [39,40]. Our method applies to any member of this class, such as [41], and covers also equations with more than one derivative of Π µν , such as [42][43][44]. Each member of this class generally gives rise to infinitely many transport coefficients in the gradient expansion (2) that are specific to it.
Results in the BRSSS model-In this Letter, we consider the restriction of the BRSSS model to conformal fluids in d spacetime dimensions (d = 4 in our numerics). We work in the Landau frame and take Π µν U ν = 0. The model is defined by promoting Π µν to a set of independent degrees of freedom subject to a relaxation equation, where we neglected terms not relevant to the flow we consider. D µ is the Weyl-covariant derivative [45], and the angle brackets instruct one to take the symmetric, transverse and traceless part of the tensor they act upon. The relaxation time τ Π , the shear viscosity η, and λ 1 are transport coefficients. Conformal invariance demands that these quantities depend on the local temperature T , defined by the relation E = E 0 T d , as where C η and C τΠ > 0. The equations of motion of BRSSS theory are given by (3) and the conservation equation ∇ µ T µν = 0, where the energy-momentum tensor is specified in terms of E, U and Π µν as in (1). The fluid flows we focus on are characterised as follows. We separate the spatial coordinates into one longitudinal direction, x, and d−2 transverse directions, x 1 , . . . , x d−2 , demanding isotropy and translational invariance in the transverse hyperplane spanned by x i . Hence, the nontrivial dynamics is confined to the longitudinal plane spanned by t and x, and our fluid flows are (1+1)-dimensional. We refer to these fluid flows as longitudinal. At the linearized level such flows would correspond to sound wave propagation. See, e.g., Ref. [46] for a study of longitudinal flows in a quark-gluon plasma context. [47] The most general fluid velocity for a longitudinal flow is parameterised by a single degree of freedom, u, as Furthermore, any tensor which is symmetric, transverse to U µ and traceless is described in terms of a single additional degree of freedom that we pick as We now consider the expansion (2) applied to the conformal BRSSS model for longitudinal flows. This is facilitated by a numerical algorithm which makes a computation of (2) to large orders tractable. Since counts derivatives it can be introduced by taking (3) and replacing ∇ µ → ∇ µ together with positing a perturbative ansatz for Π , as follows, This leads to the following recursion relation Here, E and U are not expanded in . Therefore to proceed to evaluate (8) we must first find E and U for a given choice of flow. These (as well as the exact Π ) can be obtained by numerically solving the BRSSS equations of motion as an initial value problem without invoking an expansion. Once E and U are known, the recursion relation (8) can be efficiently evaluated numerically to high orders. Careful consideration of resolution and precision is required, as this procedure involves high numbers of successive derivatives of the background solution E and U . This is further discussed in the Supplemental Material, where we show that our numerical results are convergent.
Our approach applies to the whole class of theories which build on the Israel-Stewart approach. Note that solving the recursion relation (8) analytically is prohibitively expensive due to the fast growth in the number of terms contributing at each order. In particular, for the case λ 1 = 0 we observed an exponential growth of the number of individual contributions at each order. Our method circumvents this difficulty.
We have applied our approach in the BRSSS model across a wide variety of initial conditions, transport coefficients and spacetime points. We find factorial growth in all cases considered. We illustrate this in Fig. 1 with one representative example in which we consider a strong Gaussian overdensity for our initial conditions and adopt a periodic compactification of the spatial direction. demonstrates factorial growth at these sampled points. Further details are provided in the Supplemental Material.
Momentum cutoff-In previous work [36] we showed that a momentum-space cutoff gives at most a geometrically divergent hydrodynamic expansion for linear deviations from equilibrium. This result naturally extends to strongly nonlinear scenarios, as we now demonstrate. So far in this Letter we have used a numerical grid simply as a tool to approximate the continuum, but we now push beyond the continuum picture and re-evaluate the grid in a new role as a physical lattice which naturally enforces a momentum-space cutoff. In the BRSSS model, when λ 1 = 0 the recursion relation (8) can be written as where M = −τ Π (U · ∂) − d(∂·U ) d−1 is a differential operator independent of n, depending only on the background solution E, U . For a grid of dimensions N x ×N t , each Π (j) can be written as a N x N t -sized vector, and M accordingly as a N x N t × N x N t square matrix. Thus, on a lattice the expansion is ultimately only geometrically growing at a rate set by the largest eigenvalue of M, which scales with the inverse lattice spacing. In Fig. 2 this is demonstrated by utilising a deliberately low resolution lattice to allow for evaluating the hydrodynamic expansion to order n = 8000. It shows the transition from factorial growth where the continuum approximation holds, to the geometrically divergent asymptotic behaviour governed by the aforementioned eigenvalue. We have also verified numerically that this result holds at λ 1 = 0, with the definition of M as given.
Resolving the DN model tension-Bjorken flow is a boost-invariant longitudinal flow such that the dynamics depends on the proper time τ = √ t 2 − x 2 only. In Ref. [27] the authors analysed a Knudsen number expansion for an ultrarelativistic gas of hard spheres undergoing Bjorken flow. While in all the other models such expansions have been observed to be factorially divergent, in Ref. [27] the terms grow only geometrically, with convergence ensuing for a Knudsen number smaller than a critical value. Our objective is to re-analyse this physical scenario using the expansion in . For Bjorken flow we find analogous results, namely geometric growth; however, our method allows us to explore what happens when these symmetries are relaxed. Working on a lattice gives a window of factorial growthwhere it successfully approximates the continuum-before yielding to a geometrically growing hydrodynamic expansion regulated by the lattice. The asymptotic value reached is governed by the lattice spacing as discussed in the text. Black disks are those from Fig. 1 corresponding to the disk spacetime marker point. Red circles are the same simulation and spacetime point, but on a coarse numerical grid and evaluated to hydrodynamic order n = 8000. This plot serves also as an indication of a convergence of our approach.
We work in d = 4. As in the conformal BRSSS model, the energy-momentum tensor in the DN model is traceless and decomposed as in Eq. (1), with Π µν still obeying Eq. (3) with λ 1 = 0. Hence, the recursion relations giving Π (n) still take the form (8), again with zero λ 1 . The differences start with the inclusion of a conserved current J µ = ρ U µ , where ρ is the particle density. Furthermore, τ Π and η are not fixed purely in terms of the local temperature T , but rather obey where σ T is the total cross section and a, b are positive dimensionless constants. For Bjorken flow, the conservation of the particle current J µ entails that the particle density ρ decouples from the energy-momentum tensor. One has that where ρ 0 = ρ(τ 0 ) is the initial particle density. Hence, where Kn = 1/(ρ 0 τ 0 σ T ) is the Knudsen number. In the DN model for Bjorken flow it is time independent.
To assess the large-order behavior of the expansion in , Eq. (2), we first note that one can find a closed-form expression for Π (n) , a fact that relies crucially on Eq. (12). Second, we recall that T can also be determined exactly [27] T (τ ) = T 0, where T 0,± depend on initial conditions and α ± on a, b, and Kn. Together, Eqs. (13) and (14) entail that Π (n) cannot grow factorially with n at fixed τ , since the repeated action of the differential operator τ ∂ τ on terms of the form (τ /τ 0 ) 1 3 −α± only grows geometrically. Furthermore, Eqs. (13) and (14) also imply that Π (1) (and therefore all Π (n>1) ) is a linear combination of eigenfunctions of the differential operator M defined as in Eq. (9), providing another perspective on why the gradient expansion grows geometrically in this case. We refer the reader to the Supplemental Material for further details.
The analysis above relies crucially on the symmetry restrictions of Bjorken flow. Empirically, when relaxing these symmetry restrictions in all cases studied we find that the large-order geometric growth is destroyed and the factorial divergence is restored. We illustrate this in Fig. 3 for a longitudinal flow corresponding to a small perturbation of Bjorken flow.
Summary and outlook-Understanding the behaviour of the hydrodynamic gradient expansion at large orders is a challenging question in the foundations of relativistic hydrodynamics. We have proposed a method to compute such series in a large class of models. Applying this to nonlinear longitudinal flows reveals factorially divergent series which we have illustrated with a number of examples. This shows that previously observed instances of factorial growth were not reliant on Bjorken symmetries.
It is natural to ask what are the generic conditions that lead to factorial growth. We have here established that the ability of a system to support arbitrarily large momentum is important; ways around this include working on a lattice, and appropriate restriction of initial data in the linearised case [36], both of which naturally lead to geometric growth. Second, as our analysis of the DN model shows, imposing special symmetries such as boost invariance can also lead to geometric rather than factorial growth. It would be interesting to explore this question further including other natural momentum cutoffs such as microscopic physics and turbulent cascades. eigenfunctions of M such that the hydrodynamic gradient expansion grows geometrically. This could form the basis of a rigorous mathematical formulation for investigating the genericity of factorial growth.
The picture that is emerging from this work and results in linear response [36] is that the origin of factorial growth at large n in the hydrodynamic gradient expansion is the successive action of n derivatives on the hydrodynamic variables E, U . This is intimately connected with support of a solution in momentum space. It suggests that having a factorial growth in the number of transport coefficients at each order is not necessary.
The factorial divergence of asymptotic series is not an impediment to their practical utility: such series typically provide excellent approximations as long as one does not exceed the so-called order of optimal truncation. Our work makes such investigations possible for a much wider set of flows than previously tractable.

Supplemental Material
Appendix A: The expansions in and 1/w in Bjorken flow In the studies of conformal Bjorken flow in the past decade, an important role was played by the following dimensionless clock variable [33] w ≡ τ T (τ ), whereas the conservation equation always (in all theories) reads Note that it is convenient to replace Π by a dimensionless quantity by taking its ratio with the energy density E. This is simply related to the pressure anisotropy A measuring deviations from local thermal equilibrium in the following way The 1/w-expansion has been primarily discussed for the this quantity and up to second order it takes the form Using the expansion in of Π at low orders we obtain the following expression for A where T = T (τ ). One can now use the conservation equation to replace derivatives of T (τ ) at each order in terms of powers of w, which directly leads to (A5). This shows that knowing Π in the expansion in to order n allows one to generate the expansion of A in 1/w up to and including terms O(1/w n ). Note however, that the n-th order of the expansion in 1/w contains in principle contributions from the orders 2 to n of the expansion in .
A key difference between the two expansions is that the functional expansion of Π in derivatives of E, U counted by as in Eq. (2) does not explicitly require invoking conservation laws (beyond the stated redundancies), whereas the 1/w-expansion necessarily utilizes the conservation equation. Finally, note again that the discussion above applies to conformal Bjorken flow, while in Appendix C we discuss a particular non-conformal version of Bjorken flow.

Appendix B: Further details on numerics and convergence
For the BRSSS simulations shown in Fig. 1 we use initial on a unit spatial circle (x ∼ x + 1). f was chosen to locally reproduce a Gaussian of width γ near x = 1/2, respect the spatial periodicity, and the additive constant chosen so that f has no homogeneous component in a cosine expansion. We used E = T 4 , C τΠ = 1/4, C η = 1/(4π), λ 1 = 0, γ 2 = 1/60. The initial value problem is solved using RK4 on a uniform spatial grid with N x = 1600 points for a total of N t = 3200 timesteps to reach the arbitrary choice of end time t = 1, working with 300 digits of precision. The differential operator for ∂ x used periodic fourth-order finite differences, and second-order finite differences for ∂ t one-sided at t = 0 and t = 1. Due to the large numbers of successive applications of derivative operators used in evaluating Π (n) we have placed particular importance on testing convergence of our results, both in the resolution of the numerical grid N x , N t and in the number of digits of precision used. Tests confirming convergence of |Π (n) | 1 n are shown in Fig. 4. For a given resolution and precision, when the window of factorial growth ends it can be extended by increasing one, the other, or both.
The results in the DN model had been obtained working in a curvilinear coordinate system in which the Minkowski metric reads ds 2 = −e 2α(τ,x) dτ 2 + e 2β(τ,x) dx 2 + d x 2 ⊥ , obeyed by Π , the evolution equations include one extra equation enforcing that the metric (B2) is flat. When performing numerical simulations, we worked with a compactified spatial coordinate ζ ∈ [−1, 1] instead of x. Both coordinates are related as x = γ −1 tanh −1 (ζ), with γ ∈ R + . Spatial derivatives were discretized with fourth-order centered finite-difference stencils, while the time evolution was performed with an explicit RK4 method. For the numerical results displayed in Fig. 3, we employed a spatial grid of spacing dζ = 1/1800, a time step dτ = 0.5dζ, and worked with precision 300. The time derivatives appearing in the recursion relation were computed with fourth-order finite-difference stencils.