Evolutions in first-order viscous hydrodynamics

Motivated by the physics of the quark-gluon plasma created in heavy-ion collision experiments, we use holography to study the regime of applicability of various theories of relativistic viscous hydrodynamics. Using the microscopic description provided by holography of a system that relaxes to equilibrium, we obtain initial data with which we perform real-time evolutions in 2+1 dimensional conformal fluids using the first-order viscous relativistic hydrodynamics theory of Bemfica, Disconzi, Noronha and Kovtun (BDNK), BRSSS and ideal hydrodynamics. By initializing the hydrodynamics codes at different times, we can check the constitutive relations and assess the predictive power and accuracy of each of these theories.


Introduction
Relativistic hydrodynamics is a very general theory that provides the effective description of the outof-equilibrium real-time dynamics of many physical systems of interest. In fact, hydrodynamics is better understood as a low energy effective field theory (EFT) constructed as a derivative expansion near a local thermodynamic equilibrium. As such, the fundamental variables are the basic quantities describing the equilibrium solution, such as energy density and velocities, and their gradients. The zeroth order piece in the gradient expansion corresponds to ideal hydrodynamics.
The equations of ideal hydrodynamics are hyperbolic, having a well-posed initial value problem for general initial conditions, and are amenable to numerical simulations. However, there are important physical situations where dissipative effects may play a relevant role, such as in the Quark-Gluon Plasma (QGP) [1,2] or in certain astrophysical systems [3][4][5][6][7]. However, it is wellknown that if first-order viscous terms are included in the relativistic hydrodynamic expansion (in the Landau frame), the resulting equations of motion are parabolic and hence incompatible with the causality postulates of relativity [8].
A well-known approach to address this issue of lack of hyperbolicity and causality is exemplified by the work of Müller, Israel and Stewart (MIS) [9][10][11]. The strategy consists in including second order terms and introduce new variables with the respective equations of motion, such that the overall system of equations is hyperbolic and causal. The new theory should have the same IR, i.e., hydrodynamic regime, but different UV properties. 1 The MIS formulation is not unique; there are several formulations which include extra variables and equations to recover hyperbolicity: BRSSS [12], DNMR [13], divergence type theories [14,15], etc.. As discussed in [16,17], in principle all these theories should be equivalent in the IR.
In recent years it has been realized that it is possible to write the equations of first-order viscous hydrodynamics in a way such that they are hyperbolic, without the need of adding extra variables or equations. When including first order terms in the hydrodynamic expansion, the basic hydrodynamic variables present an ambiguity when defined out of equilibrium, and fixing this ambiguity is known as choosing a 'frame'. Traditionally, it was common to use only natural frames such as the Landau frame. The insight of BDNK [18][19][20] was to notice that the choice of frame affects the hyperbolicity of the equations of motion and a certain frame can be chosen such that the latter are hyperbolic. 2 MIS-based hydrodynamic codes have been extensively used for many years to study dissipative effects. Given that there is now a new theory of viscous relativistic hydrodynamics, it is natural to explore it. The first ones to do so were the authors of [22], who studied the evolution of smooth and non-smooth initial data using BDNK in effectively 1+1 dimensions. The purpose of the present paper is to carry out the first studies of BDNK focusing on applications to the QGP. More concretely, in this article we report on studies of dynamical evolutions of the BDNK equations in a 2+1 dimensional conformal uncharged fluid. In order to acquire a better understanding of BDNK as a theory of viscous hydrodynamics, we quantitatively compare it to ideal and BRSSS [12] hydrodynamics, as well as to the microscopic theory provided by holography.
The QGP created in heavy-ion collision experiments is initially far from equilibrium and it subsequently relaxes to a hydrodynamic regime. Hydrodynamic numerical codes are used to provide effective descriptions of these evolutions, where the underlying microscopic theory is Quantum Chromodynamics (QCD). Motivated by this picture, in this paper we use numerical general relativity and holography to obtain a first principles description of the microscopic real-time dynamics of the stress tensor in a strongly coupled field theory of a system which relaxes to equilibrium. The microscopic solution is valid for all times, covering both the far-from-equilibrium and near equilibrium hydrodynamic regimes. Whilst previous studies of the applicability of hydrodynamics in holographic systems have focused on checking the constitutive relations, in this paper we go beyond the state of the art and use the microscopic data to initialize the hydrodynamic codes at different times 3 . In this way, by comparing with the UV-complete solution provided by holography, we verify both the constitutive relations and the predictivity of each hydrodynamic theory for a certain class of initial conditions. We use units c = G 4 = 1.
When our article was nearing completion, [24] appeared on arXiv, extending the numerical studies of the BDNK equations.

Hydrodynamics: the equations
We use three sets of hydrodynamic evolution equations: ideal hydrodynamics, BRSSS and BDNK. We write these equations for a 2+1 dimensional conformal fluid in Minkowski spacetime. Conformal symmetry fixes the equation of state: where is energy density and p is pressure. The temperature T is related to the energy density as = 2 3 A T 3 , where A is an arbitrary constant that in our simulations we fix to A = 4π 2 /9. 4 Considering the 3-velocity of the fluid u µ , normalized such that

2)
2 See [21] for a recent generalization to magnetohydrodynamics. 3 Previous studies in a one-dimensional setting were performed in [23]. 4 For a holographic fluid in d spacetime dimensions, this constant A is related to the bulk's Newton's constant of is symmetric, traceless and transverse to u µ . The constitutive relations of relativistic hydrodynamics up to second order in the derivative expansion in the Landau frame can be written as [12], where η is the shear viscosity and τ π , λ 1 , λ 2 , λ 3 are second order transport coefficients. The constitutive relations of ideal hydrodynamics are given by (2.3a) with Π µν = 0.
The conservation of the stress-energy tensor We now present the BDNK equations. Given a timelike unit vector u µ , a symmetric tensor can be decomposed as For a conformal fluid, the expansion in derivatives of each component to first order is where π 2 , θ 1 and η are transport coefficients. Conformal symmetry fixes the dimensions of π 2 and θ 1 , and we write them as a constant times the shear viscosity η π 2 = a 2 η , θ 1 = a 1 η . (2.9) The stress tensor, to first order in the derivative expansion, is then (2.10) The BDNK equations are obtained from plugging (2.10) into the conservation equation (2.4). The constants a 1 and a 2 specify the frame. The Landau frame corresponds to a 1 = a 2 = 0, by which we recover (2.3) up to first order terms. The BDNK equations are hyperbolic iff [25] The Landau frame lies outside this region, and there is a gap between the causal frames and the Landau frame.
The frame change to Landau frame is straightforwadly obtained from (2.10) In our simulations we use the transport coefficients of the microscopic holographic theory that will be presented in the next Section [26][27][28]: For details on the implementation of the hydrodynamic equations in the numerical codes and convergence tests see Appendix A.

Holography: microscopic evolution
Our model is Einstein's gravity with negative cosmological constant coupled to a massless scalar field in four spacetime dimensions. 5 By holography, it describes the decoupled sector of the stress 5 The scalar field is used to form a black brane by prompt collapse of the initial data. This prompt colapse happens immediately after initialization of the code and this timescale is well separated from the timescales relevant for the physics studied in this paper. Thus, after the immediate collapse the scalar field is completely negligible and the evolution will be well described by pure gravity. This pure gravity theory can be obtained as a consistent truncation in a set of top-down theories which includes the ABJM theory [29].
tensor on the conformal field theory side. The real-time quantum dynamics on the field theory maps to the bulk classical dynamics of Einstein's gravity in anti-de Sitter (AdS) space coupled to a massless scalar field. We use numerical relativity techniques to solve for the classical dynamics of the bulk fields employing the same code as in [30,31]. For more details of our implementation see these two references and Appendix B.
We consider the gravitational collapse of massless scalar field in the 3+1 dimensional Poincaré patch of AdS using time symmetric initial data. The scalar field has an initial deformed Gaussian profile along the boundary directions and is localized in the AdS radial direction. The initial data is "strong" in the sense that there is a trapped surface on the initial Cauchy surface. Therefore, one can think of our initial data as corresponding to a highly deformed black brane.
From the dual CFT point of view, the initial state corresponds to a large localized perturbation on top of a homogeneous plasma at temperature T . The initial energy density profile on the field theory is shown in Fig. 1 (top). The initial data is rotationally symmetric and has vanishing initial velocities (since the initial data is time-symmetric), even though our code does not assume rotational symmetry along the boundary directions. We have also considered initial data with small deformations that break the rotational symmetry along the boundary directions. However, the non-rotationally symmetric modes decay exponentially and by the time the system approaches the hydrodynamic regime, those modes have already decayed. We measure all quantities in units of T or E ≡ 2 3 AT 3 and we use Cartesian coordinates {t, x, y} to label the boundary directions.
The time evolution of the system is shown in Fig. 1 in three representative snapshots at times tT = 0, 0.08, 0.16. The initial peak explodes and the system expands and disperses away, and at late times it relaxes to a homogeneous solution.

Hydrodynamics: evolutions
In the following we study the applicability of hydrodynamics in this system. We do this in two ways: first we check the constitutive relations pointwise in spacetime and, second, we perform time evolution of the hydrodynamic equations using the holographic solution as initial data.
We start by evaluating the constitutive relations of hydrodynamics (2.3) pointwise in spacetime using the holographic data to calculate the various terms in this expression. In the Landau frame, gradient corrections in the hydrodynamic expansion are transverse to the velocity and, together with conformal symmetry, this implies that they vanish at the center of a spherically symmetric system. Away from the center however the gradients are non-trivial; we choose {x, y}T {0.12, 0} as a representative point where to evaluate the constitutive relations.
In Fig. 2 we compare the size of the 1st and 2nd order corrections to the ideal term in the derivative expansion of the stress tensor by plotting T 1st xx /T id xx and T 2nd xx /T id xx as functions of time. This figure indicates that, according to the constitutive relations, the system is initially far from equilibrium, as second order gradients are as large as 40% compared to the ideal terms. The size of the derivative corrections quickly decays and the system has hydrodynamized at times around tT 0.5, if our criterium is that the ratios T 1st xx /T id xx and T 2nd xx /T id xx are smaller than 2%. We now consider the evolutions of the hydrodynamic equations of the various theories using initial data from the holographic solution. We start by diagonalizing the holographic stress tensor and obtaining the local energy density and velocities (these are the energy and velocities in the Landau frame), which will be used to initialize the hydro codes. In the particular case of BRSSS, we have to also specify the dissipative part of the stress tensor Π µν from the microscopic holographic data. But notice that the use of time symmetric initial data implies that at t = 0 the dissipative part of the stress tensor vanishes identically and so in the subsequent evolution it is fairly small and hence one could initialize the code using very small initial data for the viscous tensor Π µν . However, initially the second order terms are large and therefore, alternatively, one could initialize de code using the dissipative part of the stress tensor computed with the first and second order consitutive relations (2.3), which is not small. These two possible ways of initializing the code match in the hydrodynamic regime (up to 3rd order terms), but under these far-from-equilibrium conditions they differ. We choose the second option, which captures the initial presence of large gradients. For the evolutions in the BDNK theory, we consider the following causal frames: which satisfy the hyperbolicity conditions (2.11). We obtain the energy and velocities in the causal frame from the Landau frame by inverting the expressions (2.12) and neglecting second order terms. The inverted expressions are the following: Where we have used the fact that in the first order terms we can replace = Landau , u µ = u µ Landau as the difference is second order. Recall that the evolution equations of BDNK are of second order and hence one also has to specify the time derivatives of the evolved variables in the initial data. We compute these by using the holographic data to obtain the constitutive relations in the corresponding causal frame as functions of time using (4.2), and then calculating the time derivatives in the causal frame. Other procedures might be equivalent in the hydrodynamic regime up to second order terms, but in the far from equilibrium regime they will generically differ.
We start the discussion by initializing the hydrodynamic codes using the holographic data at tT 0.019, when the system is still far from equilibrium according to the constitutive relations. Fig. 3 (top) shows the evolution of the energy density in the lab frame at a representative offcenter point of the domain {x, y}T = {0.17, 0} obtained using ideal (dashed green), BRSSS (dashed blue), BDNK in frame 2 (purple curves) and the holographic solution (solid black curve). For the BDNK evolutions, we plot two different quatities: the result of obtaining T tt from the evolved variables using (2.10) (solid line), and the result of evaluating only the 0th order terms (ideal tems) in (2.10) (dotted line). The reason for doing this is to show explicitly the size of the first order terms compared to the ideal terms. As this figure shows, the three hydrodynamic theories that we consider initialized at tT 0.019 exhibit a similar level of (large) disagreement with the microscopic theory throughout the evolution. The fact that none of the theories of hydrodynamics provides an accurate description of the evolution of the fluid confirms that at this initial time, the system is still very far from the hydrodynamic regime. This is expected since the gradients are large initially and each theory has a different UV completion. Therefore, there is no reason why these theories should agree with each other away from the hydrodynamic regime.
We continue the discussion by considering hydrodynamic evolutions with holographic initial data at a later time tT 0.16, when the size of the gradients is smaller but the system has not hydrodynamized yet according to the constitutive relations, see Fig. 3 (middle panel). Indeed, at this initial time the constitutive relations indicate that the first and second order terms are comparable, see the second vertical line in Fig. 2. We find that all hydrodynamic evolutions still exhibit large deviations from the microscopic theory which confirms that the system has not reached the hydrodynamic regime yet. Again, the deviations of BRSSS and BDNK from the microscopic theory are comparable.
Finally, we consider hydrodynamic evolutions with holographic initial data at a later time tT 0.58, when the size of the gradients is small, see Fig. 3 (bottom panel). At tT 0.58, by checking the constitutive relations (2.3) we find that the first order gradients have a maximum value in the whole domain of the order of 1%, i.e., T 1st xx /T id xx ∼ 1%, and similarly for the second order gradients, i.e., T 2nd xx /T id xx ∼ 1%. These values are small, and we expect the hydrodynamic evolutions to follow the microscopic solution when initialized at this time. Our results shown in the bottom panel of Fig. 3 confirm this expectation. This is true within the numerical errors of the holographic solution, which are inherited by the hydrodynamic evolutions through the initial data. We have checked that the energy density in the lab frame for the hydrodynamic evolutions differs by less than ∼ 2% for BRSSS, BDNK and ideal hydrodynamics, compared to the microscopic evolution in whole the domain and at all times beyond tT 0.58. These results also confirm that the various theories of hydrodynamics provide equivalent descriptions of the system, further supporting the applicability of hydrodynamics. Fig. 3 (bottom) suggests that ideal hydrodynamics provides a better description of the system at late times than viscous hydrodynamics. The reason might be the following. Our particular choice of initial data implies that the holographic system at t = 0 is exactly described by the constitutive relations of ideal hydrodynamics (2.3a) with Π µν = 0 (even if second order terms are large and the system is not within the regime of hydrodynamics) 6 . This aspect may partially survive through the far from equilibrium region explaining why ideal hydrodynamics provides a better description of the microscopic system than viscous hydrodynamics around the hydrodynamization time. We think that this might be a particularity of the specific system under consideration, and we wonder if this could be a generic feature of systems sharing similar initial conditions. See Appendix C for more details and further examples.

Discussion
We have used holography to obtain, from first principles, the real-time quantum dynamics of a large-N strongly coupled conformal field theory initially far from equilibrium which relaxes to a hydrodynamic regime. This constitutes a microscopic solution that we used to test the applicability of hydrodynamics. This solution is a very simple toy model that captures one important aspect of  the physics of the QGP created in heavy-ion collision experiments: the initial far-from-equilibrium conditions and subsequent relaxation to a hydrodynamic regime.
We have considered the applicability of three theories of hydrodynamics, namely ideal hydrodynamics, BRSSS and the newly formulated BDNK theory. For each of these theories, we evolve the equations of motion numerically using the microscopic solution as initial data at different times to explore the non-linear and far-from-equilibrium regimes. We assess the applicability of hydrodynamics by comparing both the constitutive relations, as previously done in the literature, and the evolution of the fluid as predicted by each theory. These comparisons are shown in Fig. 3, and we outline the main conclusions in the following.
The predictions of the different hydrodynamic theories differ from the microscopic solution and from each other in the far from equilibrium region, tT 0.019, 0.16 in our case. This is expected since each theory has a different a UV completion and, in this regime, the latter is relevant. On the other hand, at late times, tT 0.58 in our case, the gradients are small and hydrodynamics provides a good description of the system in the sense that the three theories that we have considered agree well with the microscopic evolution and with each other.
In this article we have studied BDNK as a theory of relativistic viscous hydrodynamics having in mind applications to the QGP. BDNK is particularly interesting because its only evolution variables are the thermodynamic quantities such as energy density and fluid velocities. To initialize the BDNK equations, we have to consider data for the fundamental thermodynamic variables and their first time derivatives in the causal frame obtained by inverting (2.12) and neglecting higher order terms. This is only justified in the hydrodynamic regime. We performed evolutions in different causal frames; only when the system is in the hydrodynamic regime we find that the physics is not affected by the choice of causal frame, up to second order terms. Away from the hydrodynamic regime, evolutions carried out in different frames in general differ.
In the QGP created in heavy-ion collisions the hydrodynamic evolutions provide a good description of the system even if gradients are not very small, which sometimes is known as 'unreasonable success of hydrodynamics'. In our case, we do not find such an unreasonable success, but instead we found that hydrodynamics applies precisely when it should. The theories that we considered and the initial data are very different from the QGP case, so there is a priori no reason why the success should happen also in our scenario.
One may wonder which causal formulation of viscous hydrodynamics, BRSSS or BDNK, provides a better description of the system. However, this may not be the right question to ask given that all theories of hydrodynamics should be equivalent in the regime of validity of hydrodynamics. Therefore, the question that one would really like to answer is which of these theories (if any!) provides a better (i.e., more in accordance to the microscopic theory) description of the system slightly outside the regime of hydrodynamics. It is likely that the answer depends on the details of the system that one wants to model and the initial conditions. In our case, where we have full control of the microscopic theory, we did not find that either BRSSS or BDNK provides a better description of the evolution outside the hydrodynamic regime. More studies are needed to address this question.
Hydrodynamics can be defined by considering a specific well-posed theory subject to suitable initial conditions. We may wonder if the constitutive relations evaluated pointwise in spacetime on the solution obtained by solving the initial value problem in hydrodynamics and those obtained from the microscopic solution may provide a different answer regarding the regime of applicability of hydrodynamics. This is an important question since in practical applications in general one does not have a theoretical control of the microscopic theory, and our solutions provide a concrete example where this question can be addressed. Above we have defined the hydrodynamization time using the second approach, obtaining tT 0.5. Alternatively, we could define the hydrodynamization time using the first approach, defined as the time at which, if we initialize the hydrodynamic code with microscopic data, the stress tensors along the rest of the evolution differ less than a 2% with the holographic one. For this case we obtain hydrodynamization times compatible with tT 0.5. Thus, both approaches provide a compatible answer. Having two different theories of relativistic viscous hydrodynamics, i.e., BDNK and BRSSS, allows us to use a third criterium, which might be useful in practical applications where the microscopic theory cannot be solved. If initialising both theories at a given time using suitable initial conditions, the evolutions differ by less than some prescribed small amount at all times, then the system is in the hydrodynamic regime. The caveat is that one has to be careful to ensure that the initial data for both theories is the "same".
With this work we initiate a program to study evolutions of the BDNK equations and to test the applicability of causal viscous hydrodynamic theories by comparing with microscopic holographic evolutions. Possible extensions include non-conformal and charged theories and initial data that models heavy-ion collisions. We hope that our work, together with [22], provide the first steps towards the implementation of the BDNK equations to describe relevant physical systems like the QGP or neutron star mergers.

A Implementation of hydrodynamic codes and tests
In our hydrodynamic codes we work directly with the primitive variables, namely energy density and velocities of the fluid, and the viscosity tensor in the case of BRSSS. This implies that our code cannot deal with shocks. That said, we point out that for the class of initial conditions that we have considered in this article, we did not observe the formation of shocks or steep features in the fluid flows that would require the use high resolution shock capturing techniques.
For the BRSSS case, we follow [32] and we impose the constraints u µ u µ = −1, u µ Π µν = 0 and Π µ µ = 0 algebraically to reduce the number of dynamical variables to 5, which we take to be U = {ρ, u x , u y , Π xx , Π xy }. Then we use Mathematica to solve the hydrodynamics equations of motion in terms of the time derivatives of U and write the them as The ideal case can be recovered from the BRSSS case by setting to zero all the transport coefficients and Π xx = Π xy = 0, as well as their derivatives. In our code, we discretize the spatial derivatives in (A.1) using 6th order finite difference stencils and we integrate them forward in time using a standard RK4 time-integrator. For simplicity, we imposed periodic boundary conditions; therefore, in our hydrodynamic simulations we have to choose a large enough domain to avoid boundary effects during the duration of the simulations. In practice this is not a problem since the hydrodynamic simulations are very cheap. The equations of motion for BDNK are second order in time and space. In this case, we implement the constraint u µ u µ = −1 algebraically to reduce the number of independent variables to 3; we choose U = {ρ, u x , u y }. Again, we use Mathematica to solve for the second time derivatives of the dynamical variables, and re-write the equations of motion as One could write (A.2) as a first order in time system of equations in the obvious way, e.g., P ≡ ∂ t U, and use a standard integrator such as RK4. However, we found that the resulting system was numerically unstable; it is not clear to us what is the origin of this instability. We did not attempt to write (A.2) as a fully first order system by further defining V i ≡ ∂ i U. As in the BRSSS and Ideal cases, in BDNK we discretise the spatial derivatives using 6th order finite differences and we impose periodic boundary conditions. To proceed, we implemented two different time integrators.
First, we implemented an implicit second order in time scheme by discretizing the time derivatives of a given variable f at the time t n = n ∆t as: where i denotes a collection of indices that labels a given spatial grid point. In this way, for known U n i and U n−1 i , the discrete equations of motion at time t n become a non-linear algebraic system for the values of the variables at the next time level, U n+1 i , which can be solved using standard techniques; in our case, for the BDNK equations we used a Newton-Raphson algorithm. This method is robust and it works well in practice, and it has been successfully used in simulations of black hole binary mergers [33] and in our holographic simulations [30,31,34]. However, for high spatial resolutions, it becomes considerably slower and memory-demanding than an explicit time integrator such as RK4. Therefore, for practical applications that one could easily run on a laptop, such as the ones presented here, it would be desirable to have a stable explicit time integrator for BDNK. We now turn to this.
We have successfully implemented an explicit time integrator for the BDNK equations of motion, written as second order in time, i.e., (A.2). The analogues of the Runge-Kutta methods for second order equations are known as Runge-Kutta-Nyström Generalized (RKNG) methods. We will review them here since they may not be as well-known as the standard first order methods. Here we follow [35] since we have implemented their RKNG34 scheme. This method is competitive in terms of both accuracy and efficiency with respect to RK methods of similar orders, see [35] for detailed comparisons.
Consider a system of N second order ordinary differential equations: The generalization of this method to PDEs is straightforward, just as in the standard RK4 case. An explicit RKNG method of s stages allows to compute the approximations y n+1 and y n+1 of the solution y(x n+1 ) and its derivative y (x n+1 ) at x n+1 ∈ I from their values in the previous steps, y n and y n , as follows: where f 1 = f (x n , y n , y n ) , (A.8) and h = x n+1 − x n . Here c i , b i and b i are s-dimensional constant vectors, and a ij and a ij are lower triangular s × s constant matrices; together, they constitute the parameters that define the method. An RKGN method is said to have order p if the local convergence order for both the approximation y n+1 that of the derivative y n+1 is p, centred at r = r 0 , with width R and amplitude A, on top of a background energy density ρ 0 . The initial values for all the remaining variables are chosen to be zero. For the runs shown below, we chose a computational domain of size L = 10, ρ 0 = 0.5, A = 0.5, R = 1 and we centred the Gaussian at the centre of our grid. For the BDNK simulations, we choose a frame with a 1 = a 2 = 10. We monitor the convergence rate Q N (t), where || · || 1 denotes the numerical L 1 -norm taken over our computational domain, and h is the grid spacing. To compute these norms, we sum over all the evolution variables.
The results of our convergence tests are displayed in Fig. 4 for the three theories that we consider in this article. As this figure shows, the convergence rate is close to 4 in all cases. Even though in our code we use 6th order spatial differences, this figure shows that the dominant error comes from the time integration. This justifies the use of 6th order Kreiss-Oliger dissipation. Furthermore, the convergence rate that we find is consistent with the order of the time integrators that we used in our code, namely RKNG34 for BDNK and RK4 for BRSSS and ideal hydrodynamics.

B Numerical relativity details
In this appendix we summarize the basic aspects of our numerical relativity simulations. We use the same code as in [30,31], so more details can be found in these references.
We solve the Einstein equations in AdS coupled to a massless scalar field in 3+1 dimensions using generalized harmonic coordinates, where n α = −∂ α t is the timelike, future-directed unit 1-form normal to slices of constant t, and κ and P are the damping parameters. Here C α ≡ H α − x α are the constraints and H α are the source functions, which we can freely prescribe. We choose both the parameters and the source functions as in [30,31]. For convenience, we couple gravity in AdS to a massless scalar field, since this allows us to easily construct far-from-equilibrium initial data by prescribing a suitable profile for the scalar field. Therefore, in addition to the metric, we also evolve a massless scalar field, with equation of motion φ = 0 , (B.1) and stress tensor To carry out the evolution, we write the general spacetime metric g µν =ĝ µν +ḡ µν , whereĝ µν is the metric of the Poincaré patch of AdS 4 (with the AdS length L = 1) written aŝ where z = (1 − ρ 2 )/ρ 2 is the usual AdS radial coordinate, and x, y are the spatial boundary directions, which we take to be Cartesian directions with infinite range. In practice we compactify  these directions as x = tan( π 2x 1 ), y = tan( π 2x 2 ).ḡ µν is a general deviation away from AdS 4 , not necessarily small, and it satisfies Dirichlet boundary conditions at the AdS boundary. Similarly, we write the scalar field as φ = (1 − ρ 2 ) 2φ , (B.4) and the evolved variableφ vanishes at the boundary of AdS 4 . Following [30,31], we write the source functions as whereĤ µ are the source functions for AdS 4 in the coordinates of equation (B.3) andH µ are the actual evolved source functions; the power of (1 − ρ 2 ) in (B.5) has been chosen so thatH µ ρ=1 = 0 at the boundary of AdS. Following the prescription of [31], we fix our gauge such that near the boundary of AdS, we haveH near the boundary. The source functions are set to zero in the interior of the spacetime. We choose the same smooth transition functions as [31] to fix the source functions everywhere in the spacetime. We construct time symmetric initial data by prescribing a Gaussian profile for the initial scalar field, where A and ∆ control the amplitude and the width of the Gaussian respectively, and Here ρ 0 allows us to localize the initial Gaussian in the AdS radial direction, while the constant parameters {a ρ , α 1 , α 2 } control the shape along the boundary directions. Having specified an initial profile for the scalar field, the Hamiltonian constraint is solved as in [31]. We choose "strong" data so that there is a trapped surface, with planar topology, in the initial data slice. We extract the stress energy tensor of the boundary CFT as in [31]. We find,    5 shows a convergence test for the trace of the boundary stress tensor of the evolution presented in the main text, which shows first order convergence with the number of points along the holographic direction N ρ .
A convergence test of the conservation equation of the stress tensor for our code was presented in [34], and we do not perform this convergence test in our case. Fig. 6 shows the temporal component of the conservation equation for the stress tensor, ∂ µ T µt , for our highest resolution case N ρ × N x × N y = 1025 × 129 × 129. We present the result at the center of the domain because the numerical error is maximal there. For the other components, ∂ µ T µx , ∂ µ T µy , the error is smaller.
We use the numerical error of the conservation of the stress tensor, Fig. 6, to provide an estimate of the numerical error in the first order terms of the hydrodynamic expansion. Schematically, the first order terms are given by the shear viscosity times derivatives of the velocity, so we estimate the error as T 1st We estimate T id xx ∼ E/2 and ∂ x v x ∼ ∂ µ T µt /(2E). In the last expression, by dividing by E we obtain the correct units, and the factor of two is because at the center the ∂ x T xt and ∂ y T yt terms contribute equally and from the actual data we check that the ∂ t T tt term is much smaller. For times around the hydrodynamization time tT 0.5, the right hand side of (B.15) gives an error of the order T 1st A numerical error of the order of ∼ 2% is of the same order of the first and second order terms around times tT 0.5, see Fig. 2. So, all the statements in this paper about the applicability of hydrodynamics for quantities beyond this time must be considered only up to numerical errors.
The change from the Landau frame to the causal frame is specially affected by the error in the conservation equation of the stress tensor. This is because the first order terms in the change of frame expression are precisely the same terms as in the ideal conservation equation of the stress tensor. For this reason, we subtract by hand the error of the conservation of the stress tensor in the change of frame to obtain the initial data for the BDNK evolutions.

C Details on the dynamical evolutions
In this appendix we provide more details about the dynamical evolutions presented in the main text. Fig. 2 shows the ratios T n xx /T ideal xx , n = 1, 2, of the constitutive relations (2.3) evaluated at {x, y}T {0.12, 0} as a function of time for the holographic solution. This shows explicitly the validity of the hydrodynamic expansion. Alternatively, the constitutive relations can also be compared with the microscopic stress tensor values, to assess if the constitutive relations describe the microscopic evolution. This is shown in Fig. 7. On general grounds, we would expect that if the system is in the hydrodynamic regime, then the gradient corrections should improve the ideal description. However, in our system, at the hydrodynamization time tT 0.5 first and second order terms do not improve the ideal description of the microscopic solution. As stated in the main text, this might be explained by our choice initial data. We expect that at later times viscous descriptions would improve the ideal result. However, at later times the numerical error (discussed in appendix B) does not allow us to make conclusive statements in this regard. Also, possibly for the same reason, ideal terms provide a good description of the system even before the hydrodynamization time tT 0.5. The microscopic T xx presented in Fig. 7 presented a numerical oscillation of a cero mode (homogeneous). This numerical oscillation converges to zero similarly to Fig. 5, and we removed it by hand in Fig. 7. Fig. 8 shows the same evolutions as in Fig. 3 but at an off center location {x, y}T = {0.17, 0}. The conclusions are similar to those obtained from Fig. 3. If the gradients are large (top and middle), the hydrodynamic evolutions provide a description that differs from the microscopic one, and are also different among them, as all theories have different UV behaviors. At late times (bottom), when gradients are small, the hydrodynamic theories provide a better description of the system, matching the microscopic theory within a maximum deviation in the domain of a ∼ 2% in the energy density. Fig. 8 also includes the evolutions of the BDNK equations in frame 1, and we discuss now the differences of evolving the BDNK equations in the different causal frames (4.1). First, recall that a 1 and a 2 are multiplying first order terms in (2.10). This means that the larger these values, the larger are those first order terms. So, for frame 1 these terms are larger than for frame 2. This is the reason why in the main text we focus on the frame 2, for which the values for a 1 and a 2 are closer to the smallest ones allowed by causality. The size of the constants a 1 and a 2 will be relevant when changing from the Landau frame to the causal frame, when evolving the equations, and also when going back to the Landau frame. We observe that the size of these first order terms (the ones proportional to a 1 and a 2 ) in the microscopic solution is typically small at all times, even far from equilibrium. On the other hand, in the BDNK simulations in the far from equilibrium regime, these terms become large, even of order one, see Fig. 8 (top and middle). Moreover, the BDNK evolutions in the different frames in the far from equilibrium region may provide considerably different results, see Fig. 8 (middle). Finally, in Fig. 8 (bottom) we see that the difference between frame 1 and 2 is compatible with second order terms, in accordance with the fact that the system has hydrodynamized.  We include the results of the hydrodynamic evolutions initialised with holographic data at tT 0.019, 0.16, 0.58, from top to bottom. For the BDNK evolutions, we include the result of using only the 0th order terms (the ideal part) in (2.10), in dotted lines. We also include the evolutions of BDNK in the frame 1.