Time evolution of entanglement for holographic steady state formation

Within gauge/gravity duality, we consider the local quench-like time evolution obtained by joining two 1+1-dimensional heat baths at different temperatures at time t=0. A steady state forms and expands in space. For the 2+1-dimensional gravity dual, we find that the shockwaves expanding the steady-state region are of spacelike nature in the bulk despite being null at the boundary. However, they do not transport information. Moreover, by adapting the time-dependent Hubeny-Rangamani-Takayanagi prescription, we holographically calculate the entanglement entropy and also the mutual information for different entangling regions. For general temperatures, we find that the entanglement entropy increase rate satisfies the same bound as in the"entanglement tsunami"setups. For small temperatures of the two baths, we derive an analytical formula for the time dependence of the entanglement entropy. This replaces the entanglement tsunami-like behaviour seen for high temperatures. Finally, we check that strong subadditivity holds in this time-dependent system, as well as further more general entanglement inequalities for five or more regions recently derived for the static case.

. At t = 0, both systems are isolated and independently at equilibrium. Evolving forward in time from t = 0, a spatially homogeneous non-equilibrium steady state develops in the middle region, carrying an energy current J E . J e ) = 0, th e steady state, develops. W ithin field theory, this was discussed by B ernard and Doyon in [35,36,38]. In their work, they showed th a t for late tim es, th e steady-state region can asym ptotically be described by a therm al d istribution a t shifted tem perature. Such a local quench-like system can be m odeled for exam ple as in [25,40], where two independently therm alised B C F T s are joined at t = 0. In this work, in contrast, we will follow [41] and utilise a different setup where equation ( 1.1) holds a t t = 0, and a steady sta te forms when evolving th e sta te b o th forwards and backwards in tim e. However, we will only consider the regime t > 0 as physically interesting.
In principle, th e tim e-dependent system described above can be set up and studied in a rb itra ry dimensions. However, in 1+1 dim ensions th e num erical analysis can be sup plem ented w ith analytical results, due to th e fact th a t in this case conform al field theory techniques can be applied. From th e holographic perspective, it is relevant th a t 2+1dim ensional gravity is non-dynam ical. We study th e 1+1-dim ensional case and its gravity dual in th e present paper. An im p o rtan t difference betw een th e 1+1-dim ensional and the higher-dim ensional case is given by th e following: in 1+1 dimensions, th e shock waves w ith which th e steady sta te region expands are dissipation-free. B oth in th e field theory and in th e gravity dual, for all tim es the tran sitio n betw een th e heat b ath s and th e steady-state region is described by a step function. On th e other hand, in th e higher-dim ensional case th e shock waves experience diffusion and th e tem p e ra tu re profile is no longer described by a step function. This m ay be referred to as a rarefaction wave.
A m ain focus of th e present paper is th e study of th e tim e dependence of entangle m ent entropy and of m utual inform ation in th e steady-state system described above. In particular, our analysis describes how these quantities change as th e shock wave moves th rough th e chosen entangling region, for instance a spacelike interval of length £ located away from x = 0.
To our knowledge, th e tim e evolution of th e entanglem ent entropy in this setup has not been studied yet. For our analysis we take th e holographic perspective, which allows us to m ake use of th e prescriptions proposed in [42][43][44]. T he original Ryu-Takayanagi prescription states th a t in a d-dim ensional C F T th e entanglem ent entropy of any region A is proportional to th e area of th e m inim al codim ension-tw o surface in th e d+1-dim ensional dual static geom etry anchored a t th e boundary of A. L ater this prescription was generalized to the tim e-dependent case, in which the entanglem ent entropy of A is proportional to the extrem ised spacelike codim ension-tw o surface in th e tim e-dependent bulk geometry. For our 1+1-dim ensional boundary setup, the extrem al surfaces we are looking for are geodesics.

JHEP10(2017)034
A d S /C F T relates therm al states to statio n ary AdS black holes. Away from equilib rium , th e tim e evolution of th e field-theory system corresponds to th e evolution of the spacetim e dynam ics subject to appropriate boundary conditions a t th e asym ptotic AdS boundary. T he holographic dual of th e initial sta te w ith th e tem p e ra tu re profile ( 1.1) is thus given by a geom etry consisting of two AdS black holes a t two different tem peratures b o th cut a t x = 0, and a half of each glued together a t t = 0. For this particu lar scenario, th e asym ptotic late-tim e geom etry is known and the steady-state region was shown to be equivalent to a boosted AdS Schwarzschild black hole a t tem p e ra tu re a / T l TR [41]. This result is in agreem ent w ith th e earlier C F T result.
Taking th e holographic approach, according to [41], th e steady-state region in the late-tim e lim it can be described by a boosted therm al sta te in th e higher-dim ensional case as well if th e system shows tim e-independent asym ptotic behaviour. T he corresponding argum ent is based on black hole uniqueness theorem s. A num erical analysis of th e 2+1 dim ensional boundary C F T [45] shows very good agreem ent w ith th e conjecture. However, while in [41] th e two wavefronts in th e higher-dim ensional case are b o th shockwaves, it was later shown in [46] and [47] th a t th e solution w ith one shockwave and one rarefaction wave is preferred.
H ydrodynam ics provides another fruitful approach to study th e tim e evolution of sys tem s subject to a local quench like ( 1.1) . Studying th e hydrodynam ic expansion to first order, th e authors of [37] conjectured th e universality of the steady sta te regime in a sense th a t its emergence a t late tim es is universal and irrespective of th e dynam ical details of th e system or details of the initial sta te configuration and th a t th e heat current can be described w ith a universal formula. T he assum ptions on which the conjecture is based are sim ilar to th e ones in [41] nam ely th a t a t late tim es th e system can be described by three regions, th e two heat reservoirs and a steady sta te regime growing in tim e as two shockwaves move tow ards spatial infinity.
A free field analysis w ithin K lein-G ordon theory in [48] showed th a t in contrast to the 1+ 1-dim ensional case in m ore dim ensions the emerging steady sta te is different from its strongly coupled analogue. A recent review [49] on quantum quenches in 1+1 dim ensional conform al theories discussed global quenches a t finite tem p eratu re and local quenches at zero tem perature.
M uch interest has also been directed tow ards th e holographic study of th e tim edependent behaviour of entanglem ent entropy after global quenches. For exam ple, in [8 14, 16, 50] it was found th a t after a global quench, an initial quadratic grow th of entangle m ent w ith tim e is followed by a universal linear grow th regime. T he special case where the final sta te is an AdS-Schwarzschild black hole is referred to as entanglem ent tsunam i [12 14]. Noteworthy, in [14] th e linear grow th is found to be independent of th e choice of entangling regions. It is also interesting to note th a t th e cases studied in [51,52] display a linear grow th independent of th e equation of state, showing m ore evidence th a t points tow ards a universal behavior. R elated work on th e evolution of entanglem ent entropy after a local quench also include [17][18][19][20][21][22][23][24][25]. In m ost of these references th e authors com plem ented th e holographic analysis w ith results from a C F T analysis. In particu lar in [21], th e authors consider a local quench w ith a small tim e w idth e and find universal features of th e tim e evolution of th e entanglem ent entropy.

JHEP10(2017)034
A related, b u t distinct line of research involves th e deform ation of strongly coupled m a tte r by tim e dependent p e rtu rb atio n s of a relevant scalar operator. N um erical investiga tions into this situation [53] involve an uncharged black brane solution which is pertu rb ed by varying th e non-norm alizable boundary m ode of a massive bulk scalar in tim e. One of th e m ost interesting results to emerge from such a study has been th e appearance of a "universal fast quench regim e" in which the change in energy density after th e quench scaled as a power law in th e quench w idth.
T he first focus of this paper (section 2) is to investigate th e tim e evolution of the steady sta te itself. We analyze the causal stru ctu re of th e dual geom etry and find th a t the hypersurfaces th a t, in th e chosen coordinate system , extend th e boundary shockwaves into th e bulk are spacelike. Therefore they ap p ear to be superlum inal. However, as we explain, causality is not violated.
In sections 3 and 4, we th en num erically com pute th e tim e evolution of th e en ta n glem ent entropy for a spacelike interval which at t = 0 is entirely enclosed in one of the heat baths, say th e left one. As we work w ith a 1+1 dim ensional boundary th e interval is one-dim ensional. T he entanglem ent entropy of such an interval quantifies th e quantum entanglem ent betw een th e interval and its com plem ent. Let us describe w hat happens when th e shockwave passes th e interval. Before th e shockwave propagating outw ards en ters th e interval on its right edge, th e entanglem ent entropy is constant in tim e. T he same is tru e once th e shockwave has left th e left edge of th e interval behind. W hile th e shock wave is passing through th e interval, we observe a universal tim e evolution: th e functional dependence on th e interval length and th e two tem p eratu res TL,R is th e same for a wide range of tem p eratu res and intervals, as long as th e tem p eratu res and th eir difference are small com pared to th e inverse of th e interval length considered. In section 5, we present an analytic proof for th e universal behaviour. For larger tem p eratu res and tem p eratu re differences there are deviations from this universality which we see b o th in our num erical result and our analytical com putation. For th e la tte r we give an estim ate.
In section 3, we also study the tim e dependence of th e m utual inform ation for which we consider equal length intervals a t equal distance from x = 0. T he m utual inform a tion quantifies th e am ount of inform ation obtained about th e degrees of freedom in the one interval from the degrees of freedom in th e other. We find th a t m utual inform ation for th e geom etry described above grows m onotonically in tim e. Furtherm ore, we look at an interval initially located w ithin th e sm aller tem p e ra tu re heat b ath . Its entanglem ent entropy increases w ith tim e. In section 5 we show analytically th a t its average increase rate is bounded. This is sim ilar to w hat is observed for th e entanglem ent tsunam i. A fu rth er analytically tra c ta b le case is when one of th e tem p eratu res is zero and th e other tem p e ra tu re is large com pared to th e inverse of th e interval length. For this case we show th a t the entanglem ent entropy grows linearly in tim e.
O ur second m ain focus, considered in section 6, concerns entanglem ent inequalities for n disconnected intervals. A fam ous exam ple is th e strong subadditivity for n = 3 inter vals. In this paper we num erically study generalized entanglem ent inequalities introduced in [54] for n = 5 intervals. These authors proved these inequalities in th e static case. We num erically verify by considering a large num ber of exam ples th a t these inequalities also T he advantage of th e m atching m ethod is th a t it allows us to study a wider spectrum of tem p eratu res com pared to th e shooting m ethod. N ote th a t in this paper we only consider intervals th a t experience at m ost two of th e three regions, th e two heat b ath s and the steady sta te region. T his paper is organised as follows. In section 2 we describe th e holographic ansatz we work w ith and discuss th e causal stru c tu re of th e geom etry considered. In section 3 and 4 we explain th e two com plem entary num erical m ethods, shooting and m atching, in detail and present th e consistent results on th e tim e evolution of th e entanglem ent entropy and th e m utual inform ation. A nalytical results are presented in section 5. In 5.2 we analytically prove th e universal behaviour of th e tim e dependence of th e entanglem ent entropy. We present analytical results for th e special case in which one of th e heat baths is a t zero tem p eratu re and discuss bounds on th e entropy increase rate in sections 5.3 and 5.4. Section 6 is devoted to th e study th e entanglem ent entropy of setups w ith m any disconnected intervals. An algorithm for dealing w ith th e large num ber of configurations is introduced and subsequently used to explore entanglem ent inequalities. In section 7 we present some analytical results on th e higher dim ensional case and com m ent on the challenges of th e num erical approach in this case. We conclude in section 8.

H olograp h ic setu p
We are interested in studying a strongly coupled C F T in d -1 = 1 spatial dimensions. T he energy flow in such a system can be qualitatively captured by pure gravity alone in holography, so the bulk action th a t we will consider is simply

JHEP10(2017)034
where A = -1 /L 2 is th e cosmological constant of AdS. A static C F T configuration at finite tem p e ra tu re T is dual to the BTZ black hole [55,56], (2 .2) where L is th e radius of AdS spacetim e, which we will set to L -1 later, and th e tem p er atu re is related to the horizon's position zH via 1 /T -2nzH . We will always assum e the spatial coordinate x to be decom pactified, such th a t -to < x < + to . In this geometry, th e calculation for th e entanglem ent entropy can be easily derived from th e fact th a t the BTZ black hole is obtained from a quotient of pure A dS3 [44]. Given a spatial interval w ith separation £ in th e C F T , the holographic result for th e entropy of th e entanglem ent betw een this region and its com plem ent is given by [42] S B T Z = 4g log (n2£2T2 sinhV £T^ , where e is a UV cut-off.
Using m inim al subtraction, this result m ay be regularized to read s BTz = 4 -g log ( n f 2 sinh2(n£T^ . (2.4) In this paper we study th e particu lar dynam ical configuration investigated in [41]. We consider two therm al reservoirs, each of them initially a t equilibrium b u t a t different tem peratures, TL and TR. A fter bringing th e two system s into therm al contact a t t = 0, a spatially homogeneous steady sta te develops, carrying a heat flow J E which transfers energy from th e heat b a th at higher tem p e ra tu re to th e other. This physical situation is presented in figure 1. T he steady sta te configuration in th e C F T is described by a L orentz-boosted stress-energy tensor, which is dual to a boosted black hole geometry, d z 2 + 0 . 0 + (dx cosh 9 -dt sinh 9)2 1 -z 2/zH (2.5) This is dual to a boosted therm al sta te w ith boost param eter 9, tem p e ra tu re T and velocity 13, which are determ ined by This is a particu lar case of equations (7.5) for d = 2. Its associated entanglem ent entropy is given by 1 Sboost = 4g log f n 2f^T^e2 sinh(n£TL) sin h (n £ T R^ . (2.7) 1By carrying out the boost, this can be obtained from th e entanglem ent entropy for a static black hole for boundary intervals w ith endpoints at different (boundary) tim es t1 = t2.
This result also gives the late-tim e lim it of our case, since th e central region expands progressively as th e shockwaves advance tow ards th e heat reservoirs located a t spatial infinity. As w ith (2.3) , it m ust be renorm alized by su btracting (L /4 G )lo g e -2 , according to our scheme of m inim al subtraction.
For th e case d = 2, th e shockwaves move w ith th e sam e speed u = 1 in b o th directions, so generically, at a tim e t > 0, th e geom etry is divided into three regions. Schematically, (2 .8) Such a dynam ical solution corresponds to th e idealized lim it in which th e initial tem p era tu re profile of th e system includes a step function of zero w idth, leading to sharp shockwaves in th e C F T . In this lim it, there are three different regions, th e central one corresponding to th e steady state, which is form ed by th e propagating shockwaves. N ote th a t this simple solution only applies to th e (1+1)-dim ensional case. In a generic num ber of dimensions, th e dynam ics is non-linear and th e n atu re of th e right and left shockwaves is very different. See section 7 and [41] for a discussion of th e higher-dim ensional case.
Given a generic sm ooth tem p eratu re profile of finite w idth, it is convenient to work in Fefferm an-G raham coordinates (z ,t, X). T he dynam ical solution in these coordinates can be found in references [41,57]. It m ay be w ritten as T he functions f L( x + t) and f R(x-1 ) are to be determ ined by th e boundary conditions. The way to do this is to calculate th e boundary stress-energy tensor. Its vacuum expectation value is given by 2 2These are expectation values of th e boundary stress-energy tensor. T he gravitational solution in th e bulk is a vacuum solution.
where c is th e central charge of th e C FT . T he initial condition ( f tx) = 0 (i.e. th e absence of a heat current a t t = 0) dem ands th a t f L(v) = f R (v). In the following we will consider a profile n 2-2 , , fL(v) = f R (v) = -( (f L + f R ) + (TR -TL) ta n h ( a v )) , (2.12) which corresponds to a step of w idth w & 1 /(2 a ). In the lim it a ^ to , it asym ptotes to a sharp step function, n 2-2 f L/R(v) ^ 2 (TL + (TR -TL) 9(v)) . (2.13) T he discontinuous m etric in this case is given by on th e left (L) and right (R) sides respectively, and in th e central region.3 T he relation betw een this solution and th e equivalent in Schwarzschild-type coordinates is discussed in section 4.
It is now w orth to point out th a t this discontinuous geom etry is formed by gluing different spacetim es together along co-dim ension one hypersurfaces which represent the extension of th e shockwaves from th e boundary into th e bulk. T he procedure by which spacetim es are m atched to one another in G R involves th e use of Israel junction condi tions [58]. Generically, each chunk of spacetim e ends in a codim ension one hypersurface, and when two of these hypersurfaces E 1, E 2 are identified, they m ust have th e sam e topol ogy and induced m etric Y ij. T he identification is generally only possible provided th a t the energy-m om entum tensor S ij , defined on th e hypersurface E 1 = S 2 which glues regions of spacetim e together, satisfies ( K + -YijK + ) -( K --YijK -) = -KSij . (2.16) Note, however, th a t Sij vanishes for our case, as th e bulk solution is supposed to be a vacuum solution everywhere. In these equations, K + , K -are extrinsic curvatures of the hypersurface, com puted from th e induced m etric on th e right and left sides respectively. They correspond to different em beddings, since this hypersurface is em bedded from both sides. We checked th a t this condition is satisfied for th e present geometry. T here is, however, a non-trivial conceptual question, since th e spacetim es th a t are being glued in clude a horizon which, apparently, is cut into three pieces. In order to visualize this, it 3This shows that when setting TL = TR, the bulk metric will just be a static BTZ black hole, and there will be no non-trivial time evolution of entanglement entropy. This distinguishes our setup from the one studied in [25,40], where two BCFTs are joined at t = 0, and non-trivial time evolution takes place even when the temperatures of both sides where equal. is useful to employ a causal K ruskal diagram of th e spacetim e [56]. This will also allow us to understand and in terpret th e fact th a t in th e bulk, th e "shockwaves" form spacelike hypersurfaces. Of course, this m eans th a t while on th e bounday, via th e A d S /C F T cor respondence, we describe actual shockwaves in th e C F T , th e spacelike hypersurfaces th a t extend these shockwaves from th e boundary into th e bulk in our dual AdS picture should not be referred to as genuine shockwaves from a bulk perspective. Nevertheless, for the sake of brevity we will from now on leave away inverted comm as and also refer to the bulk hypersurfaces as shockwaves.
T he bulk spacetim e has two spatial coordinates z,x . In order to obtain a Kruskal d iag ram ,4 we com pactify z and th e tim e t for each slice of th e x direction, thus obtaining th e K ruskal coordinates ( R ,T ), defined by Let us now analyze the resulting diagram of figure 2, in order to understand how th e space like bulk shockwaves are still in agreem ent w ith causality of the bulk. A two-dim ensional K ruskal diagram is included in th e b ottom left corner, it shows a two-dim ensional space tim e divided into four quadrants: th e external universe is captured by th e right quadrant, th e interior of th e black hole by th e q u ad ran t at th e top, and th e other two correspond to th e respective analytical continuations. T he lines of constant z correspond to hyperbolae, which get closer to th e horizon w ith increasing value of z , to the extrem e th a t the line corresponding to z = zH degenerates to th e diagonals th a t separate th e interior and the exterior of th e black hole. T he lines of constant t correspond to straight lines em anating from th e center, which at t = 0 appear as horizontal in the right q u ad ran t and as vertical in th e to p q u adrant. They get closer to th e future horizon w ith increasing value of t .
Focusing on th e exterior of the black hole in figure 2, th e shockwaves leave th e central point x = 0 a t t = 0. Therefore th e initial location of th e shockwave is m arked by the horizontal ray t = x = 0, 0 < z < zH . Any other boundary point x = xo experiences th e shockwave crossing a t t = |x0| (and this extends radially all along z), so th e location of th e shockwave is m arked by a straight line th a t increasingly separates from th e initial horizontal line as |x0| increases. G athering th e locations corresponding to a shockwave for all values of x , we obtain th e green surface in figure 2. This figure displays th a t in this causal diagram , the shockwave worldvolume does not touch th e horizon surface, except at th e line corresponding to T = R = 0, which is th e bifurcation surface. Consequently, the only p a rt of th e event horizon of the static regions on th e left and on the right th a t is retained in this construction is a half of th e bifurcation surface for each side. A part from th a t, only th e event horizon of th e steady-state region appears, it rem ains untouched by th e shockwaves on which th e gluing of spacetim es takes place.
As m entioned above, th e analysis of th e causal diagram in figure 2 reveals another im p o rtan t aspect of th e shockwaves: their spacelike nature, i.e. th e fact th a t their induced m etric will have positive determ inant. Intuitively, this m eans th a t in th e bulk, they would 4T he global extensions of m etrics of the form (2.9) have been studied in [59] in more generality.
F ig u re 2. Kruskal diagram of the bulk spacetime. The black diagonal sheets correspond to the location of the black hole's horizon. The red surfaces show the location of the singularity, and the orange surfaces show the location of the boundary at z = 0. The green surface is the worldvolume of the bulk shockwaves along which the three regions of spacetime are glued together. A steady state region forms both for t > 0 and for t < 0. However, note that the physically relevant part of this spacetime for our investigation is assumed to be t > 0 only. be perceived as being superlum inal. Of course, this raises puzzling questions concerning w hether this system should be considered to be physical or not. However, we note th a t th e present geom etry is a solution of vacuum in three dimensions, in which general rela tiv ity does not have propagating degrees of freedom. Therefore, these shockwaves do not tra n sp o rt inform ation in th e bulk, and every bulk observer will locally observe AdS space everywhere. However as we will see from th e tim e dependence of entanglem ent entropy later on, from th e dual C F T perspective th e shockwaves, which travel at th e speed of light on th e boundary, do tra n sp o rt inform ation. Considerations of energy conditions in th e bulk are also unnecessary, given th a t it is a vacuum solution (including th e fact th a t S ij = 0 in (2.16)), so m ost comm on energy conditions are trivially satisfied.
T he intuitive picture of why inform ation is not tra n sp o rte d by these shockwaves in the bulk can be understood by taking into account th a t apparent faster th a n light propagation is present in m any physical situations. In order to illustrate this point, we look a t the exam ple of two rulers, as in figure 3. T he red dot corresponds to th eir point of intersec tion. As th e rulers move in diverging directions (at speeds slower th a n light), as indicated by the arrows of th e figure, th eir point of intersection moves forward a t a higher speed.
T he velocity of this point depends on th e angle betw een th e two rulers, and it can move superlum inally if th e initial conditions conspire accordingly.
However, this point of intersection is an em ergent object and not a physical one, so it does not carry inform ation. As a consequence, causality is not violated. In o ther words, consecutive positions of th a t point are not causally related, even though th e inform ation about these events is encoded in th e initial conditions. Similarly, th e shockwaves of our system have dynam ics encoded in th e initial conditions and can develop apparent superlu m inal speeds, but since they do not carry inform ation, causality is preserved.
A nother particu lar feature present in th e 2 + 1-dimensional case is th a t spacetim e is always locally AdS, even at th e m atching surface. This m eans th a t a local observer traveling w ith th e shockwave still sees AdS everywhere. This is why th e velocity and features of the shockwave cannot be physically relevant in th e bulk.
In th e following we holographically study th e evolution of entanglem ent entropy in this setup. We find th a t it has physical behaviour in agreem ent w ith field-theory expectations. For instance, there is a well-defined velocity of propagation for entanglem ent entropy. This is in agreem ent w ith argum ents using a quasiparticle picture [60], according to which the initial condition acts as a source of pairs of quasiparticle excitations. As they propagate causally th roughout th e system , larger regions get entangled. In this picture, if a m axim um quasiparticle velocity exists, th en th e entanglem ent entropy grows linearly in tim e for cer ta in boundary regions. We will see th a t also holographically, there is a velocity associated to entanglem ent, independently of th e gluing of th e spacetim es. In fact, in th e following section we will see th a t entanglem ent entropy does evolve in a causal m anner, and obeys th e velocity bound known from th e study of entanglem ent tsunam is.

N u m erica l r esu lts I: sh o o tin g m eth o d
We now tu rn to th e num erical com putation of entanglem ent entropy in th e background introduced in th e preceding section. For this purpose, we study m inim al surfaces whose boundary at z = 0 is in x = x L and x = x R, and consider space-like intervals w ith -1 2 -

JHEP10(2017)034
t(x L ) = t( x R ) . 5 By th e HRT prescription in 2 + 1 dim ensions [44], the m inim al surface com patible w ith these boundary conditions corresponds to a geodesic in the bulk. This follows from a solution of th e geodesic equations, which read We assum e an affine param etrization of th e geodesic, which implies T he induced m etric on th e m inim al surface reads In th e following we will com pute renorm alized entropies according to this formula.   th en find th e values of these initial conditions th a t lead to th e desired boundary values at s ^ sl ,r given in (3.7) . 6 We introduce a cutoff e ^ 1 for regularization. This also induces a cutoff in th e affine param eter, i.e. sL ~ -| log e| and sR ~ | log e|. In the following we will consider space-like intervals A and B as shown in figure 4. Figures 5 and 6 display a typical solution of the geodesic equations, which satisfies th e boundary conditions of (3.7) . Once th e geodesics are obtained, th e next step is to com pute th e area of these curves and th en th e entanglem ent entropy from (3.5) . We now present results for th e tim e evolution of th e entanglem ent entropy in the system of section 2 .
6T here are in th e literatu re other num erical m ethods for th e solution of th is tw o-point b oundary value problem . A n exam ple is given by th e relaxation m ethods, in which th e solution is determ ined by startin g w ith an initial guess and im proving it iteratively, see e.g. [61]. 3.2 E n t a n g le m e n t e n tr o p y a n d u n iv e r s a l t im e e v o lu tio n For th e m om ent we consider a single interval x € [xL , x r ] denoted by A, placed in the positive semiplane, i.e. x L,R > 0. Let us study the tim e evolution of th e entanglem ent entropy S a during th e form ation of th e steady state. We are considering the m odel in d = 2, so th a t the shockwaves are at t = |x |. This m eans th a t th e shockwaves touch the two ends of th e interval at tim es t = |x L| and t = |x R |, see figure 4. We will denote these values by ti and t 2, respectively. In th e lim it a ^ to in (2.12) , th e entanglem ent entropy tu rn s out to be constant in th e regimes 0 < t < t 1 and t 2 < t, and there is a non-trivial tim e evolution only in th e interval t 1 < t < t 2, i.e.
We display in figure 7 (left) the tim e evolution of th e entanglem ent entropy of interval A of figure 4, from a num erical com putation of the geodesic equations. In this and subsequent figures we display the entanglem ent entropy renorm alized as shown in equation (3.5) . Let us focus on th e regime t 1 < t < t 2. It is convenient to define th e norm alized entanglem ent entropy / a (p) as / a (p ) = S SA(t) ) S a (? =t 0 ) 0) w ith P = (t -t 1) /^t , (3.9) Sa (t = to) -SA(t = 0) where A t = I = |x R -x L|. This corresponds to th e function SA(t) norm alized to the interval [0,1] in b o th horizontal and vertical axes for t 1 < t < t 2. It is clear from equa tions (3.8) and (3.9) th a t /a ( p ) has the values /a ( 0 ) = 0 and / a (1) = 1. We have com puted num erically th e entanglem ent entropy SA(t) in num ber of configurations w ith different tem peratures Tl , Tr and lengths ^, and find th a t for a large range of tem p e ra tu re differences, th e behavior of / a (p ) m ay be approxim ated by Specifically, this function fits extrem ely well th e num erical results of th e entanglem ent entropies as long as £TL < 1 and £TR < 1. This is illustrated in figure 7 (right) for a p articu lar case. T he result of equation (3.10) is independent of th e values of th e param eters TL, Tr and £, and so it implies th e existence of an 'almost' universal tim e-evolution of entanglem ent entropy in the theory w ith d = 2 at small tem peratures. Eq. (3.10) will be proven analytically in section 5.2 w ithin a small tem p e ra tu re expansion.
T he analysis presented above applies also to intervals in th e negative sem iplane. We show in figure 7 (left) th e entanglem ent entropy of interval B of figure 4 . N ote th a t b oth functions, S a (t) and S b (t), tend to th e sam e value when th e intervals reach th e steadysta te regime.

T im e e v o lu tio n o f m u tu a l in fo r m a tio n
A quan tity of interest related to th e entanglem ent entropy is th e mutual information. It m easures which inform ation of subsystem A is contained in subsystem B , or in other words th e am ount of inform ation th a t can be obtained from one of th e subsystem s by looking at th e other one. An advantage of this q u an tity is th a t it is finite, so th a t it does not need to be regularized. It is defined as 7An analytical guess for the time evolution of the mutual information in analogy with the universal formula of equation (3.10) turns out to be more complicated than in previous section, due to the structure of the term S (A U B).

JHEP10(2017)034
and in this case th e entanglem ent entropy for an interval [xL, x R] w ith x L < 0 and x R > 0 becomes, using m inim al subtraction, S(Tl ,xl ; Tr , x r ) (3.15) To obtain this expression we considered the length of a curve piecewise defined in th e two heat baths, consisting of two pieces a t x < 0 and x > 0 glued together at x = 0. They are parts of two geodesics, defined in a geom etry w ith a black hole a t tem p eratu re TL and Tr , respectively. T he curve is chosen such th a t its length is m inim al w ith respect to the value of th e radial coordinate at which th e two geodesics m eet a t x = 0. For a sim ilar m atching-based m ethod see section 4. If we place th e interval in ju st one semiplane, i.e. x L,R > 0 (or x L,R < 0), th e entangle m ent entropy a t t = 0 corresponds to th e one for a statio n ary black hole a t tem p eratu re T , which reads In this equation T = TL (or TR) when x L,R < 0 (or x L,R > 0). In th e other extrem e, t --to , th e system is in the steady-state regime, and th e entanglem ent entropy is th e one for a boosted black hole,8 These analytical results, equations (3.16) and (3.17) , correspond to Sa (t = 0) and S^(t = to ) in equation (3.8) , respectively. From these form ulae we easily obtain th e property where we consider intervals A w ith x A R > 0, and B w ith x f R < 0, and lengths £ = £a = £b . This property is non-trivial, as in th e left-hand side of equation (3.18) there is th e contribution of statio n ary black hole solutions a t tem p eratu res TL and TR, while in th e right-hand side there is a boosted black hole and th e corresponding energy flow contributes as well to th e entanglem ent entropy. This relation is significant as it implies th e 'conservation' of entanglem ent entropies betw een t = 0 and t = to . However, there is a non-trivial tim e evolution a t interm ediate tim es, as we discuss below. Interestingly, ( 3.18) has also been obtained in a slightly different setup in [40].
In figure 9 (left), the tim e evolution of Sa +B = Sa + S b is displayed. We see th a t our num erics confirm th e conservation law of equation (3.18) . In th e next subsection we will study this system in th e quenching regime, i.e. t i < t < t2 in equation (3.8) , and characterize the violations of th e entanglem ent entropy conservation in this case.

N o n -u n iv e r s a l e ffe c ts in tim e e v o lu tio n
As it is shown in figure 9 (left), we find from our num erics th a t Sa +b (t) = const in the quenching regime. This implies th a t the entanglem ent entropy is not conserved a t inter m ediate tim es. A straightforw ard com putation shows th a t these non-conservation effects are only possible if there are non-universal contributions in equation (3.10) , otherw ise this In the following we restrict to intervals A and B w ith the same length and placed sym m etrically w ith respect to x = 0, i.e. £a = £B and x AR = -xR L. W hile the function S a + B(t) has th e same value at t = 0 and t = to (see (3.18) ), we find from our num erics th a t it has a m axim um at t max ~ (ti + t 2)/2 . In order to characterize th e tim e evolution of S a + B (t), let us define the norm alized entanglem ent entropy where t 1 and A t are defined as in equation (3.9) . Finally, from a num erical com putation of / a + B(p) in a variety of intervals, we find th a t its behaviour is w ell-approxim ated by T his is illustrated in figure 9 (right) for several configurations. From a com bination of th e results in (3.10) and (3.20) , we conclude th a t for small tem peratures, th e norm alized entanglem ent entropy defined in equation (3.9) can be approxim ated by T he factor C (T L,T R,£) has a non-universal dependence on th e param eters of th e interval, so th a t A a (p ) is a non-universal contribution to / a (p ). Note, however, th a t C (T L,T R,£)

JHEP10(2017)034
does not appear to affect th e universal behavior of / a + B(p), see (3.20) . Some rem arks deserve to be m entioned: on th e one hand, A a (p ) is a correction of order O (p3), so th a t it does not jeopardize the early-tim e behavior S a (t) ~ t 2 which is present in a wide variety of system s, see e.g. [9-14, 16, 50, 62, 63]. O n th e o ther hand, th e effect of A a (p ) is extrem ely small in th e configurations we studied num erically.9 T he range of validity of equation ( 4.1 S e t u p a n d a s s u m p tio n s Let us take th e m etric in its piecewise form (2.14) , (2.15) . T he m etric is a piecewise sm ooth function of Fefferm an-G raham coordinates, denoted by z, t, X in this section. We use th e nam e region to refer to th e whole subset of our space on which th e m etric coef ficients are sm ooth, e.g steady state region (denoted Sb) is given by t > 0, |X| < t, z € (0, (^v'TLTR )-1 ) . In th e sam e m anner, th e left and right therm al regions will be denoted by L and R, respectively. T he dim ension two surface along which th e m etric is discontin uous will be referred to as the shockwave. O ur aim is to calculate (regularised) geodesics length betw een two points lying on the conform al boundary of th e space-tim e. If both endpoints belong to th e same region, th e answer is already known to be (2.4) in a therm al region and (2.7) in th e steady sta te region. If, however, th e boundary points belong to 9One can see from figure 9 (left) that in this case the peak in Sa + b (tmax) is a correction of order O(10-6) with respect to Sa +b (0), so that the order of magnitude of the non-universal contribution in equation

JHEP10(2017)034
F ig u re 10. A cartoon of a curve used in our procedure, projected onto a t, x hypersurface. Sq and R denote regions of spacetime, the steady state and right thermal region, respectively. The red line is a piecewise geodesic (it is a geodesic in any of both regions) connecting boundary endpoints and a point on a shockwave (X = t in our coordinates). Then the (renormalised) length of that geodesic is extremised with respect to the coordinates of the joining point, which yields the entanglement entropy.
different regions, finding th e solution is a m ore com plicated task. We therefore m ake some assum ptions ab o u t th e spacetim e. M ost im portantly, we assum e th a t th e coordinates t , t X cover all th e regions in a sm ooth way, and it is actually the m etric components as functions of t , t X th a t are discon tin u o u s.10 This assum ption can be m otivated by th e fact th a t our m etric (2.14) , (2.15) can be obtained as a lim it of a continuous m etric (2.10) , (2.12) where th e shockwave w idth (w) tends to zero (a goes to infinity). It is reasonable to assum e th a t taking th e lim it described does not influence the dom ain of our coordinates. T he agreem ent betw een num erical re sults from th e continuous m odel a t large a and the results of th is section shall confirm th a t th e assum ption m ade yields correct results. From th e assum ption follows in particu lar th a t curves which are continuous in our coordinates are also continuous on th e m anifold itself. T his will be essential in our calculation, since it is based on joining two sm ooth curves in a way th a t it is still continuous.

G e o d e s ic s a n d d is ta n c e
Given our setup, we need to calculate a spacelike distance betw een a given point on the boundary and an a rb itra ry point on th e shockwave.11 To achieve this, we shall follow the logic of [64] and utilise th e fact th a t every three-dim ensional asym ptotically AdS manifold th a t is a solution to E in stein 's equations is locally isom etric to A dS3. So, if we identify 10Note, however, that as the Israel junction conditions (2.16) are satisfied, the metric is continuous in a strict mathematical sense. Especially, the induced metrics on the shockwave both from the static side and from the steady state side agree.
11Of course on a Lorentzian manifold, not every point on the shockwave will be spatially separated from a fixed point on the boundary, as we shall directly see later. It is enough that there will always be a set of points that satisfy this condition -a situation that indeed occurs in our case.

JHEP10(2017)034
th e isometry, we m ay use th e ready formulm for geodesic distance in A dS3 space. It is w orth noting th a t th anks to this property special to three dimensions, we do not need to calculate th e geodesics explicitly, we ju st express th eir length as a function of coordinates of th e two points. This considerably simplifies th e calculation. However, we still have to find formulm for th e geodesic distance in our problem . T he m etric (2.14) , (2.15) is given in Fefferm an-G raham (FG) coordinates. To use th e results of [64], we need to change to Schwarzschild-type coordinates. T here is a technical difficulty in th a t step: in every region th e change of coordinates takes a different form. We denote F G -type coordinates as ( t, t, x) and th e Schwarzschild coordinates as (z, t, x). In Schwarzschild coordinates, th e m etric takes th e form (2.2) where T is a real, positive constant -an (effective) tem perature. Then, th e coordinate transform ations are obtained as follows. For th e steady sta te they read Eq. (4.2) is th e inverse relation to (4.3) . We quote b o th since there is a sign to be fixed. From (4.4) it follows th a t th e effective tem p e ra tu re from equation (4.1) for th e steady state can be expressed in term s of th e reservoirs' tem peratures as T = VTLTR. (4.5) For th e therm al regions it is sufficient to take (4.2) and set TL = TR, 6 = 0 to obtain T he distance can be expressed, following [64], as 12 12N ote th a t, on contrary to the setup therein, we do not use an analytical continuation of coefficients since we are interested in a single-sided black hole, not a double sided one which is th e situ atio n there.
(4.7) L w ith d th e geodesic distance. In term s of the Schwarzschild coordinates, th e functions X 1, X 2, T 1, and T 2 read , we need to connect these to th e Schwarzschild coordinates of th e respective patches. Let us note th a t th e condition for th e position of th e shockwave is identical in any of th e used coordinates, For simplicity, we regularise in Schwarzschild coordinates. Using th e asym ptotic approxi m ation of hyperbolic cosine by an exponential, we arrive a t th e conclusion th a t th e m inim al counter-term used in (3.5) is indeed th e proper one to regularise our length. At this point it is convenient to set th e A dS3 radius to unity, L = 1. Now, we are ready to w rite the full, renorm alised distance as a function of th e joining point on th e shockwave, dR ( z j,x ) = log [(1 + n 2TRzj2) cosh (2nTR (x -£)) -(1 -(n T Rt j ) 2) cosh (2nTR (t -x))] + log (1 + n 2T LTRŹ2) cosh (n (tT L -UTr + 2 T rx )) (4.12) + (n 2T LT Rzj -1) cosh (n (t(T L + t r ) -2TRx )) -^ log ( 16n §TLTRzj/) .

JHEP10(2017)034
In th e above, we used a shortened notation: £ = x max -x m;n, x = x^ -x m;n, t = Z -x m;n -th e tim e th a t has passed since th e shockwave entered th e boundary interval. Now, this q u an tity is to be extrem ised w ith respect to, Zj and x (which is th e same as extrem ising w .r.t. Z j). E xtrem al points are solutions to dzjdR = 0 , Ox dR = 0 . (4.13) These equations tu rn out to be fourth order polynom ial equations in Zj and non polynom ial13 equations in x. Therefore, we tu rn to num erical m ethods for solving non linear algebraic equations. N ote however th a t th e above-m entioned system can be solved analytically in certain simplifying cases, see section 5. U pon solving th e system (4.13) , we obtain th e coordinates of th e extrem um of dR, nam ely (z0, x 0). T hen th e desired entanglem ent entropy is given by 14 S = 4 dR(zo,xo). (4.14)

.3 N u m e r ic a l a lg o r ith m
To solve th e non-linear algebraic system (4.13) , we need to involve num erical analysis.
O ur solution was developed in W olfram M athem atica. All th e codes used in this section are available online as supplem entary m aterial to th e arX iv subm ission of th e paper. The algorithm consists of two steps: first, following th e idea of [65] in a slightly modified form (see [66]), we find a rough approxim ation of th e solution by plotting curves satisfying each of th e two equations in (4.13) . T hen we use coordinates of crossing points of those as a sta rtin g point for stan d ard N ew ton's solver built-in M ath em atica's F indR oot function.
To understand why such a tw o-step procedure is necessary, let us briefly discuss th e func tion (4.12) and equations (4.13) . As figure 11 shows, th e dom ain of th e distance function and th e full dom ain of our coordinates are not th e same. This is in agreem ent w ith the fact th a t on a L orentzian manifold, not every point is spacelike-separated from a given point. Then, the function going to zero and ceasing to be real signals th a t one of the boundary endpoints becomes null or tim e-like separated from th e joining point. However, since dR = lo g (...), b oth equations (4.13) have th e form d lo g (f ( ...) ) = 0, so looking for th eir solution is equivalent to solving d / ( ...) = 0 if / does not vanish on th e solution. This equivalent form is strongly preferable, given th e n atu re of th e num erical com putations, in which unnecessary divisions decrease num erical precision. On th e other hand, th e modified system of equations consists of two functions th a t are well-defined for the whole dom ain of our coordinates. Therefore, we begin our num erical approach w ith an algorithm capable of finding rough approxim ations of all solutions to th e system of equations in a given dom ain. From those solutions we choose th e ones satisfying our requirem ent th a t th e length evalu ated on solution is positive. Then, these solutions are refined by using N ew ton's m ethod th a t yields th e solution w ith th e desired accuracy, in our case fifteen digits. If more th a n one solution is regular in th e sense th a t the length is positive, th e final answer is taken to 13Even upon expressing hyperbolic functions in terms of exponentials and changing variable from x to £ = exp(Ax), the exponents of new the variable are non-integer for any choice of A.
14In that place we have already set Newton's constant of supergravity theory to 1. To justify th e exclusion of solutions w ith Zj > ZH , we have num erically tested th a t th e solutions lying below th e horizon, should they appear, are not th e physical ones. The argum ent behind this is based on th e analysis of th e K ruskal diagram of our space-tim e (see figure 2) . In short, we see th a t th e shockwave does not cross th e horizon except for th e bifurcation surface -it stays entirely in the o uter region of the black hole. This m eans th a t th e point where th e geodesic crosses the shockwave will generically be outside of the horizon (Z < ZH), and since this is a causality argum ent, this occurs b o th from th e point of view of th e static region and from th e point of view of th e steady sta te region. Indeed, we find th a t a solution to our m atching equations w ith these properties always seems to exist. Therefore, th e fact th a t we can find a solution beneath th e horizon (Zj > Zh ) is only an artefact of our choice of coordinates. On th e num erical side we allow, for test purposes, Zj to exceed the above m entioned bounds by large values (2 -3 tim es larger), and the solutions found in those regions were never chosen by th e algorithm . W ith the restricted dom ain of interest and given accuracy, our algorithm has an acceptable speed: com puting th e entanglem ent entropy for a given interval and a given boundary tim e takes roughly 0.2 seconds.
15Note that we always assumed TL > Tr .

R esults
In this way we obtain entanglem ent entropy for a wide range of tem peratures. To com pare (4.17) All our d a ta is generated using A dS3 radius as unit, L = 1 . We are also going to take various tem p e ra tu re differences. In

Figure 14.
Top left: sum of Sa and Sb for two different, yet close, values of TL, shifted by the asymptotic (t = to) value of that sum. The blue plot is the case studied in section 3, as previously the deviation from constancy is of order 10-7 while the sum is of order 10-1 . The yellow plot shows a similar quantity for temperature TL = 0.21, just slightly higher. The deviation is now order of magnitude bigger. Top right: normalised deviation /a + B for different temperatures TL. The blue curve has been already shown in figure 9 to follow the universal behaviour [4p(1 -p)]3. For much bigger temperature TL = 10, the shape of the curve changes considerably. Bottom : ratio of sum Sa + Sb and asymptotic value of that sum as a function of boundary time. The deviation of the sum from its asymptotic value reaches around 2.5%, but seems to be bounded even when one increases the temperature -compare left and right figures which are in the same scale. The bigger the temperature difference, the more the curve resembles a semi-circle. All lengths in units of AdS3 radius (L = 1).

JHEP10(2017)034 5 A n a ly tic a l resu lts
Some of the num erical results presented in th e two preceding sections m ay actually be obtained analytically, a t least in some lim its. This applies in particu lar to th e result ( 3.10) for the tim e evolution of th e entanglem ent entropy. Moreover, we derive a bound on the increase rate for the entanglem ent entropy. We begin this section in 5.1 by a review of velocity bounds previously obtained for global quenches. We th en obtain analytical results for th e tim e evolution of th e entanglem ent entropy for th e lim it of b o th heat b a th tem p er atures small in section 5.2 , and of one of the two tem p eratu res vanishing in section 5.3.
Finally we obtain a velocity bound on entanglem ent grow th in section 5.4.

R e v ie w o f u n iv e r sa l v e lo c ity b o u n d s
In th e past, much interest has been directed tow ards the holographic study of th e tim e dependent behaviour of entanglem ent entropy after global quenches. For exam ple, in a series of papers [8][9][10][11][12][13][14]50],16 it was found th a t after a global quench, th e entanglem ent entropy of a sufficiently large boundary region would exhibit an initial q u ad ratic grow th of entanglem ent w ith tim e, In this form ula, t is th e tim e after th e quench, A S is th e change in entanglem ent entropy, s eq is the entropy density of th e (late tim e) equilibrium therm al state, As is th e surface a velocity th a t depends on th e final equilibrium state. In th e case of an AdS-Schwarzschild black hole as final state, it was found th a t [8,[12][13][14] which is referred to as th e entanglement velocity or tsunami velocity. T he reason for this nom enclature is th a t th e behaviour (5.2) can be understood in term s of a heuristic picture in which th e entanglem ent grow th is caused by entangled quasi-particles th a t were created by th e global quench and are propagating a t th e speed (5.3) , form ing th e entanglement tsunami. See also [67][68][69][70][71][72][73][74] for furth er work on this topic. A related concept is th e so-called butterfly velocity [64,75] 16See also [15] for th e case of a background geom etry w ith a hyperscaling violating factor. 17£ is assum ed to be large com pared to th e inverse tem p eratu re of th e final equilibrium state. In th e case A S (t) a t 2 + ...,  for th e spatial propagation of chaotic behaviour in th e boundary theory. This speed is also connected to the grow th of operators in a therm al sta te [64,75]. From (5.3) and (5.4) , it is obvious th a t and th e case d = 2 is th e special case where 1 = vb = ve . Interestingly, th e velocities seem to play an im p o rtan t role in th e description of entanglem ent spreading not only for global quenches, but also for local quenches [20,22,25]. at order 50. Here, p is again th e rescaled tim e defined in (3.9) . It is interesting to note th a t th e small TL, TR expansions leading to (5.9) loose th eir analytic validity a t an order of m agnitude of TL,T R a t which our num erics are still well approxim ated by ( 5.9) . It m ight hence be interesting to do th e TL, TR expansions in a more system atic way and to study th e higher orders in more detail. This m ight also help to understand th e range of validity of equation (3.21) .
It is w orth stressing th a t th e universal dynam ics of entanglem ent entropy is not only of purely academ ic interest. It is known th a t th e low-energy spectrum of excitations of some models (i.e system s w ith ballistic conductance, possessing quasiparticle description, see [39]) are governed by effectively conform al theories. T he regime in which this approxi m ation is valid for therm al states is indeed when b o th of th e tem p eratu res are low, so the 18This means that in this section, we assume that TL and TR are both small (compared to the interval length £), but of the same order of magnitude.

JHEP10(2017)034
highest lying p arts of spectrum are not largely populated in th e therm al sta te -which is also th e range of validity of our universal form ula (3.10) . This m eans th a t in th e lim it of small tem peratures our universal evolution of entanglem ent entropy should be also valid in ballistic regimes of real, i.e electronic system s. It is therefore an interesting possible direction of investigation to com pute other quantities, as correlation functions, in this lowtem p e ra tu re lim it and com pare them against expectations from o ther theories (i.e. lattice m odels).

E n ta n g le m e n t en trop y: lim it o f z ero te m p e r a tu r e for o n e o f th e h e a t b a th s
In addition to th e case where both tem p eratu res TL and TR are small (studied in section 5.2) , there is another situation where th e m atching equations derived in section 4 can be solved analytically: th e case where TR = 0 w ith a rb itra ry TL (or th e analogous case TL = 0 w ith arb itra ry Tr , which we will not consider separately).
F irst of all, let us reassure ourselves th a t this case is actually physical. Setting TR ^ 0 in equations (7.1)-(7.5) has the consequences T ^ 0 (5.10) X ^ (5.11) P ^ + 1 (5.12) 9 ^ + to . (5.13) D espite th e divergence of th e rapidity 9, we see th a t for d = 2, th e line elem ent (2.5) of th e boosted black hole has a well-defined limit L 2 ds2 ^ ^ (dz2 + ( -1 + n 2T ' l z 2)dt2 + (1 + n 2T l z 2)dx2 -2n2T ' l z 2 dtdx) . (5.14) Similarly, instead of (4.12) , we find th e expression dR (zj ,x) = log 4n 3 (£2 -2£x -t 2 + 2tx + z 2) x (nTLz 2 cosh (ntTL) -(t -2x) sinh (ntTL)) which can be shown to be extrem ised for where we have used (5.21) w ith TR = 0. This result shows th a t for large TL, th e en ta n glem ent entropy will increase in an approxim ately linear way, sa tu ra tin g th e bound to be discussed in section 5.4. N ote th e discrepancy of a factor A^ = 2 betw een (5.18) and equation (5.2) for d = 2, where vE = 1. This m ay be because in a global quench, the entanglem ent tsunam i enters th e interval of interest from both sides, while in our local quench like setup, th e shockwave enters the interval only from one side. However, we should stress th a t in contrast to th e entanglem ent tsunam is studied in global quenches in [8, 12-14, 50, 67-74], we do not have a heuristic quasi-particle like picture explaining th e entropy increase and decrease experienced by different intervals in our inhomogeneous setup even qualitatively. We will discuss th e possible bounds on th e entropy increase or decrease rate of different intervals in our system for general TE,T R in section 5.4.

.4 B o u n d s on e n tr o p y in c re a se ra te
As shown in section 3, for th e full tim e-dependent background (2.9) we expect th a t for an interval of length L, th e tim e dependent entanglem ent entropy S(L, t) will be constant before th e shockwave enters th e interval, evolve w ith tim e t while th e shockwave passes th rough th e interval, and be constant again after th e shockwave has left th e interval. This is precisely the behaviour seen in figure 7, for exam ple. Furtherm ore, norm alising the entropy such as to obtain a dimensionless quantity, we have observed in section 3 th a t for small tem p eratu res T l ,T r , th e tim e dependence of th e entanglem ent entropy can be approxim ated by th e form ula (3.10) , which we analytically proved in section 5.2. In this section, we will have a closer look a t th e rates of entropy increase th a t we observe in the tim e dependent entanglem ent entropies S(L, t). We begin by analytically deriving some useful expressions. For a sharp shockwave moving a t th e speed of light, th e change in th e entanglem ent entropy of an interval w ith length 19 L will occur over a tim e period A t = L. From (2.7) and (2.4) (w ith T ^ T l for exam ple, as appropriate when th e interval is entirely on th e left) we th en easily find the average entropy increase rate th a t for fixed TL and TR, vav is bounded. By choosing TL and TR, however, we can make vav as large as we want.
In (5.2) , th e rate of entropy increase was norm alised by th e entropy density seq of the final state. Taking th e lim it L ^ to in (2.7) , we find th a t 19 Again, in this section we assume th a t th e interval is com pletely to th e left or to th e right of th e origin x = 0. (5.19) In particular, we find 0 < Iv av 1 < l L n | T R -Tl \ (5.20) w ith lim^o vav = 0 and l i m^^, vav = 4Gn |T R -Tl |. This is interesting, because it implies s eq = 4 G n ( T L + T b ) . Hence, m otivated by th e com parison w ith th e literatu re on entanglem ent tsunam is [8, 12 14, 50, 67-74], we can define th e norm alised average entanglem ent entropy increase rate and hence we find th e bound We consequently see th a t, when norm alising in a specific physical way, both th e average increase and decrease ra te of entanglem ent entropy in the form ation of th e steady sta te will be bounded by a sim ilar bound as observed for two dim ensional entanglem ent tsunam is or th e local quenches of [20]. It should be noted however th a t th e bound ( For this result, it can be analytically shown th a t for any param eters Tl ,£ and t the bound (5.25) is satisfied. This is in contrast to th e results of [16], where it was explic itly found in a different setup th a t th e m om entary increase rate for small regions, far away from th e tsunam i regime, can indeed violate the velocity bound (5.25) . See also [62,63] for fu rth er discussions of entanglem ent entropy grow th for small subsystem s in different setups.
A bound of th e type (5.23) is especially interesting when com pared to o ther velocities th a t are related to th e spread of entanglem ent or other disturbances on th e boundary of AdSd+I, such as th e entanglement velocity (5.3) and th e butterfly effect velocity (5.4) .
As said before, th e case d = 2 is th e special case where 1 = vB = vE , and hence the bound (5.23) can be expressed in term s of vE a n d /o r vB. As we will see in section 7.2, this m ay however not be th e case for higher dim ensions any more. Nevertheless, we can a tte m p t to interpret our findings for 2+ 1 bulk dim ensions in term s of the intuition provided by th e study of th e entanglem ent tsunam i phenom enon. As noted in section 5.3 in discussing th e result (5.18) , in th e lim it where th e entanglem ent entropy increases linearly (T » T -1 ) w ith a rate sa tu ra tin g th e bound (5.25) , th e shockwave seems to take th e role th a t th e entanglem ent tsunam i had for a global quench. As th e shockwave enters th e interval only from one side instead of from b o th sides, th e increase rate in ( 5.18) is only half of th e one calculated in a global quench according to (5.2) , where Ay = 2 and vE = 1. As pointed out in section 5.1, th e linear increase (5.2) is only valid when looking at large enough boundary regions (com pared to th e inverse of th e tem perature). O ur analytical results (5.18) and especially (3.10) th en show how this linear behaviour is modified when moving away from this limit: th e linear increase of entropy characteristic of th e entanglem ent tsunam i is replaced by a much sm oother S-shaped curve. This m ight, in analogy w ith th e tsunam i picture, be called an entanglem ent tide. It should be pointed out th a t in our m atching procedure of section 4 , th e shockwave is always assum ed to be infinitely thin, hence this m odification is n ot a result of a finite shockwave size. Also, other works where the evolution of entanglem ent entropy away from th e tsunam i regime was studied are [16,62,63], w ith som ew hat contrasting results, as explained above. W hile this is well-known and straightforw ard for th e n = 2 case ju st discussed, some interest has recently emerged [54,[80][81][82][83] for sim ilar concepts for situations involving n > 2 disconnected intervals. Here we present our study of this case, and apply it to the stead y -state spacetim e in subsection 6.4. T he W olfram M athem atica code th a t we use for this analysis is uploaded to th e arX iv as an ancillary file together w ith this paper and w ith a sam ple of th e num erical results th a t is obtained from th e m atching procedure of section 4.20 T here is some overlap betw een th e issues addressed in this section (especially subsection 6.1) and the ones investigated in [83], which was published after m ost of this section was com pleted. A lthough we are working in a covariant (tim e-dependent) setting, th e findings of [83] suggest th a t th e code used in our ancillary file m ay still be optim ized. However, it nevertheless produces th e desired results.

JHEP10(2017)034
As shown in figure 16, for two intervals (n = 2) there are two possible phases or configura tions describing th e entanglem ent entropy of the union of these intervals. Suppose we are given n € N intervals, how m any possible phases are there, and how do they look like? For a simplified situation where the lengths of all intervals are equal, as well as th e distances betw een them , this has already been studied in [80]. We are however interested in the general case here. Of course, for any given n, the above question can sim ply be answered by draw ing all possible configurations by hand. However, in this subsection we provide an algorithm th a t for any n enum erates th e possible phases in a consistent m anner, w ithout om itting any solution or counting it twice, and which can easily be im plem ented (see the corresponding ancillary file). We do not assum e a tran slatio n invariant spacetim e, however we will assum e a spacetim e w ith simple topology, such as Poincare AdS or a flat black brane, excluding possible phenom ena such as entanglem ent plateaux [85], see also [83]. O ur task th en essentially boils down to finding th e noncrossing partitions of a set w ith n elem ents, a well known com binatorial problem related to th e Catalan numbers Cn [86]. We will, however, still present our solution to this problem in detail, as this exposition also serves to explain our n otation and th e inner workings of our ancillary file.
In a 1 + 1 dim ensional C F T , th e n intervals under consideration (which we assum e to be all p a rt of a specified equal tim e slice on th e boundary) are all lined up one after the other, and we can enum erate th eir sta rt-and end-points from 1 to 2n, as was already done in figure 16. N ote th a t this is only an enum eration, and not m eant to indicate th e lengths of th e different intervals or th e coordinates of th e boundary points for example. Naively, in th e n = 2 case, we could have also draw n a configuration as th e one depicted in figure 17, w ith two curves crossing each other [78]. Such a configuration is, however, considered to be unphysical for various reasons. F irst, in th e static case where th e RT prescription holds, it can easily be shown th a t this type of configuration can never yield the lowest values for th e entanglem ent entropy, hence can be ignored in (6.1) [76,78]. Second, in a tim e dependent (HRT) case th e two curves m ay not actually cross any more. However, th e co-dim ension one surface spanned betw een them and th e boundary intervals would then become null or tim elike a t some point. As pointed out in [85], th e co-dim ension one surface required by th e homology condition has to be restricted to be spacelike everywhere in the HRT prescription. Hence th e configuration of figure 17 is also excluded in th e dynam ic case. T hird, it has been discussed in [78] th a t configurations of this intersecting type do not play a role w hen th e (regularised) entanglem ent entropy of an interval is m onotonously increasing w ith th e interval length.
W hen enum erating th e possible phases of th e entanglem ent entropy of th e union of n intervals, we therefore aim a t excluding phases w ith curves intersecting when projected into th e sam e plane, as shown in figure 17. Due to our labeling of th e boundary points, it is clear th a t each interval begins a t a point labeled by an odd num ber and ends a t an even one. In figure 16, for n = 2 we find one phase where bulk curves connect th e points 1 to 2 and 3 to 4, and one phase where th e bulk curves connect th e points 1 to 4 and 2 to 3. In th e unphysical case shown in figure 17 however, th e points 1 to 3 and 2 to 4 are connected by bulk curves. We hence realize th a t in order to avoid intersections as th e one shown in 17, each odd label has to be connected to an even label by a bulk curve. This m eans th a t each viable configuration for any n has to be a m apping of th e set of odd num bers where e.g. 1 ^ 2 stands for the curve connecting th e points 1 and 2 .

T h e p h a se s o f th e u n io n o f n in terv a ls
All th e possible phases obtained this way for n = 3 are shown in figure 18. Clearly, there are 3! = 6 of them , however we see th a t there is still one involving intersections. Of course, when sketching these six possible configurations by hand, it is easy to identify th e one involving intersections and to discard it. However, from our point of view of autom atizing this process, we need to form ulate and im plem ent th e criterion th a t distinguishes the unphysical phase 3 ^ ai(4) phase 1: 5 ^ 6 , ..., phase i: 5 ^ 0^ (6) ,...  Figure 18. 3! = 6 preliminary configurations for n = 3. Note that the 4th configuration is likely unphysical. According to the nomenclature of [80], the phase 6 is referred to as engulfed phase. (6 .8) To do so, we exclude all configurations in which there are two intervals spanned by the endpoints of bulk curves th a t intersect in such a way th a t th e intersection is an interval th a t is not either em pty or one of th e intervals spanned by th e bulk curves. For example, in the unphysical exam ple (6.7) th e curves span th e intervals [1,4], [3,6] and [2,5].21 The first and th e last of these intersect in [2,4] which is not one of the spanned intervals, hence this configuration is excluded as unphysical. In contrast, in th e exam ple ( 6.8) th e curves span th e set of intervals { [1,4], [2,3], [5,6]}, and ap a rt from th e em pty interval th e only intersection betw een these intervals is [2,3], which is an elem ent of th e above set. Hence this phase is considered physical. See th e ancillary file for a concrete im plem entation.
This approach allows us to im plem ent a general algorithm th a t gives us all possible phases for th e entanglem ent entropy of a set of n disconnected intervals, w ith any n. For th e case n = 3, we th en have to exclude phase 4 in figure 18, and are left w ith th e five physical phases already identified in [78]. For n = 4 for exam ple, th e 14 physical phases are shown in figure 19. For general n, th e num ber of these physical phases is given by the n -th C a ta la n num ber [86] Cn = -+ r ( 2n 1 , (6.9) n + 1 n which grows as Cn ~ n 3 f o r large n. However, w ith a more optim ized code, it m ay not be necessary to com pute all th e values of this num ber of phases [83].
21Rem em ber th a t the num bers 1 to 6 serve here as labels of (ordered) boundary points, and not necessarily as coordinates on th e x-axis.
V 5 ^ 6 j F ig u re 19. Physical configurations for n = 4 intervals. Note that out of the 4! = 24 configurations that our counting would naively have suggested, there are only 14 physical ones, i.e. 14 ones where the curves do not intersect when projected to the same plane.

I n e q u a litie s fo r t h e u n io n o f n in te r v a ls
O ur interest in this section is to study th e entanglem ent inequalities th a t can be form ulated when working w ith n > 2 intervals.
At th e level n = 3, th e m ost well-known inequality th a t entanglem ent entropies are expected to satisfy is th e strong subadditivity (SSA) [87], com m only stated as

JHEP10(2017)034
A different inequality often associated w ith SSA is [87] S ( A B ) + S ( B C ) -S (A ) -S ( C ) > 0, (6.11) see [78,88] for a discussion of th e relation betw een (6.10) and (6.11) in the holographic case. For th e static holographic cases in which the Ryu-Takayanagi prescription [42,43] applies, these two inequalities were proven in [76]. For th e case of tim e-dependent bulk spacetim es where th e HRT prescription [44] applies, a proof of (6.10) was given in [89] using th e null curvature condition, see also the review [90]. Similarly at n = 3, we encounter w hat is known as th e T his was proven for th e static RT case in [88], and for the tim e dependent HRT case in [89]. T here are m any m ore possible inequalities th a t entanglem ent entropies for n > 3 intervals have to satisfy [92,93], which can be seen to follow from (6.12) and th e other inequalities discussed so far for n < 3 [88]. Hence these inequalities are also proven to hold in holography, assum ing appropriate energy conditions.
T he study of entanglem ent entropy inequalities in A d S /C F T is hence of very high im portance for th e understanding of holography. On th e one hand, if it can be shown th a t holographic prescriptions satisfy certain entanglem ent inequalities th a t do not hold in general quantum theories, this would help distinguish quantum theories th a t can in principle have a simple holographic dual from those th a t cannot. On th e o ther hand, if energy conditions in th e bulk can be used to prove certain entanglem ent inequalities th a t have to hold in th e dual, th en conversely, it m ay be possible to derive novel bulk energy conditions from boundary entanglem ent entropies [94].
In th e following, we will hence use th e entanglem ent entropies th a t we have calculated in our tim e dependent background m etric (2.9) using th e m atching procedure of section 4 to test, for th e m anifestly tim e dependent HRT case, th e validity of some of th e entanglem ent inequalities derived in [54] for th e static RT case. It should however be noted th a t as the m etric (2.9) is a vacuum solution to E in stein 's equations everywhere, it trivially satisfies all com m on energy conditions, and is hence considered to be a physical spacetim e. We hence do not expect any of th e inequalities of [54] to be violated, however as th eir proof is only valid in th e static case, it is interesting to test this expectation thoroughly. At the level of n = 5 boundary intervals, these inequalities read

S ( A B C ) + S (B C D ) + S (C D E ) + S (D E A ) + S (E A B ) -S ( A B C D E ) -S ( B C ) -S (C D ) -S ( D E ) -S ( E A ) -S ( A B )
> 0 (6.13)

2S ( A B C ) + S (A B D ) + S ( A B E ) + S (A C D ) + S ( A D E ) + S (B C E ) + S ( B D E ) -S (A B) -S ( A B C D ) -S ( A B C E ) -S ( A B D E ) -S ( A C ) -S(A D )
-S ( B C ) -S (B E ) -S ( D E ) > 0 (6.14) 22See [91] for an illuminating discussion of the concept of monogamy for entanglement measures.
-39 - generalising th e concept of th re e -p a rtite inform ation introduced in (6.12) . In a holographic context, quantities such as four-and five-partite inform ation where studied for exam ple in [81,82]. In fact, based on th e exam ples studied in those papers, th e authors proposed th e entanglem ent inequalities W hile the inequalities (6.19) and (6.20) m ay be tru e for th e special cases studied in [81,82], where all intervals have th e same length and distance from th eir neighboring intervals, it was already stated in [88] th a t (6.19) and (6.20) do not hold in general holographic setups. In fact, using th e num erical d a ta for the tim e dependent backgrounds studied in this paper or sim ply d a ta valid for static backgrounds such as th e BTZ m etric (2.2) and feeding this d a ta into our ancillary file, it is possible to find explicit exam ples for sets of four or five intervals th a t will lead to violations of th e proposed inequalities (6.19) and (6.20) . However, there are different ways to assign th e labels A , B , C to th e intervals [1,2], [3,4], [5,6] in figure 18. If we impose a strict alphabetical ordering A = [1,2],B = [3,4 ],C = [5,6], th e inequality (6.10) will not be equivalent for exam ple to th e inequality In this context, it is interesting to note th a t th e degree of sym m etries uncovered this way varies from one inequality to th e other. For example, while by perm uting th e intervals A , B , C, D , E (which in this subsection are now assum ed to be ordered alphabetically on th e boundary) gives us 10 inequivalent inequalities following from (6.15) , this num ber is 60 for (6.17) .

.4 A n a ly s is a n d r e s u l ts
We now have alm ost all prerequisites th a t are needed in order to check th e validity of en tanglem ent inequalities such as (6.13)- (6.17) in th e tim e dependent system holographically described by th e bulk m etric (2.9) . As th e m atching-procedure outlined in section 4 can only be applied when th e geodesics cross th e shockwave once, we have to restrict ourselves to th e study of intervals for which all boundary points have x-coordinates either larger th a n zero or sm aller th a n zero. We generally assum e all boundary points of intervals under investigation to be located at equal boundary tim e. We have carried out this analysis for various choices of tem peratures TL and T R, for various values of th e boundary tim e t, and for intervals to th e left and to th e right of x = 0. As th e results were qualitatively sim ilar in all these cases, we will in the following only discuss th e exam ple where we chose TL = 9 ,T R = 1 (hence = 5) and th e boundary tim e slice to be a t t = 1, w ith intervals in th e range x > 0. According to our discussion -41 -

JHEP10(2017)034
in section 6.1 there will be 42 possible phases th a t th e entanglem ent entropy of n = 5 intervals can take. Furtherm ore, in a non-homogeneous system , n intervals are defined by 2n boundary points. As our calculations will be num erical, we are faced w ith a severe problem : even if we are able to check th e validity of inequalities such as (6.13)-(6.17) for any given set of five intervals, how can we make sure th a t we find all potentially interesting cases? A fter all, we have to cover a 2n dim ensional param eter space, on which phase transition s betw een 42 different phases can occur. Numerically, we will of course only be able to check a finite num ber of exam ples. Naively, th e best idea appears to be to take a finite num ber of evenly spaced points on the boundary, x = 0.1, 0.2, 0.3,..., 2.0 (6.23) and to calculate th e entanglem ent entropy for any interval formed by any two of these points. From this d ata, we m ay th en calculate th e entanglem ent entropy of the union of any possible set of n intervals th a t can be form ed from th e given boundary points, and subsequently check th e validity of all entanglem ent inequalities of interest. However, we can do b e tte r th a n this. In th e study of m utual inform ation for Poincare backgrounds, where there are only two phases as shown in figure 16, it is known th a t the phase will depend on the relation betw een th e sizes of the two intervals and th e distances betw een them . In our a tte m p t to cover th e relevant phase space for n < 5 intervals, it will hence be advantageous to allow for th e distances betw een boundary points to vary between as m any orders of m agnitude as possible. Instead of using equally spaced boundary points such as in (6.23) , th e idea is thus to use points which are positioned in a fractal-like way, e.g. where we have found th a t th e choice a = 9 /2 gives a good trade-off betw een th e orders of m agnitude of length scales covered and the overall num ber N of points, which for a reasonable runtim e of our num erics we would like to keep a t N = 20.
Using th e m atching prescription explained in section 4, we have calculated th e renor malised lengths of th e geodesics connecting any two of th e N = 20 boundary points under consideration. In our case a t hand, this requires 2 N ( N -1) = 190 calculations. As a next step, for some n < 5 we w ant to form n boundary intervals by selecting 2n boundary points out of th e N available p o in ts. 24 Obviously, there are 23Indeed, the inspiration for this comes from th e concept of fractal antennas, which in an ten n a technology can be used when attem p tin g to tran sm it in a broadband characteristic, com pared to sta n d ard dipole antennas. We thus aim a t choosing th e boundary points in such a way th a t th ey form a m etaphorical 'fractal a n ten n a' for the stru ctu re of entanglem ent entropy and n -p a rtite inform ation over m any length scales in th e quantum system th a t we are studying holographically.
24As we are selecting 2n distinct boundary points, th e intervals u nder investigation will never be adjacent, i.e. they will never share a boundary point. ways to do so. Given our N = 20 points positioned on th e boundary tim e slice t = 1 according to th e sequence (6.24) , we are hence for exam ple able to study 184756 distinct unions of n = 5 non-adjacent boundary intervals. For all these 184756 different cases, it is th en possible to calculate th e entanglem ent entropy, see figure 20. Interestingly, in th e case at hand we find th a t out of th e 184756 available unions of intervals, 100177 are in th e totally disconnected phase in which S (A 1 U A2 U ...) = S (A i) + S ( A 2) + ..., and in which hence th e entanglem ent inequalities (6.12)-(6.20) are trivially satu rated . Consequently, only the rem aining 84579 cases will be of fu rth er interest. It should also be noted th a t th e overall value of th e entanglem ent entropy for a given union of intervals is dependent on th e explicit cutoff used in our num erics. However, th e linear com binations of entanglem ent entropies appearing in th e inequalities (6.10) , (6.12) b u t also (6.13)-(6.17) , are always balanced in such a way th a t th e cut-off dependence of th e individual term s cancels, such th a t th e result is physical.25 Now, we have all th e necessary ingredients together to check the validity of various entanglem ent inequalities, as well as of th eir perm u tated versions as discussed in section 6.3. 25For (6.11) this will not be the case unless the intervals A, B and C share some of their endpoints. We will not study this inequality in this paper.

JHEP10(2017)034
T he results are as follows: • At n = 3, b o th strong subadditivity (6.10) and m onogam y of m utual inform a tion (6.12) are satisfied, as expected based on [89,90]. In contrast to [80], we find th a t even the engulfed phase (see figure 18) can be th e m inim al one in specific examples. Generically, this seems to happen when the m iddle interval is very small com pared to th e gap betw een th e other two intervals.
• At n = 4, th e only inequality th a t we are checking is th e positivity of four-partite inform ation (6.19) . In fact, in contrast to [81,82], we find a num ber of exam ples for sets of four intervals where this inequality is violated. However, as it was pointed out in [88] and was explicitly checked by us, this is not a p articular feature of the tim e dependent case, and happens already in holographic system s w ith static bulk spacetim e duals.
• At n = 5, we find num erous violations of th e negativity of five-partite inform a tion (6.20) , see th e sim ilar discussion for n = 4. In fact, out of th e 184756 to ta l and 84579 nontrivial sets of five intervals under investigation, we find a violation of ( 6.20) in 417 cases. It is also notew orthy th a t even for th e 84579 cases where five-partite inform ation does no have to vanish a priori, th e result th a t we obtain vanishes w ithin num erical accuracy for 81183 cases, see figure 21.
Furtherm ore, we check th e inequalities (6.13)-(6.17) as well as all th eir relevant per m utations, see section 6.3. T he result is th a t we find n ot a single case in which any of these inequalities is violated, neither for th e specific exam ple currently at hand (T l = 9, Tr = 1 ,t = 1) nor for any other exam ple th a t we studied. See for exam ple figure 22. We view this as a clear indication th a t th e inequalities (6.13)-(6.17) , although so far only proven in th e static case, will generally also hold in physical tim e-dependent cases. A fter investigating th e one-dim ensional case in detail, it is n atu ral to ask about a gener alisation to higher dimensions. T h a t case is, however, much subtler. It has already been addressed in various works. Let us briefly sum m arize th e current sta te of discussion about th e higher-dim ensional case. In [41], a straightforw ard generalization of th e 1+1 dim en sional m odel was suggested, nam ely a solution consisting of two shockwaves, not necessarily travelling w ith identical velocities, and a non-equilibrium steady sta te betw een th e shocks. Such a solution was num erically found in th e hydrodynam ic regime [41]. L ater, a sim ilar solution beyond th e hydrodynam ic approxim ation was found in [45], in th e fram ework of gauge/gravity duality. However, at th e hydrodynam ical level an inconsistency betw een the non-equilibrium steady-state (NESS) conjecture of [41] and therm odynam ics was pointed out in [46]. T he issue is as follows: th e setup of two heat b ath s p u t in contact at an initial

JHEP10(2017)034
tim e is essentially a classical R iem ann problem of solving p artial differential equations.26 For this type of problem , there is a so-called entropy condition. This is a requirem ent th a t characteristic lines of th e differential o perator involved, i.e. curves along which the initial condition is tran sp o rted , always end rath e r th a n begin on th e shock wave. 27 The nam e 'entropic condition' comes from th e fact th a t if characteristics end a t some point, the inform ation about th e initial sta te is lost and hence th e entropy is produced. If the char acteristics started a t th e discontinuity, th e system would require fixing an initial condition on th e shockwave, such th a t inform ation would be produced and entropy would decrease. A detailed analysis of this condition for higher dim ensions leads to th e conclusion th a t while a shockwave moving from th e region of higher tem p eratu re to the colder one is a entropically valid solution, a shock moving in th e opposite direction is not (see section 3 of [46] for details). T he results of [46] im ply th a t the solution involving two shockwaves is valid in d = 2 only, when th e velocities of th e shocks are identical and equal to one. In higher dimensions, to stay in agreem ent w ith entropic considerations, we have to replace th e unphysical shockwave by a new solution -th e rarefaction wave -which is continuous b u t not sm ooth and much wider th a n th e shockwave. Let us stress th a t the double-shock solution is not m athem atically incorrect since for com plicated non-linear PD Es, uniqueness of solutions is not always guaranteed for a rb itra ry types of boundary or initial conditions. T he double shockwave is however non-physical due to th e entropic reasons m entioned. The physical solution is unique in the sense th a t the shock-rarefaction solution is realised in n ature. Let us however em phasize th a t as shown in [46], th e double-shock solution is a valid, physically correct and unique solution to the initial value problem of our non-linear An im portant question about th e shock-rarefaction solution in higher dim ensions is w hether it does support th e existence of a NESS, defined as a region w ith constant energy current th a t can be obtained by boosting a static therm al sta te w ith some effective tem perature. T here are two possibilities: either th e rarefaction solution extends over a large enough region and reaches th e existing shock, excluding the form ation of NESS, or the rarefaction wave is relatively com pact and a NESS is formed betw een th e rarefaction and th e shock wave on th e other side. In [47], th e authors argue in favour of th e latter, based on num erical studies for hydrodynam ical setups. M oreover they discover th a t for m ost conditions, th e quan titativ e difference of observables obtained in a non-physical dual-shock solution and those obtained in th e therm odynam ically favoured rarefaction-shock solution is of order of a few percent. T he specific properties of the steady state rem ain sim ilar to th e universal behaviour of th e NESS in [41]. A fu rth er question is w hether th e dual gravity description allows for a physical rarefaction-shock solution. An exam ple of such a solution was found in [95] in th e lim it of large dim ensions d ^ rc>. However, obtaining a clear, num erical shock-rarefaction solution 26In full generality, the Riemann problem is a initial value problem for a non-linear PDE with noncontinuous, piecewise-constant initial data.
27In the characteristic formulation of a PDE, the presence of a shockwave is manifested by the intersection of characteristics. On a characteristic line, one direction is distinguished by the fact that the initial condition is evolved forward in time. So, when there is an intersection of characteristics, it is possible to distinguish whether the line 'begins' or 'ends' in that point. in th e gauge/gravity fram ework in d = 3 or 4 dim ensions and testing its properties such as stability is still an open problem . It is w orth noting th a t the au th o rs of [45] found a full num erical solution of an 'alm ost' R iem ann problem where th e initial condition is a sm ooth approxim ation of a step function, as our (2.12) in th e fram ework of A d S /C F T . Since, as discussed previously, th e values of th e observables obtained from the shock-rarefaction solution and th e double shock solution differ only by a few percent, it is not clear which of them is the holographic dual of the hydrodynam ic solution. T he au th o rs of th e previously m entioned paper them selves sta te th a t th eir solution seems to agree w ith th e double-shock conjecture, th e work however was published before th e entropic issues of th e double shock solutions were pointed out in [46].
T he argum ents m entioned here ensure th a t a qualitative analysis of th e entanglem ent entropy in higher dim ensions can be carried out, based on th e simple NESS m odel of [41]. We devote th is section to th is analysis.

A n a ly tic a l c o n s id e r a tio n s
Here we present analytical results for th e higher-dim ensional cases. These are obtained by assum ing th a t th e dual-shock solution is valid a t least approxim ately. W hile th e higher dim ensional analogue of th e tim e-dependent m etric (2.9) is not known analytically any more, we still know th e boosted black-brane line elem ent generalising ( where P is th e boost velocity. It is also im p o rtan t to note th a t in higher dimensions, the two shockwaves move w ith different velocities, [41] 28N ote th a t in contrast to [41], we are here using a n o tatio n in which th e dim ensionality of th e bulk AdS space is d + 1. Hence the case investigated so far was th e one for d = 2. A lthough the entanglem ent entropy S(£) for a strip of w idth £ and infinite extent in th e dy± directions cannot be calculated analytically in th e background ( 7.1) , we find th a t in th e lim it £ ^ to where th e entanglem ent entropy becomes extensive, th e entropy density analytically reads (7.7) Let us consider th e question w hether m eaningful statem ents, sim ilar to (5.20) , ab o u t the (average) entropy increase rates of a strip in this setup m ay also be found for higher dim ensions. Due to th e form of th e velocities (7.6) , we see th a t the tim e the shockwave takes to pass through a strip 29 of w idth £ is (7.9) J u st as in section 5.4, we m ay calculate th e average increase/decrease rate of entanglem ent entropy. We assum e for now th a t th e entanglem ent is only influenced by th e shockwaves, and not th e light cones. We find  F ig u re 23. Normalised entropy increase and decrease rates (7.11), (7.12) for d = 2 (solid), d = 3 (dashed) and d = 4 (dotted) as a function of x. Note that in the formulas (7.11), (7.12), > 0, < 0 for x > 1 and < 0, > 0 for x < 1. While 1 > > -1 for d = 2, we see that 1 > with no lower bound for d > 2.
case will produce a clearer picture of th e relation betw een th e entropy increase rate for a given intervals and other bounds or quantities such as (5.3) or (5.4) . Such a num erical solution will also allow to address th e im pact of considering a shock or a rarefaction wave in relation to the absence of a lower bound on vav in d > 2. This m ay be relevant for a general discussion of w hether choosing a gravity solution th a t decreases th e therm odynam ic entropy has unphysical consequences for th e entanglem ent entropy.

N u m e r ic a l c o n sid e r a tio n s
Refs. [96,97] give a solution of th e background equations of m otion on th e gravity side in th e case d > 2 by considering a linearization of th e system . This approach tu rn s out to be equivalent to linearized hydrodynam ics, as it is valid as long as |TL -T R < |TL + T R . Using this background, we m ay com pute num erically th e entanglem ent entropy for any num ber of dim ensions by following th e procedure of section 3.   cf. eq. (3.13) , appears to be a robust property valid also for d = 3, and th e same can be said for th e decrease of S (A U B ) w ith tim e. A more detailed study of these issues will be addressed elsewhere.

C on clu sio n and o u tlo o k
In this work we have studied a holographic m odel for far-from -equilibrium dynam ics th a t describes th e tim e-dependent properties of energy flow and inform ation flow of two therm al reservoirs initially isolated. In this system , a universal steady sta te develops, described by a boosted black brane. A relevant observable th a t provides physical insight into th e evolution of th e system is th e entanglem ent entropy, which m easures th e inform ation flow between two subsystem s. By using th e exact solution for d = 2 provided in [41], we have studied the tim e evolution of th e entanglem ent entropy, and characterized some universal properties of th e quenching process. We also studied the tim e evolution of m utual inform ation and found it to m onotonically grow in tim e.
In section 5 , after a brief overview of velocity bounds for entropy spread and increase, we have investigated th e m atching procedure outlined in section 4 in m ore detail, showing th a t in certain circum stances an analytical solution is possible. This allowed us to prove -50 -

JHEP10(2017)034
th e validity of th e universal form ula (3.10) in th e appropriate low tem p e ra tu re lim its. In subsection 5.4, we th en investigated th e increase rates of entanglem ent entropy obtained using th e num erical and analytical results of th e previous sections. We find th a t both averaged and m om entary entanglem ent entropy increase and decrease rates are bounded by th e speed of light (5.25) . W hile this bound is close to being sa tu ra te d for intervals th a t are large com pared to th e scale set by th e tem p eratu re, this is not th e case for sm aller intervals, where th e universal form ula (3.10) becomes a good approxim ation, see again figure 15. This indicates th a t th e shockwave in our setup, which in m any ways is sim ilar to a local quench, mimics an entanglem ent tsunam i for large interval sizes t, leading to a linear entropy increase w ith th e appropriate rate. For small t however, th e universal behaviour (3.10) w ith its characteristic S-shape takes over. We refer to this as an 'entanglem ent tid e '. As discussed in section 7, it will be very interesting to study these questions for analogous system s in higher dimensions, where the speed of light, the entanglem ent velocity vE and th e butterfly velocity vB are not equivalent any more. This m ay help to get a b e tte r understanding of th e m echanism s related to entanglem ent tsunam is.
A part from th e m onogam y of m utual inform ation and strong subadditivity, other in equalities involving a large num ber of subsystem s have been proven in th e static case, see [54]. In section 6, we have studied various entanglem ent entropy inequalities, which were proposed for up to n = 5 intervals, in th e present tim e-dependent system . W hat we found was th a t th e inequalities proven in [54] also hold in th e tim e-dependent system under consideration in this paper, at least in all cases th a t we num erically checked. How ever, we found th a t th e signs of four-and five-partite inform ation are not definite in this holographic system , in contrast to th e results of [81,82]. As th e bulk m etric investigated in this paper is a vacuum solution everywhere, and hence trivially satisfies th e m ost common energy conditions, we did not have any a priori reason to expect encountering a violation of th e entanglem ent entropy inequalities of [54]. It m ay hence be an interesting possibil ity for future research to check th e validity of these inequalities for tim e-dependent bulk spacetim es th a t violate, for exam ple, th e null energy condition (N EC), sim ilarly to w hat was done for strong subadditivity in [78,90,98,99]. W ith this paper, we also upload the num erical code used to obtain th e results of section 6 to th e arXiv. We hope th a t this will facilitate future research in this direction.
One of th e possible furth er directions of investigation is suggested by th e elegant ana lytical behaviour of the entanglem ent entropy in th e small tem p eratu re lim it. It is known th a t low-energy behaviour of ballistic, quantum -m echanical models is well described by conform al field theories. For a therm al sta te this m eans th a t in th e low -tem perature regime of lattice m odel m ay be approxim ated by a therm al sta te of a C F T since th a t is a situ a tion in which lower p art of energy spectrum determ ines properties of th e theory, as more excited states are not occupied. Therefore, we presum e th a t th e simple universal evolu tion of th e entanglem ent entropy we observe should be as well visible in lattice (i.e. tensor netw ork or exact diagonalisation) calculations. It will be interesting to com pare to th a t kind of models, as local quenches in such system s have recently draw n some attention, see for exam ple [100]. Moreover, it is conceivable th a t fu rth er physically observables can be com puted in th a t low -tem perature limit.

JHEP10(2017)034
Finally, com parisons to non-equilibrium hydrodynam ics m ay provide fu rth er useful inform ation. Recent work on this includes [101]. Universal structures in a holographic m odel of non-equilibrium steady states, which are spatial analogues of quasinorm al modes, have recently been considered in [102]. O p e n A c c e ss. This article is distrib u ted under th e term s of th e Creative Commons A ttrib u tio n License (CC-BY 4.0) , which perm its any use, d istribution and reproduction in any m edium , provided the original author(s) and source are credited.