Evolution of holographic entanglement entropy in an anisotropic system
Abstract
We determine holographically 2point correlators of gauge invariant operators with large conformal weights and entanglement entropy of strips for a timedependent anisotropic 5dimensional asymptotically antide Sitter spacetime. At the early stage of evolution where geodesics and extremal surfaces can extend beyond the apparent horizon all observables vary substantially from their thermal value, but thermalize rapidly. At late times we recover quasinormal ringing of correlators and holographic entanglement entropy around their thermal values, as expected on general grounds. We check the behaviour of holographic entanglement entropy and correlators as function of the separation length of the strip and find agreement with the exact expressions derived in the small and large temperature limits.
Wiedner Hauptstr. 810, A1040 Vienna, Austria
List of Figures
 I Anisotropy function and pressures.
 II Geodesics and their relative position to event and apparent horizon.
 III Renormalized length of geodesics.
 IV Comparison of longitudinal and transverse geodesic lengths.
 V Longitudinal and transverse HEE for different separations.
 VI Longitudinal and transverse HEE for same separation.
 VII Correlators and HEE multiplied by lowest QNM.
 VIII Late time correlators and HEE as function of separation.
 IX Comparison between AdS and conformal geodesics.
1 Introduction
Far from equilibrium dynamics and thermalization of strongly coupled systems has attracted a lot of attention in the past decade. The reason for this is twofold. On the one hand experiments at RHIC and the LHC revealed that the quark gluon plasma created in heavy ion collisions behaves as a strongly coupled liquid that thermalizes extremely fast Teaney:2000cw (); Huovinen:2001cy (); Hirano:2002ds (); Romatschke:2007mq (). On the other hand, in condensed matter experiments it is now possible to drive an isolated system to a far from equilibrium state by a quantum quench, i.e. a control parameter of the system is varied rapidly Polkovnikov:2010yn ().
One powerful tool to study strongly coupled systems out of equilibrium is the gauge gravity duality where classical supergravity in dimensions is dual to a dimensional strongly coupled large field theory that loosely speaking lives on the boundary of the gravity theory Maldacena:1997re (); Witten:1998qj (); Gubser:1998bc (). The thermalization process of field theories with a conformal field theory fixed point in the ultraviolet is then mapped to black hole (or black brane) formation in asymptotically Antide Sitter (AdS) space Witten:1998zw ().
Particularly useful theoretical observables to monitor the thermalization process are the vacuum expectation value of the stressenergy tensor, correlation functions of gauge invariant operators and entanglement entropy (EE). The former two are easy to define on the field theory side and, given some standard assumptions that we recall below, easy to calculate on the gravity side. We recall now relevant aspects of EE, see Eisert:2008ur () for more details and references.
On the field theory side EE is defined in the following way. Dividing a system into two subsystems and the EE of the subsystem is defined as the von Neumann entropy of the reduced density matrix obtained by tracing out the degrees of freedom of the subsystem .
On the gravity side, the holographic entanglement entropy (HEE) is defined as the area of a minimal surface extending from some predefined surface on the boundary into the bulk Ryu:2006bv (); Nishioka:2009un (). For timedependent backgrounds the minimal surface has to be replaced by an extremal surface Hubeny:2007xt (). It is noteworthy — both for conceptual and for technical reasons — that the relevant surfaces used to determine HEE can extend beyond the apparent horizon. Conceptually, this implies that HEE provides a ‘classical’ way (on the gravity side) to extract (quantum) information from the region beyond the horizon. Technically, this requires the numerical determination of (part of) the spacetime region beyond the apparent horizon, which can be a challenge since that region ends in a singularity.
EE is notoriously hard to compute in quantum field theories; so far this is only possible in highly symmetric theories such as (relativistic) conformal field theories Calabrese:2004eu (); Holzhey:1994we (); Vidal:2002rm () or Galilean conformal field theories Bagchi:2014iea (). In a seminal work, Calabrese and Cardy Calabrese:2005in () were able to compute the time evolution of the EE after a quench in a two dimensional conformal field theory and in the Ising spin chain model in a transverse magnetic field. In both cases they find that for an entangling interval of length , the EE increases linearly with time until after which it saturates. The linear scaling with time and the crossover at can be understood in terms of entangled quasiparticles pairs emitted from the initial state and is therefore expected to hold for a wider class of systems.
The holographic duality offers a playground where one can study various different systems and search for universal and novel properties of HEE in dynamical situations.
The simplest example where one can study the analog of quenches in the holographic setup are spacetimes where thin shells collapse to form a black hole, first utilised in AbajoArrastia:2010yt (). The behaviour of the EE in these setups has been studied extensively AbajoArrastia:2010yt (); Balasubramanian:2011ur (); Albash:2010mv (); Baron:2012fv (); Galante:2012pv (); Keranen:2011xs (); Liu:2013qca () and indeed shows universal behaviour consistent with the findings of Cardy and Calabrese. Namely the initial short early time epoch is followed by a long period where the EE grows linearly with time and is independent of the entangling boundary region. This was first worked out in the Vaidya spacetime AbajoArrastia:2010yt (); Liu:2013iza (); Liu:2013qca () where the shell is composed out of null dust and later generalized to matter with arbitrary equation of states Keranen:2014zoa (); Keranen:2015fqa (). The linear scaling is even present in geometries with Lifshitz scaling and hyper scaling violation Alishahiha:2014cwa (); Fonda:2014ula (). In addition, in the above works the HEE is a monotonically increasing function that saturates to its equilibrium value from below. Note that these quenches are thermal quenches because the end state is a thermal state.
Another way of studying thermal quenches is by turning on a radially collapsing scalar field that forms a black hole. This situation is more complicated because in order to obtain the geometry from which one can then obtain the HEE, Einsteins equations have to be solved numerically Bizon:2011gg (); AbajoArrastia:2014fma (); Buchel:2014gta (). In AbajoArrastia:2014fma (); daSilva:2014zva () this was done for a radial collapsing massless scalar field in global AdS where the scalar field can have many bounces between the boundary and the center of AdS before a black hole forms resulting in a periodic behaviour of the HEE.
In Buchel:2014gta () a massive scalar field dual to a massive fermionic operator was turned on, treating the quench as a perturbation on the static spacetime. In this setup the HEE is also not monotonic and in some cases approaches the equilibrium value from above. This reveals a qualitative difference of EE to the thermal entropy which, on general grounds, must be monotonically growing in a closed system.
An additional motivation to study HEE comes from the question how to measure entropy production in (holographic models for) heavy ion collisions Muller:2011ra (). Within the gauge/gravity duality the entropy of the (stationary) black hole corresponds to the entropy of the field theory. However, in timedependent backgrounds entropy as defined from the area of the apparent horizon is ambiguous because it depends on the choice of time slicing. By contrast the definition of the HEE is unique and therefore may serve as an alternative measure for entropy production.
So far nearly all studies of HEE have relied on the simplifying assumption of a spherically symmetric spacetime (see Narayan:2012ks () for a notable exception). In the paper at hand we drop this assumption and investigate the effect of anisotropies on the evolution of the HEE.
The paper is organized at follows. In section 2 we introduce the anisotropic spacetime we will use, recall general aspects of 2point correlators and simplify the determination of HEE to a geodesic problem in an auxiliary spacetime. In section 3 we discuss the numerical strategy that we used, relegating details to the appendices. In section 4 we present results for the background geometry, holographic stress tensor, 2point correlators and HEE. In section 5 we analyze the late time behaviour of correlators and HEE and relate it to the lowest lying quasinormal mode. In section 6 we conclude with a brief summary of our results and some possible next steps. In appendix A we comment on the spectral method and other numerical routines used to solve the 5dimensional vacuum Einstein equations. In appendix B we explain the relaxation code used to determine geodesics and extremal surfaces.
Before starting we mention some of our conventions. We use mostly signature and set the AdS radius and final black brane mass to unity.
2 Theoretical setup
2.1 Anisotropic asymptotically AdS spacetimes
In this section we review the most important details of the model first introduced in Chesler:2008hg () and studied further in Chesler:2009cy (); Heller:2012km (); Heller:2013oxa (); Chesler:2013lia ().
The 5dimensional bulk metric that introduces anisotropy between the longitudinal and transverse directions with an rotational invariance in the transverse plane can be written conveniently in Eddington–Finkelstein coordinates,
(1) 
where the functions and only depend on the holographic coordinate and (advanced) time . In this coordinate system the vacuum Einstein equations read
(2a)  
(2b)  
(2c)  
(2d)  
(2e) 
where prime denotes radial derivative and dot time derivative, viz.
(3) 
for any function .
The Einstein equations (2) have to be solved for special initial conditions and appropriate boundary conditions. There are, at least, two ways to create a far from equilibrium state. On the one hand one can turn on a time dependent anisotropy function at the boundary as in the original works of Chesler:2008hg (); Chesler:2009cy () and let the system evolve. In this case the boundary metric is curved and the conformal anomaly is present Henningson:1998gx (). On the other hand one can specify the initial state in the absence of external sources by specifying the metric in the bulk on the initial time slice Heller:2012km () with a flat boundary geometry. For simplicity, in the following we will study the setup where the boundary metric is flat and time independent.
Requiring that the spacetime is asymptotically AdS, Einstein’s equations in the near boundary expansion are solved by^{1}^{1}1We fix a residual gauge freedom that preserves the form (1), , by demanding that the first subleading term in the function falls off like . Note that in several numerical simulations in the literature the same freedom is fixed by placing the apparent horizon at a specific value of the radial coordinate .
(4a)  
(4b)  
(4c) 
The coefficients in the asymptotic expansion determine the expectation value of the stress energy tensor in the dual field theory deHaro:2000xn ()
(5) 
where
(6) 
In order to determine the fourth order coefficients one needs to solve Einstein’s equations numerically for some initial conditions. When doing the actual calculation we work with the inverse radial coordinate so that the boundary is located at .
For our initial data we follow Heller:2013oxa (); Chesler:2013lia () and choose for the anisotropy function on the initial time slice
(7) 
with , and . In addition the initial conditions have to be supplemented with a value for the coefficient which sets the energy density of the initial state, for which we take , corresponding to an equilibrium temperature , as we now recall.
At late times we expect isotropization, . In that case we recover the usual static AdS black brane solution as follows. Solving (2e) and using residual gauge transformations yields . This implies and . Solving (2a) then yields , where we fixed the integration constant such that . [The other equations are either trivial, (2b) and (2d), or redundant, (2c).] The result for is the usual Killing norm for the static AdS black brane. Surface gravity is given by so that the Hawking temperature is .
In the generic anisotropic case, , we solve the Einstein equations (2) numerically for the initial conditions (7). In this background we then study the evolution of 2point correlation functions for operators of large conformal weights and the HEE. This in turn requires us to determine the background sufficiently far beyond the apparent horizon. In section 3.1 below we discuss the numerical implementation.
2.2 2point correlators
The equal time 2point function for an operator of large conformal weight can be computed via a path integral as Balasubramanian:1999zv (); Festuccia:2005pi ()
(8) 
where the integral is a sum over all possible paths with endpoints at and and is the proper length of the path. The first approximation neglects perturbative corrections and is the so called geodesic approximation, which holds in the limit when the conformal weight of the operator is large. The conformal weight effectively plays the role of in usual perturbative expansions of path integrals. Then it can be shown that the sum over all paths reduces to a sum over all geodesics where denotes the length of the corresponding geodesic. To leading order only the geodesic with the smallest value of contributes, whose length we denote by , which explains the second approximation^{2}^{2}2 For a comparison of the 2point correlation function obtained by using the “extrapolate” dictionary and the geodesic approximation in AdS Vaidya spacetime see Keranen:2014lna ().. It neglects instanton corrections.
However, the length of the geodesic has a divergence originating from the asymptotically AdS boundary and therefore needs to be renormalized. We choose to subtract the length of a geodesic in the static black brane background, which we denote by . In terms of the renormalized length the 2point function becomes
(9) 
This means that we can obtain the time evolution of 2point functions by looking at spacelike geodesics that are anchored at the boundary at fixed separation and calculating their length at different times. Due to the anisotropy in the system we only solve for the subset of correlation functions that are either separated in the longitudinal direction or in the transverse directions.
To this end we let all the coordinates depend on one parameter , which lies in the interval . To obtain the lengths of the geodesics we have to solve the geodesic equation for the two subspaces given by the line elements
(10)  
(11) 
For the separation in the transverse direction the geodesics end at , where is the boundary time. Similarly, for the longitudinal separation we take . With this choice of boundary conditions the lengths of the geodesics in the background (1) are given by
(12)  
(13) 
where prime denotes the derivative with respect to .
It is important to point out that we can only study geodesics after some advanced time with boundary separations below a maximal separation . This comes from the fact that by solving Einstein’s equations numerically we have to choose a finite computational domain. Also, by specifying the initial state in the entire bulk on the initial time slice the advanced time interval at our disposal is . As the geodesics reach into the bulk they bend back in advanced time leaving the computational domain for advanced times as well as extending too far into the bulk for separations .
In section 3.2 below we discuss how to solve the geodesic equations numerically.
2.3 Holographic entanglement entropy
In time dependent systems the covariant HEE Hubeny:2007xt () for some boundary region is obtained by extremizing the 3surface functional
(14) 
that ends on the boundary surface . In the dual field theory the EE is then conjectured to be given by Ryu:2006bv (); Ryu:2006ef (); Hubeny:2007xt ()
(15) 
Usually the boundary regions of interest are either a sphere or a strip that has finite extent in one direction and infinite extent in the other two directions. In spacetimes with spherical symmetry in the three spatial dimensions the problem of finding the extremal area functional (14) effectively reduces to finding geodesics. In our case where spherical symmetry is broken this is not the case anymore. For example, finding the extremal area for a spherical boundary region would require to solve nonlinear coupled partial differential equations. However, in the case of a strip with finite extent either in the transverse or longitudinal direction it is possible to reduce the problem to finding geodesics in a suitable auxiliary spacetime, as we now demonstrate.
We introduce two scalar fields and write the line element as
(16) 
where is a 3dimensional metric with coordinates where represents the coordinate we choose to have finite spatial extent, i.e. either or one of . The remaining (noncompact) coordinates are then denoted by , , which we choose to be two of our three worldvolume coordinates; the third one is denoted by . Parametrizing the 3dimensional coordinates as , the area functional (14) can be written as
(17) 
Performing the integration over the Killing coordinates and yields a (possibly infinite) constant volume factor through which we are going to divide. Thus, instead of calculating HEE we calculate a HEE density per Killing volume. The problem of extremizing the 3surface corresponding to a boundary region of striptopology is then reduced to a 1dimensional problem.
In fact, from the expression (17) on can see that the problem of finding the extremal 3surfaces reduces to finding geodesics of the conformal metric
(18) 
The 3dimensional conformal metrics for separation in the transverse and longitudinal directions for which we have to solve the geodesic equation in our case are given by
(19a)  
(19b) 
3 Numerical implementation
3.1 Einstein equations
We solve the Einstein equations (2) using pseudo spectral methods as described in detail in Chesler:2013lia (), with the only difference that we do not fix the location of the apparent horizon, see appendix A for some details. Not fixing the apparent horizon facilitates the study of geodesics and extremal surfaces that reach behind the apparent horizon, which is of relevance for 2point functions and HEE. For that reason we want a large computational domain in the holographic coordinate . In all the computations we took with the final position of the horizon located at . For the time evolution it is sufficient to use a fourth order Runge Kutta method with time steps . All the computations were done with the open source software GNU Octave octave:2014 ().
3.2 Geodesics
To compute 2point functions we need to find curves of extremal length in a curved spacetime whose endpoints reside on fixed positions on the boundary of that spacetime. These curves are solutions to the geodesic equation subject to boundary conditions at the endpoints. For numerical reasons it turns out to be convenient to use a nonaffine parameter , where is the usual affine parametrization, . In terms of the geodesic equation reads
(20) 
where and denotes the Jacobian which originates from the change in parametrization. This form of the geodesic equation gives us the freedom to choose parametrizations resulting in better convergence behaviour of the relaxation algorithm than the affine parametrization does. In physical terms the right hand side in (20) introduces a fictitious viscous force that enhances numerical convergence.
In our case the geodesic equation is given by a set of three coupled nonlinear ODEs of second order for the geodesic coordinates . This set of equations can be reduced to a set of six first order equations in terms of the geodesic coordinates and their first derivatives:
(21a)  
(21b)  
(21c)  
(21d)  
(21e)  
(21f) 
This set of equations is a twopoint boundary value problem, which is usually either solved with shooting methods or relaxation methods Press:2007:NRE:1403886 (). We do not shoot but relax. The geodesics obtained in appendix B serve as the initial guess for the relaxation algorithm. Since we are interested in oneparameter families of geodesics (evaluated at different constant time slices) we can take the solution for the family member as initial guess for the family member. More details of our implementation of the relaxation method are described in appendix B.
3.3 Extremal surfaces
The same method as above works also for HEE. Namely, for our problem at hand the evaluation of extremal surfaces is reduced to the evaluation of geodesics in some auxiliary spacetime, as we have shown in section 2.3.
4 Results
In this section we display and discuss our main results. In all figures where a timeaxis is plotted we measure the boundary time or the bulk advanced time in units of the temperature of the final black brane, . The separation length of 2point functions and HEE and the corresponding boundary coordinates are given in units of as well. To make the approach to thermal equilibrium most transparent we use normalized quantities for the geodesic length and HEE defined by
(22a)  
(22b) 
where () is the unrenormalized length (HEE) and () is the corresponding thermal value.
4.1 Background geometry and holographic stress tensor
Figure I displays the most salient features of the background geometry. The left figure plots the anisotropy function and displays the regions outside and inside the apparent horizon, as well as the event horizon. The black lines depict a null congruence of geodesics close to the event horizon to exhibit their ingoing/outgoing nature. The right figure plots transversal and longitudinal pressures as function of boundary time . Note the quick thermalization of the pressure components.
4.2 2point correlators
As we have noted before geodesics can extend beyond the apparent horizon. This is made explicit in Fig. II where the blue curve serves as our initial guess for the relaxation code. The red geodesics at late times approach the apparent horizon without crossing it. At sufficiently early times (and sufficiently large separation) the geodesics cross the apparent horizon, an example of which is depicted by the cyan curve.
The evolution of the renormalized lengths in the transverse and longitudinal directions for different separations are depicted in Fig. III. Depending on the separation the 2point functions start at which is the time when the geodesics extend beyond the computational domain.
The first observation is that the transverse and longitudinal directions oscillate out of phase as shown in Fig. IV. The same feature is seen in the transverse and longitudinal pressure. By comparing the thermalization times of the one point functions, i.e. the expectation value of the stress energy tensor with the 2point functions we see that the 2point functions thermalize later as expected. Also, the thermalization time increases if the boundary separation is increased.
4.3 Holographic entanglement entropy
The extremal surface equations — which we mapped to geodesic equations in an auxiliary spacetime — are solved again by a relaxation method. We observe the same qualitative features as for geodesics in Fig. II above: at early times extremal surfaces can extend beyond the apparent horizon, while at sufficiently late times they approach it from the outside without crossing. However, there are also notable differences to geodesics, which we discuss now.
As can be seen from Fig. IX in appendix B.2 conformal geodesics reach much further into the bulk compared to the pure AdS case. Therefore the boundary separations we can study for the HEE are smaller compared to the 2point functions. This is also the reason why for the same boundary separation the HEE reaches equilibrium later as the 2point functions. For the same boundary separation and at the same boundary time conformal geodesics reach deeper into the bulk and further back in time and therefore are more sensitive to out of equilibrium effects which are most pronounced at early times. In addition the shape of the curves differ from the 2point functions with the oscillations less pronounced. We exhibit these features now in some plots.
5 Late time behaviour and quasinormal modes
After the early far from equilibrium phase the geometry relaxes to the static Schwarzschild black brane solution. As noted in Heller:2013oxa (); Chesler:2013lia () the anisotropy of the system is exponentially damped and at sufficiently late times one enters the linearized regime. In this regime the approach to equilibrium is accurately described by the lowest lying quasinormal mode (QNM) which characterizes the response of the system to infinitesimal metric perturbation. In the case at hand the relevant channel for the gravitational fluctuations is the spin two symmetry channel which coincides with the fluctuations of a massless scalar field in the static black brane geometry. The asymptotic response of the pressure anisotropy then takes the form
(23) 
with the lowest QNM given by Starinets:2002br (); Kovtun:2005ev ()
(24) 
On the field theory side QNMs appear as poles in the retarded Green function Birmingham:2001pj (); Son:2002sd (); Kovtun:2005ev (); Berti:2009kk (). It is therefore expected that also the late time behaviour of the correlation functions obtained in the previous section is described by the lowest QNM. We now show that this is indeed the case.
In figure VII (left) we plot the renormalized geodesic length multiplied with the imaginary part of the lowest QNM for transverse and longitudinal separations. One clearly sees that after a short period of time the evolution of the correlator is accurately described by the ringdown of the black brane with constant amplitude and frequency. The connection between the late time behaviour of correlation functions and QNMs was previously also observed in Balasubramanian:2012tu (); Ishii:2015gia (); David:2015xqa ().
It turns out that HEE also follows this pattern. In Bhattacharya:2013bna () a connection between QNMs and the behaviour of the HEE was found. From linearised Einstein equations one can derive a differential equation for the first order correction of the HEE describing its change when a given ground state is excited. By imposing infalling boundary conditions at the horizon one obtains a QNM dispersion relation putting a constraint on HEE. With our numerical solution we can demonstrate that the late time behaviour of HEE indeed follows the QNM ringdown even without imposing infalling boundary conditions. In Fig. VII (right) we show the HEE multiplied with for the infinite strip with finite separation in longitudinal and transverse direction. As for the correlation function, at late times, the HEE shows quasinormal ringing with constant amplitude and frequency. These oscillations show that HEE must not approach its thermal value from below but rather shows oscillatory behaviour around its thermal value.
To conclude this section we finally study the departure of the length of the geodesics and HEE from equilibrium for different times as a function of the boundary separation. This time we normalized the length of geodesics by subtracting a cutoff dependent piece. The case for the 2point function with separation in the transverse direction is displayed in Fig. VIII. Out of equilibrium effects manifest themselves as oscillations around the thermal value. The curves terminate when the geodesics leave the computational domain, so for early times we only have access to rather small boundary separations. The same effect is seen for HEE.
In the thermal limit the scaling for the geodesic length at small and large boundary separation is dictated by conformal symmetry and is proportional to and respectively. At large separation the HEE also scales linearly with the separation length, whereas at small separation it is proportional to . All these results agree precisely with the perturbative expressions derived in the limits of small and large temperatures Ryu:2006ef (); Fischler:2012ca (). Our numerical results are shown in Fig. VIII where we also plot the corresponding zero temperature results, which coincide with the thermal curves for small separations.
6 Conclusions
In the paper at hand we have studied the time evolution of 2point correlation functions in the geodesic approximation and holographic entanglement entropy in an anisotropic SYM plasma.
To obtain the 2point correlation function for separation in transverse and longitudinal direction we solved for the length of the geodesics in the anisotropic background (1) by using a relaxation method with the zero temperature geodesic as an initial guess. At the th step the solution of the st step is used as initial guess.
Choosing an infinite strip with finite separation either in the transverse or longitudinal direction and using the symmetries of the system the problem of finding minimal surfaces was reduced to finding geodesics in an auxiliary conformally related spacetime and the same strategy as for the 2point correlation function goes through.
The 2point correlation functions as well as holographic entanglement entropy show oscillatory behaviour around their thermal value, where transverse and longitudinal directions oscillate out of phase. Both quantities approach their equilibrium value by exponentially damped oscillations and after the initial far from equilibrium epoch the thermalization process is accurately described by the lowest quasinormal mode. That the late time behaviour of holographic entanglement entropy is captured by the ring down of the black brane geometry is one of the main results of this work and to our knowledge the first explicit verification of the connection between QNMs and holographic entanglement entropy found in Bhattacharya:2013bna ().
The methods developed and applied in this work can be generalized to other timedependent asymptotically antide Sitter backgrounds of interest, such as colliding shockwave backgrounds Grumiller:2008va (); Gubser:2008pc (), which have been constructed numerically in the past few years Chesler:2010bi (); Wu:2011yd (); CasalderreySolana:2013aba (). Another interesting generalization would be the consideration of compact entangling regions, e.g. with spherical topology to check the (in)dependence of our results on the form of the entangling region.
Acknowledgments
We thank Max Attems, Arjun Bagchi, Rudranil Basu, Paul Chesler, Jan de Boer, Ville Keränen, Esperanza Lopez, Masahiro Nozaki, Florian Preis, Max Riegler, Paul Romatschke and Wilke van der Schee for discussions.
This work was supported by the following projects of the Austrian Science Fund (FWF): Y435N16, I952N16, P27182N27, DKW1252N27 and P26328.
Appendix A Spectral method
In this appendix we give some details on how we numerically solve the Einstein equations (2). We work with the inverse radial coordinate such that the boundary is located at . Due to asymptotic AdS boundary conditions, the metric functions and diverge as . It is numerically favourable to define new functions with the known divergent pieces removed and rescale them with appropriate powers of so that the resulting functions are finite or vanish as ; the precise boundary conditions on the new functions are presented in (30) below. This leads us to the following field redefinitions
(25a)  
(25b)  
(25c) 
The redefined anisotropy function allows to extract simply from the boundary value of
(26) 
where . In terms of the redefined fields the first four Einstein equations (2) can be rewritten in the form
(27a)  
(27b)  
(27c)  
(27d) 
with the source functions given by
(28a)  
(28b)  
(28c)  
(28d) 
The relation between dotderivative and timederivative, originally given in (3), turns into
(29) 
The boundary conditions for the redefined fields read
(30a)  
(30b)  
(30c)  
(30d) 
where we set . On the initial time slice the boundary condition (30c) is computed from the initial conditions
(31) 
where we set , and .
For a given set of initial and boundary data the system of equations (27) allows for the following solution strategy:

Finally we integrate (29) to get the new on the next time slice and repeat the whole procedure all over again.
For constant we solve each of the equations (27) with a pseudospectral method boyd2001chebyshev () where the grid points on an intervall are located at
(32) 
In our computations we choose and . These points can be used to construct the entries of the spectral differentiation matrix Trefethen:2000:SMM:343368 ():
(33)  
(34)  
(35) 
where and otherwise . The derivative of a function on the spectral grid points is obtained by multiplication with this diffentiation matrix
(36) 
This allows to turn each of the equations (27) into a system of linear equations. For example the second equation in (27) translates to
(37) 
where the matrix is given by
(38) 
and the source vector by
(39) 
The boundary condition is implemented by setting and . The solution vector is then obtained by multiplying the inverse of with the source vector
(40) 
We do this with the OCTAVE routine linsolve. The equations for , and can be solved in the same way.
To advance the solution for to the next time slice it is sufficient to use a simple fourth order Runge Kutta method Press:2007:NRE:1403886 ()
(41) 
where the koefficients are given by
(42)  
(43)  
(44)  
(45) 
and is computed from (29). In our simulations we use a stepsize of which is sufficient to achieve a stable time evolution.
Appendix B Relaxation method
In relaxation methods differential equations are replaced by approximate finite difference equations (FDEs) on a discrete set of points. The solution is determined by starting with an inital guess and improving it iteratively. In this iterative procedure the result is said to relax to the true solution.
First we define a grid of equidistant spacing with . We use a grid with points in all our simulations. The upper and lower bound of this grid are given by and respectively, where denotes a UV cutoff. The discretized version of our trial solution on this grid is written as
(46) 
The explicit form of the trial solutions used in the computations of 2point functions and HEE are given in appendix B.1 and appendix B.2, respectively. In all our simulations we evaluate the cutoff such that the position of the corresponding cutoff surface is fixed at , i.e. at . The boundary conditions are imposed at this cutoff surface
(47)  
(48)  
(49) 
where denotes the time and the separation on the cutoff surface.
The finite difference representation of the geodesic equation (21) is given by
(50a)  
(50b)  
(50c)  
(50d)  
(50e)  
(50f) 
where is the residual at point in equation ; further denotes the first derivatives of the geodesic coordinates; quantitities with bar are averaged like ; the Christoffel symbols are evaluated from the averaged metric functions; the explicit form of the Jacobian for the case of 2point functions and HEE are given in appendix B.1 and appendix B.2.
The initial guess will in general not satisfy these FDEs very well, i.e. the residua will be rather large. To quantify the deviation of a given trial solution to the true solution we use the following measure
(51) 
The strategy is to compute increments for the geodesic coordinates such that is an improved approximation to the previous solution . This we do iteratively until we reach
(52) 
Equations for the increments are obtained by demanding the first order Taylor expansion of the finite difference equations with respect to small changes in the coordinates to vanish. We do this following exactly the procedure described in the section on relaxation methods in reference Press:2007:NRE:1403886 ().
The correction generated from the first order Taylor expansion is in general only an improvement close to the true solution. We account for this by introducing a weight that modifies the correction in each relaxation step
(53) 
We choose the weight such that the full correction is used only close to the true solution.
(54) 
In the time evolution we use as ansatz geodesic at time the geodesic form the previous time step . With a step size of usually less than 5 relaxation steps turned out to be sufficient to reach our accuracy goal of .
b.1 Initialization with Poincaré patch AdS geodesics
The line element in 3dimensional Poincaré patch AdS space is given by
(55) 
We parametrize the geodesics by an affine parameter and denote derivatives with respect to it by dot. The geodesic equations of motion allow first integrals
(56) 
where and are constants of motion and the third equation comes from the spacelike condition . The last equation shows that the geodesics have two branches. It is convenient to reparametrize both branches of the geodesics in terms of the holographic coordinate , i.e., to find expressions and . We want to initialize the relaxation code with geodesics that have a symmetric advanced time with respect to the exchange of the branches, . This is achieved by setting . Then integrating the above system yields
(57) 
where is the boundary time and we have set to zero the additive integration constant in . In the solution for the advanced time coordinate one clearly sees the ‘bending back in time’effect as one goes deeper into the bulk. By this we mean simply that becomes more negative as the holographic coordinate increases.
With the choices above the constant of motion is related to the separation of the two endpoints on the boundary
(58) 
For our numerical algorithm it turns out to be useful to choose a nonaffine parametrization that lies in a fixed interval and covers both branches at the same time, namely
(59) 
where . We realize a UVcutoff at a given value by truncating the nonaffine parater with given by
(60) 
To get the Jacobian we need in the nonaffine geodesic equation we first integrate the third equation in (57) to express the affine parameter in terms of
(61) 
Next we use the first equation in (59) and express the affine parameter in terms of the nonaffine parameter
(62) 
The Jacobian is then given by
(63) 
b.2 Initialization with conformal AdS geodesics
We proceed here similarly to appendix B.1. Starting with either of the line elements (19) in the isotropic static limit (, , ) and introducing again yields
(64) 
where we called the third coordinate , regardless of the specific case. Integrating the geodesic equations, using and setting yields
(65) 
With help of the identity
(66) 
the third equation in (65) can be used to express the affine parameter in terms of the holographic coordinate