Holographic constraints on Bjorken hydrodynamics at finite coupling

In large-$N_c$ conformal field theories with classical holographic duals, inverse coupling constant corrections are obtained by considering higher-derivative terms in the corresponding gravity theory. In this work, we use type IIB supergravity and bottom-up Gauss-Bonnet gravity to study the dynamics of boost-invariant Bjorken hydrodynamics at finite coupling. We analyze the time-dependent decay properties of non-local observables (scalar two-point functions and Wilson loops) probing the different models of Bjorken flow and show that they can be expressed generically in terms of a few field theory parameters. In addition, our computations provide an analytically quantifiable probe of the coupling-dependent validity of hydrodynamics at early times in a simple model of heavy-ion collisions, which is an observable closely analogous to the hydrodynamization time of a quark-gluon plasma. We find that to third order in the hydrodynamic expansion, the convergence of hydrodynamics is improved and that generically, as expected from field theory considerations and recent holographic results, the applicability of hydrodynamics is delayed as the field theory coupling decreases.


Introduction
Hydrodynamics is an effective theory [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15] of collective long-range excitations in liquids, gases and plasmas.Its applicability across energy scales has made it a popular and fruitful field of research for over a century.A particularly powerful aspect of hydrodynamics is the fact that it provides a good effective description over a vast range of coupling constant strengths of the underlying microscopic constituents.This is true so long as the mean-free-time between microscopic collisions t mft is smaller than the typical time scale (of observations) over which hydrodynamics is applicable, t mft t hyd .At weak coupling, the underlying microscopic dynamics can be described in terms of kinetic theory [16][17][18][19][20][21][22][23][24], which relies on the concept of quasiparticles.On the other hand, at very strong coupling, the applicability of hydrodynamics to the infrared (IR) dynamics of various systems without quasiparticles has been firmly established much more recently through the advent of gauge-gravity duality (holography) [25][26][27][28].In infinitely strongly coupled CFTs with a simple holographic dual, the mean-free-time is set by the Hawking temperature of the dual black hole, t mft ∼ /k B T . 1 In a CFT in which temperature is the only energy scale, this implies that hydrodynamics universally applies to the IR regime of strongly coupled systems for ω/T 1, where the frequency scales as ω ∼ 1/t hyd (and similarly for momenta, q/T 1).A natural question that then emerges is as follows: how does the range of applicability of hydrodynamics depend on the coupling strength of the underlying microscopic quantum field theory?Qualitatively, using simple perturbative kinetic theory arguments (see e.g. a recent work by Romatschke [29] or Ref. [30]), one expects the reliability of hydrodynamics to decrease (at some fixed ω/T and q/T ) with decreasing coupling constant λ.The reason is that, typically, the mean-free-time increases with decreasing λ.From the strongly coupled, non-perturbative side, the same picture recently emerged in holographic studies of (inverse) coupling constant corrections to infinitely strongly coupled systems in [31][32][33][34],2 which we will further investigate in this work.
In holography, in the limit of infinite number of colors N c of the dual gauge theory, inverse 't Hooft coupling constant corrections correspond to higher derivative gravity α corrections to the classical bulk supergravity.In maximally supersymmetric N = 4 Yang-Mills (SYM) theory, dual to the IR limit of ten-dimensional type IIB string theory, the leading-order corrections to the gravitational sector (including the five-form flux and the dilaton), are given by the action [37][38][39][40][41] compactified on S 5 , where γ = α 3 ζ(3)/8, κ 10 ∼ 1/N c and the term W is proportional to fourth-power (eight derivatives of the metric) contractions of the Weyl tensor 2) The 't Hooft coupling of the dual N = 4 CFT is related to γ by the following expression: , where L is the anti-de Sitter (AdS) length scale.For this reason, perturbative corrections in γ ∼ α 3 are dual to perturbative corrections in 1/λ 3/2 .Another family of theories, which have been proven to be a useful laboratory for the studies of coupling constant dependence in holography, are curvature-squared theories [31-34, 42, 43] with the action given by Although the dual(s) of (1.3) are generically unknown, 3 one can treat curvature-squared theories as invaluable bottom-up constructions for investigations of coupling constant corrections on dual observables of hypothetical CFTs. 4 From this point of view, it is natural to interpret the α n coefficients as proportional to α .Since the action (1.3) results in higher-derivative equations of motion, the α n need to be treated perturbatively, i.e. on the same footing as the γ ∼ α 3 corrections in N = 4 SYM.The latter restriction can be lifted if one instead considers a curvature-squared action with the α n coefficients chosen such that α 1 = −4α 2 = α 3 .The resulting theory, known as the Gauss-Bonnet theory results in second-derivative equations of motions, therefore enabling one to treat the Gauss-Bonnet coupling, λ GB ∈ (−∞, 1/4], at least formally, non-perturbatively. 5Even though this theory is known to suffer from various UV causality problems and instabilities [47][48][49][50][51][52][53][54][55][56][57][58][59][60][61][62][63][64], one may still treat Eq. (1.4) as an effective theory which can, for sufficiently low energy and momentum, provide a well-behaved window into non-perturbative coupling constant corrections to the low-energy part of the spectrum.This point of view was advocated and investigated in [31,34,42,43] where it was found that a variety of weakly coupled properties of field theories, including the emergence of quasiparticles, were successfully recovered not only from the type IIB supergravity action (1.1) but also from the Gauss-Bonnet theory (1.4). 6An important fact to note is that these weakly coupled predictions follow from the theory with a negative λ GB coupling (increasing |λ GB |).
We can now return to the question of how coupling dependence influences the validity of hydrodynamics as a description of IR dynamics by using the above two classes of top-down and bottom-up higher derivative theories.The first concrete holographic demonstration of the failure of hydrodynamics at reduced (intermediate) coupling was presented in [31].The same qualitative behaviour was observed in both N = 4 and (non-perturbative) Gauss-Bonnet theory.Namely, as one increases the size of higher derivative gravitational couplings (decreases the coupling in a dual CFT), there is an inflow of new (quasinormal) modes along the negative imaginary ω axis from −i∞.Note that at infinite 't Hooft coupling λ, these modes are not present in the quasinormal spectrum.However, as λ decreases, the leading new mode on the imaginary ω axis monotonically approaches the regime of small ω/T .In the shear channel, 7 which contains the diffusive hydrodynamic mode, the new mode collides with the hydrodynamic mode after which point both modes acquire real parts in their dispersion relations.Before the modes collide, to leading order in q, the diffusive and the new mode have dispersion relations [31,34] ) where the imaginary gap ω g , the shear viscosity η and energy density ε, and pressure P depend on the details of the theory [31,34].Note also that both the IIB coupling γ and the Gauss-Bonnet coupling −λ GB have to be taken sufficiently large in order for this effect to be well described by the small-q expansion (see Ref. [34]).In the sound channel, ) where c s = 1/ √ 3 is the conformal speed of sound and Γ = 2η/3 (ε + P ).In both channels, it is clear that the IR is no longer described by hydrodynamics.To quantify this, it is natural to define a critical coupling dependent momentum q c (λ) at which Im |ω 1 (q c )| = Im |ω 2 (q c )| in the shear channel, and Im |ω 1,2 (q c )| = Im |ω 3 (q c )| in the sound channel.With this definition, hydrodynamic modes dominate the IR spectrum for frequencies ω(q), so long as q < q c (λ).To leading order in the hydrodynamic approximation, in N = 4 theory, q c scales as q c ∼ 0.04 T /γ ∼ 0.28 λ 3/2 T , while in the Gauss-Bonnet theory, q c ∼ −3.14 T /λ GB .Even though these scalings are approximate, they nevertheless reveal what one expects from kinetic theory: the applicability of hydrodynamics is limited at weaker coupling by a coupling dependent scaling whereas at strong coupling, hydrodynamics is only limited to the region of small q/T , independent of λ 1.8 Understanding of hydrodynamics has been important for not only the description of everyday fluids and gases, but also a nuclear state of matter known as the quark-gluon plasma that is formed after collisions of heavy ions at RHIC and the LHC.Hydrodynamics becomes a good description of the plasma after a remarkably short hydrodynamization time t hyd ∼ 1 − 2 fm/c measured from the moment of the collision [66][67][68][69][70][71].In holography, heavy ion collisions have been successfully modelled by collisions of gravitational shock waves [72][73][74][75][76][77][78][79], including the correct order of magnitude result for the hydrodynamization time (at infinite coupling).Coupling constant corrections to holographic heavy ion collisions were studied in perturbative curvature-squared theories (Gauss-Bonnet) in [32], which found that for narrow and wide gravitational shocks, respectively, the hydrodynamization time is where T hyd is the temperature of the plasma at the time of hydrodynamization.For λ GB = −0.2,which corresponds to an 80% increase in the ratio of shear viscosity to entropy density, we thus find a 25% and 290% increase in the hydrodynamization time [32].Thus, t hyd was found to increase for negative values of λ GB , which is consistent with expectations of the behavior of hydrodynamization at decreased field theory coupling.Consistent with these findings, the investigation of [33,80] further revealed that for negative λ GB , the isotropization time of a plasma also increased, again reproducing the expected trend of transitioning from infinite to intermediate coupling.
In this paper, we continue the investigation of coupling constant dependent physics by studying the simplest hydrodynamic model of heavy ions-the boost-invariant Bjorken flow [81]-in higher derivative bulk theories of gravity.The Bjorken flow has widely been used to study the evolution of a plasma (in the mid-rapidity regime) after the collision.While the velocity profile of the solution is completely fixed by symmetries, relativistic Navier-Stokes equations need to be used to find the energy density, which is expressed as a series in inverse powers of the proper time τ .The details of the solution will be described in Section 2.
In N = 4 SYM at infinite coupling, the energy density of the Bjorken flow to third order in the hydrodynamic expansion (ideal hydrodynamics and three orders of gradient corrections) takes the following form [82][83][84][85][86][87][88]: where w is a dimensionful constant. 9Physically, the energy density of the Bjorken flow must be a positive and monotonically decreasing function of the proper time τ , capturing the late-time expansion and cooling of the fluid.For a conformal, boost-invariant system, the energy density (1.10) uniquely determines all the components of the stress-energy tensor.Energy conditions then imply that the solution becomes unphysical at sufficiently early times, when (1.10) is negative.For instance, by considering the first two terms in (1.10), it is clear that the solution becomes problematic at times τ < τ 1st hyd , where Physically, the reason is that for τ < τ hyd , the first viscous correction becomes large and the hydrodynamic expansion breaks down, making the Bjorken flow unphysical. 10Ref. [90] further analyzed the evolution of non-local observables in a boost-invariant Bjorken plasma, finding stronger constraints on the value of initial τ for the Bjorken solution.For instance, equal-time two-point functions and space-like Wilson loops are expected to relax at late times as for some f and g such that f (τ w 3/2 ) → 0 and g(τ w 3/2 ) → 0 as τ → ∞.In the hydrodynamic regime, both f and g must be positive and monotonically decreasing functions of τ , implying that, as the plasma cools down, these non-local observables relax smoothly from above to the corresponding vacuum values.Such exponential decays have indeed been observed from the full numerical evolution in shock wave collisions [91,92].The interesting point here is that, if we were to truncate the hydrodynamic expansion to include only the first few viscous corrections, then f and g may become negative or non-monotonic at some τ crit > τ hyd , imposing further constraints on the regime of validity of hydrodynamics.In [90], it was found that a much stronger constraint (approximately 15 times stronger than (1.11)) for first-order hydrodynamics comes from the longitudinal two-point function: while for Wilson loops, the constraint was weaker: In addition, Ref. [90] also studied the evolution of entanglement (or von Neumann) entropy in a Bjorken flow, but found that the bound obtained in that case was equal to τ 1st hyd given by Eq. (1.11), i.e. weaker than the two constraints above.The reason for this equality is that in the late-time and slow-varying limit considered for the computation, the entanglement entropy satisfies the so-called first law of entanglement, where V A is the volume of the subsystem and T A is a constant that depends on its shape.Such a law holds for arbitrary time-dependent excited states provided the evolution of the system is adiabatic with respect to a reference state [93].
In this paper, we ask how higher-order hydrodynamic and coupling constant corrections affect the critical time τ crit after which the Bjorken flow yields physically sensible observables.In particular, we extend the analysis of [90] focusing on equal-time two-point functions and expectation values of Wilson loops.From the point of view of our discussion regarding viscous corrections and their role in keeping ε(τ ) positive, it seems clear that at decreased coupling, when the viscosity η becomes larger, the applicability of the Bjorken solution should become relevant at larger τ .Our calculations provide further details regarding the applicability of hydrodynamics.As a result, we will be computing an observable that is related to a coupling-dependent hydrodynamization time [32], but is analytically-tractable and therefore significantly simpler to analyze, albeit for realistic applications limited to the applicability of the Bjorken flow model.In this way, we obtain new holographic coupling-dependent estimates for the validity of hydrodynamics, analogous to the statement of Eq. (1.9), which allow us to compare top-down and bottom-up higher derivative corrections.
We will consider both the effects of higher-order (up to third order [94]) hydrodynamics and coupling constant corrections.Up to third order in the gradient expansion, we find no surprises as the Bjorken flow observables become well defined in higher-order hydrodynamics at earlier times.In other words, no effects of asymptotic expansion divergences [95] are found to third order.As for coupling dependence, what we find is that the most stringent constraints arise from the calculations of a longitudinal equal-time two-point function, i.e. with spatial insertions along the boost-invariant flow direction.For the two higherderivative theories, to first order in the coupling and to second order in the hydrodynamic expansion, ) where τ crit is the initial critical proper-time.At γ = 6.67 × 10 −3 (λ = 7.98, having set L = 1) and at λ GB = −0.2(each increasing η/s by 80%), we find that τ 2nd crit w 3/2 increases by 92.3% and by 150% in N = 4 and a linearized dual of Gauss-Bonnet theory, respectively (see Tables 1 and 2 for other numerical estimates).In a fully non-perturbative Gauss-Bonnet calculation, the increase is instead found to be 145%, which shows a rather quick convergence of the perturbative Gauss-Bonnet series for this observable to the full result at λ GB = −0.2(see also [32]).Thus, our results lie inside the interval of increased hydrodynamization time found in narrow and wide shocks obtained from non-linear shock wave simulations [32].
The paper is structured as follows: In Section 2, we discuss higher-order hydrodynamics and details of the hydrodynamic Bjorken flow solution, including all necessary holographic transport coefficients that enter into the solution.In Section 3, we discuss the construction of holographic dual geometries to Bjorken flow.We focus in particular on the case of the Gauss-Bonnet theory which, to our understanding, has not been considered in previous literature. 11In Section 4, we analyze the relaxation properties of two-point functions and Wilson loops, extracting the relevant critical times at which the hydrodynamic approximation breaks down.Finally, Section 5 is devoted the discussion of our results.

Hydrodynamics and Bjorken flow
We begin by expressing the equations that describe the boost-invariant evolution of chargeneutral, conformal relativistic fluids, which will be studied in this work.In the absence of any external sources, the equations of motion (relativistic Navier-Stokes equations) follow from the conservation of stress-energy The constitutive relations for the stress-energy tensor of a neutral, conformal (Weylcovariant) relativistic fluid can be written as (see e.g [97]) where we have chosen to work in the Landau frame.The transverse projector ∆ ab is defined as ∆ ab ≡ g ab + u a u b , with u a being the velocity field of the fluid flow.In four spacetime dimensions, the pressure P and energy density ε are related by the conformal relation P = ε/3.The transverse, symmetric and traceless tensor Π ab can be expanded in a gradient expansion (in gradients of u a and a scalar temperature field).To third order in derivatives [94,98,99], where we have used the longitudinal derivative D ≡ u a ∇ a and a short-hand notation which ensures that any tensor A ab is by construction transverse, u a A ab = 0, symmetric and traceless, g ab A ab = 0.The tensor σ ab is a one-derivative shear tensor The vorticity Ω µν is defined as the anti-symmetric, transverse and traceless tensor The transport coefficients appearing in (2.3) are the shear viscosity η, 5 second order coefficients ητ Π , κ, λ 1 , λ 2 , λ 3 , and 20 (subject to potential entropy constraints) conformal third order transport coefficients λ i , which multiply 20 linearly independent, third order Weyl-covariant tensors O ab i that can be found in [94].The boost-invariant Bjorken flow [81] is a solution to the hydrodynamic equations (Eq.(2.1)), and has been widely used as a simple model of relativistic heavy ion collisions (see [77]).Choosing the direction of the beam to be the z axis, the Bjorken flow is boost-invariant along z, as well as rotationally and translationally invariant in the plane perpendicular to z (denoted by x ⊥ ).By introducing the proper time τ = √ t 2 − z 2 and the rapidity parameter y = arctanh(z/t), the velocity field, which is completely fixed by symmetries, and the flat metric can be written as ) Note that the solution is also invariant under discrete reflections y → −y.What remains is for us to find the solution for the additional scalar degree of freedom that is required to fully characterize the flow.In this case, it is convenient to work with a proper time-dependent energy density ε(τ ) and write Eq. (2.1) as in [98]: (2.9) By using the conformal relation P = ε/3 and the fact that the only non-zero component of ∇ a u b is ∇ y u y = ∇ ⊥y u y = τ , Eq. (2.9) then gives with Π yy from Eq. (2.3) expanded as (3) 2

+ 3λ
(3) 8 (2.11) Each transport coefficient appearing in (2.11) can only be a function of the single scalar degree of freedom-the energy density-with dependence on ε determined uniquely by its conformal dimension under local Weyl transformations [94,98]: where C, η, τΠ and λ(3) n are constants.Finally, the Bjorken solution to Eq. (2.1) for the energy density, expanded in powers of τ , becomes with ν = 2/3.Terms at order O τ −2−3ν are controlled by the hydrodynamic expansion to fourth order, which is presently unknown.
In this work, we will not look beyond third-order hydrodynamics.What is important to note is that the gradient expansion is believed to be an asymptotic expansion, similar to perturbative expansions.As a result, the Bjorken expansion in proper time formally has a zero radius of convergence [95].In practice, this means that at some order, the expansion in inverse powers of τ breaks down and techniques of resurgence are required for analyzing long-distance transport (see e.g.[95,[107][108][109][110][111][112]).

Gravitational background in Gauss-Bonnet gravity
In this section, we begin our analysis of holographic duals to Bjorken flow.Throughout this paper, we will be interested in three separate cases: • Einstein gravity.Bjorken flow in N = 4 SYM at infinite coupling, expanded to third order in the hydrodynamic series.
• λ GB -corrections.Bjorken flow in a hypothetical dual of Gauss-Bonnet theory with λ GB coupling corrections, expanded to second order in the hydrodynamic series.
In the first case, the holographic dual geometry is well known (see Refs. [82][83][84][85][86][87][88]).What one finds is that in the near-boundary region, which is the only region relevant for computing the non-local observables studied in this paper (two-point correlators of operators with large dimensions and Wilson loops), the geometries are specified by symmetry and (relevant order) hydrodynamic transport coefficients. 12As we will see, the same conclusions can also be drawn in higher-derivative theories.As a check, we derive here the full geometric Bjorken background in non-perturbative Gauss-Bonnet theory.All details of the perturbative calculations in Type IIB supergravity with α corrections will be omitted, but we refer the reader to [96] for the explicit derivation.

Static background
Equations of motion for Gauss-Bonnet gravity in five dimensions can be derived from the action (1.4) and take the following form: where This set of differential equations admits a well-known (static) asymptotically AdS black brane solution: with the emblackening factor In the near-boundary limit, the asymptotically AdS region exhibits the following scaling: where η ab is the flat metric and the AdS curvature scale, L, is related to the length scale set by the cosmological constant, L, via The Hawking temperature, entropy density and energy density of the dual theory are then given by13 ) In what follows, we will set L = 1 unless otherwise stated.
To make the metric manifestly boost-invariant along the spatial coordinate z, we transform (3.2) by introducing a proper time coordinate τ = √ t 2 − z 2 .Next, we perform an additional coordinate transformation to write the metric in terms of ingoing Eddington-Finkelstein (EF + ) coordinates with which gives the metric It should be noted that the EF + time, τ + , mixes the proper time, τ , and r in the bulk.At the boundary, however, lim A static black brane with a constant temperature cannot be dual to an expanding Bjorken fluid, which has a temperature that decreases with the proper time, T fluid ∼ τ −1/3 .As in the fluid-gravity correspondence [99], where the black brane is boosted along spatial directions, here, one may make an informed guess and allow for the horizon to become time-dependent by substituting r h → wτ where w is a constant and τ + is the fluid's proper time at the boundary.The Hawking temperature is then and the static black brane metric (3.10) takes the form Of course, as in the fluid-gravity correspondence, Eq. (3.14) is not a solution to the Gauss-Bonnet equations of motion.As will be shown below, however, the background solution asymptotes to (3.14) at late times, i.e.Eq. (3.14) is (approximately) dual to Bjorken flow in the regime dominated by ideal hydrodynamics.

Bjorken flow geometry
The full (late-time) geometry is systematically constructed following the procedure outlined in Ref. [113] (see also [114]).In EF + coordinates, the most general metric respecting the symmetries of Bjorken flow is where a, b, c are functions of r and τ + and our boundary geometry is expressed in proper time-rapidity coordinates (see the discussion above Eq.(2.8)).
At late times, the equations of motion (3.1) can be solved order-by-order in powers of τ + fixed.To perform the late time expansion, we will change coordinates from {τ + , r} → {v, u}, where and assume the metric functions a, b and c can be expanded as We then solve the equations order-by-order in powers of u and impose Dirichlet boundary conditions (at the boundary) at every order: At a given order, i, the equations of motion form a system of second-order differential equations for a i , b i and c i along with two constraint equations.We therefore have six integration constants at each order.One integration constant is related to a residual diffeomorphism invariance of our metric under the coordinate transformation [113] r → r + f (τ + ) , (3.20) and can be freely specified without affecting the physics of our boundary field theorya feature that will be exploited to simplify the solutions.Three of the five remaining integration constants can fixed by requiring the bulk geometry to be free of singularities (apart from at v = 0) and imposing the asymptotic AdS boundary conditions above.In practice, to the order considered, we find that the integration constant which ensures bulk regularity can be set by requiring ∂ v c i to be regular at a particular value 14 of v.The remaining integration constants are specified by the two constraint equations.For i > 0, one of the constraint equations can specify a constant at order i, while the other specifies a constant at order i − 1.

Solutions
We now present the full zeroth-and first-order solutions in the late-time (hydrodynamic gradient) expansion.At second order, we were unable to find closed-form solutions analytically that would extend throughout the entire bulk.However, sufficiently complete solutions for the purposes of this work can be found non-perturbatively in λ GB near the boundary, or perturbatively in the full bulk.

Zeroth Order
At zeroth order in the hydrodynamic expansion (ideal fluid order), the equations of motion are solved by 15 One can see immediately that the zeroth-order solution is the boosted black brane metric given by Eq. (3.14).Near the boundary we find

First Order
At first (dissipative) order, our equations of motion are solved by 14 With the next section in mind, we require lim v→w + ∂vci < ∞. 15 We note that this is not the most general solution to the equations of motion at this order-there is an additional nonphysical integration constant corresponding to a gauge degree of freedom.A simple coordinate transformation [113] brings the solution into the form presented here.Similar remarks apply for the first-order solution. where (3.24) For simplicity, here we have presented c 1 in an integral representation.An explicit evaluation of the integral would result in an Appell hypergeometric function (see Ref. [34]). 16ear the boundary,

Second Order
As in Gauss-Bonnet fluid-gravity calculations [34], at second order in the hydrodynamic expansion, one is required to solve non-homogeneous differential equations with sources depending on complicated expressions involving Appell hypergeometric functions.For this reason, we were only able to find non-perturbative solutions (in λ GB ) near the boundary and solve the full equations perturbatively.
Near the boundary we find where A 2 and C 2 are, as yet, unspecified constants.To determine them, we would need to know the full bulk solutions and the constants would then follow from horizon regularity.Instead, as will be shown below, we will use known properties of the dual field theory (the transport coefficients and energy conservation) to show that they must take the following values: Full perturbative first-order (in λ GB ) solutions are presented in Appendix A. Here, we only state their near-boundary forms:

Stress-energy tensor and transport coefficients
We can now compute the boundary stress-energy tensor by following the well-known holographic procedure (see e.g.[34,115,116]), which we review here.First, we introduce a regularized boundary located at r = r 0 = const.The induced metric on the regularized boundary is given by γ µν ≡ g µν − n µ n ν , where n µ ≡ δ µ r / √ g rr is the outward-pointing unit vector normal to the r = r 0 hypersurface.The boundary stress-energy tensor is then where µν is the induced Einstein tensor on the regularized boundary, K µν is the extrinsic curvature17 and J = g µν J µν .The constants δ 1 and δ 2 , fixed by holographic renormalization, are given by For the background derived in Section 3.3, the non-zero components of the four dimensional boundary stress-energy tensor, T ab , are found to be where we identify τ + with the proper time, τ , at the boundary.
Before analyzing T ab , we note three immediate observations: 1. T ab is traceless: with η ab given by Eq. (2.8).
2. Conservation implies a relationship between A 2 and C 2 : 3. The stress-energy tensor is completely specified by a single time-dependent function, ε(τ ) ≡ T τ + τ + : The three properties above are the defining properties of the hydrodynamic description of a relativistic, conformal Bjorken fluid.The only thing that remains to be specified is a single integration constant A 2 (see discussion below Eq. (3.26)).Now, the energy density of a Bjorken fluid, given by Eq. (2.13), can be written to second order in the hydrodynamic gradient expansion as where Σ(2) represents the relevant linear combination of second-order transport coefficients: By comparing the energy density of the Gauss-Bonnet fluid derived in the previous section with that of the Bjorken fluid, we identify At zeroth order in the hydrodynamic expansion, the energy density of our plasma is, as required, where we have used Eq.(3.13) to express our answer in terms of T .The shear viscosity is then which agrees with Eq. (2.20).At second order, we find Collecting our results, the energy density, as a function of proper time, takes the final form:

Breakdown of non-local observables
In this section we study various non-local observables in the boost-invariant backgrounds described above.As advertised in the Introduction, we will see that requiring a physically sensible behavior for the observables leads to several constraints on the regime of validity of hydrodynamic gradient expansions at a given order.

Two-point functions
According to the holographic dictionary [117,118], bulk fields φ are dual to gauge-invariant operators O with conformal dimension ∆, specified by their spin s, the mass m and the number of dimensions d.For scalar fields, the relation is given by ∆(∆ − d) = m 2 .The equivalence between the two sides of the correspondence can be made more precise by the identification: The left-hand-side of the above equation is the bulk partition function, where we impose the boundary condition φ → d−∆ φ .The right-hand-side is the generating functional of correlation functions of the CFT, where the boundary value φ acts as a source of the dual operator O.The equivalence (4.1) becomes handy by treating the bulk path integral in the saddle point approximation.In this regime, the above relation becomes where on the left-hand side we have the bulk action evaluated on-shell and the righthand side is the generating functional of connected correlation functions of the CFT.For instance, two-point functions can be computed by differentiating two times with respect to the source: For operators with large conformal dimension ∆ (or equivalently, bulk fields with large mass m) the above problem simplifies even further.It can be shown that, in this limit, the relevant two-point functions reduce to the computation of geodesics in the given background geometry [119,120], i.e.
where S reg is the regularized length of a geodesic connecting the boundary points x and x .

Perturbative expansion: Eddington-Finkelstein vs. Fefferman-Graham
We can now compute the late-time behavior of scalar two-point functions probing the out-of-equilibrium Bjorken flow.In order to do so we will follow the approach of [90]. 18onsider the functional L[φ(y); α] for the geodesic length, i.e. S ≡ dy L[φ(y); α].Here, φ(y) denotes collectively all of the embedding functions, y is the affine parameter and α is a small parameter related to the hydrodynamic gradient expansion in which the perturbation is carried out.Its precise definition will be given below.We can expand both L and φ(y) as: L[φ(y); α] = L (0) [φ(y)] + αL (1) [φ(y)] + O(α 2 ) , φ(y) = φ (0) (y) + αφ (1) (y) + O(α 2 ) . (4.5) The functions φ (n) (y) can in principle be found by solving the geodesic equation order-byorder in α.However, the embedding equations are in most cases highly non-linear making closed form solutions difficult to find.The key point here is that at first order in α, so we only need φ (0) (y) to obtain the first correction to the geodesic length.
Let us now discuss the expansion parameter α in more detail.In particular, what we will see is that there is a natural choice for α depending on whether we work in Eddington-Finkelstein or Fefferman-Graham coordinates, so we must proceed with some care before we interpret our results. 19Let us start with the Fefferman-Graham expansion, which was first considered in [90].In this case, the metric coefficients can be expanded as in Eq. (B.8) so each hydrodynamic order is suppressed by a factor of the dimensionless quantity ũ = τ −2/3 w −1 , where w is the same dimensionful parameter that appears in the energy density.On the other hand, the near-boundary expansion stipulates that we can alternatively expand all metric coefficients in powers of ṽ ≡ zτ −1/3 w.This is the expansion that will be relevant for our perturbative calculation (4.6).Notice that when ṽ → 0, we recover pure AdS, for which the embedding function φ (0) (y) is analytically known.The first correction in this expansion enters at order O(ṽ 4 ) so we can identify α ∼ ṽ4 .Now, according to the UV/IR connection [122][123][124], the bulk coordinate z can roughly be mapped to the length scale z ∼ in the boundary theory.In our setup, the only length scale of the problem is given by the separation the two points (x, x ) so ∼ ∆x ≡ |x − x |. 20 Therefore, in terms of CFT data, our expansion parameter in Fefferman-Graham coordinates is given by α = 4 τ −4/3 w 4 (Fefferman-Graham).( As mentioned already in Appendix B, the leading correction to the metric in the nearboundary expansion receives contributions at all orders in hydrodynamics, so one can obtain non-trivial results by studying contributions to the two-point correlators to only first order.For instance, as found in Ref. [90], in order to have a well behaved late-time relaxation of longitudinal two-point functions, first-order hydrodynamics puts a constraint on the regime of validity of ũ.Namely, the approximation breaks down when21 In this work, we are interested in studying both i) higher order hydrodynamic corrections and ii) (inverse) coupling constant corrections in the N = 4 plasma and a hypothetical dual of Gauss-Bonnet theory.
In Eddington-Finkelstein coordinates, the hydrodynamic expansion is performed in terms of u, and the near boundary expansion in terms of v, both given in Eq. (B.3).However, notice that these definitions involve τ + instead of τ , which at the leading order becomes Eq.(B.6).If we perform a similar analysis in Fefferman-Graham coordinates, we find that in Eddington-Finkelstein coordinates the expansion parameter is given by α ∼ v −4 , or equivalently, α = 4 (τ − ) −4/3 w 4 (Eddington-Finkelstein). (4.9) Notice that in this case, truncating the expansion (4.6) at the leading order in α is problematic for τ < .Furthermore, if we expand (4.9) for τ , even the first subleading term is not complete since, due to the coordinate mixing, we would require higher order terms in the near-boundary expansion to have a full result at the given order in /τ .Thus, in Eddington-Finkelstein coordinates the results can only be trusted in the limit /τ → 0. 22  To avoid this issue we will convert first to Fefferman-Graham coordinates and perform our calculations in that chart. 23Explicit expressions for the metric functions are given in Appendix B.1.

Transverse correlator
In Fefferman-Graham coordinates, a generic bulk metric dual to Bjorken hydrodynamics can be written as follows: where {ã, b, c} are functions of (τ, z) that can be expanded in terms of ũ = τ −2/3 w −1 and ṽ = zτ −1/3 w 1 as in (B.10), i.e., ã(ṽ, ũ) = ã4 (ũ)ṽ 4 + . . ., and similarly for b and c.Notice that we have set the AdS radius to unity L = 1.The AdS radius generally depends on the cosmological constant Λ as well as all higher derivative couplings of the gravity theory that we consider.Since L is just an overall factor of our metric, it will only appear as an overall factor in the various observables we study, and can be easily restored via dimensional analysis.
Let us begin by considering space-like geodesics connecting two boundary points separated in the transverse plane: (τ 0 , x) and (τ 0 , x ), where x ≡ x 1 and all other spatial directions are identical.Because the metric (4.10) is invariant under translations in x, we can parameterize the geodesic by two functions τ (z) and x(z), satisfying the following UV boundary conditions:

.11)
22 For longitudinal correlators, this would imply that only the ∆y → 0 limit is valid.Fortunately, this is exactly the limit for which the constraint (4.8) was found. 23We explicitly checked that the results in both coordinate systems agree at the leading order in /τ .
At the end of the calculation, we can shift our coordinate x → x+x 0 , where x 0 = 1 2 (x+x ), and express the results in terms of ∆x = |x − x |, for any x and x .The length of such a geodesic is given by: (4.12) We can now use (B.10) and expand the above as: S = S (0) + S (1) + • • • , where The first term is just the pure AdS contribution, which is UV divergent.To see this, we can use the zeroth order embeddings: with z * = ∆x 2 .Integrating from → 0 to z * and subtracting the divergence S div = −2 ln , we obtain: which is the expected result for a two-point correlator in the vacuum of a CFT.At next order, the correlator can be written as follows: where S (1) is given in (4.13).The functions {ã 4 (ũ), b4 (ũ), c4 (ũ)} are generically theorydependent (see Appendix B.1 for explicit expressions) and contain information about all orders in hydrodynamics.On general grounds, we expect S (1) to be positive definite at late times, so the correlator relaxes from above as the plasma cools down.Below, we will use the explicit form of {ã 4 (ũ), b4 (ũ), c4 (ũ)} to put constraints on the regime of validity of hydrodynamics, at each order in the derivative expansion.
For the transverse correlator, there is a very drastic simplification: once we evaluate S (1) using the zeroth order embeddings (4.14), we have: where x = z/z * and ũ0 = τ −2/3 0 w −1 .Therefore, the positivity of S (1) follows directly from the positivity of c4 (ũ).Let us specialize to the particular cases of interest: Einstein gravity (which is dual to a Bjorken flow at infinite coupling), and higher derivative gravities with αand λ GB -corrections (two different models of Bjorken flow with finite coupling corrections).
• Einstein gravity.The function c4 (ũ) is known up to third order in hydrodynamics and is given by equation (B.12).Up to first order in hydrodynamics c4 (ũ) is positive definite but it becomes negative for τ < τ 2nd crit and τ < τ 3rd crit in second-and third-order hydrodynamics, respectively, where It is interesting to note that for this particular obsevable, the above criterion would naively imply that third-order hydrodynamics is more constraining than second-order hydrodynamics.However, as we will see below, the most stringent bound on the applicability of hydrodynamics will come from the longitudinal correlator, which decreases at each order in hydrodynamics (up to third order), as expected.
• α -corrections.The function c4 (ũ) is known to linear order in γ = α 3 ζ(3)/8 = λ −3/2 ζ(3)L 6 /8, and up to second order in hydrodynamics, and is given by equation (B.13).The coefficient c4 (ũ) is positive definite for first-order hydrodynamics, but becomes negative for τ < τ 2nd crit (γ) in second order hydrodynamics, where Finite coupling corrections (γ > 0) are shown to increase τ 2nd crit , which is in accordance with our expectations that they should reduce the regime of validity of hydrodynamics.As we will see below, the most stringent bound will again come from the longitudinal correlator.
• λ GB -corrections.The function c4 (ũ) is known non-perturbatively in λ GB and up to second order in hydrodynamics, and is given by equation (B.14).c4 (ũ) is positive definite for first-order hydrodynamics, but becomes negative for τ < τ 2nd crit (λ GB ) in second-order hydrodynamics, where Negative values of λ GB tend to increase τ 2nd crit so they reduce the regime of validity of hydrodynamics.This is indeed the expected behavior as we flow from strong to weak coupling.It is also interesting to study the full dependence of τ 2nd crit on λ GB ∈ (−∞, 1/4], which we plot in Figure 1.For negative λ GB , we observe that τ 2nd crit increases monotonically.However, for positive λ GB , τ 2nd crit is non-monotonic.We note that, also for this case, the true bound will come from the longitudinal correlator. Finally, it is worth noting that the results above can be expressed generically in terms of a few theory-specific constants { Σ, Σ(γ) , Σ(λ GB ) , Λ}, which can be found in Appendix C. At second order in the hydrodynamic expansion, the critical time is given by Expressing our coupling constants γ and λ GB collectively as β, first-order corrections to τ 2nd crit then take the form The expressions for τ 3rd crit are complicated, but correspond to the smallest real root of the equation 1 − Σξ 4/3 − 2 Λξ 2 = 0 , ( where ξ = τ −1 0 w −3/2 .

Longitudinal correlator
We are now interested in a space-like geodesic connecting two boundary points in the longitudinal plane: (τ 0 , y) and (τ 0 , y ) for any y and y .We can make use of the invariance under translations in y and parameterize the geodesic by functions τ (z) and y(z) with boundary conditions At the end, if desired, we can simply shift our rapidity coordinate y → y + y 0 , where y 0 = 1 2 (y+y ), and express our results in terms of x 3 = τ 0 sinh(y 0 + ∆y 2 ) and x 3 = τ 0 sinh(y 0 − ∆y 2 ).The length of such a geodesic is given by: S = 2 We can now use (B.10) and expand the above as: S = S (0) + S (1) + • • • , where Again, the first term gives the pure AdS contribution.To see this, we can use the zeroth order embeddings, which in this case are given by:

.27)
Integrating from → 0 up to z * = ∆x 3 2 = τ 0 sinh( ∆y 2 ) and subtracting the divergent part S div = −2 ln , we obtain: (4.28) At zeroth order, the longitudinal correlator depends only on |x 3 − x 3 |.This is expected because this is the result for a two-point correlator in the vacuum of a CFT, which is translationally invariant.At next order, the correlator can be written as follows: where S (1) is given in (4.26).Again, we expect S (1) to be positive definite at late times, so the correlator relaxes from above as the plasma cools down.However, we will see below that there are crucial differences with respect to the transverse case, which will ultimately lead to stricter bounds on the regime of validity of the hydrodynamic expansion.The next step is to evaluate S (1) using the zeroth-order embeddings (4.27) and then use the explicit forms of {ã 4 (ũ), b4 (ũ), c4 (ũ)} which are theory-dependent.Defining a dimensionless variable x = z/z * , we arrive at the following expression: where ũ(x) = 1 for some numbers {ã Different values of k correspond to contributions from different orders in hydrodynamics; for example, k = 0 corresponds to the perfect fluid approximation, k = 1 corresponds to first-order hydrodynamics, and so on.Therefore, we can rewrite S (1) as follows: where Both integrals can be performed analytically for any value of k, although we refrain from writing them out here, since they are not particularly illuminating.Nevertheless, it is interesting to study the ∆y → 0 limit, from which we can extract τ crit at different orders in hydrodynamics [90].A simple observation is that both of I (k) ± are positive definite and decrease monotonically as ∆y increases.In the limit ∆y → 0, both integrals are finite and independent of k: However, it is clear that the first term of (4.33) dominates since in this limit cosh( ∆y 2 ) → 1, while sinh( ∆y 2 ) → O(∆y).Putting everything together, we find that for ∆y → 0: where ũ0 = τ −2/3 0 w −1 .Therefore, in this limit the positivity of S (1) follows directly from the positivity of b4 (ũ).In the cases we considered, this criterion was enough to guarantee the positivity of S (1) for any other value of ∆y.However this does not trivially follow from (4.33): at finite ∆y, the value of S (1) will generally depend on the interplay between the coefficients {ã In the following, we will study in more detail the behavior of S (1) as a function of ∆y and τ 0 w 3/2 , specializing to the particular cases of interest: Einstein gravity and higher derivative gravities with α -and λ GB -corrections.
• Einstein gravity.The functions ã4 (ũ) and b4 (ũ) are known up to third order in hydrodynamics and are given in (B.12).With these functions at hand we can extract the numbers ã(k) 4 and b(k) 4 and then use formula (4.33).Figure 2 (left) shows some representative curves for S(1) ≡ S (1) 3 as a function of ∆y for various values of ξ = τ −1 0 w −3/2 = {0, 0.15, 0.3, 0.45, 0.6} depicted in blue, orange, green, red and purple, respectively.The solid lines correspond to third-order hydrodynamics; the dashed and dotted lines correspond to second-and first-order hydrodynamics, respectively.For ξ = 0.45 the dotted curve becomes negative for small ∆y, indicating that first-order hydrodynamics is no longer valid.For ξ = 0.6 both the dotted and dashed curves are negative for small ∆y.This indicates that second-order hydrodynamics is also invalid at this time.Finally, for all values of ξ that were plotted, the solid lines are always positive, so third-order hydrodynamics is valid for these values.However, if we keep on increasing ξ, the solid lines will become unphysical for small ∆y at some point.We observe the following behavior for any finite value of ξ (in the range of parameters that we plotted): the value of S(1) increases up to a maximum S(1) max > 0 and then decreases monotonically to zero as ∆y → ∞.This implies that the positivity of S(1) at ∆y = 0 is enough to guarantee a good physical behavior for any ∆y.In Figure 2 (right) we show the behavior of S(1) (0) as a function of ξ for first-, secondand third-order hydrodynamics, depicted in blue, orange and green, respectively, and we indicate the times at which it becomes negative.From the ∆y → 0 limit of the correlator (4.36) we obtain the critical times: These bounds are stricter than the ones derived from the transverse correlator (4.18), and decrease at each order in hydrodynamics, as expected.
• α -corrections.The functions ã4 (ũ) and b4 (ũ) are known to linear order in γ = for various values of ξ = τ −1 0 w −3/2 = {0, 0.15, 0.3, 0.45, 0.6} depicted in blue, orange, green, red and purple, respectively.The solid lines correspond to third-order hydrodynamics; the dashed and dotted lines correspond to second-and first-order hydrodynamics, respectively.Right: Plots for S(1) (0) for first-, second-and third-order hydrodynamics, depicted in blue, orange and green, respectively.The dashed vertical lines correspond to the critical times at each order in hydrodynamics.

and b(k)
4 and then use the formula (4.33).Figure 3 (left) shows some representative curves for S(1) ≡ S (1)  4  3 as a function of ∆y for various values of ξ = τ −1 0 w −3/2 = {0, 0.12, 0.25, 0.28, 0.5} depicted in blue, orange, green, red and purple, respectively.The solid lines correspond to γ = 0 (Einstein gravity) while the dashed lines correspond to γ = 10 −3 , both for second-order hydrodynamics.For all the ξ that were plotted the solid lines are well behaved because we have chosen ξ < ξ 2nd crit (γ = 0) = 0.503.For ξ = 0.5 the dashed curve becomes negative for small ∆y, indicating that second-order hydrodynamics becomes invalid faster at finite coupling.We observe the same behavior as in Einstein gravity, namely that the positivity of S(1) at ∆y = 0 is enough to guarantee a good physical behavior for any ∆y.In Figure 3 (right) we show the behavior of S(1) (0) both for γ = 0 and γ = 10 −3 as a function of ξ for first-and second-order hydrodynamics, depicted in blue and orange, respectively, and we indicates the times at which it becomes negative.From the ∆y → 0 limit of the correlator (4.36) we obtain the following critical times: These bounds increase as we increase the value of γ and are stricter than the ones derived from the transverse correlator (4.19).Based on this, we can conclude that finite coupling corrections indeed tend to reduce the regime of validity of hydrodynamics.
• λ GB -corrections.The functions ã4 (ũ) and b4 (ũ) are known non-perturvatively in λ GB and up to second order in hydrodynamics, and are given in (B.14  for various values of ξ = τ −1 0 w −3/2 = {0, 0.12, 0.25, 0.38, 0.5} depicted in blue, orange, green, red and purple, respectively.Solid lines correspond to γ = 0 (Einstein gravity) while the dashed lines correspond to γ = 10 −3 (α -corrections), in both cases for second-order hydrodynamics.Right: Plots for S(1) (0) for first-and second-order hydrodynamics, depicted in blue and orange, respectively.Solid lines correspond to γ = 0 while dashed lines correspond to γ = 10 −3 .The dashed vertical lines correspond to the critical times at each order in hydrodynamics, including the leading α -corrections.(4.33).For small and negative values of λ GB we observe qualitatively the same behavior as for the γ−corrections: the critical time below which first-and second-order hydrodynamics break down increases, which is the expected behavior for a theory that flows from strong to weak coupling.On the other hand, positive values of λ GB behave in the opposite way, and thus appear unphysical for λ GB interpreted as a coupling constant.From the ∆y → 0 limit of the correlator (4.36) we obtain the following critical times: It is interesting to consider the behavior of the correlator for negative values of λ GB in the non-perturbative regime.Figure 4 (left) shows S(1) ≡ S (1) τ 4/3 0 /w 4 ∆x 4 3 plotted as a function of ∆y for a few representative values of ξ = τ −1 0 w −3/2 = {0, 0.1, 0.2}, depicted in blue, orange and green, respectively.The solid lines correspond to λ GB = 0 (infinite coupling limit) while the dashed and dotted lines correspond to λ GB = −0.5 and λ GB = −2, respectively, all for second-order hydrodynamics.For all the ξ that were plotted the solid lines are well behaved because we have chosen ξ < ξ 2nd crit (λ GB = 0) = 0.503.For ξ = 0.2 the dashed curve becomes negative for small ∆y, indicating that second-order hydrodynamics becomes invalid faster for λ GB = −0.5.As mentioned earlier, this is what is indeed expected as the theory flows to weak coupling.However, the dotted curves are always positive in this range of ξ, which means that something qualitatively different is happening for sufficiently negative values of λ GB .In Figure 4 (right) we investigate this behavior in more 3 for some representative values of ξ = τ −1 0 w −3/2 = {0, 0.1, 0.2} depicted in blue, orange and green, respectively.Solid lines correspond to λ GB = 0 (infinite coupling result) while the dashed and dotted lines correspond to λ GB = −0.5 and λ GB = −2, respectively, all cases for second-order hydrodynamics.Right: Plot of ξ 1st crit (blue) and the two branches of ξ 2nd crit (orange and green) as a function of λ GB .In the ranges of λ GB ∈ (−∞, −1.657) and λ GB ∈ (0.073, 1/4], the correlator is positive but non-monotonic as a function of ξ.Here, ξ 2nd crit is found instead by requiring a monotonic decay at late times and is depicted in red.The dashed blue and orange lines correspond to the perturbative results to leading order in λ GB .The vertical line indicates the maximum allowed value for λ GB = 1/4.The behavior observed for negative values of λ GB in the range λ GB ∈ (−1.583, 0) is what is expected for a theory that flows from strong to weak coupling, i.e. ξ 2nd crit decreases as the coupling decreases.However, ξ 2nd crit increases in the range λ GB ∈ (−1.657, −1.583).The small square on top of the figure is a zoomed-in version of the same around this region.The dashed vertical line there signals the value of λ GB = −1.583for which dξ 2nd crit /dλ GB = 0.The discontinuous jump in the derivative of ξ 2nd crit at λ GB = −1.657 is likely to be an artifact of a truncated hydrodynamic gradient expansion or a truncated gravitational derivative expansion.
detail.In this plot we show the behavior of ξ 1st crit and ξ 2nd crit as a function of λ GB .The blue curve corresponds to ξ 1st crit and has precisely the expected behavior: it decreases monotonically as we decrease the value of λ GB .However, we observe something different for ξ 2nd crit : it has two branches for each value of λ GB , depicted in orange and green, respectively, which merge at two values of the coupling, λ GB = −1.657and λ GB = 0.073.For values of the coupling within the ranges λ GB ∈ (−∞, −1.657] and λ GB ∈ [0.073, 1/4] the correlator is always positive, however, non-monotonic with respect to ξ.In these ranges of λ GB we can find ξ 2nd crit by requiring monotonicity of the late-time correlator.The result of applying the latter criterion is depicted in red in Figure 4 (right).Combining these two criteria, we find that ξ 2nd crit decreases monotonically as λ GB varies from 0 to −1.583, but then increases again as λ GB goes from −1.583 to −1.657.Moreover, the derivative of ξ 2nd crit is discontinuous at λ GB = −1.657.Such behavior does not match the expectations for a theory that flows from infinite to zero coupling.It is likely that the inclusion of higher-than-secondderivative terms in the gravity action (beyond R 2 Gauss-Bonnet terms) or a higherorder hydrodynamic expansion would cure these problems.As a result, we conclude that the qualitative resemblance between non-perturbative λ GB -corrections and (nonperturbative) finite coupling corrections to the longitudinal two-point correlator, to second order in the hydrodynamic gradient expansion, is restricted to the range of λ GB ∈ (−1.583, 0].The critical times found for the longitudinal correlator can also be expressed generically in terms of a few theory-specific constants {η, Appendix C, and take the form: ) Expressing our coupling constants γ and λ GB collectively as β, first order corrections to τ 1st crit and τ 2nd crit take the form: The expression for τ 3rd crit now corresponds to the smallest real root of the equation where ξ = τ −1 0 w −3/2 .

Wilson loops
Wilson loops are another phenomenologically relevant non-local observable that can be studied within the framework explored in this work.The Wilson loop operator is a pathordered integral of the gauge field, defined as where the trace runs over the fundamental representation and C is a closed loop in spacetime.In AdS/CFT, the recipe for computing the expectation value of a Wilson loop, in the strong-coupling limit, is given by [125] W (C) = e −S NG (Σ) , ( where S NG = (2πα ) −1 × Area(Σ) is the Nambu-Goto action and Σ is an extremal surface with boundary condition ∂Σ = C.
Here, we consider two separate cases.The first case consists of a rectangular loop in the plane transverse to the boost-invariant direction of the Bjorken flow, where and → ∞.In the second case, we consider a rectangular loop with two sides extended along the longitudinal (beam) direction, y ∈ [− ∆y 2 , ∆y 2 ], x 1 ∈ [− 2 , 2 ] and → ∞.The calculation of the Wilson loop is qualitatively similar to that of the two-point function, so we will omit some of the redundant details below.

Transverse Wilson loop
The Nambu-Goto action for the transverse Wilson loop in the Fefferman-Graham chart is Using Eq. (B.10), we can expand this expression as NG + . . ., where

S
(1) and we used α = λ −1/2 .The first term is the pure AdS contribution, which we can see by using the zeroth-order embeddings: with z * = ∆x Γ[1/4] 2 /(2π) 3/2 .Integrating from → 0 to z * , and subtracting the divergent part, S div = √ λ/π , we obtain which gives the vacuum expectation value of the Wilson loop, (4.54) At next order, after using the zeroth-order embeddings and defining a dimensionless variable x = z/z * , we find where ũ0 = τ −2/3 0 w −1 .We observe that S NG depends linearly on c4 (ũ), similarly to S (1) for the transverse two-point function.Therefore the resulting values of τ 2nd crit and τ 3rd crit will be the same as those obtained in that case, for both Einstein gravity and the higher derivative gravities with α and λ GB corrections.As a result, the transverse Wilson loop provides no new bounds on the validity of the hydrodynamic description.

Longitudinal Wilson loop
The Nambu-Goto action for the longitudinal Wilson loop is which gives via (B.10)

S
(1) Again, the first expression gives the pure AdS embedding when we use the zeroth-order embeddings: with Integrating from → 0 to z * , and subtracting the divergent part S div = √ λ/π , we find i.e. the same result as in the transverse case.The next step is to evaluate S NG using the zeroth-order embeddings (4.59)-(4.60)along with the explicit forms of {ã 4 (ũ), b4 (ũ), c4 (ũ)}.Defining the dimensionless variable x = z/z * and expanding {ã 4 (ũ), b4 (ũ), c4 (ũ)} as in (4.32), we find where  We can extract τ crit at different orders in the hydrodynamic expansion by studying the ∆y → 0 behavior of S (1) NG .In this limit, the sinh 2 (∆y/2) term of (4.62) vanishes and the relevant I (k) integrals are finite and independent of k: Collecting our results, we find that for ∆y → 0: NG (0) itself does not provide a useful criterion for establishing the regime of validity of the hydrodynamic description at all orders in the hydrodynamic expansion, so we have to also impose monotonicity.The positivity criterion is enough only at first order, however S (1) NG (0) is strictly positive at second and third order in the backgrounds we consider.In these cases, we find that S (1) NG (0) decreases with decreasing τ until it reaches some minimum value, S NG,min (0, τ = τ min ), and then turns around and grows without bound (this behavior is demonstrated in Figure 5 for Einstein gravity).Therefore, for τ < τ min , the longitudinal Wilson loops are unphysical.This will be our criterion for establishing τ crit for the higher order hydrodynamic descriptions.
In the following, we will study the full behavior of S NG as a function of ∆y and ξ = τ −1 0 w −3/2 for our three cases of interest.In each case, the bounds on the validity of the hydrodynamic description are less constraining than those coming from the longitudinal correlator.
• Einstein gravity.Using the expansions of ã4 (ũ), b4 (ũ) and c4 (ũ) up to third order in hydrodynamics in (B.12), we evaluate S NG via (4.62), and plot the results for some representative values of ξ in Figure 5. From the ∆y → 0 limit of S (1) NG , we find: • α -corrections.Using the expansions of ã4 (ũ), b4 (ũ) and c4 (ũ) up to second order in hydrodynamics in (B.13) we evaluate S NG via (4.62), the results of which are shown in Figure 6.The solid lines correspond to γ = 0 (Einstein gravity) while the dashed NG τ 4/3 0 /w 4 √ λ∆x 3  3 for various values of ξ = τ −1 0 w −3/2 = {0, 0.15, 0.3, 0.45, 0.6} depicted in blue, orange, green, red and purple, respectively.The solid lines correspond to 3 rd order hydrodynamics; the dashed and dotted lines correspond to 2 nd and 1 st order hydrodynamics, respectively.Right: Plots for S(1) NG (0) for 1 st , 2 nd and 3 rd order hydrodynamics, depicted in green, orange and blue, respectively.The dashed vertical lines correspond to the critical times at each order in hydrodynamics.
7. The solid lines correspond to λ GB = 0 (Einstein gravity) while the dashed lines correspond to λ GB = −0.2.Following the same line of reasoning as in the previous two cases, we find: Finally, the critical times found above can be expressed generically in terms of the theory-specific constants defined in Appendix C, and take the form: Expressing our coupling constants γ and λ GB collectively as β, first order corrections to τ 1st crit and τ 2nd crit take the form:

Discussion
This work provides a new tile in the mosaic of recent developments on coupling-dependent thermal physics from the point of view of holography.With a view towards a better understanding of heavy ion collisions, the goal of this program has been to uncover qualitative and quantitative features of physical phenomena across a wide range of coupling constants-an understanding of which will likely require an interpolation between weaklycoupled perturbative field theory and strongly-coupled holographic techniques.Non-linear shock wave collisions were recently analyzed in perturbative Gauss-Bonnet theory to, for the first time, numerically model coupling-dependent heavy ion collisions [32] and, for example, compute the corrected hydrodynamization time.The extension of those results to either non-perturbative Gauss-Bonnet gravity or to type IIB supergravity is technically demanding.Therefore, it is useful to also study other, simpler models and probes of phenomena related to hydrodynamization.In this paper, we studied the gravity backgrounds dual to a boost-invariant Bjorken flow, which are good models for the late time dynamics of heavy ion collisions, at least in the regime of mid-rapidities.We considered non-perturbative Gauss-Bonnet gravity, studied in the present context for the first time, and type IIB supergravity (to leading order in α ), both to second order in hydrodynamics.Following up on [90], we provided an example of an analytically-tractable computation of a critical time defined through relaxation properties of non-local observables (equal-time correlators and Wilson loops), after which hydrodynamics becomes a good description.
Numerical estimates of the critical times obtained for second-order hydrodynamicscomputed to leading order in inverse 't Hooft coupling corrections in N = 4 theory and nonperturbatively in λ GB in Gauss-Bonnet theory-are summarized in Tables 1 and 2, where we show the increase of the critical time at decreased field theory coupling corresponding to a 10% and an 80% increase of η/s compared to its infinitely strongly coupled value of η/s = 1/4π.In both theories, the most stringent critical time is set by the longitudinal two-point correlator, φφ .and in a dual of Gauss-Bonnet theory at λ GB = −0.025.Both choices of the coupling correspond to a 10% increase of η/s.We use ⊥ and subscripts to denote transverse and longitudinal operators, respectively.
Several interesting features can be extracted from our analysis.One is the possibility of direct comparison between the size of effects of the 't Hooft coupling in N = 4 SYM and λ GB in the hypothetical dual of Gauss-Bonnet gravity.Such results should come in handy when using Gauss-Bonnet theory for phenomenologically relevant studies.The sec- and in a dual of Gauss-Bonnet theory at λ GB = −0.2.In this case, the choices of the coupling correspond to an 80% increase of η/s.
ond is the comparison between the sizes of perturbative and non-perturbative corrections in Gauss-Bonnet theory.As noted before, in both N = 4 SYM and Gauss-Bonnet gravity, the strictest bound on the regime of validity of hydrodynamics comes from the longitudinal two-point correlator.Since all other bounds are weaker, their non-convergent behavior in terms of the gradient expansion (third-order hydrodynamics giving a stricter bound than second-order hydrodynamics for φφ ⊥ , W (C) ⊥ and W (C) ) and in the perturbative λ GB expansion should not be taken seriously: at their respective critical times, the hydrodynamic description assumed in the derivation is no longer valid.What is important, however, is that for the critical time derived from the longitudinal φφ , the perturbative λ GB corrections converge remarkably quickly to the non-perturbative results, even for the increase of η/s by 80%.While perhaps surprising at first, this observation is compatible with the results of [32].
Another interesting consequence of our analysis is the emergent restriction on the range of the (non-perturbative) Gauss-Bonnet coupling for the second-order hydrodynamic approximation to a boost-invariant flow.While the Gauss-Bonnet theory with negative λ GB very well reproduces the expected behavior of a thermal CFT with finite coupling [31][32][33][34]43], it is also known that the theory suffers from instabilities and UV problems for large (or finite) values of λ GB .For the non-linear setup studied in this work, our computations suggest that the range of the non-perturbative coupling needs to be restricted to the interval λ GB ∈ (−1.583, 0].If we continue to decrease the Gauss-Bonnet coupling, then the bound on hydrodynamics becomes weaker, which is incompatible with the expectations for the behavior of a theory that flows from infinite to zero coupling.As is usual in holographic higher-derivative theories, we expect that in order to (reliably) flow from an infinitely coupled theory dual to Einstein gravity to a free thermal CFT, one would need to include an infinite tower of higher-order curvature corrections, beyond the R 2 terms considered in the Gauss-Bonnet theory, or the R 4 terms derived from type IIB string theory.We leave the investigation of these, and issues pertaining to finding phenomenologically relevant applications of non-local observables and the validity of hydrodynamics investigated in this work for the future.

A Second order solutions in perturbative Gauss-Bonnet gravity
As discussed in Section 3, the Gauss-Bonnet equations of motion can be solved at second order in the late-time expansion to first order in λ GB by writing the metric functions a 2 , b 2 and c 2 as a 2 = a 0 2 + λ GB â2 (and similarly for the other two functions) and expanding the equations of motion to first order in λ GB .The resulting system of equations is solved by where we have presented the solutions for c 0 2 and ĉ2 as first order derivatives due to the complexity of their integrated forms.Upon integration, the resulting integration constants are set by imposing AdS boundary conditions (see Eq. (3.19)).

B.1 Explicit expansions in Fefferman-Graham coordinates
We will consider three gravity solutions dual to Bjorken flow: Einstein gravity including 3 rd order hydrodynamics, perturbative α -corrections up to second order in hydrodynamics and non-perturbative λ GB -corrections up to second order in hydrodynamics: • Einstein gravity.The full gravity solution is known analytically only up to second order in hydrodynamics.However, the near-boundary metrics can be easily obtained for 3 rd order hydrodynamics from the expected stress-energy tensor and the corresponding transport coefficients [94].In particular, we find that at this order: where γ = α 3 ζ(3)/8 = λ −3/2 ζ(3)L 6 /8.As we can see, in the limit of infinite 't Hooft coupling λ → ∞ (or γ → 0) we recover the coefficients for second-order hydrodynamics in Einstein gravity (B.12).
• λ GB -corrections.The full gravity solution including non-perturbative λ GB -corrections and first-order hydrodynamics was obtained for the first time in the present paper.
Since the transport coefficients are known non-perturbatively up to second order in hydrodynamics [34], we can reconstruct the near-boundary coefficients explicitly.We find that: where γ GB = √ 1 − 4λ GB .For λ GB → 0 (or γ GB → 1) we recover the coefficients for second-order hydrodynamics in Einstein gravity (B.12).

C Useful definitions
We can express the critical times found in the previous sections generically in terms of a few theory-specific constants {η, Σ, Λ}, which correspond to contributions from first-, secondand third-order hydrodynamics, respectively.

Figure 1 :
Figure 1: Behavior of τ 2nd crit (λ GB ), non-perturbative in λ GB , coming from the transverse correlator.Negative values of λ GB resemble qualitatively the expected behavior as we flow from strong to weak coupling.

Figure 4 : 3 0
Figure 4: Left: Plots for S(1) ≡ S (1) τ 4/3 0 /w 4 ∆x4  3 for some representative values of ξ = τ −1 0 w −3/2 = {0, 0.1, 0.2} depicted in blue, orange and green, respectively.Solid lines correspond to λ GB = 0 (infinite coupling result) while the dashed and dotted lines correspond to λ GB = −0.5 and λ GB = −2, respectively, all cases for second-order hydrodynamics.Right: Plot of ξ 1st crit (blue) and the two branches of ξ 2nd crit (orange and green) as a function of λ GB .In the ranges of λ GB ∈ (−∞, −1.657) and λ GB ∈ (0.073, 1/4], the correlator is positive but non-monotonic as a function of ξ.Here, ξ 2nd crit is found instead by requiring a monotonic decay at late times and is depicted in red.The dashed blue and orange lines correspond to the perturbative results to leading order in λ GB .The vertical line indicates the maximum allowed value for λ GB = 1/4.The behavior observed for negative values of λ GB in the range λ GB ∈ (−1.583, 0) is what is expected for a theory that flows from strong to weak coupling, i.e. ξ 2nd crit decreases as the coupling decreases.However, ξ 2nd crit increases in the range λ GB ∈ (−1.657, −1.583).The small square on top of the figure is a zoomed-in version of the same around this region.The dashed vertical line there signals the value of λ GB = −1.583for which dξ 2nd crit /dλ GB = 0.The discontinuous jump in the derivative of ξ 2nd