Evolution of holographic entanglement entropy in an anisotropic system

We determine holographically 2-point correlators of gauge invariant operators with large conformal weights and entanglement entropy of strips for a time-dependent anisotropic 5-dimensional asymptotically anti-de 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 quasi-normal 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.


List of Figures
14 IX Comparison between AdS and conformal geodesics. 22

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 [1][2][3][4].
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 [5].
One powerful tool to study strongly coupled systems out of equilibrium is the gauge gravity duality where classical supergravity in d + 1 dimensions is dual to a d-dimensional strongly coupled large-N field theory that loosely speaking lives on the boundary of the gravity theory [6][7][8]. 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 Anti-de Sitter (AdS) space [9].
Particularly useful theoretical observables to monitor the thermalization process are the vacuum expectation value of the stress-energy 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 [10] for more details and references.
On the field theory side EE is defined in the following way. Dividing a system into two subsystems A and B the EE S A of the subsystem A is defined as the von Neumann entropy of the reduced density matrix obtained by tracing out the degrees of freedom of the subsystem B.
On the gravity side, the holographic entanglement entropy (HEE) is defined as the area of a minimal surface extending from some predefined surface A on the boundary into the bulk [11,12]. For time-dependent backgrounds the minimal surface has to be replaced by an extremal surface [13]. 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 [14][15][16] or Galilean conformal field theories [17]. In a seminal work, Calabrese and Cardy [18] 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 l, the EE increases linearly with time until t ∼ l/2 after which it saturates. The linear scaling with time and the crossover at t ∼ l/2 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 [19]. The behaviour of the EE in these setups has been studied extensively [19][20][21][22][23][24][25] 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 [19,25,26] where the shell is composed out of null dust and later generalized to matter with arbitrary equation of states [27,28]. The linear scaling is even present in geometries with Lifshitz scaling and hyper scaling violation [29,30]. 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 [31][32][33]. In [32,34] 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 [33] 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 [35]. Within the gauge/gravity duality the entropy of the (stationary) black hole corresponds to the entropy of the field theory. However, in time-dependent 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 [36] 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 2-point 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, 2-point 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 5-dimensional 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.

Anisotropic asymptotically AdS 5 spacetimes
In this section we review the most important details of the model first introduced in [37] and studied further in [38][39][40][41].
The 5-dimensional bulk metric that introduces anisotropy between the longitudinal and transverse directions with an O(2) rotational invariance in the transverse plane can be written conveniently in Eddington-Finkelstein coordinates, where the functions A, B and Σ only depend on the holographic coordinate r and (advanced) time v. In this coordinate system the vacuum Einstein equations read where prime denotes radial derivative and dot time derivative, viz.
for any function h(r, v). The Einstein equations (2.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 B(r = ∞, v) = B 0 (t) as in the original works of [37,38] and let the system evolve. In this case the boundary metric is curved and the conformal anomaly is present [42]. 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 [39] 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 5 , Einstein's equations in the near boundary expansion (r → ∞) are solved by 1 The coefficients in the asymptotic expansion determine the expectation value of the stress energy tensor in the dual field theory [43] 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 z = 1/r so that the boundary is located at z = 0. For our initial data we follow [40,41] and choose for the anisotropy function on the initial time slice with β = 6.6, r 0 = 4 and w = 1. In addition the initial conditions have to be supplemented with a value for the coefficient a 4 which sets the energy density of the initial state, for which we take a 4 = −1, corresponding to an equilibrium temperature T = 1/π, as we now recall. At late times we expect isotropization, B = 0. In that case we recover the usual static AdS black brane solution as follows. Solving (2.2e) and using residual gauge transformations yields Σ = r. This implies Σ = 1 andΣ = 1 2 A. Solving (2.2a) then yields A = r 2 (1−1/r 4 ), where we fixed the integration constant such that a 4 = −1. [The other equations are either trivial, (2.2b) and (2.2d), or redundant, (2.2c).] The result for A is the usual Killing norm for the static AdS black brane. Surface gravity is given by κ = 1 2 A r=1 = 2 so that the Hawking temperature is T = κ/(2π) = 1/π.
In the generic anisotropic case, B = 0, we solve the Einstein equations (2.2) numerically for the initial conditions (2.7). In this background we then study the evolution of 2-point 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-point correlators
The equal time 2-point function for an operator of large conformal weight ∆ can be computed via a path integral as [44,45] where the integral is a sum over all possible paths with endpoints at (t, x ) and (t, x) and L(P) 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 1/ 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 L g denotes the length of the corresponding geodesic. To leading order only the geodesic with the smallest value of L g contributes, whose length we denote by L, which explains the second approximation 2 . 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 L therm . In terms of the renormalized length δL = L − L therm the 2-point function becomes (2.9) This means that we can obtain the time evolution of 2-point functions by looking at spacelike geodesics that are anchored at the boundary at fixed separation l 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 σ ∈ [−σ m , σ m ]. To obtain the lengths of the geodesics we have to solve the geodesic equation for the two subspaces given by the line elements For the separation in the transverse direction the geodesics end at (v(±σ m ) = t, . With this choice of boundary conditions the lengths of the geodesics in the background (2.1) are given by where prime denotes the derivative with respect to σ.
It is important to point out that we can only study geodesics after some advanced time v > v min with boundary separations below a maximal separation l < l max . 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 v ∈ [v 0 , ∞]. As the geodesics reach into the bulk they bend back in advanced time leaving the computational domain for advanced times v < v min as well as extending too far into the bulk for separations l > l max .
In section 3.2 below we discuss how to solve the geodesic equations numerically.

Holographic entanglement entropy
In time dependent systems the covariant HEE [13] for some boundary region A is obtained by extremizing the 3-surface functional that ends on the boundary surface A. In the dual field theory the EE is then conjectured to be given by [11,13,47] 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 (2.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 φ i (x α ) and write the line element as where h αβ is a 3-dimensional metric with coordinates (v, r, x 1 ) where x 1 represents the coordinate we choose to have finite spatial extent, i.e. either x or one of x ⊥ . The remaining (non-compact) coordinates are then denoted by x 2 , x 3 , which we choose to be two of our three world-volume coordinates; the third one is denoted by σ. Parametrizing the 3-dimensional coordinates as x α = (v(σ), r(σ), x 1 (σ)), the area functional (2.14) can be written as Performing the integration over the Killing coordinates x 2 and x 3 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 3-surface corresponding to a boundary region A of strip-topology is then reduced to a 1-dimensional problem. In fact, from the expression (2.17) on can see that the problem of finding the extremal 3-surfaces reduces to finding geodesics of the conformal metric The 3-dimensional 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 3 Numerical implementation

Einstein equations
We solve the Einstein equations (2.2) using pseudo spectral methods as described in detail in [41], 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 2-point functions and HEE. For that reason we want a large computational domain in the holographic coordinate z = 1/r. In all the computations we took z ∈ [0, 1.6] with the final position of the horizon located at z = 1. For the time evolution it is sufficient to use a fourth order Runge Kutta method with time steps δt = 10 −3 . All the computations were done with the open source software GNU Octave [48].

Geodesics
To compute 2-point 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 non-affine parameter σ , where τ = τ (σ) is the usual affine parametrization, dX µ dτ dX ν dτ g µν = 1. In terms of σ the geodesic equation readsẌ whereẊ µ = dX µ dσ and J = d 2 τ dσ 2 / dτ dσ 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 (3.1) 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 V, Z, X. 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: This set of equations is a two-point boundary value problem, which is usually either solved with shooting methods or relaxation methods [49]. 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 one-parameter families of geodesics (evaluated at different constant time slices) we can take the solution for the n th family member as initial guess for the (n + 1) st family member. More details of our implementation of the relaxation method are described in appendix B.

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.

Results
In this section we display and discuss our main results. In all figures where a time-axis is plotted we measure the boundary time t or the bulk advanced time v in units of the temperature of the final black brane, T = 1/π. The separation length l of 2-point functions and HEE and the corresponding boundary coordinates are given in units of T as well. To make the approach to thermal equilibrium most transparent we use normalized quantities for the geodesic length L ren and HEE S ren defined by where L (S) is the unrenormalized length (HEE) and L th (S th ) is the corresponding thermal value.

2-point 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.

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 2-point functions. This is also the reason why for the same boundary separation the HEE reaches equilibrium later as the 2-point 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 2-point functions with the oscillations less pronounced. We exhibit these features now in some plots.

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 [40,41] 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 with the lowest QNM given by [50,51] On the field theory side QNMs appear as poles in the retarded Green function [51][52][53][54]. 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 e −Im[ω 1 t] L ren 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 [55][56][57].
It turns out that HEE also follows this pattern. In [58] 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 ∆S A 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 e −Im[ω 1 t] S ren 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 2-point 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 2 log(l/2) and 2l respectively. At large separation the HEE also scales linearly with the separation length, whereas at small separation it is proportional to 1/l 2 . All these results agree precisely with the perturbative expressions derived in the limits of small and large temperatures [47,59]. 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.

Conclusions
In the paper at hand we have studied the time evolution of 2-point correlation functions in the geodesic approximation and holographic entanglement entropy in an anisotropic N = 4 SYM plasma.
To obtain the 2-point correlation function for separation in transverse and longitudinal direction we solved for the length of the geodesics in the anisotropic background (2.1) by using a relaxation method with the zero temperature geodesic as an initial guess. At the n-th step the solution of the (n − 1)-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 2-point correlation function goes through.
The 2-point 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 [58].
The methods developed and applied in this work can be generalized to other timedependent asymptotically anti-de Sitter backgrounds of interest, such as colliding shockwave backgrounds [60,61], which have been constructed numerically in the past few years [62][63][64]. 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.
In this appendix we give some details on how we numerically solve the Einstein equations (2.2). We work with the inverse radial coordinate z ≡ 1/r such that the boundary is located at z = 0. Due to asymptotic AdS boundary conditions, the metric functions A and Σ diverge as z → 0. It is numerically favourable to define new functions with the known divergent pieces removed and rescale them with appropriate powers of z so that the resulting functions are finite or vanish as z → 0; the precise boundary conditions on the new functions are presented in (A.6) below. This leads us to the following field redefinitions The redefined anisotropy function B allows to extract b 4 (t) simply from the boundary value where B = ∂ z B. In terms of the redefined fields the first four Einstein equations (2.2) can be rewritten in the form with the source functions given by The relation between dot-derivative and time-derivative, originally given in (2.3), turns intȯ The boundary conditions for the redefined fields read where we set a 4 = −1. On the initial time slice the boundary condition (A.6c) is computed from the initial conditions where we set β = 6.6, z 0 = 1/4 and ω = 1. In our computations we choose z i ∈ [0, 1.6] and N = 60. These points can be used to construct the entries of the spectral differentiation matrix D ij [66]: where the matrix L ij is given by .14) and the source vector (jΣ) i by The boundary conditionΣ 1 = a 4 is implemented by setting (jΣ) 1 = a 4 and L 1j = δ 1j . The solution vectorΣ i is then obtained by multiplying the inverse of L ij with the source vectoṙ We do this with the OCTAVE routine linsolve. The equations for Σ i , A i andḂ i can be solved in the same way.
To advance the solution for B to the next time slice it is sufficient to use a simple fourth order Runge Kutta method [49] where the koefficients k i are given by and ∂ v B is computed from (A.5). In our simulations we use a stepsize of δv ≈ 0.01 which is sufficient to achieve a stable time evolution.

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 σ i = h i of equidistant spacing h = σ N −σ 1 N with i = 1, . . . , N . We use a grid with N = 500 points in all our simulations. The upper and lower bound of this grid are given by σ 1 = −1+ and σ N = 1− respectively, where denotes a UV cutoff. The discretized version of our trial solution on this grid is written as The explicit form of the trial solutions used in the computations of 2-point 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 z-position of the corresponding cutoff surface is fixed at Z U V = 0.05, i.e. at r U V = 20. The boundary conditions are imposed at this cutoff surface 3) where t denotes the time and l the separation on the cutoff surface. The finite difference representation of the geodesic equation (3.2) is given by where E k i is the residual at point i in equation k; further (p V ) i =V i denotes the first derivatives of the geodesic coordinates; quantitities with bar are averaged like ; the Christoffel symbols (Γ X XZ ) i are evaluated from the averaged metric functions; the explicit form of the JacobianJ i for the case of 2-point 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 E k i will be rather large. To quantify the deviation of a given trial solution to the true solution we use the following measure The strategy is to compute increments ∆X µ for the geodesic coordinates X µ such that X µ new = X µ old + ∆X µ is an improved approximation to the previous solution X µ old . This we do iteratively until we reach δ < 10 −15 . (B.7) 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 [49]. The correction ∆X µ 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 We choose the weight α such that the full correction is used only close to the true solution.
In the time evolution we use as ansatz geodesic at time t i = t i−1 + δt the geodesic form the previous time step t i−1 . With a step size of δt ≈ 0.05 usually less than 5 relaxation steps turned out to be sufficient to reach our accuracy goal of δ < 10 −15 .

B.1 Initialization with Poincaré patch AdS geodesics
The line element in 3-dimensional Poincaré patch AdS space is given by We parametrize the geodesics by an affine parameter τ and denote derivatives with respect to it by dot. The geodesic equations of motion allow first integralṡ where L and E are constants of motion and the third equation comes from the spacelike condition ds 2 = 1. 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 z, i.e., to find expressions x ± (z) and v ± (z). We want to initialize the relaxation code with geodesics that have a symmetric advanced time with respect to the exchange of the branches, v + (z) = v − (z). This is achieved by setting E = 0. Then integrating the above system yields where v 0 is the boundary time and we have set to zero the additive integration constant in x ± . 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 v(z) becomes more negative as the holographic coordinate z increases. With the choices above the constant of motion L is related to the separation l of the two endpoints on the boundary For our numerical algorithm it turns out to be useful to choose a non-affine parametrization that lies in a fixed interval and covers both branches at the same time, namely where σ ∈ [−1, 1]. We realize a UV-cutoff at a given value Z U V by truncating the non-affine parater σ ∈ [−1+ , 1− ] with given by To get the Jacobian we need in the non-affine geodesic equation we first integrate the third equation in (B.12) to express the affine parameter τ in terms of z τ (z) = ± dz Next we use the first equation in (B.14) and express the affine parameter τ in terms of the non-affine parameter σ τ (σ) = ∓artanh σ 2 − σ 2 . (B.17) The Jacobian is then given by (B.18)

B.2 Initialization with conformal AdS geodesics
We proceed here similarly to appendix B.1. Starting with either of the line elements (2.19) in the isotropic static limit (B = 0, Σ = r, A = r 2 ) and introducing again z = 1/r yields where we called the third coordinate x, regardless of the specific case. Integrating the geodesic equations, using ds 2 = 1 and setting E = 0 yieldṡ In order to obtain the solution for the geodesic it turns out to be convenient to solve for x in terms of z x ± = ± dz Lz 3 √ 1 − L 2 z 6 = ∓ l 2 ± Lz 4 4 2 F 1 1 2 , 2 3 , 5 3 ; L 2 z 6 . (B.23) The subscript ± corresponds to the two branches of the geodesic and l/2 is an integration constant corresponding to the boundary value of the geodesic since the hypergeometric function is zero at the boundary z = 0. The endpoint of the geodesic is given by L 2 z 6 = 1, meaning that the two branches do not join in general. In order to get a smoothly joining geodesic we need to adjust L and express it in terms of the boundary separation l L = π where σ ∈ [−1, 1] and Z max = 2l √ π Γ( 7 6 ) Γ( 5 3 ) denotes the z-position at which the two branches join. We realize a UV-cutoff at a given value Z U V by truncating the non-affine parater σ ∈ [−1+ , 1− ] with given by The affine parameter τ in terms of the non-affine parameter σ reads The Jacobian is then given by As can be seen from Fig. IX for the same boundary separation the conformal geodesic extends much further into the bulk compared to the Poincaré patch AdS case. As a consequence, the boundary separation we can choose for the HEE is much smaller than for the 2-point functions for any given value of the boundary time.