Exploring nonlocal observables in shock wave collisions

We study the time evolution of 2-point functions and entanglement entropy in strongly anisotropic, inhomogeneous and time-dependent N=4 super Yang-Mills theory in the large N and large 't Hooft coupling limit using AdS/CFT. On the gravity side this amounts to calculating the length of geodesics and area of extremal surfaces in the dynamical background of two colliding gravitational shockwaves, which we do numerically. We discriminate between three classes of initial conditions corresponding to wide, intermediate and narrow shocks, and show that they exhibit different phenomenology with respect to the nonlocal observables that we determine. Our results permit to use (holographic) entanglement entropy as an order parameter to distinguish between the two phases of the cross-over from the transparency to the full-stopping scenario in dynamical Yang-Mills plasma formation, which is frequently used as a toy model for heavy ion collisions. The time evolution of entanglement entropy allows to discern four regimes: highly efficient initial growth of entanglement, linear growth, (post) collisional drama and late time (polynomial) fall off. Surprisingly, we found that 2-point functions can be sensitive to the geometry inside the black hole apparent horizon, while we did not find such cases for the entanglement entropy.


Introduction
The gauge/gravity duality has established itself as a valuable tool in the quest for a better understanding of strongly coupled systems. In particular it is used to gain insights into the thermalization process of non-abelian plasmas (such as the quark gluon plasma generated in heavy ion collisions at RHIC and LHC) by studying the gravitational dual of N = 4 super Yang-Mills (SYM) theory, a maximally supersymmetric conformal field theory (CFT) in four spacetime dimensions. The equilibration of the field theory is then mapped to black hole formation on the gravity side. In the last decade there has been considerable progress in setting up collisions of SYM matter in various scenarios and studying its evolution. One of the starting points was the study of perfect fluid dynamics in a boost invariant setting [1,2]. In [3] it was possible to study far-from-equilibrium dynamics by numerically solving the full Einstein equations in an anisotropic but otherwise completely homogeneous system. Trying to come closer to mimic a heavy ion collision led to the idea [1] of colliding delta like gravitational shock waves [4,5], which are dual to lumps of energy in the SYM theory moving at the speed of light. The next step was to make the system anisotropic and inhomogeneous by the collision of gravitational shock waves which are homogeneous in the transverse direction and have finite width in the longitudinal direction [6]. It was found that a hydrodynamic description of the plasma is valid even when the anisotropy is still large [7]. This onset of hydrodynamic behavior is now termed hydrodynamization. Further advances include radial flow [8], the effect of different initial conditions [9], the collision of two black holes [10], and more [11][12][13].
Now it is even possible to simulate the collision of two localized lumps of matter to mimic off-central nucleus-nucleus [14,15] and proton-nucleus collisions [16].
Despite all the advances one has to keep in mind that in heavy ion collisions there are many energy scales involved and to get an accurate understanding of the thermalization mechanisms involved strong and weak coupling phenomena must be combined. One step towards this direction is the combination of different effective descriptions [17] or by constructing a semi-holographic framework where the weakly and strongly coupled sector can interact with each other [18,19].
So far, in most colliding shock wave studies the quantities of interest are local quantities, i.e. the components of the energy momentum tensor, such as the energy density or the pressures. This allows to determine if local equilibrium is reached, here understood as the local applicability of hydrodynamics. In order to gain further insight into the thermalization process the time evolution of nonlocal quantities, such as various correlation functions (e.g. Wightman function or Feynman propagator), in coordinate space need to be considered. This is still a complicated task but two such nonlocal quantities can be obtained relatively easily with the help of the gauge/gravity duality, namely the equal time 2-point function for scalar operators of large conformal weight and entanglement entropy (EE). In the context of thermalization these quantities where first computed to study the analog of quenches in conformal field theories [20] via the collapse of thin shells [21,22] in AdS space, where the EE shows universal behavior. After the initial short early time epoch the EE grows linearly with time, which is independent of the entangling regions [23] or the equation of state [24,25]. In these works the EE is a monotonically increasing function that approaches the final equilibrium value from below. However, this universal behavior disappears in more complicated setups. For example, when a radially collapsing scalar field forms a black hole the EE can be non-monotonic [26][27][28][29][30][31]. In anisotropic N = 4 SYM the EE and equal time 2-point functions show oscillatory behavior with exponential damping at late times which is given by the lowest quasinormal mode [30]. Analytic progress has been made in [32] where the late-time behavior of two-point functions, Wilson loops and entanglement entropy has been studied perturbatively in a boost-invariant system.
The equal time 2-point function can be obtained from the length of space like geodesics which are anchored to the boundary of anti-de Sitter (AdS) and extending into the bulk. Although the geodesic approximation is only valid for operators of large conformal weight, a comparison of the Feynman propagator for a scalar field with conformal dimension ∆ = 3/2 with the geodesic approximation revealed that qualitatively they show the same behavior [33]. Similarly the holographic entanglement entropy (HEE) can be obtained from the area of extremal surfaces [34,35].
In this work we extend the existing studies by investigating the time evolution of equal time 2-point functions and HEE in the colliding shock wave geometry for different initial conditions, carefully differentiating between wide, intermediate and narrow shocks, which turn out to have quite different phenomenology. Our results allow to use HEE to distinguish between the phases corresponding to wide or narrow shocks, in a sense that we shall make precise.
This paper is organized as follows. In Section 2 we introduce the geometry and the different initial conditions. The results for the equal time 2-point function and EE are discussed in Sections 3 and 4, respectively. In Section 5 we conclude.
2 Gravitational shock waves in asymptotically AdS 5 The holographic setup we consider describes the collision of two sheets of energy having Gaussian shape in their direction of motion and which are homogeneous in the remaining two spatial directions. These shocks serve as caricatures of two highly Lorentz contracted nuclei which approach each other at the speed of light and induce non-abelian plasma formation as they collide.
On the gravity side the corresponding 5-dimensional bulk metric is rotationally invariant and homogeneous in the transverse plane (x 1 , x 2 ) but inhomogeneous in the longitudinal direction y, which is the direction of motion of the shocks. The metric ansatz in Eddington-Finkelstein coordinates reads where the functions A, S, B and F depend on the holographic coordinate r, (advanced) time v and longitudinal coordinate y, but are independent from the transversal coordinates x. The equations of motion can be found e.g. in [6] and are solved near the boundary by where ξ(v, y) encodes the residual diffeomorphism freedom r → r + ξ(v, y). It is possible, though not necessarily numerically convenient, to choose ξ = 0. As usual the normalizable modes a 4 (v, y), b 4 (v, y) and f 4 (v, y) are undetermined by the near-boundary expansion and require a solution of the full bulk dynamics. These coefficients in the asymptotic expansion determine the expectation value of the conserved and traceless stress energy tensor in the dual field theory [36]

Initial conditions
The pre-collision geometry describing two shocks moving in ±ỹ-direction can be written down in Fefferman-Graham coordinates (r,t,ỹ, x) as follows [1] ds 2 =r 2 η µν dx µ dx ν + 1 where η µν denotes the usual 4-dimensional Minkowski boundary metric and h(t ±ỹ) is an arbitrary function for which we choose a Gaussian of width ω and amplitude µ 3 In this gauge the non-vanishing components of the energy momentum tensor read and describe two lumps of energy with maximum overlap att = 0. At early timest −w, when the shocks h(t ±ỹ) have negligible overlap, the line-element (2.5) is close to an exact solution to Einstein's equations, but aroundt = 0 their dynamics can only be computed numerically.
We do this for three different initial conditions h n,i,w (ỹ) describing qualitatively different situations that we shall refer to as narrow, intermediate and wide shocks, where in all cases the initial position of the shocks is atỹ 0 = ±7/4. For the width of the shocks we take ω n,i,w = 0.1, 0.25, 0.5 and we will display all our results in units of µ.
For the numerical evolution scheme the initial data needs to be transformed to Eddington-Finkelstein coordinates (r, v, y, x 1 , x 2 ) by solving for radially infalling null geodesics in the background (2.5), leading to ordinary differential equations, which are solved for appropriate boundary conditions at the boundaries of the radial domain. We omit a discussion of the numerical details concerning this coordinate transformation and the subsequent evolution and refer the reader to [37,38], where the full procedure is explained.

Evolution of the energy momentum tensor
The time evolution of the energy momentum tensor for colliding shocks has been studied extensively in [6,9,38,39]. In Fig. I we show the evolution of the energy density E(t, y) extracted from the numerical evolution for the different initial conditions stated above. As discussed in [9] the energy density behaves qualitatively different in collisions of narrow shocks and in those of wide shocks. This cross-over is not only of academic interest, but also for applications, since it was argued that the narrow shocks describe more adequately the situation at LHC, while the wide shocks are more suitable for RHIC [9] (see also [15]). We list below some relevant properties that differ between wide and narrow shocks: • Narrow shocks exhibit transparency, in the sense that they pass through each other and, even though their shape gets altered and they decay, they continue to move at the speed of light after the collision. By contrast, wide shocks realize a full-stopping scenario, in the sense that the energy density is localized mostly in the central region after the collision, and the shocks themselves not only change their shape but also get slowed down. Wide shocks then lead to initial conditions for hydrodynamics where all velocities are close to zero, i.e. there is a hydrodynamical explosion in close similarity to the Landau model of heavy ion collisions [40].
• Narrow shocks can yield regions of negative energy density after the collision right behind the original shocks on the lightcone. Curiously, this region does not admit a local restframe [41], but also does not violate general principles of quantum field theory, such as the averaged null energy condition [42]. At y = 0 after the shocks pass through each other, the energy density grows at early times as E = 2µ 6 t 2 + O(t 5 ), which implies pressures equal to P /E = −3 and P ⊥ /E = 2. This feature was first observed for δ-like shockwaves analytically [4] and then numerically for sufficiently narrow Gaussian profiles [9]. By contrast, for the wide shocks the energy density and pressures remain positive everywhere.
Given the substantial differences in local observables one may expect that the characteristic features for narrow and wide shocks also show up in nonlocal observables, like 2-point functions and HEE. In the remainder of this work we verify this expectation by explicit computations, starting with the 2-point functions in the next section.

Two-point functions
Within AdS/CFT the equal time 2-point function of operators O with large conformal weight ∆ can be computed from the length L of spacelike geodesics in the bulk geometry [43,44] via In asymptotically AdS the length of a geodesic which is attached to the boundary is infinite and a regularization scheme must be adopted. A natural way to regularize is to subtract the length L 0 of a geodesic in AdS corresponding to the vacuum value of the correlator For illustrative purposes we set ∆ = 1 when we display our results which is the same as interpreting L reg to be given in units of ∆. Thus, the two point functions we compute are defined as follows In order to obtain the geodesic length we solve the geodesic equation numerically with a relaxation algorithm which iteratively relaxes an initial guess to the true solution. For a detailed description of the relaxation algorithm we refer the interested reader to [30].

Geodesics in the shock wave geometry
For simplicity we restrict our attention to geodesics that only extend along the y-direction and not along the transverse directions (x 1 , x 2 ), i.e. we consider geodesics in the three dimensional bulk-subspace where z = 1/r. To find these geodesics we solve the (non-affine) geodesic equation subject to the following boundary conditions at z = 0 where X µ (σ) are the embedding functions of the geodesic and dots denote derivatives with respect to the non-affine parameter σ ∈ [−1, 1]. The quantity J = d 2 τ dσ 2 / dτ dσ denotes the Jacobian of the reparametrization from the affine parameter τ , defined by ( dX dτ ) 2 = 1, to σ. The boundary time and separation for which the geodesics are computed are denoted by t and l respectively. The fictitious viscous force provided by the Jacobian J helps with the numerics, resulting in better convergence of the relaxation algorithm.
Working in asymptotically AdS makes it natural to choose as an initial guess a geodesic in pure AdS which can be written as In this parametrization the affine parameter is given by τ (σ) = ∓arctanh σ √ 2 − σ 2 from which the Jacobian needed in (3.5) can be computed We assume the boundary separation to be centered around y = 0. Describing off-central geodesics requires some straightforward modifications of our formulas. The bulk part of the geodesic length, which is the contribution from z > z cut , follows from integrating the line elements (3.4) and (3.7) where the metric functions (A, B, S, F ) have to be evaluated along the geodesic X µ (σ). In order to realize an IR-cutoff at a given value z cut the range of the non-affine parameter σ ∈ [σ − , σ + ] has to be bounded by The near boundary part of the geodesic length, which is the contribution from 0 ≤ z ≤ z cut , can be extracted form the near boundary solution of the geodesic equation. Near z = 0 the embedding functions and the Jacobian can be expressed in terms of a power series in z In Appendix A we give the explicit expressions for the expansion of the metric that we have used. The coefficients (t n , y n , j n ) in Eq. (3.12) can be computed by solving the geodesic equation order by order in z, which leads to the following expressions Here we fixed the leading coefficients by the boundary conditions (3.6), but the coefficients v 2 and y 2 cannot be determined by a near boundary expansion. This is analogous to the normalizable modes of the metric, which are also sensitive to the full bulk geometry. The pure AdS solution is given by which hence has v 2 = 0 and y 2 = ∓1/l. We can now compute the near boundary expansion of the geodesic length, which for one branch is given by where the leading AdS divergent 1 z term nicely cancels out. The regularized geodesic length L reg , which we need to evaluate Eq. (3.3), is the sum of the bulk contribution and the near boundary contribution 1 When using Eq. (3.16) to evaluate Eq. (3.3) numerically one has to ensure that the results are, to some required accuracy, independent of the discretization and the cutoff. We require this accuracy to be of the same order as the maximal residual (= 10 −5 ) we allow in the geodesic equation and below which we stop to iterate the relaxation procedure. We checked the convergence of the 2-point function with the gridsize in the range from 50 up to 400 gridpoints and find that for more than 200 gridpoints the change is smaller than O(10 −5 ) which is the same order as the allowed residual (see Appendix B). Sample checks are presented in Appendix B, where only a mild cutoff dependence of O(10 −5 ) is obtained for a range z cut = [0.01, 0.1], which is again of the same order as the allowed residual. Based on this analysis we choose 200 gridpoints to discretize our geodesics and set z cut = 0.075 in all our calculations.

Evolution of two-point functions
In this section we present our numerical results for 2-point functions in holographic shock wave collisions. Before we discuss the actual results let us start with some remarks regarding the computational domain used in these simulations. As input for the relaxation algorithm we provide numerical results of the shock wave metric in a finite domain in (t, y, z). This computational domain, in which we can solve the geodesic problem, is bounded by µt ∈ [−1.5, 6], µy ∈ [−5, 5], where in the radial coordinate we have chosen the apparent horizon as a natural bound z ∈ [0, 1.08z AH ]. That means whenever we display geodesics which reach beyond this radial domain, which can happen as we discuss below, an extrapolated version of the metric is used 2 . For a given choice of boundary conditions (µt, µl) the final shape of the geodesic in the bulk is a priori unknown, i.e. initially we do not know if the geodesic resides entirely within or extends beyond the computational domain in which the metric is known. Therefore finding a feasible set of parameters (µt, µl) for a given computational domain requires some trial and error. The geodesics bend back in advanced time as they reach into the bulk, leaving the computational domain for too early boundary times. Therefore we can display our results only in a finite time near the collision time t = 0 where all geodesics with different boundary separation lie in the computational domain. All these points apply accordingly to the EE simulation. For the time evolution it is of advantage, after using the pure AdS geodesic at the initial time, to use the previous solution to initialize the next time step. This approach turns out to be numerically extremely efficient and the relaxation algorithm reveals its full power, since in most cases the result at a given time is an excellent estimate at the next time step. A time step of ∆t = 0.1 allows to resolve nicely the shape of the 2-point function and reduces the required number of iterations almost to a minimum. Usually two iterations are sufficient to achieve relative residuals in the geodesic equation which are < 10 −5 and in many cases even one or two orders smaller.
We follow the same logic when we compute the evolution in the boundary separation, where this approach is not only numerically efficient but also crucial to reach large separations. Undeformed ansatz geodesics of large separation typically reach far beyond the radial domain and finding the true solution using such geodesics to initialize the relaxation inevitably fails. We circumvent this problem by initializing with an ansatz of small separation (µl = 0.2), which comfortably resides within the computational domain. Then we increase step by step the boundary separation and use the result for a given separation as ansatz for the next separation step. By using a step size of ∆l = 0.1 we can nicely resolve the shape of the 2-point function and the relaxation usually converges again after two iterations. Since the relaxed geodesics are typically strongly deformed in direction away from the apparent horizon, i.e. the upper bound of the radial domain, we can reach separations which were inaccessible by simply relaxing the corresponding ansatz geodesic. We like to discuss first the results from the time evolution before we go to the evolution in the separation. In Fig. II (left) the whole setup for wide, intermediate and narrow shocks is displayed. The dark surface represents the radial position of the apparent horizon z AH (t, y). The evolution of the energy density of the boundary conformal field theory is shown by a contour plot located at the boundary z = 0. The green lines are geodesics at various time steps for a given separation. For narrower shocks the minimum of the apparent horizon is closer to the boundary and the influence on the shape of the geodesics is bigger. One can see that the tips of the geodesics tend to avoid the apparent horizon and the evolution of the tips show a similar shape as the apparent horizon. Once the profile of the geodesics is found the evolution of the 2-point functions can be extracted by computing their length. On the right hand side of Fig. II the evolution of the 2-point functions for various boundary separations for the different geometries are displayed.
Let us now summarize the most salient features in the time evolution of the 2-point function during a holographic shock wave collision.
• rapid onset of linear de-correlation: The system starts in some correlated state.
As the shocks are getting closer more and more short range correlations are destroyed and the system rapidly starts to de-correlate in a linear fashion until a local minimum is reached. The rapid onset of the linear regime is clearly visible for the narrow shocks in Fig. II, where for intermediate and wide shocks the onset lies outside our computational domain for larger separations, but the linear regime is still visible. For intermediate and narrow shocks the minimum is located close to t = 0 where the energy density is maximal. For wide shocks this minimum is reached significantly before t = 0.
• premature de-correlation: A careful tracking of the position of the minimum as a function of the boundary separation reveals that it is shifted to earlier times as the separation increases. This effect, which is very small and therefore hardly visible in Fig. II, is a robust feature of all three kinds of shocks that we have studied.
• linear correlation restoration: During the collision, when the shocks interact, new correlations are formed in the system. As the shocks move outwards (t > 0), the correlations are linearly restored for all three kinds of shocks.
• correlation overshooting of narrow shocks: After the linear restoration regime, the correlations in wide and narrow shocks approach their final values in very different ways. For intermediate and narrow shocks the correlations significantly overshoot their final values before they finally approach them from above. In the case of wide shocks this effect is strongly damped and the correlations approach their initial value almost monotonically from below.
We switch now to the scaling of the 2-point function with the separation. The holographic setup and the results for the evolution of the 2-point function are displayed in shocks. The position of this additional maximum is centered around the position of the outgoing shocks. It is suggestive that narrow shocks which pass through each other almost transparently remain correlated for some time after the collision while wide shocks stop each other before they explode hydrodynamically and the correlations are completely lost. This motivated us to study the correlations between the shocks themselves, which we do systematically in Section 3.3. There we find that the correlations between intermediate and narrow shocks significantly grow after the collision before they start to decay, where the correlations between wide shocks decay immediately. Interestingly, for larger separations the geodesics remain outside the horizon for early times, but they cross the horizon after a time of around µt = 1.5. This can be seen from the blue curves in Fig. III and is displayed more clearly in Fig. IV where we plot the tip of the geodesic located at y = 0, for different separations and the position of the apparent horizon at y = 0. This happens for all the initial conditions (wide, intermediate, narrow) we have studied and is in strong contrast to the EE case where we do not find extremal surfaces which cross the horizon. The crossing after a time of µt = 1.5 is perhaps counterintuitive since geodesics are expected to remain outside the horizon when the system is close to equilibrium. Indeed, hydrodynamics applies after a time µt = 0.89 [9], which is well before the crossing of the geodesics. At later times presumably the geodesics indeed remain outside again, though our numerics did not allow to determine the precise time at which this is the case.

Correlations of colliding shocks
Instead of studying the time evolution of the 2-point function between two fixed points in space, in the context of heavy ion collisions it might be more interesting to actually study the correlation between the two shocks itself. In order to do so, the endpoints of the geodesics follow the maxima of the energy density.
When the separation of the endpoints becomes smaller than three times the cutoff we fix the endpoints to this value until the distance between the two maxima after the collision As already discussed in Section 3.2, for wide shocks the behavior is qualitatively differ-ent than for intermediate and narrow shocks. As the two wide shocks approach each other their correlation increases almost linearly until it reaches a plateau, which is the point when the separation of the endpoints is smaller than three times the cutoff. Once the shocks separate again from each other their correlation decreases.
As the shocks get narrower the initial growth slows down because the shocks start to overlap later. After the fixed separation period a local minimum appears after which the correlations continue to grow to reach another maximum which appears later for narrow shocks. In addition, the maximum correlation is highest for narrow shocks.
This behavior is reminiscent of the full stopping and transparency scenario for wide and narrow shocks considered in [9]. As the wide shocks start to interact the energy density starts to pile up and all the energy density is contained in a small region after which hydrodynamical explosion occurs. This behavior is also encoded in the 2-point function which reaches a maximum and can only decrease when hydrodynamic explosion occurs.
For narrower shocks the situation is different. The shocks almost move through each other. Their shape gets altered but no hydrodynamic explosion occurs. The shocks separate from each other and plasma between them forms resulting in a growth of the correlation also after the collision. At sufficiently late times, when the shocks are separated far enough and a hydrodynamical description is applicable the 2-point function decreases rapidly.
To summarize, there is a general pattern appearing. As the shocks become narrower the initial growth slows down, the maximum correlation increases and occurs later.

Entanglement entropy
In this section we monitor the evolution of EE. In time dependent systems the covariant HEE [35] 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 [34,35,45] Under certain circumstances the problem of finding extremal surfaces can be reduced to finding geodesics in an auxiliary space-time and the problem of solving nonlinear partial differential equation can be circumvented [30]. In the case at hand this can be achieved by considering a stripe entangling region with finite extent in the longitudinal direction y and infinite extent in the homogeneous transverse directions (x 1 , x 2 ) for which (4.1) simplifies to The surface functional (4.3) suffers from two kinds of infinities, one from the integral V = dx 1 dx 2 over the homogeneous directions and another one from the infinite geodesic lengthL in the auxiliary spacetime Ω 2 h µν . Since the infinite volume factor V contains no dynamical information these singularities are avoided by considering EE densities S EE V . Analogous to the 2-point function we regularize the geodesic lengthL by subtracting the corresponding auxiliary vacuum contributionL 0 . The observable we compute is the regularized EE density per Killing volume in units of 4G N

Geodesics in the auxiliary spacetime
Our aim is to compute the EE for a stripe region with finite extent in y-direction and infinite extent in (x 1 , x 2 ) using formula (4.4). Therefore we have to find geodesic lengthsL andL 0 in the corresponding auxiliary spacetimes. The auxiliary spacetime, which is related to the metric (2.1) by a conformal factor Ω 2 = S 4 e 2B , reads This time we initialize the relaxation algorithm with a geodesic in Poincaré patch AdS (3.7) times a conformal factor Ω 2 0 = 1 z 4 Like for the Poincaré patch AdS geodesics we choose a non-affine parametrization 8l 3 Γ[7/6] 3 ensures that the two branches, discriminated by sgn(σ), join smoothly at Z max = 2l √ π Γ( 7 6 ) Γ( 5 3 ). The affine parameter τ in terms of σ reads and the Jacobian evaluates to Using the ansatz (4.7) and the corresponding Jacobian (4.11) in the relaxation algorithm allows us to compute geodesics in the auxiliary spacetime (4.5). The bulk parts of the geodesic lengths in Eq. (4.4), which are the contributions from z > z cut , follow from integrating the line elements (4.5) and (4.6) where in this case the bounds of the integral σ ± , implementing the IR-cutoff at z = z cut , are given by We build the near boundary part (0 ≤ z ≤ z cut ), like for the 2-point function, from the asymptotic solution of the geodesic equation in the conformal spacetime, which leads to the following near-boundary expansion where the normalizable modes a 4 (v, y), b 4 (v, y) and f 4 (v, y) are evaluated at v = t 0 and y = ± l 2 . We again have two undetermined constants v 4 and y 4 , which now appear two orders higher than for the case of the 2-point function. Again we also have the analytic solution in the auxiliary pure AdS space time The near boundary contribution to the geodesic length for both endpoints evaluates tõ where the divergent term cancels out again. Now this formula is clearly more useful, as the two leading contributions do not depend on the unknown coefficients v 4 and y 4 , which hence allows to reduce the cut-off dependence significantly. The regularized EE of Eq. As for the 2-point function we checked the convergence of S reg with the gridsize in the range from 50 up to 400 gridpoints and find again that for more than 200 gridpoints the change in S reg is smaller than O(10 −5 ) which is the same order as the allowed residual we choose in the relaxation algorithm.
To achieve cutoff independence of S reg turns out to be more delicate than for the 2-point function. Now for a range z cut = [0.05, 0.1] we obtain a slightly worse cutoff dependence of O(10 −3 ) which is however sufficient for our qualitative studies where S reg = O(10 −1 ) and the influence of the cutoff can be estimated to be ≈ 1% (see Appendix B). Again we choose 200 gridpoints to discretize our geodesics and set z cut = 0.075 in all the calculations we present below.

Evolution of entanglement entropy
In this section we present our numerical results for the EE. The shape of EE as a function of time originates from a complicated interplay between the different metric functions appearing in the energy momentum tensor. However, most features can be understood in terms of energy density and pressures. In Fig. VI we display the time evolution of HEE for various separations in the two different scenarios. It can be characterized by four distinct regions: 1. rapid initial growth: Once some energy density enters the entangling region the rapid initial growth starts. The narrower the shocks the more rapidly the initial growth happens, because the rate at which the energy density enters the entangling region is bigger than for wider shocks.
2. linear growth: The linear growth starts when the two shocks start to overlap and the energy piles up, with a steeper slope for larger separations. This is the same behavior as the post-local equilibration growth after a global quench [46]. The maximum occurs with a short delay compared to the maximum energy deposited in the entangling region, with a more pronounced delay for wider shocks.
3. post collisional regime: The post collisional regime is quite different for the three cases considered. For wide shocks the EE falls off without any additional features.
In the case of intermediate shocks a small shoulder appears. In the case of narrow shocks this shoulder turns into a new feature, where an additional minimum appears where the coefficient a w,i,n depends on the initial conditions and the separation. In Table 1  where the function S is evaluated at the position of the apparent horizon and integrated over the same intervals as for the EE. The late time behavior is displayed in Table 2 and barely depends on the separation. It is expected on general grounds that at very late times and large separations, far beyond our computational domain, the effective entropy density and EE show the same fall off behavior.
Let us now discuss the results from the evolution in the separation. The geometrical setup and the evolution in the separation at different times are shown in This is again perhaps counter-intuitive, since one would usually think about the EE as a more 'nonlocal' quantity than the 2-point functions, and hence probing deeper in the bulk. Indeed, this is the case for pure AdS and also for thermal AdS, but in this case for large  enough separations the 2-point function at the same time and length probes deeper in the bulk than the EE. Of course our simulations only probed a limited set of times and lengths for our extremal surfaces and hence we cannot make a general statement if the EE never probes beyond the apparent horizon in geometries produced by shock wave collisions. Nevertheless, we think we have strong evidence that this is so, mainly since increasing the lengths at our chosen times clearly moves the tip of the surface along the horizon. We furthermore checked that extremal surfaces centered around y = 0 behave similarly, so that the property is not due to our symmetric set-up.

Conclusions
In the paper at hand we studied the time evolution of equal time 2-point functions and HEE in strongly coupled anisotropic and inhomogeneous N = 4 super Yang Mills theory via its dual description. In the dual description this amounts to finding geodesics and extremal surfaces in the gravitational background of two colliding gravitational shock waves. We used three different initial conditions, corresponding to wide, intermediate and narrow shocks.
When the separation is held fixed the 2-point functions decrease before and increase after the collision. During the collision new correlations form such that the system becomes more correlated than in the beginning. The narrower the shocks the higher the gain in correlations before they reach their final value.
We also studied the correlation between the two shocks itself by following the maximum of the energy density. In this case the correlation between the two shocks increases linearly before the collision. After the collision correlations decrease for wide shocks, whereas for the narrower shocks they continue to grow before they fall off again.
The time evolution of the EE can be divided into four regimes, namely highly efficient rapid initial growth, linear growth, post collisional regime and late time fall off. The smaller the shocks the more rapid the initial growth, reflecting the fact that the rate at which the energy density enters the entangling region is larger for smaller shocks. The post collisional regime is qualitatively different for the different initial conditions. As the shocks get smaller an additional minimum appears which we attribute to the fact that the longitudinal pressure becomes negative. The existence or absence of a minimum in EE in the post collisional regime thus serves as an order parameter to discriminate between the transparency (narrow shocks) and full-stopping (wide shocks) scenarios. At late times we observe polynomial fall off behavior where the exponent depends on the initial conditions. Surprisingly, we found that 2-point functions can probe behind the horizon, but only after the system has hydrodynamized. In contrast, the EE surface did not probe behind the horizon in our simulations, which is perhaps counter-intuitive.
This finding has to be contrasted to the observations made in [21], where the authors studied the holographic entanglement entropy in Vaidya AdS 3 and found geodesics which cross the apparent horizon. In AdS 3 /CF T 2 , however, the holographic entanglement entropy and the two-point function are equivalent, whereas in our AdS 5 they have manifestly different behavior.
An interesting application of our results is to check numerically the quantum null energy condition [47][48][49] in a regime where the classical null energy condition breaks down. Namely, for the narrow shock-waves shortly after the collision there are regions where the classical null energy condition fails. We intend to perform this check in future work using the results for HEE established in the present work.
An interesting generalization of our results could be the consideration of shock-wave collisions in non-conformal theories, holographically modeled by the addition of a scalar field with judiciously chosen self-interactions [50,51].

B Numerical checks
In any numerical analysis it is important to check the underlying algorithm for programming mistakes and to track numerical errors. In order to check the correctness of our numerical results two completely independent relaxation codes were developed, one by the Vienna group and another one by Wilke van der Schee. The first algorithm employs first order finite differences, the second one a spectral method. We find excellent agreement (see Fig. IX).
In both computer codes the embedding functions of the geodesics are represented on a finite number of grid points. The numerical result must converge to the true solution when the number of gridpoints is increased. Table 3     Our numerical scheme employs a cutoff z cut in the holographic coordinate. The final result for our observables should not depend on this cutoff which purely serves numerical purposes. In Table 4 we show the scaling of the 2-point function of separation µl = 1 and the EE of separation µl = 0.5 evaluated at two different times (µt = 0, 2) for the narrow shocks. The results for the 2-point function are nicely independent of the cutoff in the range z cut ∈ [0.01, 0.1]. In case of the EE the cutoff dependence turns to be ≈ 1% in the range z cut ∈ [0.05, 0.1] which is sufficient for our qualitative studies. In all our simulations presented in this work we have set the cutoff to z cut = 0.075.