Holographic Isotropisation in Gauss-Bonnet Gravity

We study holographic isotropisation of homogeneous, strongly coupled, non-Abelian plasmas in Gauss-Bonnet gravity with a negative cosmological constant. We focus on small values of the Gauss-Bonnet coupling parameter $\lambda_{GB}$ and linearise the equations of motion around a time-dependent background solution with $\lambda_{GB}=0$. We numerically solve the linearised equations and show that the entire time evolution of the pressure anisotropy can be well approximated by the linear in $\lambda_{GB}$ corrections to the quasinormal mode expansion, even in the cases of high anisotropy. We finally show that, quite generally, the time evolution of the pressure anisotropy with the Gauss-Bonnet term is approximately {\it shifted} with respect to the evolution without it, with the sign of the shift being directly related to the sign of the $\lambda_{GB}$ parameter. This suggests that finite coupling corrections generically {\it increase} the isotropisation time of strongly coupled plasmas.


Introduction
The gauge/string duality is a powerful theoretical laboratory in which to study the dynamics of strongly coupled, non-Abelian gauge theories. In the limit of large number of colours (N c ) and large 't Hooft coupling (λ), the duality maps complicated, fully quantum computations on the gauge theory side into simple and tractable classical gravity problems. This access to the strongly coupled sector of a large class of non-Abelian theories has led to many insights into the physics of strongly coupled matter at very different energy scales, from deconfined QCD matter to the behaviour of non-Fermi liquids, superconductors and cold atom systems (see [1][2][3][4] for recent reviews).
One of the areas of physics in which the gauge/string duality is most powerful is the analysis of far-from-equilibrium dynamics. Gravity computations have been employed to study, among many other subjects, the reaction of non-Abelian gauge theories to quenches [5][6][7][8][9], the formation of turbulence [10][11][12][13], the approach towards hydrodynamic behaviour of large disturbances of non-Abelian theories [14][15][16][17][18][19], as well as the characterisation of the debris of the collision of energetic projectiles [20][21][22][23][24]. Nevertheless, most of the analyses performed up to date focus on the infinitely strongly coupled limit and very little is known about finite coupling effects to those dynamics. In this work we present the first step towards understanding finite coupling corrections to off-equilibrium dynamics.
According to the holographic dictionary, finite coupling corrections in the gauge theory correspond to high curvature corrections on the gravity side. For the particular case of N = 4 supersymmetric Yang-Mills (SYM) theory, the complete set of high curvature, R 4 , terms responsible for leading order corrections in the inverse 't Hooft coupling are known [25]. These have been used in the past to determine the corrections to equilibrium properties of N = 4 SYM, such as the free energy [25], to the transport properties [26][27][28], as well as the equilibrium photon emission rate of the plasma [29]. Corrections to nearequilibrium dynamics have been also addressed [30][31][32], by computing the relaxation rates of small deviations from equilibrium, which on the gravity side are controlled by characteristic relaxation modes of black branes, known as quasinormal modes (QNM) [33]. Higher curvature corrections to thermalisation dynamics within a specific holographic construction in which AdS-Vadya black holes model the thermalisation after a sudden injection of energy, have also been studied [30,31,[34][35][36]. Quadratic curvature effects in the off-equilibrium dynamics of homogenous matter in an expanding universe have been also recently considered in [37]. For a compilation of finite coupling corrections to infinitely strongly coupled N = 4 SYM see [38].
In generic holographic constructions, one expects the leading higher curvature correction to be quadratic [39]. A particular model of this kind is Gauss-Bonnet gravity, which is based on a particular set of R 2 corrections to the Einstein-Hilbert action with negative cosmological constant controlled by a single parameter λ GB [40][41][42][43][44][45][46]. While on the field theory side the holographic dual of Gauss-Bonnet gravity is unknown, this setup is appealing since the equations of motion remain of second order, so in principle we expect the theory to be free of the pathologies induced by generic higher derivative terms. 1 In addition, unlike the higher curvature corrections of N = 4, the black hole solution of the Gauss-Bonnet theory with negative cosmological constant can be found analytically [42], so some calculations are amenable to treatment beyond perturbation theory in λ GB in a convenient way. Recently, the analysis of the relaxation of small out-of-equilibrium perturbations for arbitrary values of λ GB has been performed in [32]. This revealed a very interesting behaviour of the relaxation for large values of the λ GB parameter, which qualitatively resembled the expected behaviour of N = 4 SYM plasma in the small 't Hooft coupling limit. Furthermore, for small negative values of λ GB , the authors of [32] found that the structure of small corrections to the QNM spectrum is qualitatively similar to that of infinitesimal 1/ √ λ corrections to the relaxation dynamics in strongly coupled N = 4 SYM computed in [30,32,34]. These facts indicate that Gauss-Bonnet gravity provides a good testing ground for understanding the effects of finite coupling corrections on dynamical processes in strongly coupled N = 4 SYM.
In this paper, we study the problem of isotropisation in Gauss-Bonnet gravity. At an initial time, we prepare a homogeneous, anisotropic off-equilibrium state of the field theory dual to Gauss-Bonnet gravity, and study the process by which the system becomes isotropic. The initial disturbances are large, meaning that the difference in pressure between 1 The inclusion of higher curvature corrections in the holographic context has been called into question by [47], which found that certain causality pathologies arise in the graviton three point functions of such theories unless one includes the full tower of corrections coming from a stringy model; see however [48]. the anisotropic direction and the transverse directions is large in units of the energy density, which is constant as a consequence of the homogeneity of the state. On the gravity side, this setup corresponds to specifying the dual metric on the initial time-slice and using Einstein's equation to evolve the different metric fields. In the limit of small λ GB , we will consider the effect of higher curvature corrections as a small (linear) perturbation on top the non-linear evolution of the metric at λ GB = 0. In this way, we extract the leading order correction in λ GB to the isotropisation dynamics of the dual field theory.
The problem of non-linear isotropisation in N = 4 SYM has been studied in the past [15,[17][18][19]24]. Quite remarkably, those analyses have shown that, in spite of the intrinsically non-linear setup, the full out-of-equilibirum dynamics of this strongly coupled theory can be described as a linear superposition of the relaxation modes of the static black brane dual to the thermal system for a large class of initial conditions. This observation provides a tremendous simplification, since it renders the complicated dynamics of isotropisation into a linear problem. In this paper we will find that in the presence of higher curvature corrections this simplification is preserved for moderately large (O(1)) anisotropies, although not for arbitrarily large ones. As we will show, for the initial configurations described above, the evolution of the perturbed plasma can be approximated as a linear combination of the QNM of the black brane solution in Gauss-Bonnet gravity.
We also investigate the effect of the Gauss-Bonnet term on isotropisation time, defined as the time at which the ratio of the pressure anisotropy to the average pressure relaxes beyond a given threshold criterion. As observed in previous studies [17][18][19]24], the isotropisation time is not unique, but depends on the initial configuration. Similarly, the corrections induced by the Gauss-Bonnet term depend on the initial perturbation. Nevertheless, by studying many different initial configurations, we have observed that for (small) negative values of λ GB , the isotropisation time is always larger than for λ GB = 0, while the opposite is true for positive λ GB values. In the linearised regime of moderately large anisotropies, we demonstrate that this correlation of the change of isotropisation time with the sign of λ GB always holds, irrespective of initial conditions. Quite satisfactorily, Gauss-Bonnet gravities with negative (positive) λ GB values are dual to gauge theories with shear viscosity to entropy density ratio η/s, larger (smaller) than 1/4π [45], the value of the ratio for infinitely strongly coupled SYM. This correlation further supports the interpretation that Gauss-Bonnet gravity represents a good tool to understand the finite coupling corrections of N = 4 SYM. This paper is organised as follows: in Section 2 we will review the isotropisation of large initial configurations in N = 4 SYM, and also gather the key ingredients of Gauss-Bonnet gravity that we will need for our analysis. In Section 3 we introduce the numerical procedure we employ to study small λ GB corrections to the isotropisation process, and discuss the main systematics of our numerical results over many different initial configurations. In Section 4 we focus on the effective linear regime, and show how the systematics observed in Section 3 follow from the analysis of the associated quasinormal mode expansion. In Section 5 we perform a systematic exploration over initial conditions to determine the validity of our findings in the effectively linear regime. Finally, in Section 6 we summarise our main results and discuss their implications for finite coupling corrections of holographic theories.

Isotropisation in strongly coupled N = 4 SYM
In this section we briefly review the procedure to study the isotropisation of far-fromequilibrium initial states in strongly coupled N = 4 SYM, following [15,17,18,24]. This section will also serve to establish the notation and the general strategy to solve the isotropisation problem with higher derivative corrections, which we will perform in the next section.
We prepare an anisotropic yet homogeneous state of N = 4 by specifying an initial value for the bulk metric in the gravity dual. We will parametrise the (4+1)-dimensional space dual to the field theory state by the following metric ansatz: where A 0 , Σ 0 and B 0 are all functions of t and r only, and the boundary is at r → ∞. This ansatz enjoys and residual gauge freedom which arises from reparametrizations of the holographic coordinate, r → r + λ(t). Introducing this ansatz into Einstein's equations and imposing that the space is asymptotically AdS, the near boundary expansion of the different metric fields is given by where λ(t) is arbitrary, andĝ (4) 0,11 and a (4) 0 are unknown functions of time which cannot be determined from a power series expansion. Here L is the AdS radius defined in terms of the cosmological constant Λ via The asymptotic expansion (2.2) determines the stress tensor of the dual gauge theory via holographic renormalisation [49,50] to bê with T ab = 2T ab L 3 /κ 2 5 and κ 2 5 the five dimensional Newton constant, which is related to the number of colours of the dual gauge theory by L 3 /κ 2 5 = N 2 c /4π 2 . We work in units of L = 1 henceforth.
With the gauge choice in Eq. (2.1), Einstein's equations contain only 5 non-vanishing components. After defining the modified time derivative as Einstein's equations take the following nested form where primes denote r−derivatives. One can check that Eq. (2.10) yields da (4) 0 (t)/dt = 0, which, via Eq. (2.4), implies that the energy density of the system is conserved. Otherwise, it does not participate in the dynamics and can be dropped as long as a (4) 0 is held constant. This nested form provides a convenient setup in which the metric functions Σ 0 , d + Σ 0 , d + B 0 and A 0 can be determined sequentially at every time slice by solving linear ODE's, once the r-dependence of B 0 is known. Therefore, the time evolution of the metric functions can be obtained after determining ∂ t B 0 from Eq. (2.5) and knowledge of d + B 0 and A 0 at each time slice. The above nested pattern allows us to specify a full set of initial states of the dual gauge theory by specifying the functional form of the metric field B 0 (t 0 , r) at some initial time t 0 . These states provide a convenient framework in which to study far-from-equilibrium dynamics with a variety of initial conditions.
To solve the dynamical equations (2.6)-(2.9), it is convenient to redefine the metric fields to facilitate the imposition of the boundary conditions (2.2). Following [24], after the coordinate transformation u ≡ 1/r, we define: (2.11) Note that in these redefinitions the functionsσ andḃ are not the time derivatives of σ and b, but rather independent functions in the same way in which the metric functions d + B and d + Σ are independent from B and Σ at each time slice in the nested procedure to solve Einstein's equation. Imposing Eq. (2.2) and the boundary expansion of d + Σ 0 and d + B 0 , these redefined fields satisfy the following boundary conditions in the u → 0 limit: 0 ,ḃ 0 → −2ĝ For a given gauge choice λ(t), the nested system of equations (2.6) together with the redefinitions (2.11) and the boundary conditions (2.12), fully specify the time evolution of the system once the initial data b 0, init = b 0 (0, u) is specified. For the computations of this work, in the λ = 0 gauge, we consider manny different arbitrary initial data b 0, init . These are constructed as the ratio of two 10th order polynomials in u whose coefficients are chosen randomly in the range [0,1]. We will furthermore multiply this ratio by a random amplitude (in the range [0,10]), and also subtract the constant term, so that the boundary conditions (2.12) are obeyed 2 . We have also tested the Gaussian initial conditions studied in [24] with varying amplitude.
At each time slice, we solve for the metric components via pseudospectral methods on the Chebyshev grid in the holographic coordinate u. Since generic choices of b 0, init will lead to the formation of an apparent horizon in the bulk metric [18], we will choose our grid to be u ∈ [0, u H ], where u H is the location of the apparent horizon, given by the condition (2.13) In the λ = 0 gauge this location changes with time which complicates the application of pseudospetral methods. For this reason, after locating the apparent horizon in the initial time slice, we reparametrise the holographic coordinate such that the position of the horizon is fixed. Without loss of generality we choose λ such that u H = 1 in our numerics. This, in turn, implies that λ becomes a dynamical variable that must be updated at each step on the time integration. The time derivative of λ may be found via the horizon stationarity condition [24], i.e. the time derivative of Eq. (2.13), which, after using the equations of motion, boils down to 3 (2.14) To obtain ∂ t λ at a given time slice, we solve the first three nested equations (2.6)-(2.8) using the explicit boundary conditions (2.12). After solving these equations, Eq. (2.14) may be viewed as imposing a boundary condition for Eq. (2.9) at the apparent horizon. We can then solve for a 0 (t, u) and from Eq. (2.12), ∂ t λ is then given by−a 0 (t, 0). Combining this derivative with d + B 0 , obtained from Eq. (2.8) and with the redefinition of the fields Eq. (2.11), we can determine ∂ t b 0 . This allows us to determine b 0 at the next time slice and, subsequently, all other metric functions. Iterating the procedure at every time step determines the metric at all times. We perform the time evolution using a fourth order Runge-Kutta algorithm. From the evolution of the metric we can extract the stress tensor of the anisotropic state in the dual field theory at later times. Given that the energy density remains constant throughout the evolution and that the trace of the stress tensor vanishes identically, the only non-trivial component is the pressure anisotropy, ∆p 0 ≡ T 0,zz − (T 0,xx + T 0,yy )/2.
Combining Eqs. (2.4), (2.11) and (2.12), the pressure anisotropy is given by where the hatted quantities have the same definition as in Eq. (2.4). For very anisotropic initial estates, the pressure anisotropy ∆p 0 can be large at initial time t 0 . As time passes, the anisotropy decays until, at sufficiently late time the geometry becomes that of the AdS-Schwarzschild black hole, given by Eq. (2.1) with In terms of the field theory dual, the stress tensor relaxes to the equilibrium stress tensor in which all pressures are the same and equal to p eq = /3, with the energy density which is determined by a This relaxation process is called isotropisation and has been studied in detailed in previous works [15,17,18,24]. These studies constitute the basis we will employ to study isotropisation in Gauss-Bonnet theories.

Gauss-Bonnet gravity
In this section we collect some useful results about Gauss-Bonnet gravity with a negative cosmological constant [42]. The five-dimensional action we consider takes the form where λ GB is a dimensionless number, constrained (at least for holographic purposes) by causality [45] and positive definiteness of the boundary energy density [46] to be 4 Einstein equations coming from Eq. (2.18) can be written as where E 0 are the zeroth order (λ GB = 0) Einstein equations 21) and the contribution from the Gauss-Bonnet term is given by The near boundary behaviour of the metric fields can be easily found to be where This follows from the fact that the solutions of Eq. (2.18) are asymptotically AdS with the effective radius L c [52], and can be of course checked explicitly. The holographic renormalisation for Gauss-Bonnet gravity has been performed in [52,53], where the covariant counterterms were computed. The resulting boundary stress tensor turns out to be (2.25) where Roman indices go over the boundary directions, γ ab is the induced metric on the boundary, K ab is the extrinsic curvature of the boundary (and K is its trace), and Q ab is a tensor (whose explicit form can be found in [52]) given in terms of the extrinsic curvature, the Ricci scalar R, Ricci tensor R ab and the Riemann tensor R abcd associated with the boundary metric γ ab .
The expectation valueT ab of the dual theory stress tensor is then given by where hat denotes the same rescaled definition of the stress tensor as in Eq. (2.4), and where h ab is the background metric on which the dual theory lives, given by which, using the boundary conditions (2.2) with L → L c , evaluates to diag(−1, 1, 1, 1), as it should. In the end, we get: Note that setting λ GB = 0 or equivalently L c = 1, we arrive at the same expression as in Eq. (2.4).
Similarly to the asymptotically AdS case, we expect the final state of the isotropisation in Gauss-Bonnet to be the corresponding black brane solution in this theory. The line element takes the form (2.1) with [42] The event horizon is located at r = r H . The choice of the minus sign in front of the square root in Eq. (2.29) makes the limit λ GB = 0 well defined so that this geometry is smoothly connected to that of the AdS-Schwarzschild black hole. 6 Note that the boundary speed of light associated with Eq. (2.29) is unity. The Hawking temperature of this solution is given by and the energy density, as given byˆ =T 00 in Eq. (2.28), turns out to bê

Time evolution in Gauss-Bonnet gravity
Our strategy to study time evolution in Gauss-Bonnet gravity is to linearise the problem in the coupling λ GB . In addition to the obvious technical advantage that the evolution equations are linear, this approach ensures that the dynamics remain hyperbolic, so that we can evolve in time for a given initial condition specified on a null slice, as we shall see explicitly below. For finite values of λ GB , characteristic surfaces -i.e. the places on which we wish to specify initial data -in general do not coincide with metric null cones, which in certain cases leads to the initial value problem being ill-defined [51,54] 7 . Working in perturbation theory around λ GB = 0 allows us to overcome this difficulty because the dynamical structure of the perturbations is that of the underlying Einstein-Hilbert equations of motion, thus, much of the formalism of [15] which we employ to determine the background solution carries over to the Gauss-Bonnet case.

Linearised time evolution
Having found a (time-dependent) background solution g 0 (t, r) of the non-linear equations E 0 , we wish to linearise Einstein equations around it by postulating a solution of the form Plugging this into Eq. (2.20) and keeping terms up to first order in λ GB , we get where E 0,lin are the zeroth order Einstein equations linearised around the background solution g 0 (t, r), for which the Gauss-Bonnet contribution E GB evaluated on the background acts as a source. In addition to the terms determined explicitly by the background evolution, the source receives contributions which involve can be determined from the zeroth order equations E 0 , these do not contain d + d + B 0 nor d + A 0 . This does not present a problem because we can obtain these extra terms (or simply ∂ 2 t B 0 and ∂ t A 0 if one wishes) numerically since we know the background solution for all t and r.
The form of the resulting linear equations (3.2) makes it clear that the dynamics of the perturbations δg inherit the structure of the background equations of motion so that, in particular, we can cast the evolution problem as a set of nested ODE's. 8 To do so, it suffices to appropriately linearise the definition of d + by letting, for any perturbation δF , The steps to solve for the linear evolution problem then closely parallel those of the background: after specifying the initial condition for the metric anisotropy on an initial time slice and the boundary condition that fixes the perturbed energy density, we solve the set of ODE's for δΣ,d + δΣ,d + δB and δA, which allows us to calculate ∂ t δB and time evolve to the next slice. Notice that the small perturbation problem we are considering can be formulated in the same coordinates as the background evolution and no additional λ GBdependent shift in the gauge parameter λ(t) needs to be performed. This is a consequence of demanding regularity of the solution at the event horizon, which, in the Eddington Finkelstein coordinates we use, implies in-falling boundary conditions for the perturbation. Because of this regularity, infinitesimally small shifts in the position of the event horizon do not lead to information loss in the boundary theory. We should also note that, since we are consistently linearising in λ GB , we do not have to specify its value at any step in the algorithm, since this parameter appears outside of δg (Eq. (3.1)), and hence the solution In order to impose the boundary conditions, we record the asymptotics satisfied by the linearised fields, which directly follow from Eq. (2.23) after choosing a (4) = a (4) for the perturbations, we find it convenient to introduce the redefinitions which resemble those of Eq. (2.11). Notice that, as in Eq. (2.11), δσ and δḃ are not the time derivatives of δσ and δb, but independent fields. The boundary conditions as u → 0 for the redefined fields are (3.7) Since we want to study the change in isotropisation between the same state at different values of λ GB , we impose boundary conditions for the perturbations such that the energy density and initial pressure anisotropy of the Gauss-Bonnet solution coincide with that of the λ GB = 0 solution. From the expression for the stress tensor, equating the energy densities at zero and non-zero λ GB , taking into account Eq. (3.5) and linearising in λ GB , we arrive at Similarly, up to linear order in λ GB , the pressure anisotropy is given by where we have defined the perturbed pressure anisotropy The requirement that the initial expectation value of the stress tensor is the same independently of the value of λ GB does not completely fix the state since it only constraints its near-boundary behaviour. However, our approach requires to specify the perturbation δb(u, t = 0) for all u. The simplest possible initial choice consistent with the former requirement is that Eq. (3.10) vanishes at all u, i.e.
All the numerical results presented in this Section will be performed with this initial condition. In Sec. 4 we will come back to this point and discover that for a large class of initial perturbations, the isotropisation dynamics depend very weakly on the particular choice of initial condition for δb. The specification of the boundary condition (3.8) and the initial condition (3.11) determines the linearised solutions uniquely. Implementing the algorithm described above for a variety of choices of background initial conditions b 0,init , we obtain the numerical results depicted in Fig. 1 (all our numerics are produced with a (4) 0 = −1/2). The 4 rows in the figure correspond to 4 different families of initial conditions. Each of the families was generated by specifying a given functional form of the initial condition F(u) (as described in Sec. 2), and different members of a family correspond to choosing different amplitudes such that for values of A ranging from 1 to 20.
In the left column of Fig. 1 we show the evolution for the background pressure anisotropy ∆p 0 , scaled by the amplitude A. The functional forms of b 0,init (u) are shown in the insets in each of the panels. We can see that, in agreement with the results of [17,18], this evolution is, to a remarkably good approximation, linear for the whole range of amplitudes, as evident from the fact that all the curves ∆p 0 /A lie approximately on top of each other. This is the case even for large anisotropies where one would be expect to have a non-linear evolution.
In the right column of Fig. 1 we show the perturbed anisotropy δ (∆p) for the same family of solutions studied in the left panels. Unlike the λ GB = 0 case, the time evolution does not scale with A for all values of A. While for moderate values of A ∼ 2, the perturbed anisotropy does exhibit a clear linear behaviour, the more extreme values of A ∼ 20 show deviations from linearity. Therefore, while the isotropisation dynamics at λ GB = 0 are effectively linear for the whole range of initial conditions we have explored, the higher derivative corrections only exhibit the approximately linear behaviour for large, but not arbitrarily large, anisotropies. Note that the apparent linear behaviour still occurs in a regime of anisotropies where full non-linear evolution would be expected, since the initial conditions that exhibit scaling still have rather large values of anisotropy δ (∆p) /p eq > 1.
In Fig. 2 we show the direct effect of higher curvature corrections on the isotropisation process for the same family of initial conditions studied in Fig. 1. In the left column, we show the full pressure anisotropy ∆p/p eq assuming a fixed value of λ GB = −0.1 (solid lines) in comparison with the N = 4 SYM case (dashed lines) for different families of initial conditions. As in Fig. 1, the anisotropy is scaled by A for each configuration. As we can see, for all families, the effect of the higher curvature corrections is to approximately shift the pressure anisotropy profile towards later times. We can quantify the approximate shift between ∆p and ∆p 0 by comparing the time derivative of the background anisotropy, ∂ t ∆p 0 (t), with the perturbed anisotropy δ (∆p) (t), which we do in the right column of Fig. 2. To leading order in λ GB , the full evolution may be understood as a shift of the background anisotropy, ∆p(t) = ∆p 0 (t + λ GB ∆t(t)) with ∆t a slowly varying function of t, as long as (3.13) The inspection of Fig. 2 shows that Eq. (3.13) holds for almost the entire time evolution, at least as long as the anisotropy is not too large. Note that this statement is independent of value of λ GB and hence it holds for all (small) values of λ GB . This observations leads us to conclude that for all initial conditions we have studied, the effect of the higher curvature corrections is to delay or advance the isotropisation, depending solely on the sign of λ GB .
In the next Section we will show that for those initial conditions in which the full time evolution is approximately linear, this correlation between the signs of λ GB and the shift always holds, irrespective of the initial conditions.

Analysis of the results with QNM expansions
The results from the previous section indicate that, at linear order in λ GB , the evolution of the pressure anisotropy in Gauss-Bonnet gravity behaves in an approximately linear fashion, similarly to the case of λ GB = 0 described in [17,18], although for a smaller range of initial pressure anisotropies. This suggests that the QNM of the final state can capture the full time evolution well, even at early times. In this section we check this hypothesis and conclude that, indeed, for amplitudes for which the linear behaviour in δ (∆p) is observed, its time evolution is well-described by the QNM's of the Gauss-Bonnet black hole. Therefore, in the regime in which the QNM description is applicable, the dynamics of the system are largely simplified and, to a good approximation, the results can be obtained without resorting to full numerical time evolution. Armed with this simplification, we argue that the introduction of the Gauss-Bonnet term induces a time-shift of ∆p with a definite sign, consistent with our numerical findings in Sec. 3. Our argument applies to configurations which satisfy initial conditions more general than Eq. (3.11), suggesting that the observed time shift in the pressure anisotropy is a generic feature of finite coupling corrections.

QNM in AdS/CFT
Generically, QNM's are linearised fluctuations around a black hole background which satisfy ingoing boundary conditions at the horizon 9 and an appropriate boundary condition in the UV. In the context of AdS/CFT, the interpretation of the slow/fast fall-offs of the fields at infinity as sources/vevs for the dual operators indicates what this boundary condition should be if one is to interpret the QNM as poles of the corresponding 2-point function: we simply set to zero the coefficient which corresponds to the source. The frequency of these excitations is generally complex, and their imaginary part naturally provides a time scale for the decay of the perturbation. The longest time-scale is then controlled by the lowest lying QNM's and these govern the late-time dynamics. See [33] for a review.
In the context of holographic isotropisation, the relevant QNM are the linearised fluctuations of the field B which quantifies the anisotropy of the brane. It is easy to check that at the linear level these fluctuations decouple and can be studied on their own. In Fourier space, the perturbation equation is an ODE which defines an eigenvalue problem for the (complex) frequency ω. The frequency spectrum is fixed once the boundary conditions at 9 Note that the ingoing boundary condition translates simply into regularity in coordinates of the form the horizon and the boundary are provided. As mentioned above, the boundary condition at the horizon is simply regularity there, which translates into B having a regular power series expansion in (r − r H ). On the boundary, we demand the induced metric to be fixed B = 0, since this gives the QNM the interpretation of poles of the stress-tensor correlator. The QNM's relevant for isotropisation have been studied in the case of AdS-Schwarzschild planar black hole in [55], and more recently, their Gauss-Bonnet counterparts have been obtained at finite λ GB [32]. We refer the reader to these references for details of their calculation and the specific values of their frequencies.

Matching of QNM in AdS
Our approach is based on the work of [17,18], where it was shown that the time evolution of the pressure anisotropy in N = 4 SYM, even in highly anisotropic cases, can be well captured by a truncated expansion in QNMs, i.e. by linearising Einstein equations around the final equilibrium state of the AdS-Schwarzschild black hole. The idea is to decompose the solution for b 0 (t, u) in QNM's as follows: where N QNM is the number of QNMs we will use in practice, C

Matching of QNM in Gauss-Bonnet
Following the ideas from the previous section, we would like to decompose the time evolution of b(t, u) in QNMs in the Gauss-Bonnet case as well: where now the lack of (0) subscripts indicates that these quantities are the full expressions in the Gauss-Bonnet case. We will expand those linearly in λ GB : Here we note that, even though we have a term linear in t, this expansion is still valid up to arbitrarily large times, since δb(t, u) is multiplied by an infinitesimally small λ GB . In order to compute δω i and δφ i (u), we have numerically solved the Gauss-Bonnet QNM equations at finite λ GB 11 and evaluated the derivatives of ω i and φ i with respect to λ GB , at λ GB = 0. Care must be taken when computing δφ i and δω i , as the former is a function of a dimensionful quantity, and the latter is a dimensionful quantity itself. Since our theories are finite temperature CFT's, these quantities are typically computed as dimensionless numbers, using the temperature to make them such. However, because of our requirement that the two systems, at zero-and nonzero-λ GB , have the same energy densities (Eq. (3.8)), their final temperatures are not the same. Equating Eq. (2.17) and Eq. (2.31) we find the relation between these two temperatures to be We will use this relation to relate the perturbations of the quasinormal modes at fixed energy density with those at fixed temperature. We start with the explicit expression for δω i . The numerically generated dimensionless frequencyω i = ω i /(π T ) as a function of λ GB , can be expanded around λ GB = 0, yielding: After multiplying this expression by πT , with T given by Eq. (4.5), we obtain Note that with the choice of a 0 in our numeric simulations, T 0 = 1/π. For δφ i , the procedure is analogous. From the numerical computation of the QNM, we extract the variation of the wave function evaluated at the same argument whereũ = πT u. Using Eq. (4.5) and expanding to leading order we obtain where we have explicitly used that in our choice of units T 0 = 1/π. 11 We performed the numerics by discretising the eigenvalue problem using pseudospectral methods with a Chebyshev lattice. This allows us to go up to NQNM = 10.
To determine δC i we impose that at t = 0 Eq. (4.4) matches the initial condition Eq. (3.11) for δb as closely as possible. As in the zeroth order case, we use multilinear regression to determine the best-fit coefficients. It is worth noting that even if we start with an initial condition b 0,init such that only certain zeroth order QNM is excited (and so the zeroth order time evolution of the pressure anisotropy is governed by that mode only), in the time evolution of the perturbations, more than one Gauss-Bonnet QNM will be excited.
Having fixed all the components of Eq. (4.4) for δb and using the QNM decomposition Eq. (4.1) for b 0 , we use Eq. (3.9) for the change in the pressure anisotropy to arrive at (4.10) where the u−derivatives are evaluated at u = 0. Similarly, and for later reference, the background pressure anisotropy is given by In Fig. 3 we compare the time dependence of δ (∆p) as predicted by the QNM fit, Eq. (4.10), with our numerical simulations for the four families of solutions studied in Sec. 3. We see that, as long as the initial conditions are such that the system behaves effectively linearly, Eq. (4.10) provides a good description of the full time evolution. However, as the anisotropy increases to very large values, clear deviations from Eq. (4.10) are observed, consistent with the absence of linear behaviour observed in Fig. 1.

Time shift from QNM
We will now argue that the introduction of the Gauss-Bonnet term induces a time shift in the pressure anisotropy in the sense of Eq. (3.13). Our derivation will be valid for those initial conditions in which the full evolution can be described via QNM, which, as indicated in the previous section, can be valid even for large initial pressure anisotropies, δ∆p/p eq ∼ O(1). We will come back to the validity of the QNM approximation in Sec. 5.
As we have already mentioned, the expansion of the initial conditions in terms of a (finite) set of N modes is performed via a multilinear regression on a given grid {u m }, m = 1, ..., M . In this section we will provide explicit expressions for this procedure and use them to determine a relation between the coefficients δC i and C . This relation will allow us to express the change of the pressure anisotropy δ (∆p) directly in terms of initial conditions, which will be the basis to argue that the relation (3.13) is true. Althought it is possible to find an equivalent formulation in terms of a basis of functions in the u-interval, our analysis in the discretised grid is closer to the actual numerical procedure we employ to approximate the time evolution via QNM, and, for this reason, we will focus on that approach.  3.12)). The quasinormal mode evolution describes the numerical results well for those amplitudes that exhibit effective linear behaviour, as shown in Fig. 1.
For convenience, let us first introduce some notation: where δb init (u) is a particular initial condition for δb at t = 0. Note that the quantities x  With this notation, our multilinear regression problem of approximating b 0,init (u) with (4.1) at t = 0 and δb init (u) with (4.4) at t = 0 is equivalent to determining the coefficients C n and δC n by minimising the square errors on our grid: where we have introduced b The analytical solution to this problem is provided by the normal equation, which gives the following maximum likelihood estimate of the coefficients: where we packaged b where we have rewritten x (m) n as a matrix X mn , with the first index being the row index and the second one being the column one. The inverse in Eq. (4.15) is meant to be the Moore-Penrose pseudoinverse in case X T X is non-invertible.
The first equation in (4.14) provides an explicit expression for the coefficients C (0) i in terms of the background metric function b 0 (0, u); therefore the vector C is fully specified by the initial conditions. We will now use the second equation in (4.14) to express δC as a function of the initial conditions as well. Examining the definition of y in Eq. (4.12) we identify two distinct terms: the first term depends solely on the initial condition for perturbations δb init (u), while the second term depends on the initial conditions for the background anisotropy (through C) and on the modifications of the QNM wave functions, which are a property of the asymptotically late stable state and are common to all perturbations. Therefore, after packaging {δb init (u m )} into an M -vector δb, we reorganise the second equation in (4.14) as δC = ρ · δb +ρ · C , (4.16) with δX mn an (M × 2N )-dimensional matrix with components δx (m) n . Since δX only depends on the leading λ GB perturbation of the QNM wave functions, the matrixρ does not depend on initial conditions. As announced, expression (4.16) shows how coefficients δC are directly determined via initial conditions. However, they do not depend solely on the initial conditions for the background anisotropy, i.e. coefficients C n , but also on the initial conditions for the perturbation of the metric, δb init . As we will show in a moment, the influence of a particular choice of δb init on the pressure anisotropy is small. With explicit expressions (4.14), we may now directly relate the pressure anisotropy to the initial conditions for the metric. Before doing so, let us introduce some more notation to Eqs. (4.10) and (4.11): where the definitions of K n (t) and D n follow directly from Eqs. (4.10) and (4.11). Note that these two sets of quantities are independent of initial conditions and are completely determined by the QNM of the late time solution, and that, since ∆p 0 and δ (∆p) only depend on the near boundary behaviour of the metric function, these variables do not depend on the holographic coordinate u. The final bit of bookkeeping consists of defining: T n = e −Γ (0) n t cos Ω (0) n t Re K n + sin Ω (0) n t Im K n , T N +n = e −Γ (0) n t − cos Ω (0) n t Im K n + sin Ω (0) n t Re K n , δT n = e −Γ (0) n t cos Ω (0) n t Re D n + sin Ω (0) n t Im D n , where n = 1, ..., N and where we have explicitly separated the real and imaginary parts of the QNM frequencies, ω n . Finally, we plug Eqs. (4.14) and (4.16) into the QNM decompositions for the time evolution of the background pressure anisotropy and its perturbation, Eqs. (4.11) and (4.10), and use the notation (4.19) to write them as: where H 0 (t) ≡ δT (t) · ρ , are M -dimensional vectors with values at the grid points given by components H (m) 0 (t) and H (m) (t) 12 , and T (t) and δT (t) are 2N -dimensional vectors with components given in Eq. (4.19). We should point out that both H 0 (t) and H(t) are collections of universal functions of time, independent of the initial conditions. In this way, Eq. (4.20) provides a direct link between initial conditions on the gravity side and the pressure anisotropy on the field theory side.
In Eq. (4.20), the perturbed pressure anisotropy is expressed in terms of the initial conditions for the background metric field b 0 and the perturbation δb via two sets of functions H (m) (t) and H (m) 0 (t), which encode the time evolution of the anisotropy. Examining these two sets of functions numerically, we observe that the magnitude of the components of H(t) is at least one order of magnitude larger than those of H 0 (t). To illustrate this fact, in the right panel of Fig. 4 we show the ratio of norms of these two functions, defined simply as the square sum of the components over grid points, |Z| 2 ≡ Z (m) 2 . As claimed, after a short transient time, the magnitude of H 0 is one order of magnitude smaller than H over the whole time range. This, in turn, means that as long as the magnitudes of initial conditions for δb and b 0 are comparable, the influence of the exact form of the initial perturbation δb init on the final perturbation of the pressure anisotropy δ(∆p) is small.
To show this more explicitly, in the right panel of Fig. 4 we compare the perturbations of the pressure anisotropy (as given by the QNM expansion) evaluated for different choices of δb init . For the background initial condition b 0,init we chose the same one as in the first row of Fig. 1, while the different choices of δb init include our canonical choice Eq. (3.11), and four sets of control initial conditions δb init = 0 and δb init = b 0,rand , where b 0,rand are the background initial conditions b 0,init from the rest of the panels in Fig. 1. As we can see, the evolution of the perturbed pressure anisotropy becomes rather insensitive to the particular choice of δb init already at early times.
The observation above allows us to effectively neglect the second term in the expression for δ(∆p) in Eq. (4.20), and describe the evolution of both the background pressure anisotropy ∆p 0 and its perturbation δ(∆p) as linear combinations of values of only the background initial condition b 0,init . Because of this convenient structure, we are now in position to demonstrate that the striking feature of our numerical analysis in Sec. 3 -that the time evolution in the presence of Gauss-Bonnet terms can be understood as a time shift of the time evolution without it -holds true for any initial condition b 0,init , as long as the resulting time evolution of the pressure anistropy is approximately linear, so that the QNM expansion is valid.
As noted earlier, a necessary condition for the (nonconstant) time shift is that signs of ∆p 0 and δ(∆p) are correlated, Eq. the second term in the expression for δ(∆p)), this will hold true for any initial condition, as long as the signs of H (m) and ∂ t (H For all values of u, after a short transient time these two initial-condition independent sets of functions become very similar, and satisfy Eq. (3.13) for most of the time range. This observation implies that for all those large (but not too large) pressure anisotropies which enjoy effectively linearised dynamics, the effect of the Gauss-Bonnet corrections may be understood as a time-shift of the isotropisation dynamics. As already stated in Sec. 3, although the magnitude of this time-shift depends on the initial conditions, whether such shift is positive or negative is completely determined by the sign of λ GB .

Exploring initial conditions
As we have seen, both from the of numerical simulations considered so far and the analysis of the QNM expansions, the effect of a (negative) λ GB correction leads to a delay in the isotropisation time. In this section, we will inspect how general this observation is by exploring the perturbed dynamics of 460 different random numerical configurations (generated using the procedure explained in Sec. 2), which we present in Fig. 6.
For a large portion of those configurations, we have observed that for times t > t iso , with the isotropisation time t iso defined as the time after which ∆p 0 /p eq < 0.1, the effect of λ GB -corrections can be understood as an approximately constant shift of the dynamics. To see this, we have, following the discussion around Eq. (3.13), fitted the numerical results for the perturbed anisotropy to the time derivative of the background anisotropy times a constant, α: and used linear regression to extract the best-fit value for the constant time shift ∆t = λ GB α ≡ λ GB ∆t. The left plot in Fig. 6 shows the ratio of ∆t/t iso as a function of the maximum background pressure anisotropy for all the numerical configurations we considered. The points are color coded with the value of the coefficient of determination R 2 , a measure of how good the simple linear fit (5.1) is. As we can observe, for a large portion of the simulations, Eq. (5.1) is a rather good fit: 83% of the simulations have R 2 > 0.8, which signals a satisfactory fit. We should note that even the simulations with lower R 2 , irrespective of the value of the maximum anisotropy, still display a time shift, albeit not a constant one. Another striking feature of Fig. 6 is that, apart from a handful of low-R 2 outliers, the vast majority of points are spread over a relatively narrow time interval with ∆t/t iso mean = 0.98 and a standard deviation of σ ∆t/t iso = 0.32.
Another key result of this work is the apparent effective linearisation of the far-fromequilibrium dynamics of the perturbation of the pressure anisotropy for large (but not arbitrarily large) background anisotropies, which admits a description in terms of a QNM expansion. To quantify when the non-linearities become important we have inspected the effect of describing the time evolution of the pressure anisotropy (both the background and the perturbed one) with QNM expansions on our observable of interest, the shift in the isotropisation time. In the right panel of Fig. 6 we show the difference between the relative time shift ∆t num /t iso obtained from the numerics and the one obtained from the corresponding QNM expansions, ∆t QNM /t iso . This plot shows that for values of the maximum background pressure anisotropy ∆p 0 /p eq < 3, the difference between the numerical result and the linearised approximation is less than 20%. For larger values, however, the spread of this difference is larger and, even though there are some sets of initial conditions which can be nicely described by QNM, the quality of the description becomes worse for more general initial conditions. 13 Even though the values of maximum background pressure anisotropy for which the perturbed dynamics is well described by the QNM are generally lower than the values for which the background dynamics is well described, they are still rather high and in the regime in which we may not normally expect the linearised approximation to hold.

Conclusions
We have performed the first analysis of higher curvature effects on far-from-equilibrium anisotropic dynamics in holography. Although it is yet unknown what is the precise gauge theory dual to Gauss-Bonnet gravity (if any), in this first exploratory work we have focused on this model since it provides a simple setup in which to address the effect of quadratic curvature corrections, which are the leading order corrections expected in a generic holographic construction. In the particular case of finite coupling corrections of N = 4, the leading order quadratic corrections identically vanish and the first non-vanishing correction is the quartic one. Nevertheless, as recently pointed out in [32], when the higher derivative 13 We should note that the isotropisation time tiso is a rather sensitive observable, and what may seem like large discrepancies between the numerics and the QNM results in the right plot in Fig. 6, often are not. This was noted already in [18], and has to do with the fact that around times of the order of tiso, the pressure anisotropies become rather oscillatory and small differences between the numerics and the corresponding QNM approximation may sometimes lead to large differences between tiso.
corrections are small, the features of quasinormal modes of Gauss-Bonnet and N = 4 SYM at large but finite coupling are qualitatively similar. Since QNMs control the far-fromequilibrium dynamics of large (although not arbitrarily large) anisotropies, we expect that many of our results apply to the case of finite coupling corrections to N = 4 SYM, at least qualitatively.
One of the main results of this work is the determination of the effect of the Gauss-Bonnet terms on the isotropisation dynamics of far-from-equilibrium plasmas. In all the numerical simulations we have studied in Sec. 5, with pressure anisotropies spanning over several orders of magnitude, the inclusion of (negative) λ GB corrections leads to a λ GBproportional shift of the time evolution of the full pressure anisotropy, for the most part of the evolution in which the anisotropy remains large. This shift then implies a delay in the isotropisation time (see Figs. 2 and 6). In Subsec. 4.4, we have demonstrated that, as long as the description of the perturbed dynamics is effectively linear, this observation holds true for any choice of initial conditions (Fig. 5).
Our finding that, for all the numerical configurations considered, ∆t is always positive, and hence the sign of the physical time shift ∆t is alway opposite to the sign of λ GB nicely resonates with the corrections to transport properties in this theory, since negative λ GB implies larger η/s ratio and vice versa. Quite intuitively, gauge theories with higher viscosities, which may be viewed as less strongly coupled, possess longer isotropisation times. It would be interesting to analyse other examples of large curvature corrections to verify these findings.
Another main result of this work is the apparent effective linearisation of the far-fromequilibrium dynamics at strong coupling for large, but not arbitrarily large anisotropies. Such effective linearisation is one of the most remarkable features of the relaxation of farfrom-equilibrium anisotropic states in holography. This linearisation has been observed in configurations with strong anisotropies [17,18] with or without the magnetic field and chemical potential [19]; in the dynamics of baryon charge in the collisions of shocks [23]; and in non-relativistic holography [56]. Those studies indicate that, at least in the large coupling limit, the dynamics of those far-from-equilibrium processes are tremendously simplified.
In this paper we have performed the first steps towards understanding whether this simplification survives once corrections to those extreme limits are considered. Within the context of Gauss-Bonnet gravity we have found that linearised dynamics provide a good description of the full far-from-equilibrium processes even when the initial anisotropies are large (∆p/p eq ∼ O(1)). This regime is of particular phenomenological relevance, since the initial anisotropies of the far-from-equilibrium dynamics in heavy ion collisions at RHIC and the LHC are of that order. For those anisotropies, we have shown that the dynamics of the system can be predicted by a QNM analysis, much like in the unperturbed λ GB = 0 case [17][18][19]. However, such simplification does not extend to arbitrarily large values of the initial anisotropy, as illustrated in Fig. 1.
From the gravity point of view, the onset of the non-linear behaviour is controlled by the right hand side of Eq. (3.2), which is sensitive to the higher derivatives of the background metric. As a consequence, when the pressure anisotropy is pushed to very large values, the corrections induced by the higher curvature terms exhibit clear non-linearities, even though the unperturbed solutions still behave surprisingly linear. On the gauge theory since, our results indicate that the dynamics of gauge theories at finite but large coupling become more non-linear than in the infinitely strongly coupled limit.
A clear limitation of our approach is the lack of control over the small λ GB corrections to the initial out-of-equilibrium state. As we have discussed, in our setup we compare states that, independently of the value of λ GB , have the same initial anisotropy in units of the background energy density. However, this is not enough to fully specify the state. On the gravity side this ambiguity translates into the choice of initial δb(u). One possible way to constrain this function is to demand the higher point correlators of certain operators are identical; however, to fully constrain the initial data, one would need an infinite set of those. Another possibility is to prepare the initial state as the result of a rapid quench of some anisotropic deformation of the theory independently of the value of λ GB . This procedure fixes the initial state uniquely. Nevertheless, our study of different functional forms of the initial conditions shows that our results remain independent of the details of the initial disturbances, provided that these are comparable to the initial anisotropy profile b 0 . For anisotropies such that the system behaves effectively linearly, the independence on the initial disturbance can be expressed by the smallness of the initial condition independent function H 0 (t) as compared to H(t), see Subsec. 4.4. For these reasons, we expect our main conclusions to remain unchanged for a generic far-from-equilibrium excitation of the system.
The method we have employed in this paper is easily extendable to other sets of higher curvature corrections, as long as these corrections are small. Such small higher-curvature corrections will then lead to small deviations on top of the non-linear relaxation of the initial far-from-equilibrium state without corrections. In particular, the left hand side of Eq. (3.2) is common to all gravities with small higher derivative corrections, since it only depends on the details of the time evolution of the far-from-equilibrium background we are perturbing. The difference between the various sets of corrections resides in the explicit form of the right hand side of Eq. (3.2), which may be understood as a source for the equation of motion for small perturbations on top of the background solution. This source term depends on different curvature tensors of the dynamical background. It would be interesting to generalise our study to other sets of corrections, such as the quartic curvature terms dual to finite coupling corrections of N = 4 SYM. Note, however, that the higher the order of curvature corrections is, the higher the derivatives of the background contained in the source term are. Therefore, to implement the quartic curvature terms, the dynamical background should be determined to a higher precision than what our current implementation allows. For this reason, we postpone this interesting analysis for future work.
Note added: As this paper was being finalised, we became aware of a related work by S. Grozdanov and W. van der Schee in which they analyse the effect of the Gauss-Bonnet term on the far-from-equilibrium dynamics of the debris of two shock-wave collisions in holography.