The holographic dual of a Riemann problem in a large number of dimensions

We study properties of a non equilibrium steady state generated when two heat baths are initially in contact with one another. The dynamics of the system we study are governed by holographic duality in a large number of dimensions. We discuss the “phase diagram” associated with the steady state, the dual, dynamical, black hole description of this problem, and its relation to the fluid/gravity correspondence.


Introduction
The Riemann problem may provide a relatively simple setting in which to study the nonequilibrium physics of quantum field theory. The problem asks for the time evolution of piece wise constant initial conditions with a single discontinuity in the presence of some number of conservation laws, for example of energy, momentum, mass, or charge. In our case, we consider a fluid phase of a conformal field theory (CFT) with an initial planar interface, where the energy density jumps from e L on the left of the interface to e R on its right. We also allow for a discontinuity in the center of mass velocity of the fluid across the interface.
For simplicity, we will make a number of further restrictions. We assume a conformal field theory that has a dual gravity description via the AdS/CFT correspondence. A priori, this will allow us to study the system beyond the hydrodynamic limit. We also take the JHEP08(2016)120 Figure 1. A phase diagram for the solution to the Riemann problem in a large d limit. Given a pair (e L , 0) and (e R , j R ), the selection of shock and rarefaction waves is determined by the value of e R /e L and j R /e L . The dashed and solid lines are "critical": the dashed line indicates the values of (e R , j R ) connected to (e L , 0) by a single rarefaction wave while the solid line indicates the values of (e R , j R ) connected to (e L , 0) by a single shock wave.
limit that the number of spatial dimensions d is very large. In this limit, we find that the system is described by two conservation equations where e is, up to gradient corrections, the energy density and j the energy current. These equations are a special case of equations derived in ref. [1]. In these variables the Riemann problem amounts to a determination of e and j given an initial configuration of the form (e, j) = (e L , j L ) z < 0 (e R , j R ) z > 0 . (1.2) By choosing an appropriate reference frame, we may set j L = 0 without loss of generality. As it happens, there are extensive treatments of this type of Riemann problem in hydrodynamics textbooks. See for example ref. [2]. Typically, a pair of rarefaction and/or shock waves form and move away from each other, creating in their wake a region with almost constant e and j. In recent literature, this intermediate region has been called a nonequilibrium steady state (NESS) [3,4]. One of the main results of this paper is a "phase" diagram valid in a large d limit (see figure 1) that describes, given the conservation equations (1.1) and initial conditions (1.2), which pair of waves are formed: rarefaction-shock (RS), shock-shock (SS), shock-rarefaction (SR), or rarefaction-rarefaction (RR). A physical reason for the preference of a rarefaction wave to a shock wave is entropy production.
Recent interest in this type of Riemann problem was spurred by a study of the problem in 1+1 dimensional conformal field theory [3] where the evolution is completely determined by the conformal symmetry and a hydrodynamic limit need not be taken. Conservation JHEP08(2016)120 and tracelessness of the stress tensor imply that the stress tensor is a sum of right moving and left moving parts. When j R = j L = 0 one finds a NESS in between the two asymptotic regions, characterized by an energy density (e R + e L )/2 and an energy current proportional to e R − e L . The NESS is separated from the asymptotic regions by outward moving shock waves traveling at the speed of light. (An extension of the analysis of [3] which includes a discontinuity in the center of mass velocity, holomorphic currents and chiral anomalies can be found in [5]. An analysis of shock waves and their relation to two dimensional turbulence was carried out in [6].) In more than two space-time dimensions, conformal symmetry alone is not enough to specify the evolution completely and one needs additional assumptions about the structure of the conserved currents. Recent work appealed to the gauge/gravity duality [7][8][9][10], an analogy with 1 + 1 dimensions [5], and hydrodynamics [7,[11][12][13]. These papers focused on the case j R = j L = 0 and e L > e R such that from a hydrodynamic perspective a left moving rarefaction wave and a right moving shock wave are expected to emerge.
The distinction between rarefaction and shock waves was ignored in some of these papers [5,7,11]. Indeed, when working with 2 + 1 or 3 + 1 dimensional conformal field theories, the difference between, say, an SS solution to the Riemann problem and an RS solution to the Riemann problem is very small for all but extreme initial energy differences. As the spacetime dimension d increases however, the difference between a rarefaction wave type of solution and a shock wave solution becomes significant [13]. This amplification of the difference between the two solutions serves as a motivator for studying this Riemann problem in a large number of dimensions.
Interestingly, a large d limit has independently been a topic of recent interest [1,[14][15][16][17][18][19][20][21][22][23][24][25] in the study of black hole solutions to Einstein's equations. Of particular relevance to our work is the connection between black holes in asymptotically AdS spaces and hydrodynamics [26]. Certain strongly interacting conformal field theories are known to have dual classical gravitational descriptions. In the limit where these conformal field theories admit a hydrodynamic description, a solution to the relevant hydrodynamic equations can be mapped to a solution of Einstein's equations, in a gradient expansion where physical quantities change slowly in space and time. Transport coefficients such as shear viscosity are fixed by the form of Einstein's equations. Thus, one may study the Riemann problem in conformal field theories with a large number of dimensions by studying an equivalent Riemann-like problem involving an initially discontinuous metric of a black hole in an asymptotically AdS background.
Given that extensive analyses of conservation equations like (1.1) can be found in many hydrodynamics textbooks and papers, one can legitimately ask why we bother to redo the analysis here. The reason is that when working in a large number of dimensions, one can solve for the black hole metric exactly, independent of the derivative expansion (which is naturally truncated), thus obtaining an exact solution to the Riemann problem which includes possible viscous terms and is in general valid even when gradients of thermodynamic quantities are large (as is the case with discontinuous initial conditions).
Our work is organized as follows. In section 2, we rederive the equations (1.1) by taking a large d limit of Einstein's equations. We show how to rewrite them as the conservation JHEP08(2016)120 condition on a stress-tensor, ∂ µ T µν = 0. In section 3, we compare the large d stress tensor and equations of motion to those arising from the fluid-gravity correspondence [26]. We find that both eqs. (1.1) and the stress tensor T µν are equivalent to the hydrodynamic equations that come from the fluid-gravity correspondence at large d, at least up to and including second order gradient corrections. In the same section we also construct an entropy current J µ S using an area element of the black hole horizon and show that the divergence of the entropy current is positive ∂ µ J µ S ≥ 0 in this large d limit. In section 4, we solve the Riemann problem for eqs. (1.1) and derive the phase diagram given in figure 1. Finally, we conclude in section 5 with some directions for future research. Appendix A contains a short calculation of the entropy produced across a shock, while appendix B contains plots of auxiliary numerical results.
2 The holographic dual of the Riemann problem for large d We wish to construct a holographic dual of the Riemann problem. Consider the Einstein Hilbert action A canonical stationary solution of the resulting equations of motion is the black brane solution where T is an integration constant which denotes the Hawking temperature. The solution (2.2) is dual to a thermal state of a conformal field theory with temperature T . For instance, the thermal expectation value of the stress tensor in such a state is given by is the pressure with p 0 a theory dependent dimensionless parameter. (The indices µ and ν run over the d − 1 dimensions of the (d − 1)-dimensional CFT.) As discussed in [8] a dual description of the Riemann problem necessitates an initial black hole configuration which is held at some fixed temperature T L for all z < 0 and at a different temperature T R for z > 0. This would correspond to a configuration where the expectation value of the stress tensor is given by (2.3) with T = T L for z < 0 and by (2.3) with T = T R for z > 0. Since the initial black hole is out of equilibrium it will evolve in time. Its dual description will provide a solution for the time evolution of the stress tensor which we are after. Thus, our goal is to solve the equations of motion following from (2.1) and use them to construct the dual stress tensor.

JHEP08(2016)120
An ansatz for the metric which is compatible with the symmetries and our initial conditions is given by where the metric components are functions only of t, r, and z. (A more general ansatz which involves a transverse velocity can be found in [1].) A numerical solution of the equations of motion for g tt , g tz and g ii (i = x ⊥ or z) with smoothened initial conditions has been obtained for d = 4 in [8] for relatively small initial temperature differences, A solution for finite d > 4 and for large temperature differences, In this work we use the methods developed in [1,14] (see also [15][16][17][18][19][20][21][22][23]) to address the Riemann problem in the limit that d is very large. Such a limit can be understood as follows. In an appropriate gauge, the near boundary expansion of the metric gives Thus, in the large d limit at any finite value of r, the spacetime looks like the AdS vacuum. Only by keeping R = r n finite with n ≡ d − 1 will the O(r −n ) corrections to the metric remain observable. Our strategy is to solve the equations of motion in the finite R region subject to the boundary conditions (2.6). Following [1], we also use the scaling x ⊥ = χ/ √ n and z = ζ/ √ n so that in this coordinate system the line element takes the form (2.8) (In a slight abuse of notation i is now either χ ⊥ or ζ.) We have used the letters E and J to emphasize these quantities' (soon to be seen) close connection with an energy density and energy current in the dual hydrodynamic description.
One can now solve the equations of motion order by order in 1/n. The equations of motion are simply Einstein's equations in the presence of a negative cosmological constant: (2.9) JHEP08(2016)120 setting L = 1 for convenience. Let a and b index the t, r, and ζ directions only, while i and j index the remaining perpendicular directions. Furthermore, letR ab be the Ricci tensor with respect to the three dimensional metric in the t, r, and ζ directions. Then Imposing that the boundary metric is Minkowski and choosing a near boundary expansion of the form (2.6) we find where the O(n −2 ) correction to g tt and the O(n −3 ) contributions to g ζζ are too long to write explicitly. The functions e and j are functions of t and ζ only and must satisfy the additional constraints (1.1). Equations (1.1) are identical to those obtained in [1,14]. We can rewrite them in terms of a conservation law (2.14) where g is an arbitrary function. Likewise, the functions e 2 and j 2 must also satisfy a set of equations which can be obtained from the conservation of (2.16) We will use and ∂ ζ interchangeably in what follows.

JHEP08(2016)120 3 Comparison with hydrodynamics
Let us pause to understand (2.14). Within the context of the gauge-gravity duality it is possible to construct a solution to the Einstein equations which is perturbative in t, ζ and χ ⊥ derivatives of the metric components [26]. Such a perturbative solution to the equations of motion, which is available for any dimension d [27,28], allows for a dual description of the theory in terms of fluid dynamical degrees of freedom.

Stress tensor from fluid-gravity correspondence
To construct the dual hydrodynamic description of a slowly varying black hole, we boost the black hole solution (2.2) by a constant velocity u µ in the t, z, x ⊥ directions. The resulting line element is given by (3.1) Allowing for u µ and T to become spacetime dependent implies that (3.1) will get corrected. By setting gradients of u µ and T to to be small, one can solve for the corrections to (3.1) order by order in derivatives so that the line element will take the schematic form where ds 2 (i) denotes the ith order gradient corrections to the line element. The stress tensor T µν which is dual to (3.1) takes the form also expanded in gradients. One finds [27,28] T µν which is nothing but a boosted version of (2.3) and then, in the Landau frame, (Note that our definition of σ µν is somewhat unconventional.) An initial analysis of third order gradient corrections has been carried out in [29] for d = 5. A full analysis of all third order transport terms for arbitrary dimension d is currently unavailable. Since (2.14) has been obtained from a large d limit of a gravitational dual theory, we expect that (2.14) coincides with (3.3) when the former is expanded in derivatives and the latter is expanded around large n = d − 1. In short, we expect that taking a gradient expansion commutes with taking a large d limit. To make a direct comparison let us consider the hydrodynamic stress tensor (3.3) in the t, ζ, χ ⊥ coordinate system where the metric tensor takes the form One important effect of this rescaling is to keep the sound speed to be an order one quantity. Scaling the spatial component of the velocity field by 1/ √ n, viz., and maintaining that = (d − 2)P is finite in the large d limit, we find, and thus, and O ∂ 3 denotes third order and higher derivative corrections. Note that this constitutive relation for the stress tensor includes and encodes the large d limit of the transport coefficients (3.7). Now, we insert the redefinitions

JHEP08(2016)120
into the large d constitutive relation for the stress tensor (2.14), use the large d stress tensor conservation equations (1.1), and throw out terms that have three or more derivatives. We claim that in this fashion, we recover the stress tensor (3.11) in the gradient expansion. Thus, the large d limit and the gradient expansion seem to commute. Note that while the conservation equations (1.1) are of second order in gradients of ζ and t, the stress tensor includes at least second order gradients. The implications of (3.13) are worth emphasizing. The equations of motion (1.1) are equivalent to the standard equations of motion of relativistic hydrodynamics when the latter are expanded in a large d limit. When working with the e and j variables one obtains equations of motion which are second order in derivatives and therefore include dissipative effects. When carrying out a frame transformation to the more traditional Landau frame, more derivatives will appear. When considering the stress tensor associated with the equations of motion (1.1) one obtains more terms with higher gradients which do not contribute to the equations of motion. It would be interesting to see if one can construct an alternative to the Israel-Stewart theory using a "large d-frame" where gradients naturally truncate.

Entropy from gravity
Within the context of our forthcoming analysis, it is instructive to compute the dual entropy production rate which is associated with the evolution of the horizon. Due to its teleological nature, it is usually difficult to identify the location of the event horizon. However, in the large d limit the analysis is somewhat simplified. Let us look for a null surface of the form R = r h (t, ζ). The normal to such a surface is (3.14) Demanding that Ξ 2 R=r h = 0 implies, to leading order in the large d limit, that The spacetime singularity which exists in our solution implies that an event horizon must be present. Since the only null surface available is (3.15), it must be the location of the event horizon. Subleading corrections to the location of the event horizon are given by To compute the change in the black hole entropy over time we compute the area form of the event horizon. Following the prescription of [30], we find that

JHEP08(2016)120
where h is the spatial (t = constant) part of the induced metric on the horizon and N µ is defined via (3.21) Thus, where we have normalized the entropy density so that it is compatible with our conventions for the energy density. The second law of black hole thermodynamics amounts to In our large d limit we find that (3.24) The expectation from hydrodynamics, to second order in derivatives, is that the divergence of the entropy current is given by (See for example (8) of ref. [31].) This expectation matches (3.24) on the nose. Note that to leading order in the large d limit the entropy current vanishes. This somewhat surprising feature of the large d limit follows from the fact that entropy production terms are suppressed by inverse powers of the dimension in the large d limit. Another way of understanding this suppression comes from thinking about the temperature T ∼ e 1/(d−1) . In the large d limit, T is constant to leading order in d. From the thermodynamic relation de = T ds, it then follows that changes in energy are proportional to changes in entropy, and entropy conservation follows from energy conservation at leading order in a large d expansion. 1

JHEP08(2016)120 4 Near equilibrium steady states
We now analyze the dynamics controlled by the partial differential equations (1.1) which encode the dynamics of an out of equilibrium black hole (2.5) and its dual stress tensor (2.14). Various related holographic analyses can be found in [32][33][34][35][36][37][38][39][40][41]. As discussed in the introduction, the particular question we would like to address is a Riemann problem: what is the time evolution following from an initial condition (1.2)? We are particularly interested in the steady state solution which will emerge at late times. For convenience we will consider a reference frame for which j L = 0. Indeed, if e(x, t) and j(x, t) satisfy the conservation equations (1.1), then so do e(x − vt, t) and j(x − vt, t) + ve(x − vt, t). Thus, for constant values of e and j, we can choose a v such that j will be set to zero. The non-relativistic nature of the boost symmetry reflects the fact that the large d limit we have taken is effectively a non-relativistic limit where the speed of light c ∼ √ d has been pushed off to infinity.

Rarefaction waves vs. shock waves
Before addressing the Riemann problem in its entirety let us consider a simplified system which is less constrained. Consider (2.14) with gradient terms neglected. The resulting expression is the large d limit of the energy momentum tensor of an inviscid fluid which is known to support (discontinuous) shock waves [2] for any finite value of d. While the solution to the full Riemann problem will consist of a pair of shock and/or rarefaction waves, we begin in this section with a single discontinuous shock wave moving with velocity s.

Conservation of energy and momentum imply
where [Q] = Q l − Q r and Q r/l specify the value of Q to the left or right of the shock respectively. 2 The conservation conditions (4.1) are very general and are often referred to as the Rankine-Hugoniot (RH) relations. In our setup they reduce to where e r/l and j r/l are the energy density and current immediately to the right or left of the shock. While these Rankine-Hugoniot relations hold for an arbitrary, piece-wise continuous fluid profile, in what follows, we are interested in the much simpler situation where e and j are constant functions away from the shocks. Amusingly, e r satisfies a cubic equation, 3 2 In this section we use subscripts r and l to denote values of quantities to the right or left of the shock.
In other sections we use subscripts R and L to denote quantities in the right and left asymptotic regions.
In the latter case there is generally an interpolating region which we denote with a 0 subscript. 3 In general d, one finds the relation where β = tanh α is the fluid velocity.

JHEP08(2016)120
a plot of which as a function of j r resembles a fish: fixing (e l , j l ), each value of s is mapped to a point on the (e r , j r ) plane. The collection of such points is given by a fish-like curve, an example of which is given in the left panel of figure 2.
We make two observations about the fish. The vacuum (e r , j r ) = (0, 0) always lies on the cubic (4.3), corresponding to the fact that a shock can interpolate between any value of (e l , j l ) and the vacuum. Also (e r , j r ) = (e l , j l ) is the point of self-intersection of the cubic and has s = ±1 + j l /e l . The physical content of this observation is that when (e r , j r ) is close to (e l , j l ) but still lies on the cubic, we can find a close approximation to the fluid profile by linearizing the equations of motion. As we will describe in greater detail below, linearized fluctuations correspond to damped sound modes, and indeed the two regions can be connected by sound waves propagating at the local sound speed s = ±1 + j l /e l .
The shock solutions we found all solve the conservation equations (4.2). However, some of these solutions are unphysical in the following sense. Let us boost to a frame where the shock speed vanishes, s = 0. In half of the shock solutions, a quickly moving fluid at low temperature is moving into a more slowly moving fluid at higher temperature, converting kinetic energy into heat and producing entropy. We will refer to these shocks as "good" shocks. The other half of the solutions correspond to the time reversed process where a slowly moving fluid at high temperature moves into a rapidly moving but cooler fluid, turning heat into kinetic energy. This second solution, as we shall see shortly, should be discarded.
Strictly speaking, entropy is conserved in the large d limit (see the discussion following equation (3.25)). A more formal way of understanding why one should discard the bad shocks is to restore the gradient corrections but take a limit where these are small. Let us assume that in the frame where the shock velocity is zero there is an approximately stationary configuration such that time derivatives are much smaller than spatial derivatives. Boosting back to a shock with velocity s, we expect that e and j depend only on the combination ζ − st, i.e., j(t, ζ) = j(ζ − st) and likewise, e(t, ζ) = e(ζ − st). The equations of motion (1.1) become ordinary differential equations which can be integrated once to obtain e = − s(e − e l ) + (j − j l ) , (4.4) We have picked the two integration constants such that e and j vanish in the left asymptotic region. The Rankine-Hugoniot conditions (4.2) imply that e and j also vanish in the right asymptotic region. As e and j themselves vanish in the left and right asymptotic regions, we can describe e and j well near these points by looking at a gradient expansion.
Near the left asymptotic region

JHEP08(2016)120
There is a similar looking equation for e and j near the right asymptotic region The solutions near (e l , j l ) and near (e r , j r ) have an exponential nature with the sign of the exponents depending on the eigenvalues of M l and M r appearing on the right hand side of (4.5) and (4.6) given by We now observe that the signs of the eigenvalues of M l and M r determine whether the shock is a viable solution to the equations of motion.
• If both eigenvalues of M l are negative, then e and j will not vanish as x → −∞.
Thus we require that at least one eigenvalue of M l is positive in order for a shock solution to exist.
• If we assume there is exactly one positive eigenvalue, then 1 + j l /e l > s and −1 + j l /e l < s. Note that the value 1 + j l /e l corresponds to the slope of one of the characteristics (i.e. the local speed of one of the sound waves), and this condition implies that this characteristic will end on the shock. Since λ l − is assumed to be negative, we have to tune one of the two integration constants of the system of differential equations to zero. This tuning means that generically the solution to the right of the shock will be a linear combination of both of the solutions near (e r , j r ). If both solutions are to be used, then it had better be that both eigenvalues of M r are negative. (Otherwise, it will not be true that e and j vanish in the limit x → ∞.) In particular, the larger of the two eigenvalues must be negative, which implies that 1 + j r /e r < s. (In terms of characteristics, both will end on the shock.) Thus, we find the constraint 1 + j r /e r < s < 1 + j l /e l . (4.8a) • If both eigenvalues of M l are positive, we still need at least one negative eigenvalue of M r to be able to connect the solutions in the left and right asymptotic regions. Moreover, for M r to have two negative eigenvalues would be inconsistent with momentum conservation (4.2). An analysis similar to the previous one yields The constraints (4.8) choose the good shocks over the bad ones. 4 4 In appendix A, we discuss a third RH relation one can write down for the entropy current. If the RH relations for energy and momentum are satisfied, the RH relation for the entropy current will typically be violated due to entropy production associated with viscous effects. In the weak shock limit, we demonstrate that gradient corrections produce the entropy that leads to this violation of the third RH relation. Since bad shocks are not allowed, one may inquire as to the time evolution of a discontinuity with initial conditions which would have generated a bad shock. As it turns out, bad shocks can be replaced by the more physical rarefaction solutions [2]. The rarefaction solution assumes that between the asymptotic regions specified by (e l , j l ) and (e r , j r ), there is an interpolating solution where e and j are functions of ξ = ζ/t. As was the case for the shock wave, given e l and j l , there is a one parameter family of allowed values of e r and j r . These are given by e r =e l exp (±j l /e l − 1 ∓ ξ r ) , j r =e l (±1 + ξ r ) exp (±j l /e l − 1 ∓ ξ r ) . (4.9) The curve traced by (e r , j r ) also resembles a fish, and for moderate values of the shock parameters e r and j r it closely follows the cubic curve corresponding to a shock solution.
(See the central panel of figure 2.) The vacuum (0, 0) = (e r , j r ) solution can always be connected to (e l , j l ) through a rarefaction wave. The self-intersection point (e r , j r ) = (e l , j l ) has ξ = ∓1 + j l /e l , again corresponding to a sound wave type interpolation between the two regions (e r , j r ) ≈ (e l , j l ). Given that bad shocks are replaced by rarefaction waves, one should remove from the fish diagram (left panel of figure 2) the portion of the curve which corresponds to bad shocks and replace it with a curve corresponding to a rarefaction solution (central panel of figure 2). The resulting curve can be found on the right panel of figure 2: the belly of the JHEP08(2016)120 Figure 3. A graphical determination of the "good shocks" and "bad shocks". The red fish corresponds to (e r , j r ) while the blue fish is built from (e l , 0). See the main text for a discussion. fish and the lower part of its tail corresponds to a good shock and its back and upper tail to a rarefaction solution. One may compute the curve explicitly by imposing (4.8), but it can also be understood from a graphical viewpoint as we now explain.
Recall that the self intersection point of the shock wave fish (solid curve on the left panel of figure 2) corresponds to a shock velocity, s, which takes the values of the local speed of sound, ±1+j l /e l . On the tail, s is either larger than 1+j l /e l (upper tail) or smaller than −1 + j/e (lower tail). Thus, on the tails, the eigenvalues are either both positive or both negative. The top portion of the tail has λ ±l < 0 while the bottom portion of the tail has λ ±l > 0. As a result, the top portion of the tail must be replaced by a rarefaction wave while the bottom portion can be a shock. To decide which portion of the body of the shock fish to replace by a rarefaction wave, one must study λ ±r .
Consider a second fish which exhibits the solution to the cubic (4.3) for a given value of (e r , j r ). We will call this second fish an r-fish and the first an l-fish. Similar to the analysis of the tail of the l-fish, we find that the bottom portion of the tail of the r-fish should be constructed from a rarefaction solution while the top portion from a shock.
Consider an r-fish whose point of self intersection lies somewhere on the body of the l-fish. When the r-fish is drawn so that it intersects the back of the l-fish, the bottom portion of the r-fish's tail will go through the point of self-intersection of the l-fish (see the left panel of figure 3). As the bottom portion of the tail of the r-fish is a rarefaction, the region (e r , l r ) can be connected to (e l , j l ) by a rarefaction. Reciprocally, since we're describing a single shock or rarefaction interface between two regions, the back of the l-fish should be replaced by a rarefaction wave. We can run the argument again for an r-fish drawn to intersect the belly of the l-fish. We conclude that the belly of the l-fish must be a shock (see the right panel of figure 3).

Solving the Riemann problem using ideal hydrodynamics
Armed with our understanding of shock waves and rarefaction solutions, let us now tackle the Riemann problem we set out to solve. At t = 0, we consider a pair (e L , 0) which describes the fluid for z < 0 and another pair (e R , j R ) describing the fluid for z > 0. For a single interpolating shock or rarefaction, we have seen that given (e L , 0) there is a one parameter family of solutions that determine (e R , j R ). Thus, generically, there will not be a single shock or rarefaction solution that joins (e L , 0) to an arbitrary (e R , j R ). However, we can connect the two regions using a pair of shock and/or rarefaction waves. That is, we could connect (e L , 0) to an intermediate regime with values of e and j given by (e 0 , j 0 ) using a shock or rarefaction wave and another shock wave or rarefaction wave to connect the intermediate regime to the right asymptotic region (e R , j R ). In all cases, given the initial conditions, the pair of rarefaction and/or shock waves should be such that they move away from each other.
The strategy for determining which type of solution is allowed is to prefer good shocks over rarefaction solutions and rarefaction solutions over bad shocks. Thus, given a pair (e L , 0) and (e R , j R ) we need to establish which of the four possibilities for the time evolution of the initial state is allowed: two shocks (SS), a rarefaction wave followed by a shock (RS), or the remaining two configurations which we will denote by SR and RR.
To understand the possible solutions to the Riemann problem, let us first consider two fish diagrams: one associated with (e l , j l ) = (e L , 0) (the l-fish) and another with (e r , j r ) = (e R , j R ) (the r-fish). The points of overlap of the diagrams will give us the possible value of e 0 and j 0 . We will always choose a point where the two disturbances are moving away from each other. See, for example, figure 4.
Instead of plotting the r-and l-fishes, we can obtain closed form expressions for the various types of solutions by solving (4.8) and (4.9) on a case by case basis. In the following we provide some simple examples of such expressions.

JHEP08(2016)120
• RS configurations. As an example of the RS case, we take (e L , 0) and (e R , 0) as the asymptotic regions with e L > e R . The SR case is a left-right reflection of the RS case and therefore does not warrant further discussion.
To estimate the values of e 0 and j 0 we can follow the strategy laid out in [12,13].
For the left region we use the solution (4.9) with e l = e L , j l = 0, e r = e 0 and j r = j 0 . For the right region we use (4.2) with e l = e 0 , j l = j 0 , e r = e R and j r = 0. We find which, unsurprisingly, coincides with the large d limit of the hydrodynamic analysis of [12,13].
As pointed out in [12] the rarefaction solution will cover the location of the original shock discontinuity whenever At the point ζ = 0 in the rarefaction wave, the values of e and j are time independent (since any function of ζ/t will have a fixed point at ζ = 0). Moreover for a conserved stress tensor T µν = T µν ζ t , the first spatial derivative of T tζ and the first and second spatial derivatives of T ζζ vanish at this fixed point. Thus, one may think of the pressure at the fixed point as a "short" steady state for long enough times. "Short" implies that the region is of small spatial extent. From this perspective one has split steady states for large enough initial temperature differences. The values of e and j at the short steady state are given by e s = j s = e L exp(−1) . (4.12) • SS configurations. A simple example of the SS case has (e L , 0) on the left and (e L , j R ) on the right with j R < 0. We compute the NESS by gluing two shock waves to an intermediate region with (e, j) = (e 0 , j 0 ), similar to the RS case. Setting β = j R /e L , the intermediate NESS is given by 13) and the shock velocities for the left and right moving shocks, s L and s R respectively, are given by

JHEP08(2016)120
• RR configurations. Using e L = e R and j R > 0, we can find simple solutions that involve two rarefaction waves. 5 In this case, the NESS is characterized by where the left moving rarefaction wave extends from ξ = −1 to ξ = ξ − while the right moving rarefaction wave extends from ξ = ξ + to ξ = 1 with Similar to the RS case we find that there is a fixed point associated with the left moving wave whenever j R 2e L ≥ 1 , We claim that given (e L , 0), the "phase diagram" of figure 1 immediately allows us to choose the correct configuration of shocks and rarefaction waves for any (e R , j R ). Indeed, following figure 4, the location of the self intersection point of the r-fish will determine the nature of the intersection of the r-and l-fish: if the intersection point of the r-fish lies above the l-fish we will always get an RR solution; if the intersection point of the r-fish is below the l-fish we get an SS solution; and RS and SR solutions will correspond to an intersection point of the r-fish in the body or tail of the l-fish respectively. Conformal invariance dictates that the phase diagram can depend on the only two dimensionless parameters of this problem, and we obtain the phase diagram in figure 1.
Note that even though the r-fish and the l-fish intersect at (0, 0), we can always rule out an intermediate point that corresponds to a vacuum. The vacuum intersection point is always along the bodies of the two fish where we have λ −,l/r < 0 < λ +,l/r . As discussed, we can not in general connect the two asymptotic solutions if we do not have two eigenvalues of the same sign (positive for l and negative for r) in one of the regions. 5 As it turns out in the RR phase, there is a simple expression for the steady state for all values of eL, eR, jL and jR, where A fixed point associated with a left moving rarefaction solution occurs whenever and a fixed point associated with the right moving rarefaction solution occurs whenever eR eL ≥ exp jL eL + jR eR + 2 with es = −js = eR exp −1 + jR eR . JHEP08(2016)120

A numerical solution to the Riemann problem
In the previous sections we have obtained predictions for the evolution of e and j starting from an initial configuration (1.2) and assuming that gradient corrections to the equations of motion are small. It is somewhat unfortunate that this assumption stands in stark contrast to the discontinuous jump in the initial state and one may inquire whether the analysis of the previous section is relevant for the problem at hand. In order to resolve this issue we solve the full equations of motion (1.1) numerically. We give numerical examples of the RR, SS, and RS phases described above. To our numerical accuracy, the difference in e 0 and j 0 between the ideal case which we have studied analytically and the case with gradients included which has been obtained numerically appears to disappear in the long time limit.
As it turns out, the equations (1.1) are easy to evolve numerically with canned PDE solvers, such as Mathematica's NDSolve routine [42]. To obtain various solutions one can evolve the initial condition e = e (1 + δe tanh(c sin(2πx/L))) , (4.20) in a periodic box of length L. (In appendix B, we use a more elaborate piecewise continuous initial condition.) For c sufficiently large, the initial condition approaches a square wave. As long as the disturbance has not travelled a distance of order L, causality ensures that the behaviour of e and j are very close to that of an infinite system where the values of e and j in the asymptotic region are fixed at some constant value. If we denote these asymptotic values as e L and e R then δe = e L − e R e L + e R and e = 1 2 (e L + e R ) . We can similarly define j and δj.
In figures 5, 6, and 7, we have plotted typical results for numerical solutions to (1.1), corresponding to RS, SS, and RR configurations. The resulting values of e and j seem to approach the predicted values of e 0 and j 0 at long times -at least as far as our numerical precision can be trusted (see appendix B). In particular, in the RS case, we approach the steady state value (4.10); in the SS case, we approach (4.13); and in the RR case, we approach (4.16). As we discuss in greater detail in the next section, one place where gradient effects show up and do not disappear as a function of time is in the shock width.
One may speculate that the agreement between the predicted steady state in the absence of gradient corrections and the numerical results is associated to the fact that the gradient corrections, even though order one in our system of units, come with dimensionful coefficients. In the language of the renormalization group, they conform to irrelevant couplings. Perhaps it is for this reason that at long enough time and in a large enough box, we may be able to ignore these corrections for the most part. JHEP08(2016)120 Figure 5. A numerical solution to the Riemann problem. The plots were obtained starting with an initial condition (B.5) with L = 8000, c = 300 and j = 0. Only one half of the box, centered around the origin, is depicted. The dashed curve corresponds to values of e and j at t = 0 while the solid curve corresponds to values of e and j at t = 800. The black, red and blue horizontal lines correspond to the predicted near equilibrium steady state associated with a rarefaction wave and shock pair (cf., equation (4.10)), a bad shock and good shock pair (cf., references [5,7]), and a non thermodynamic shock pair (cf., reference [5]) respectively. The fixed point associated with a rarefaction solution which exists for δe ≥ 0.7536 . . . is represented by a black dot.

Restoring gradient corrections
In this section, we try to gain a better handle over the gradient corrections and their affect on the predicted steady state values. The analysis here is incomplete and approximate. To overcome the deficiencies of paper and pencil estimates, we include some numerical solutions to the conservation equations (1.1) that provide support for the estimates. We will consider separately corrections to each of the features we found in the idealized limit: the steady state and asymptotic regions with constant e and j, a shock wave, a rarefaction wave, and the discontinuity at the edge of the rarefaction.
Corrections to constant regions. Corrections to a constant e and j region are easiest to analyze. Assuming the fluctuations are small, we look for linearized solutions of the form e = e 0 +δe exp(−iωt+ikζ) and j = j 0 +δj exp(−iωt+ikζ). We find two propagating modes

JHEP08(2016)120
These two modes are damped sound modes whose speed is shifted by the fluid velocity β = j/e. The gradient corrections appear here in the form of the damping term ik 2 in the dispersion relation. Given this result, we anticipate that we will be able to correct a constant e and j region by taking an appropriate linear superposition of sound waves. The damping suggests that at long times the solution can only involve constant e and constant j. As a side comment, an odd thing about these mode relations is that they are exact. Recall that in first order viscous hydro, we would typically solve an equation of the form ω 2 + iΓk 2 ω − k 2 = 0 for ω, in the case of vanishing background fluid velocity. If this equation were treated as exact, the solutions for ω would be non linear in k and therefore have higher order contributions, i.e. O(k 3 ), O(k 4 ), etc., when expanded around small k.
Corrections to shocks. The gradient corrections should act to smooth a shock and give it some characteristic width. We estimate this width in a frame in which the shock is not moving, i.e. s = 0. In this frame, j r = j l and e r e l = j 2 l . We can find a solution for the shock profile in the case where the shock is weak e r ∼ e l : where we have defined e ≡ e r + e l 2 , δe ≡ e r − e l e r + e l , and j ≡ j r + j l 2 .
We can see in figure 8 that even for values of δe ∼ 1/2, that e δe 2 /2 appears to be a good estimate for the slope of the shock. 6 In appendix A, we show that this shock profile produces, at the correct subleading order in a large d expansion, the correct (positive) amount of entropy predicted by the RH relations.
Corrections to a rarefaction. We will perform two estimates of gradient corrections to the rarefaction wave. The first estimate is a correction to the interior of the wave far from the edges where it joins onto constant e and j regions. The second estimate is a correction to the discontinuity where the rarefaction joins a constant region. For the first estimate, we assume an ansatz for the long time behavior of the rarefaction wave: 6 We found that when δe = 0.8 the relative error between (4.24) and the numerical solution grew to ∼ 13%. As δe gets closer to one numerical error is more difficult to control. JHEP08(2016)120 Figure 8. A numerical simulation of stationary shocks. We start from an initial condition e = e (1 + δe tanh(c sin(2πx/L))), j = 1 with parameters L = 8000 and c = 1.2(Lδe/4π). We chose e r and e l to produce a stationary shock (e l = √ 1−δe √ 1+δe , e r = √ 1+δe √ 1−δe ) using the RH relations. We then plot the value of the slope of the shock after the system has settled into a steady state. This is compared with the weak shock solution (4. With an appropriate choice for the integration constant c 1 , the expressions for e 0 and j 0 become the same as we had before (4.9). There are subleading corrections that scale as 1/t and log(t)/t that depend on a second integration constant c 2 and an arbitrary function e 1 (ξ), both presumably set by the initial conditions. Note that the combination ξe − j is independent of the arbitrary function e 1 (ξ) at order 1/t. In figure 9, the numerics confirm that the corrections to ξe − j do indeed scale as 1/t. Last, we would like to heal the discontinuity at the edge of a rarefaction wave. The tanh function we found above heals the discontinuity in the shock case, making the question of what happens at the edge of a shock less pressing. Consider a case where the rarefaction wave meets a steady state at ζ = 0, with the rarefaction region to the right and the steady state to the left. (We can always move the meeting point away from ζ = 0 by boosting the solution ζ → ζ +vt.) With the intuition that the second order gradients in the conservation equations are dominant and render the behavior similar to that of a heat equation with JHEP08(2016)120 Figure 9. A plot of δ(ξe − j) vs. time at three different points in a single rarefaction wave. The quantity δ(ξe − j) is the difference between the zeroth order prediction (4.9) and numerics. The rarefaction wave spreads from ξ l = −1 to ξ r = 1. The three points correspond to ξ = −1/2 (red), ξ = 0 (purple) and ξ = 1/2 (green). The dashed line 1/(2t) is a guide to the eye. Inset: the rarefaction profile at t = 3000. Dashed lines correspond to e while the solid lines correspond to j. The blue curve is numeric, while the red curve is the ideal result (4.9).
Note that the relation j 0 = ±e 0 is consistent with a rarefaction meeting a steady state region at ζ = 0. These relations for the j i lead to a second order, nonlinear differential equation for e 1 : Remarkably, this equation can be written as a total derivative and integrated to yield where c 1 is another integration constant. The integration constants reflect a translation symmetry of both e 1 and χ. We can shift χ → χ + j 1 /e 0 and e 1 (χ) → e 1 (χ − j 1 /e 0 ) ± j 1 /2. JHEP08(2016)120 The shifts send j 1 → 0 and c 1 → c 1 ∓ 3j 2 1 /8e 0 in the equation (4.32). If we apply the boundary condition that both e 1 (χ) and e 1 (χ) vanish in the steady state region χ → −∞, then we must set c 1 = 0, and the resulting first order differential equation becomes separable. To match onto the rarefaction region, we require that e 1 → ±e 0 as χ → ∞. This boundary condition fixes the remaining integration constant associated with the first order equation (4.32), and the solution for e 1 is then 2e 0 e −χ 2 /4 √ π erfc(χ/2) . (4.33) As we choose the rarefaction region to match onto the steady state at χ = 0, we conclude that the integration constant j 1 in the original differential equation must be zero as well. We can check numerically that a 1/ √ t scaling is consistent with the behavior at the endpoints of a rarefaction solution. See figure 10.

Discussion
We presented a solution to the Riemann problem for the conservation equations (1.1). Through fluid-gravity and the AdS/CFT correspondence, these equations describe, in a large d limit, both the dynamics of a black hole horizon and also the dynamics of a strongly interacting conformal field theory.
There are a number of possible future directions for research. The simplest is perhaps to include a transverse velocity. With a transverse velocity, in addition to the shock and rarefaction waves, there will in general be a contact discontinuity [13,[43][44][45]. It is known JHEP08(2016)120 (and perhaps intuitive given the similarity to a counter flow experiment), that the contact discontinuity is in general unstable to the development of turbulence [46]. It would be interesting to see what precisely happens in our large d limit. Another more complicated extension is the inclusion of a conserved charge. The large d equations of motion in the presence of a conserved charge are available from ref. [14]. Once again, a contact discontinuity is expected (see for example [13]) although whether such a discontinuity is stable or unstable to turbulence is unclear. More ambitiously, one could consider what happens for the holographic dual of a superfluid or superconductor [19,25,[47][48][49][50][51].
Another possible direction is the addition of higher curvature terms to the dual gravitational description. One could presumably tune the d dependence of these terms such that higher order gradient corrections appear in the conservation equations (1.1) and also such that the first and second order transport coefficients are tuned away from the values examined in this paper.
Perhaps the most interesting direction for future study is the connection to black hole dynamics. What can we learn about black holes through the connection to hydrodynamics in a large d limit?

JHEP08(2016)120
Equation (A.3) can be obtained by using a large d expression for the entropy current (3.22) along with the Rankine-Hugoniot relations for energy and momentum, (4.1) supplemented by (2.14) and (2.15). Note that in the asymptotic regions, the gradient terms will all vanish. (It is also possible to start with a finite d result, using for example refs. [12] or [13], and then take a large d limit directly.) The non-conservation of entropy (A.3) can be captured by the leading viscous corrections to the shock width (4.24) when the energy difference is small. Indeed, using (3.24) Integrating this divergence over the ζ direction leads to

B A bestiary of plots
In section 4.3 we studied the numerical solutions to the Riemann problem for various initial energy and velocity profiles associated with RR, RS and SS type solutions. In what follows we provide additional evidence that at late times the full numerical solution to the Riemann problem approaches the appropriate predicted steady state values e 0 and j 0 and fixed point values e s and j s .

B.1 RR configurations
To generate an RR configuration we used the initial data Once j * ≥ 2 one should find a fixed point with e s = j s = exp(−1). We find that the numerical solution approaches the predicted states via power law behavior, see figure 11. JHEP08(2016)120

B.2 SS configurations
To generate an SS configuration we used the initial data (B.1) with j * < 0. The analysis of section 4.2 predicts a steady state of the form

B.3 RS configurations
To generate an RS configuration we used the initial data exp(1) we will obtain a fixed point at the origin with e s = j s = exp(−1). An analysis of the late time behavior of the numerical solution can be found in figure 13.

B.4 Error analysis
In sections B.1 and B.3 we have fit the late time approach of the data to the predicted steady state and (or) fixed point values to a power law behavior. The fit was done using Mathematica's NonLinearModelFit routine [42]. In detail, the late time data was discretized into order 1 time steps which were then fit to a a/t α curve with a and α as parameters. The standard errors for the fit were usually of order 10 −3 to 10 −4 . Fits involving very JHEP08(2016)120 Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.