Phase Transitions, Inhomogeneous Horizons and Second-Order Hydrodynamics

We use holography to study the spinodal instability of a four-dimensional, strongly-coupled gauge theory with a first-order thermal phase transition. We place the theory on a cylinder in a set of homogeneous, unstable initial states. The dual gravity configurations are black branes afflicted by a Gregory-Laflamme instability. We numerically evolve Einstein's equations to follow the instability until the system settles down to a stationary, inhomogeneous black brane. The dual gauge theory states have constant temperature but non-constant energy density. We show that the time evolution of the instability and the final states are accurately described by second-order hydrodynamics. In the static limit, the latter reduces to a single, second-order, non-linear differential equation from which the inhomogeneous final states can be derived.


I. INTRODUCTION
Hydrodynamics is one of the most successful theories in physics, capable of describing the dynamics of very different systems over an enormous range of length scales. Traditionally, it is thought of as an effective theory for the conserved charges of a system, constructed as a derivative expansion around local thermal equilibrium. From this perspective hydrodynamics is only expected to be valid when these gradient corrections are small compared to all the microscopic scales. However, in recent years it has been discovered that the regime of validity is actually much broader.
In this paper we use the gauge/string duality to test the validity of hydrodynamics in a theory with a first-order thermal phase transition. We place the theory on a cylinder in a variety of homogeneous, unstable initial states. In the gauge theory the instability is a spinodal instability associated to the presence of a first-order phase transition. On the gravity side the instability is a long-wave length, Gregory-Laflamme-type instability [30]. We use the gravity description to follow the evolution of these states until they settle down to a stationary, inhomogeneous final state. We then compare both the time evolution and the final configurations to the corresponding hydrodynamic predictions. Recent related work includes [31][32][33].
To the best of our knowledge, from the boundary theory viewpoint our results provide the first example of a first-principle, non-perturbative, complete dynamical evolution of a spinodal instability in a strongly coupled gauge theory. From the gravity viewpoint, they provide the first example of the full dynamical evolution of a Gregory-Laflamme-type instability in an asymptotically anti-de Sitter space.

II. MODEL
Motivated by simplicity, we study the Einstein-plus-scalar model with action The potential is with the asymptotic curvature radius and φ M = 2.3. The dual gauge theory is a Conformal Field Theory (CFT) deformed with a dimension-three scalar operator with source Λ, which appears as a boundary condition for the scalar. This is exactly the potential of [34] except for the sign reversal of the fourth term in V . This difference has dramatic implications for the thermodynamics of the gauge theory, which we extract from the homogeneous black brane solutions of the gravity model. In particular, the gauge theory possess a first-order phase transition at a critical temperature T c = 0.247Λ, as illustrated by the multivalued plot of the energy density as a function of the temperature in Fig. 1. States on the dashed red curve are locally thermodynamically unstable since the specific heat is negative, c v = dE/dT < 0. In this region the speed of sound, c 2 s = s/c v , with s the entropy density, becomes imaginary. This leads to a dynamical, spinodal instability (see e.g. [35]) whereby the amplitude of small sound excitations grows exponentially with a momentumdependent growth rate dictated by the sound dispersion relation: where η and ζ are the shear and bulk viscosities. In our model η/s = 1/4π [36] and we compute ζ numerically following [37]. The corresponding statement on the gravity side is that the black branes dual to the states on the dashed red curve are afflicted by a long-wave length instability. Although this is similar to the Gregory-Laflamme instability of Ref. [30], we will see below that there are also crucial differences between the two.

III. INHOMOGENEOUS HORIZON
To investigate the fate of the spinodal instability we compactify the gauge theory direction z on a circle of length L 57/Λ. This infrared cut-off reduces the number of unstable sound modes to a finite number. We then consider a set of homogeneous, unstable initial states with energy densities in the range E/Λ 4 (0.002, 0.016) 1 , as indicated by the grey, dashed horizontal lines in Fig. 1. For concreteness we will show results for the state with E/Λ 4 0.0096, whose temperature and entropy density are T i 0.251Λ and s i = 0.037Λ 3 . To trigger the instability, we introduce a small z-dependent perturbation in the energy density corresponding to a specific Fourier mode on the circle. For concreteness we will show results for the case with k = 3(2π/L) 1.3T i . This mode is unstable with positive growth rate Γ = 0.0247 Λ according to (3). For numerical simplicity we impose homogeneity along the transverse directions. 2 On the gravity side the sound-mode instability may be viewed [39][40][41] as a Gregory-Laflamme instability [30]. We follow it by numerically evolving the Einstein-plus-scalar equations as in [27,38] until the system settles down to a state with an inhomogeneous Killing horizon with constant temperature T f = 0.250Λ. Note that this is close but not identical to T i or T c . From the dynamical metric we extract the boundary stress tensor. The result for the energy density is shown in Fig. 2. The time dependence of the amplitudes of several modes of the energy density is shown in Fig. 3. The n = 3 mode grows with a rate that agrees with (3) within 4%. Resonant behavior makes the modes with n = 6, 9, . . . grow at a rate that is roughly the corresponding multiple of the n = 3 rate. Numerical noise makes some non-multiples of the n = 3 mode (of which only three are shown in the figure) grow too. At late times the non-linear dynamics stops the growth of all these modes and the system settles down to a stationary, inhomogeneous configuration consisting of three identical domains. The fact that their size is comparable to the inverse temperature, ∆z = L/3 0.4/T f , is our first indication that spatial gradients are large.
In Fig. 4 we show the final entropy density as extracted from the area of the horizon. The fact that s is not constant proves that the horizon itself (not just the boundary energy density) is inhomogeneous. We see in the figure that s is very well estimated by a point-wise application of the equation of state to the final energy density of Fig. 2, suggesting that the evolution is quasiadiabatic. This is confirmed by the fact the final average entropy density is only 1% larger than the initial one.
In the paragraphs above we have been careful to use the term "stationary" as opposed to "static" to refer to the final, time-independent state to which the system settles down. The reason is that the final metric is stationary but not static. The physical reason for this is that the dtdz off-diagonal terms in the final metric do not vanish. As a consequence, the final state is time-independent but not time-reversal invariant (although it is invariant under the simultaneous sign reversal of both t and z). The fact that this is a true property of the solution and not just a consequence of our metric ansatz follows from the fact that the timelike Killing vector of our solution is not hypersurface-orthogonal. Interestingly, the g tz component of the final-state metric falls off quickly enough near the AdS boundary so that there is no energy flux in the z-direction. In other words, the one-point function of the T t z component of the boundary stress tensor vanishes identically. As a consequence, at the level of one-point functions the final state on the gauge theory side looks not only stationary but static. Since hydrodynamics is an effective theory for one-point functions, in the next section we will speak of a (hydro)static final state. However, it must be kept in mind that the full microscopic state in the gauge theory is only stationary, since n-point functions with n > 1 would be sensitive to the properties of the dual metric deep in the bulk and hence would reveal the lack of staticity.

IV. HYDROSTATIC FINAL STATE
Previous holographic studies have established that hydrodynamics can describe the evolution of the stress tensor in the presence of large gradients. Here we investigate whether it can also describe the static inhomogeneous configuration at the endpoint of the spinodal instability.
In the hydrodynamic approximation the stress tensor is expanded in spatial gradients in the local fluid rest frame. Up to second order, the constitutive relations are where, in the fluid rest frame, T ideal µν = Diag{E, P eq (E)}, P eq (E) is the equilibrium pressure of the homogeneous states shown in Fig. 1, σ µν and Π are the shear and bulk stresses, and ∆ µν is the projector onto the spatial directions. The tensor Π (2) µν contains all the second-order terms. In a non-conformal theory this tensor contains thirteen E-dependent second-order transport coefficients and its explicit expression may be found in [42]. A subset of these coefficients for the model of [34] has been computed in [43].
In a static configuration the fluid three-velocity field is identically zero, the stress tensor is diagonal, both σ µν and Π vanish, and the leading gradient corrections are those in Π (2) µν . This tensor also simplifies since only two independent second-order terms survive in this limit. In the Landau frame, the constitutive relations reduce to The statement that hydrodynamics describes the final states is the statement that the corresponding pressures are given by these equations with second-order transport coefficients c L, T (E) and f L, T (E) that depend on the local energy density 3 but not on any other details of the final states. We have verified this state independence by varying both the average energy densities within the range shown in Fig. 1 and the length of the circle L. Equivalently, a hydrodynamic description of the final states means that, once c L, T (E) and f L, T (E) have been extracted from a given state (or computed microscopically), they can be used to predict the pressures of a different state given its energy density profile. This is illustrated in Fig. 5, where we compare the true pressures obtained from the gravity side for the end state of Fig. 2 with the hydrodynamic predictions (5) and (6) based on coefficients extracted from a different state. It is remarkable that hydrodynamics works despite the fact that the second-order terms in these equations are as large as the equilibrium pressure, as can be seen in the figure from the difference between the continuous gray curve and the true pressures. This is particularly dramatic for the longitudinal pressure, P L , which must be z-independent since, in a static state, conservation of the stress tensor reduces to ∂ z T zz = 0. While the equilibrium contribution to P L inherits z dependence from the energy density, this modulation is precisely compensated by the second-order contribution. We conclude that the second-order gradients sustain the inhomogeneous state.
Taking the applicability of hydrodynamics one step further, we can use it to predict all other static, inhomogeneous configurations. 4 Time independence implies that P L must be constant, which in the hydrodynamic approximation reduces to a second-order, non-linear differential equation for E(z) via Eq. (5). This depends on two integration constants that can be traded for the length of the circle (more precisely, the size of the domains) and the average energy density. Once E(z) is known P T is predicted by Eq. (6).

V. HYDRODYNAMIC EVOLUTION
The second-order coefficients that we have extracted also control the dynamical evolution of the instability. However, unlike in the stationary final configuration, during the dynamical evolution there are non-zero momentum fluxes, leading to a small but non-vanishing three-velocity field. To approximate the full time evolution by hydrodynamics it is thus necessary to include the firstorder shear and bulk tensor contributions in Eq. (4). Similarly, while the systems evolves there are additional second-order gradients beyond those considered in Eqs. (5) and (6). However, these additional second-order terms are quadratic in the velocity field and can be consistently neglected.
In Fig. 6 we show the time evolution of the pressures at z Λ = 0. The very early time behavior is the exponential decay of the quasi-normal modes (QNM) that are excited by the perturbation that we introduce to trigger the instability. After this short time, the predictions of second-order hydrodynamics match the true pressures. Second-order terms become increasingly important as the instability saturates, where the first-order approximation, let alone the equilibrium pressure, fails to predict the true pressures while the second-order approximation continues to do so accurately. We conclude that hydrodynamics with large second-order gradients describes the evolution and the saturation of the spinodal instability.

VI. DISCUSSION
We have uncovered a new example of the applicability of hydrodynamics to systems with large gradients. As in other known examples, part of the reason behind this success is the relaxation of the non-hydrodynamic QNMs at the relevant time scales of the evolution, as in the early transient behaviour observed in Fig. 6. The analysis of the exponential decay of these modes reveals that this process occurs over a time Γ QNM ∼ c πT with c 3.4, consistent with [44]. In contrast, for the unstable hydrodynamic mode, assuming ζ/η ∼ 1 and η/s = 1/4π, Eq. (3) yields the typical growth rate Γ ∼ |c s | 2 πT , which in our case is suppressed with respect to Γ QNM by the small value of |c s | 2 0.03. Although this argument may explain why hydrodynamics provides a good description at intermediate times, once the QNMs have decayed, it remains surprising that it also describes the late-time evolution and the final state, where the spatial gradients are large. Moreover, it suggests that it would be interesting to explore other values of the parameter φ M for which the first-order transition becomes stronger and |c s | can become of order unity. Our static inhomogeneous configurations do not describe the separation of the system into two stable phases. For example, the maxima of the end state of Fig. 2 do not lie on the green stable branch but on the red unstable branch of Fig. 1. Presumably the reason is that, if the available average energy density is small, the cost of the necessary gradients for the peaks to reach the green stable branch makes the corresponding configuration disfavored. Similar considerations seem to be responsible for the fact that the final state of Fig. 2 exhibits three identical domains. We have found another static configuration with five peaks which is also well described by hydrodynamics but whose entropy is smaller than that of the three-peak state. Moreover, we have studied the dynamics of the system when perturbed by an n = 1 mode instead of by an n = 3 mode. Despite the fact that we have followed the evolution for very long times we have found that this configuration does not seem to settle down to a stationary state; nevertheless, we have observed that throughout its non-linear evolution it develops a second peak. This suggest that a phase separated configuration that naively one would expect to reach by starting with the n = 1 mode may simply not exist for the average energy density and the box size that we have considered. A logical possibility compatible with all the above is therefore that, under these conditions, the three-peak configuration that we have presented is the thermodynamically preferred state. We emphasize, however, that establishing this definitively would require comparing the entropies of all possible inhomogeneous states of the system with the same average energy density and box size, which is certainly beyond the scope of this paper.
We suspect that the reason why the phase-separated configuration is seemingly not allowed for the box size that we have considered is the fact that the two stable phases in our phase diagram in Fig. 1 differ by several orders of magnitude. This suggests that, in order to realise the phaseseparated configuration, one may need a much larger box.
As pointed out above, from the gravity viewpoint the gauge theory spinodal instability translates into an instability of the dual black brane. This is similar to the Gregory-Laflamme instability of a black string in five-dimensional flat space [30] in that it is a long wavelength instability associated to the sound mode of the system that appears below a certain mass or energy density. However, there is also one crucial dissimilarity that makes the time evolution and the final state of our system and those of a black string [45] completely different. This is the fact that the black string is unstable for any mass density below a certain critical value, whereas in our case the black brane is only unstable in a certain energy band (the red dashed curve in Fig. 1). In the case of the black string, at intermediate stages in the evolution the horizon can be described as a sequence of three-dimensional spherical black holes joined by black string segments. Since these segments are themselves subject to a Gregory-Laflamme instability, the evolution results in a self-similar cascade, where ever-smaller satellite black holes form connected by ever-thinner string segments. In contrast, in our case this evolution stops (roughly speaking) once the energy density at a certain point is low enough to lie in the stable region of the phase diagram, i.e. below the dashed red curve in Fig. 1.