Crossing a large-$N$ phase transition at finite volume

The existence of phase-separated states is an essential feature of infinite-volume systems with a thermal, first-order phase transition. At energies between those at which the phase transition takes place, equilibrium homogeneous states are either metastable or suffer from a spinodal instability. In this range the stable states are inhomogeneous, phase-separated states. We use holography to investigate how this picture is modified at finite volume in a strongly coupled, four-dimensional gauge theory. We work in the planar limit, $N \to \infty$, which ensures that we remain in the thermodynamic limit. We uncover a rich set of inhomogeneous states dual to lumpy black branes on the gravity side, as well as first- and second-order phase transitions between them. We establish their local (in)stability properties and show that fully non-linear time evolution in the bulk takes unstable states to stable ones.


Introduction
Phase coexistence is an essential feature of systems with a first-order phase transition. Consider for example Fig. 1. This shows the energy density as a function of the temperature, in the infinite-volume limit, for the four-dimensional gauge theory that we will study in this paper. The blue curve indicates homogeneous states with energy density E, which we measure in units of the microscopic scale in the gauge theory Λ. In the canonical ensemble there is a first-order phase transition at a critical temperature T c indicated by the dashed, vertical line. The thermodynamically preferred, lowest-free energy states at T > T c lie on the upper branch and have energies above that of point D. Similarly, at T < T c the preferred states are on the lower branch with energies below that of point E. States between points A and D, and between B and E, are locally but not globally thermodynamically stable. Finally, states between points A and B are locally thermodynamically unstable. The region between A and B is known as the "spinodal region".
In the canonical ensemble, setting T = T c does not select a unique state. For this reason, it is convenient to work in the microcanonical ensemble, in which the control The blue curve corresponds to homogeneous states. The red line corresponds to inhomogeneous, phase-separated states. Solid segments indicate locally dynamically stable states; dashed segments indicate unstable ones. The black curves with arrows indicate the sequence of maximum-entropy states as the average energy density E decreases. In the canonical ensemble there is one first-order phase transition at which the system jumps betweens points D and E. In the microcanonical ensemble there are two second-order phase transitions at points D and E between homogeneous and inhomogeneous states.
parameter is the energy instead of the temperature. In this case the preferred, maximumentropy configuration for energy densities between points D and E is well understood in the infinite-volume limit: it is a phase-separated state in which part of the volume is in the phase associated to point D and the other part is in the phase associated to point E (see Sec. 2.4.1). The fraction of volume occupied by each phase is determined by the average energy density E, which lies between D and E. The two phases are separated by a universal interface, i.e. by an interface whose spatial profile is independent of the way in which the phase-separated configuration is reached. Since the temperature is constant and equal to T c across the entire volume, these states lie on the red, vertical segment DE in Fig. 1. We conclude that, at infinite volume in the microcanonical ensemble, the sequence of preferred states as the energy density decreases is that indicated by the black arrows in Fig. 1. The thermodynamic statements above have dynamical counterparts. Since the total energy is conserved under time evolution, it is again convenient to think of the system in the microcanonical ensemble. Imagine preparing the system in a homogeneous state. If the energy density lies above point D or below point E then this state is dynamically stable against small or large perturbations. If instead the energy density is between points A and D or between B and E then we expect the system to be dynamically stable against small perturbations but not against large ones. This means that, if subjected to large enough a perturbation, the system will dynamically evolve to a phase-separated configuration. The average energy density in this inhomogeneous configuration will be the same as in the initial, homogenous state, but the entropy will be higher. Finally, if the initial energy density is between A and B then the state is dynamically unstable even against small perturbations. This instability, known as "spinodal instability", implies that the slightest perturbation will trigger an evolution towards a phase-separated configuration of equal average energy but higher entropy.
If the system of interest is an interacting, four-dimensional quantum field theory then following the real-time evolution from an unstable homogeneous state to a phase-separated configuration can be extremely challenging with conventional methods. For this reason, in [1,2] holography was used to study this evolution in the case of a four-dimensional gauge theory with a gravity dual (see also [3,4] for a case in which the gauge theory is three-dimensional). In order to regularise the problem, Refs. [1,2] considered the gauge theory formulated on R 1,2 × S 1 with periodic boundary conditions on a circle of size L. For simplicity, translational invariance along the non-compact spatial directions was imposed, thus effectively reducing the dynamics to a 1+1 dimensional problem along time and the compact direction. The compactness of the circle makes the spectrum of perturbations discrete and simplifies the technical treatment of the problem. Ref. [1] provided a first example of the time evolution from a homogeneous state to an inhomogeneous one. A systematic study was then performed in [2]. In this reference the focus was on the infinitevolume limit, understood as the limit in which L is much larger than any other scale in the problem such as the microscopic gauge theory scale Λ, the size of the interface, etc. It was shown that, if slightly perturbed, an initial homogeneous state with energy density between A and B always evolves towards a phase-separated configuration, and that the latter is dynamically stable.
In addition to its implications for gauge theory dynamics, the spinodal instability of states between A and B is interesting also on the gravity side, where it implies that the corresponding black branes are afflicted by a long-wavelength dynamical instability. Although this is similar [5][6][7] to the Gregory-Laflamme (GL) instability of black strings in spacetimes with vanishing cosmological constant [8], there is an important difference: In the GL case all strings below a certain mass density are unstable, whereas in our case only states between points A and B are unstable. Having clarified this, since the term "GL-instability" is familiar within part of the gravitational community, in this paper we will use the terms "spinodal instability" and "GL-instability" interchangeably to refer to the dynamical instability between points A and B.
The purpose of this paper is to extend the analysis of the equilibrium states summarised in Fig. 1, as well as the systematic analysis of their dynamical stability properties of [2], to the case of finite volume. In particular, we would like to: (i) classify all possible states, homogeneous or inhomogeneous, available to the system; (ii) determine which ones are thermodynamically preferred; (iii) establish the local dynamical stability or instability of each state; and (iv) investigate the time evolution from unstable states to stable ones. For this purpose we will place the system in a box, impose translational invariance along two of its directions and vary the size L of the third direction. We will then see that the results depend on the value of L compared to a hierarchy of length scales (1.1) These three scales are an intrinsic property of the system at finite volume that cannot be determined through an infinite-volume analysis. Depending on the ratio of L to these scales we will uncover: (i) a large configuration space of inhomogeneous states, both stable and unstable; (ii) a rich set of first-and second-order thermodynamic phase transitions between them; and (iii) the possible time evolutions from dynamically unstable to dynamically stable states. Note that the existence of phase transitions is not in contradiction with the finite volume of the gauge theory because we work in the planar limit, N → ∞, which effectively acts as a thermodynamic limit. Our results are summarised in Sec. 4. The reader who is only interested in this summary can go directly to this section.
2 Nonconformal lumpy branes: nonlinear static solutions 2.1 Setup of the physical problem and general properties of the system We consider the AdS-Einstein-scalar model with action where κ 2 = 8πG 5 with G 5 Newton's constant, g is the determinant of the metric g µν , R is the associated Ricci scalar and φ is a real scalar field. The potential V (φ) can be derived from the superpotential through the usual relation The positivity of energy theorem for the AdS-Einstein-Scalar model is subtle [9]. For a given potential, one can find up to two superpotentials W ± that satisfy (2.3). One way of distinguish these two possible solutions is to inspect the small φ behaviour of the superpotential, namely where ∆ + = 3 and ∆ − = 1 for our potential. It is not a coincidence that ∆ ± correspond to the conformal dimensions of the operator dual to φ in standard, and alternative quantisation, respectively. One can show that if W − exists globally, then so does W + [9]. Our superpotential (2.2) is of the W − form irrespectively of the value of the dimensionless parameters φ M and φ Q . For a sourced solution such as ours, the existence of W − ensures that all solutions of our model have positive energy [9]. For any value of φ M and φ Q the potential V (φ) has a maximum at φ = 0, corresponding to an ultraviolet (UV) fixed point of the dual gauge theory. We will choose the values φ Q = 10 and φ M = 1, for which W and V take the form shown in Fig. 2. In this case both functions have a minimum at φ min ≈ 1.54, corresponding to an infrared (IR) fixed point of the gauge theory. The potential has an additional maximum at φ max ≈ 3.65 and diverges negatively, i.e. V (φ) → −∞, as φ → +∞. However, values of φ larger than φ min will play no role in our analysis. Our motivation to choose this model is simplicity. The superpotential (2.2) is the same as in [1,2,10] except for the φ 6 term, which was absent in those references but was introduced in [11]. As in [1,2,10,11], the dual gauge theory is a Conformal Field Theory (CFT) deformed by a dimension-three scalar operator with source Λ. On the gravity side this scale appears as a boundary condition for the scalar φ. The first two terms in the superpotential are fixed by the asymptotic AdS radius and by the dimension of the dual scalar operator. As in [1,2,10], the present model also possesses a first-order phase transition. In those references this was achieved by choosing the value of φ M appropriately, with no need to include a φ 6 term. However, this leads to a phase diagram in which the energy densities at points D and E differ from one another by three orders of magnitude. This huge ratio makes the numerical treatment of the system extremely challenging. In contrast, by including the φ 6 term as in [11] and choosing the values of φ M and φ Q as quoted above, this ratio is of order unity, as is clear from Fig. 3.
We are interested in finding static, "lumpy" black brane solutions of (2.1) that can break translational invariance along a spatial gauge theory directionx while being isometric along the remaining two spatial directions x 2 and x 3 directions. The most general ansatz compatible with such symmetries is where Z is the holographic coordinate. We shall be interested in solutions which are isotropic in x 2 and x 3 , so we take Q 5 = Q 6 and Q 7 = 0 . To fix the gauge completely we demand Q 3 = 0, and together with the condition q 1 (x, Z + ) 2 q 2 (x, Z + ) = α Λ , (2.6f) where α Λ > 0 is a positive constant whose physical significance will be discussed later; see (2.27). The coordinatex is periodic with period L, and we takex ∈ [−L/2, L/2] and Z ∈ [0, Z + ]. The translationally invariant directions x 2,3 have arbitrary periods L 2 and L 3 . These will play essentially no role in our discussion, since only densities per unit area in the 23-plane will matter. Thus we will take them to be the same, i.e. L 2 = L 3 ≡ L and We shall also be interested in solutions which are Z 2 -symmetric aroundx = 0, which means that we can restrict our domain of integration tox ∈ [0, L/2], at the expense of imposing ∂xq j x=0 = 0, for j = 1, 2, 3, 4. In order to vary L in a numerically efficient manner, we further change to a new coordinate and take all functions to take values in x ∈ [0, 1]. Note that our ansatz (2.8) together with our periodicity conditions further imply that ∂ x q j x=1 = 0 (j = 1, 2, 3, 4).
Putting everything together brings (2.5) to the following simplified form Our gauge choice is such that the determinant of the metric along the Killing directions (t, x 2 , x 3 ) is fixed and defines the radial (holographic) direction Z. The conformal boundary is located at Z = 0 where we demand In this sense we can denote this gauge choice as the "double Wick rotation Schwarzschild gauge" and, as far as we are aware, this is the first time it is introduced. The advantage of this gauge choice (at least in the present system) is that the fields q j (j = 1, 2, 3, 4) have an asymptotic power law decay without irrational powers (nor logarithmic terms; more below), unlike e.g. the DeTurck gauge. 1 Our gauge choice -condition (2.6f) -reveals that Z = Z + is a null hypersurface, where the norm of the Killing vector field ∂/∂t vanishes. 2 Thus, Z = Z + is a Killing horizon, and α Λ controls its associated surface gravity or, equivalently, temperature. In fact, we find (2.10) Our solutions have two important scaling symmetries. The first one is {t, Z, x i } → {λ 1 t, λ 1 Z, λ 1 x i }, {q 1,2,3 , q 4 } → {q 1,2,3 , λ −1 This leaves the equations of motion and scalar field invariant and rescales the line element as ds 2 → λ 2 1 ds 2 , namely g µν → λ 2 1 g µν . It follows that we can use this scaling symmetry to fix the AdS radius to ≡ 1. In other words, under the scaling g µν → λ 2 1 g µν , the affine connection Γ γ µν , and the Riemann (R α βµν ) and Ricci (R µν ) tensors are left invariant. It follows from the trace-reversed equations of motion that the AdS radius must scale as → λ 1 and we can use this scaling symmetry to set ≡ 1.
The second scaling symmetry (known as a dilatation transformation, one of the conformal transformations) is This leaves the metric, scalar field and equations of motion invariant. It follows that we can use this symmetry to set the horizon radius at Z = Z + ≡ 1 or, equivalently, the temperature (2.10) to (2.13) We will see below that this is just a convenient choice of units with no effect on the physics. Let us turn now our attention to the scalar field. It follows from (2.3) and (2.2) that the scalar field potential has the Taylor expansion about φ = 0 and thus it describes a scalar field with mass µ 2 = V (0) = −3/ 2 . According to AdS/CFT, the conformal dimension of the dual operator is simply given by 15) and these give the two independent asymptotic decays Z ∆ ± of the scalar field. Actually, since ∆ ± are integers, the nonlinear equations of motion might also generate logarithmic decays of the form n=3 c n Z n ln Z where the coefficients c n depend exclusively on the amplitude of the two independent terms. When this is the case, the conserved charges depend on c 3 . However, for the potential we use (only with even powers of φ and double Wick rotation Schwarzschild gauge choice) it turns out that logarithmic terms are not generated by the equations of motion. So, for our system, the scalar field decays asymptotically as where Λ and φ 2 are two arbitrary integration constants and · · · represent higher order powers of Z (with no logarithms) whose coefficients are fixed in terms of Λ and φ 2 by the equations of motion. The fact that ∆ − = 1 motivates our choice of ansatz for the scalar field in (2.8). The Breitenlöhner-Freedman (BF) bound of the system is µ 2 BF 2 = −4 and thus µ 2 2 = −3 coincides precisely with the unitarity bound µ 2 BF 2 + 1. It follows that only the mode Z ∆ + with the faster fall-off is normalizable. In the AdS/CFT correspondence, the non-normalizable mode Λ is the source of a boundary operator O φ since it determines the deformation of the boundary theory action. On the other hand, the normalizable modes φ 2 are identified with states of the theory with φ 2 being proportional to the expectation value O φ of the boundary operator (in the presence of the source Λ). ∆ + = 3 is then the (mass) conformal dimension of the boundary operator O φ dual to φ. Under the scaling symmetries (2.11)-(2.12) the scalar source transforms as Λ → λ 1 λ 2 Λ. As a consequence, the ratio T /Λ is left invariant by these scalings. Since the physics only depends on this ratio, setting T = √ α Λ /π as we did above is just a convenient choice of units with no effect on the physics. In general, throughout this paper we will measure all dimensionful physical quantities in units of Λ. The undeformed boundary theory -a CFT -corresponds thus to the Dirichlet boundary condition Λ = 0, and we have a pure normalizable solution. For this reason the planar AdS 5 Schwarzschild solution (2.8) with q 1,2,3 = 1, q 4 = 0 is often denoted as the (uniform) "conformal" brane of the theory (2.1). In contrast, if we turn-on the source the dual gauge theory is no longer conformal. In particular, there are such solutions (2.8) with q 1,2,3,4 (Z, x) = q 1,2,3,4 (Z) (and q 4 (0) = Λ) that are translationally invariant along x, x 2 , x 3 . These are often denoted as the "uniform nonconformal" branes of the theory. This is one family of solutions that we will construct in this manuscript. Still with Λ = 0, we can then have solutions that break translational invariance along x (while keeping the isometries along the other two planar directions). Our main aim is to construct these nonuniform solutions, which we denote as "lumpy nonconformal branes", and study their thermal competition with the uniform nonconformal branes in a phase diagram of static solutions of (2.1), both in the micro-canonical and canonical emsembles. Of course there are also nonconformal branes that break translation invariance along the other two directions x 2,3 . These are cohomogeneity-4 solutions that we will not attempt to construct. Fortunately, the cohomogeneity-2 lumpy branes that we will find seem to already allow us to understand the key properties of the most general system.

Setup of the boundary-value problem
Finding the nonconformal brane solutions necessarily requires resorting to numerical methods. For that, we find it convenient to change our radial coordinate into y = Z 2 Z + , and define y + ≡ 1 Z + , (2.17) so that the ansatz (2.8) now reads with compact coordinates y ∈ [0, 1] and x ∈ [0, 1]. The horizon is located at y = 1 and the asymptotic boundary at y = 0. Note, that we used the two scaling symmetries (2.11)-(2.12) to set ≡ 1 and y + ≡ 1.
In these conditions we now need to find the minimal set of Einstein-scalar equationsthe equations of motion (EoM) -that allows us to solve for all q j (x, y) while closing the full system of equations, g µν = T ab and 2φ = 0, of the action (2.1). This is a nested structure of PDEs. This structure motivates also in part our original choice of gauge in the ansätze (2.8) and (2.18). For this reason, rather than presenting the final EoM, it is instructive to explain their origin and nature. Prior to any gauge choice and symmetry assumptions, the differential equations that solve (2.1) are second order for all the fields. The symmetry requirements we made fix some of these fields. Additionally, the fact that we have chosen to fix the gauge freedom of the system using the 'double Wick rotation Schwarzschild' gauge means that our system of equations should have the structure of an ADM-like system but with the spacelike coordinate x playing the role of "time". Of course, our problem is ultimately an elliptic problem. However, it is instructive to analyse our EoM adopting the above "time-dependent" viewpoint. In doing so, one expects that the equations of motion take a nested structure of PDEs that include a subset of "evolution" equations (to be understood as evolution in the x-direction) but also a subset of non-dynamical ("slicing" and "constraint") equations. We now describe in detail this nested structure.
We have a total of five equations of motion. Two of these EoM are dynamical evolution equations for q 1 and q 4 . The main building block of these equations is the Laplacian operator ∂ 2 x + ∂ 2 y (acting either on q 1 or q 4 ) that describes the spatial dynamics as the system evolves in x. With respect to the familiar ADM time evolution in the Schwarzschild gauge, this Laplacian replaces the wave operator ∂ 2 t − ∂ 2 y . Additionally, we have two (firstorder) slicing EoM for ∂ y q 2 and ∂ y q 3 that can be solved at each constant-x spatial slice for q 2 and q 3 . Besides depending on ∂ y q 2,3 and q 2,3 -but, quite importantly, not on ∂ x q 2,3 -these equations also depend on q 1,4 (that are determined "previously" by the evolution equations) and on their first derivatives (both along x and y). This means that we can integrate the slicing equations of motion along the radial direction to find q 2,3 at a particular constant-x slice. Finally, the EoM still includes a fifth PDE that expresses ∂ x q 2 as a function of (q 1,2,3,4 , ∂ x q 1,4 , ∂ y q 1,4 ). Let us schematically denote it as C(x, y) = 0. This is a constraint equation. To see this note that, after using the evolution and slicing EoM and their derivatives, the Bianchi identity ∇ µ (R µν − g µν R/2) = 0 implies the constraint evolution relation 3 where F (x, y) is a function that is regular at the horizon whose further details are not relevant. 4 It follows from this constraint evolution relation that if the constraint equation is obeyed at a given y, say at the horizon C(x, y) y=1 = 0, then it is obeyed at any other y ∈ [0, 1]. In practice, this means that we just need to impose the constraint as a boundary condition at y = 1, say. It is then preserved into the rest of the domain. Altogether, the strategy to solve the EoM is thus the following. There are effectively four EoM, two second order PDEs for q 1,4 (x, y) and two second order PDEs for q 2,3 (x, y). We need to solve these equations for q 1,2,3,4 (x, y) as a boundary-value problem. One of the 3 Note that in standard ADM time-evolution problems the constraint relation that must vanish involves the time derivative, i.e. it is schematically of the form ∂t √ −gC +F (t, y) √ −gC = 0, where y is a radial coordinate. Interestingly, in our 'double Wick rotation of the ADM gauge' -where x is the evolution coordinate -it is the radial derivative ∂y (and not ∂x) that appears in the vanishing constraint relation. 4 LetĈ ≡ √ −g C. Solving ∂yĈ + F (x, y)Ĉ = 0 yieldsĈ(x, y) =Ĉ(x, 1) exp y 1 F (x, Y )dY which converges if F (x, y) is regular at the horizon y = 1. boundary conditions is imposed at the horizon and takes the form C(x, y) y=1 = 0, whereas the others are the physically motivated boundary conditions discussed next.
Our integration domain is a square bounded in the radial direction by y = 0 (the asymptotic boundary) and y = 1 (the horizon). Along the x-direction the boundaries are at x = 0 and x = 1. The angular coordinate x is periodic in the interval [0, 1]. We can thus use this symmetry to impose Neumann boundary conditions for all q j at x = 0 and x = 1: Consider now the asymptotic UV boundary at y = 0. The invariance of the EoM under dilatations (2.12) guarantees that asymptotically our PDEs are of the Euler type and thus y = 0 is a regular singular point. The order of our PDE system (two second-order and two first-order PDEs) is 6. It follows that we have a total of 6 free UV independent parameters. We have explicitly checked that, for our gauge choice and scalar potential, all our functions q j admit a Taylor expansion in integer powers of y (in particular, without logarithmic terms) that contains precisely 6 independent parameters. In this Taylor expansion, at a certain order (as expected due to the fact that we fixed the gauge and our system is cohomogeneity-2) we need to use a differential relation that is ultimately enforced by the Bianchi identity: The requirement that our nonconformal branes asymptote to AdS fixes two of the six UV integration constants to unity, namely q 1 (x, 0) and q 3 (x, 0) at the boundary (it then follows directly from the EoM that q 2 (x, 0) = 1). We complement these boundary conditions with a Dirichlet boundary condition for the scalar field function q 4 which introduces the source Λ. This will be a running parameter in our search for solutions. Altogether we thus impose the boundary conditions at the UV boundary: Finally, we discuss the boundary conditions imposed at the horizon, y = 1. It follows directly from the four equations of motion that q 1,2,3,4 (x, 1) are free independent parameters and q 1,2,3,4 must obey a set of four mixed conditions that fix their first radial derivative as a function of q 1,2,3,4 (x, 1) and ∂ x q 1,2,3,4 (x, 1) that is not enlightening to display. When assuming a power-law Taylor expansion for we are already imposing boundary conditions that discard two integration constants that would describe contributions that diverge at the horizon.
Having imposed the boundary conditions, we must certify that we have a well-defined elliptic (boundary-value) problem. For that, we need to confirm that the number of free parameters at the UV boundary matches the IR number of free parameters. Recall that, before imposing boundary conditions, we have 6 integration constants at the UV and another 6 in the IR. A power-law Taylor expansion about the asymptotic boundary, concludes that, after imposing the boundary conditions (2.23) that fix a 1(0) = 1, a 3(0) = 1 and a 4(0) = Λ (a 2(0) is not a free parameter since it is fixed by the equations of motion), we are left with three free UV parameters, namely a 1(2) (x), a 2(2) (x), a 4(1) (x). Note that we will give Λ as an input parameter so the boundary-value problem will not have to determine it. As we shall find later, the energy density depends on these three parameters, whereas the expectation value of the dual operator sourced by Λ is proportional to a 4(1) (x). On the other hand, at the horizon, after imposing the aforementioned boundary conditions that eliminate two integration constants, one finds that there are 4 free IR parameters. So we have 3 UV free parameters but 4 IR free parameters. In order to have a welldefined boundary value problem (BVP) the number of free UV parameters must match the IR number. Note however that in our discussion of the boundary conditions we have not yet imposed the constraint equation C(x, y) = 0. As described above we just need to impose it at the horizon y = 1, where C(x, 1) = 0 simply reads This is solved by q 1 (x, 1) = √ α Λ / q 2 (x, 1) where α Λ is a constant to be fixed below. It follows that, if we impose this Dirichlet condition together with the three aforementioned mixed conditions for q 2,3,4 at the horizon (2.27) then we have just three free IR parameters, namely q 2,3,4 (x, 1). We fix the value of α Λ as follows. For a given source Λ, the lumpy nonconformal branes that we seek merge with the uniform nonconformal branes at the onset of the Gregory-Laflamme-type instability of the latter. We use this merger (where the fields of the two solutions must match) to fix the constant α Λ to be the value of q 1 (y) 2 q 2 (y) y=1 (no x-dependence) of the nonuniform solution when it merges with the uniform brane.
This discussion can be complemented as follows (which also allows us to set the IR boundary condition we impose to search for the uniform branes). Note that when looking for uniform branes, the PDE system of EoM reduces to an ODE system without any xdependence. In particular, the constraint equation C(x, 1) = 0 reduces to C(1) = 0 and is trivially obeyed, since all of its terms involve terms with partial derivatives in x. So when searching for uniform branes we use the IR boundary conditions (2.27) but with the first condition replaced by the Dirichlet condition q 1 (y) y=1 = 1: (if uniform branes). (2.28) This is a choice of normalization that does not change physical thermodynamic quantities. Solving the EoM for a given source Λ we then find, in particular, the value of q 2 (y) at y = 1. We can then read the constant α Λ = q 1 (y) 2 q 2 (y) y=1 = q 2 (y) y=1 of the uniform brane with source Λ. This value of α Λ (Λ) is then the one we plug in the boundary condition (2.27) to find the non-uniform branes with the same source Λ. The UV boundary conditions for the uniform branes is still given by (2.23) (with the replacement q j (x, y)| y=0 → q j (y)| y=0 ).

Thermodynamic quantities
Having found the nonconformal solutions (2.18) that obey the boundary conditions (2.20)-(2.27), we will now implement the holographic renormalization procedure in order to obtain the relevant thermodynamic quantities (recall that ≡ 1 and y + ≡ 1). Our (non)uniform branes are asymptotically AdS 5 solutions with a scalar field with mass µ 2 = −3. In these conditions, the holographic renormalization procedure to find the holographic stress tensor T ab and expectation value O φ of the operator dual to the scalar field φ was developed in [14,15]. We apply it to our system. We first need to introduce the Fefferman-Graham (FG) coordinates (z, χ) that are such that the asymptotic boundary is at z = 0 and g zz = 1/z 2 and g za = 0 (with a = t, χ, w 2,3 ) at all orders in a Taylor expansion about z = 0. In terms of the radial and planar coordinates (y, x) of (2.18) the FG coordinates are The expansion of the gravitational and scalar fields around the boundary up to the order that contributes to the thermodynamic quantities is then The holographic quantities can now be computed using the holographic renormalization procedure of Bianchi-Freedman-Skenderis [14,15]. 5 At the end of the day, for our system, the expectation value of the holographic stress tensor is given by where the metric components g (0) ab , g (4) ab and the scalar decay φ 2 can be read directly from (2.30)-(2.33), and we recall that φ M is a parameter of the superpotential (2.2) that we will eventually set to φ M = 1. Similarly, the expectation value of the dual operator sourced by Λ is The trace of the expectation value yields the expected Ward identity associated to the conformal anomaly which reflects the fact that our branes are not conformal. 6 Furthermore, after using the Bianchi relation (2.22), we confirm that the expectation value of the holographic stress tensor is conserved, i.e.
In (2.34) and (2.35) we have reinstated the appropriate power of in order to remind the reader that, for an SU (N ) gauge theory, the prefactor in these expressions typically scales as In the rest of the paper we will work with rescaled quantities obtained by multiplying the stress tensor and the scalar operator by the inverse of this factor. Evaluating (2.34) explicitly we find the following expressions for the energy density E, the longitudinal pressure P L (along the inhomogenous direction x) and the transverse pressure P T (along the homogenous directions x 2 and x 3 ): 5 Note however that we use different conventions for the Riemann curvature, that is to say, with respect to [14,15] our action (2.1) has the opposite relative sign between the Ricci scalar R and the scalar field kinetic term (∇φ) 2 . 6 Note that the holographic gravitational conformal anomaly contribution Agrav [14,15] and the scalar conformal anomaly contribution A scalar vanish for our system.
Note that P L is the pressure conjugate to the dimensionful coordinatex, not to the dimensionless coordinate x. Moreover, for static configurations, conservation of the stress tensor implies that P L is constant along the inhomogeneous direction, i.e. independent of χ. The temperature T and the entropy density s of the nonconformal branes can be read simply from the surface gravity and the horizon area density of the solutions (2.18), respectively: where have already used the boundary condition (2.27) that introduces the constant α Λ (that we read from the uniform solutions; see discussion below (2.27)). The Helmoltz free energy density is F = E − T s. The total energy E, entropy S and free energy F are obtained by integrating over the total volume: where we have made use of the fact that the system is homogeneous in the transverse directions. It will also be useful to define average densities by dividing the integrated quantities by the total volume: For uniform branes these averages coincide with the corresponding densities, since the latter are constant, but for nonuniform branes they do not. A quantity that will play a role below is an analogous integral for the expectation value of the scalar operator: Because of the translational invariance in the transverse directions it will also be convenient to work with densities in the transverse plane, namely with quantities that are only integrated along the inhomogeneous direction. Thus we define the energy, the entropy, the free energy and the expectation value densities per unit area in the transverse plane as We will refer to these type of quantities as "area densities" or "Killing densities". In order to write the first law we will also need the integral of the transverse pressure along the inhomogeneous direction. We therefore define Note that p L and p T have mass dimension 4 and 3, respectively. Finally, we will choose to measure all dimensionful quantities in units of the only gauge theory microscopic scale Λ.
We will use a "ˆ" symbol to denote the corresponding dimensionelss quantity obtained by multiplying or dividing a dimensionful quantity by the appropriate power of Λ, thus: (2.48) We are now ready to write down the first law. In order to do this, we first note that the extensive thermodynamic variables of the system are the total energy E, the total entropy S, the scalar source Λ, and the three lengths L, L 2 ≡ L, L 3 ≡ L of the planar directions. It follows that the first law for the total charges of the system is: We see that T , O, p L L 2 and p T L are the potentials (intensive variables) conjugate to S, Λ, L and L, respectively. Since the system is translationally invariant along the x 2 and x 3 directions, under the associated scale transformation x 2,3 → λ 0 x 2,3 the energy transforms as E(x 2,3 ) → λ 2 0 E(x 2,3 ) and thus it is a homogeneous function of λ 0 of degree 2. This means that for any value of λ 0 one has: We can now apply Euler's theorem for homogeneous functions to write the energy as a function of its partial derivatives: 7 The relevant partial derivatives in (2.51) can be read from (2.49) and, recalling that we are taking L 2 = L 3 ≡ L, this yields the Smarr relation for the total charges of the system Dividing by L 2 we obtain a Smarr relation for the area densities along the transverse plane: Rewriting the first law (2.49) in terms of these densities and using (2.53) we find the first law for the area densities: Since we will measure all dimensionful quantities in units of Λ, it will be useful to find a first law and a Smarr relation for the dimensionless densitiesρ,σ, etc. In order to do this we first use the dilatation transformation (2.12). Under this scale transformation x µ → λ 2 x µ the energy density ρ transforms as ρ(x µ ) → λ 3 2 ρ(x µ ) and thus it is a homogeneous function of λ 2 of degree 3, i.e. ρ λ 2 2 σ, λ 2 Λ, L/λ 2 = λ 3 2 ρ (σ, Λ, L) .  Reading the associated derivatives from the first law (2.54) we find (another) Smarr relation for the dimensionful densities: that nonconformal branes with two Killing planar directions x 2,3 must obey. In a traditional thermodynamic language the first law (2.59) and the Smarr relation (2.58) are also known as the Gibbs-Duhem and Euler relations, respectively. In our case they provide valuable tests of our numerical results. Moreover, they will be useful to discuss the dominant thermal phases in the microcanonical and canonical ensembles. Indeed, in the microcanonical ensemble the dominant phase will be the one that maximisesσ for fixed values ofρ andL. Similarly, in the canonical ensemble the dominant phase will be the one that minimisesf for fixedT andL.

Perturbative construction of lumpy branes
In the previous sections we have setup the BVP that will allow us to find the uniform and nonuniform nonconformal branes (2.18) that obey the boundary conditions (2.20)-(2.27). This nonlinear BVP can be solved in full generality using numerical methods. In the uniform case we have a system of coupled quasilinear ODEs that can be solved without much effort. However, in the nonuniform case the ODEs are replaced by PDEs and it is harder to solve the system. We will do this numerically in Sec. 2.5. In the present section we will complement this full numerical analysis with a perturbative nonlinear analysis that finds lumpy branes in the region of the phase diagram where they merge with the uniform branes. This perturbative analysis will already provide valuable physical properties of the system. Additionally, these perturbative results will also be important to test the numerical results of Sec. 2.5. We solve the BVP in perturbation theory up to an order in the expansion parameter where we can distinguish the thermodynamics of the uniform and nonuniform branes. We follow a perturbative approach that was developed in [16] (to find vacuum lattice branes) and that has its roots in [17][18][19] (to explore the existence of vacuum nonuniform black strings). More concretely, our strategy to find perturbatively the lumpy branes has three main steps: 1. The first step is to construct the uniform branes.
2. Then, at linear (n = 1) order in perturbation theory, we find the locus in the space of uniform branes where a zero-mode, namely a mode that is marginally stable, exists.
We will refer to this mode as the GL-mode. In practice, we will identify this locus by finding the critical length L = L GL (wavenumber k GL = 2π/L GL ) above (below) which uniform branes become locally unstable (stable). As expected from the discussion in Sec. 1, this critical length only exists for energy densities between points A and B.
3. The third step is to extend perturbation theory to higher orders, n ≥ 2, and construct the nonuniform (lumpy) branes that bifurcate (in a phase diagram of solutions) from the GL merger curve of uniform branes.
We describe in detail and complete these three steps in the next three subsections.

Uniform branes: O(0) solution
The first step is to construct the uniform branes. We solve the system of four coupled ODEs for q j (x, y) ≡ Q j (y) (here and below, j = 1, . . . , 4) subject to the boundary conditions (2.23) and (2.28), as described in Sec. 2.2. This can be done only numerically: we use the numerical methods detailed in the review [20].
There is a 1-parameter family of uniform nonconformal branes. We can take this parameter to be the scalar field source Λ. This is actually how we construct these solutions since Λ is an injective parameter: we give the source Λ via the boundary condition (2.23) and find the associated brane; then we repeat this for many other values of Λ. The dimensionless energy densityÊ = E/Λ 4 decreases monotonically as Λ grows, so this procedure maps out all possible uniform branes. Recall that, once we have found q j (x, y) ≡ Q j (y), the thermodynamic quantities of the solution follow straightforwardly from Sec. 2.2.
The properties of uniform branes are summarized in Figs. 3 and 4. In the left panel of Figs. 3 we plot the dimensionless energy densityÊ ≡ E/Λ 4 as a function of the dimensionless temperatureT ≡ T /Λ. We see the familiar S-shape associated to the multivaluedness of a first-order phase transition. Specifically, for a given temperatureT in the window of temperaturesT A ≤ T /Λ ≤T B there are three distinct families or branches of uniform branes with different values ofÊ. We will refer to these families as the "heavy", "intermediate" and "light" branches. The heavy branch (with higher energy density) starts in the conformal T /Λ → ∞ limit and then extends through point D all the way down to point A as the temperatureT decreases. The intermediate branch extends from point A, passes though point C, towards point B. This branch has negative specific heat and is both thermodynamically and dynamically locally unstable. A general discussion of these features can be found in Sec. 2 of Ref. [2]. In the present paper we will analyse the zero-mode properties of this instability in Sec. 2.4.2, and its timescale in Sec. 2.7. Finally, the light branch (with lower energy density) starts at point B, passes thought point E and extends all the way down towards T /Λ → 0. We do not show the plots ofŝ(T ) andp L,T (T ) because they are qualitatively similar to the plot ofÊ(T ). The relevant phase diagram for the canonical ensemble, namely the dimensionless free energyF ≡ F/Λ 4 as a function of the dimensionless temperatureT , is displayed in the right In the microcanonical ensemble, the relevant phase diagram is the average entropy densityŝ ≡ s/Λ 3 as a function of the average energy densityÊ ≡ E/Λ 4 . It is important to consider averaged quantities (which involve integration along the x direction) because inhomogeneous state will play a role. The qualitative form of the functionŝ(Ê) is shown in Fig. 5. The key features are as follows.ŝ is convex (ŝ > 0) in the region between A and B. This indicates local thermodynamical instability, since the system can increase its total entropy by rising the energy slightly in part of its volume and lowering in another so as to keep the total energy fixed. In the regions EB and AD the entropy function is concave (ŝ < 0) but there are states with the same total energy and higher total entropy, namely phase-separated configurations in which the phases E and D coexist at the critical temperature. These states are characterised by the fractions 0 ≤ ν, (1 − ν) ≤ 1 of the total volume occupied by each phase, so their total entropy is of the formŝ E + (ŝ D −ŝ E )ν, as indicated by the red segment in Fig. 5. Therefore the regions EB and AD are locally but not globally thermodynamically stable. Finally, all states outside the region ED are globally stable. For our system, these qualitative features are difficult to appreciate directly on a plot ofŝ versusÊ because the curveŝ(Ê) is very close to a straight line. For this reason we show the convexity/concavity property (the second derivative) in Fig. 6(left) and the difference between the phase-separated configurations and the homogeneous solutions in Fig. 6(right).

Gregory-Laflamme physics: O(1) solution and the spinodal zero-mode
The intermediate uniform branes withÊ B <Ê <Ê A (see left panel of Fig. 3), and only these, can be Gregory-Laflamme (GL) unstable. Roughly speaking, we expect this to happen if their dimensionless length L Λ (along the x direction) is bigger than the dimensionless thermal scale Λ/T of the system. This linear instability is ultimately responsible for the nonlinear existence of the lumpy solutions. Therefore, our second step is to consider static perturbations about the uniform branes, q j (x, y) = Q j (y) + q (1) j (x, y), that break the U (1) symmetry along x (see Sec. 2.7 for time-dependent perturbations). Here, 1 is the amplitude of the linear perturbation and, ultimately, it will be the expansion parameter of our perturbation theory to higher order.
We adopt a perturbation scheme that is consistent with our nonlinear ansatz (2.18)where we recall that x ∈ [0, 1] -since we want to simply linearize the nonlinear equations of motion that we already have (Sec. 2.2) to get the perturbative EoM. In this perturbation scheme we assume an ansatz for the perturbation of the form 8 (2.60) This means that the length L of the periodic coordinate x is given in terms of the wavenumber k of the perturbation by L = 2π/k, and it will change as we climb the perturbation ladder (this is because k, and thus L, will be corrected at each order; see Sec. 2.4.3). Since the EoM depend on L, this relation L = 2π/k introduces the zero mode wavenumber k in the problem. 9 Under these circumstances the linearized EoM become a simple eigenvalue problem in k 2 of four coupled ODEs. Henceforth, we denote this leading-order wavenumber by k GL . So we need to solve our eigenvalue problem to find the eigenvalue k GL as well as the associated four eigenfunctions q (1) j (y). Note however that we "just" need to solve an ODE system of four coupled equations (not PDEs) subject to the linearized versions of the boundary conditions (2.23)-(2.27). For example, when we linearize (2.23) using j | y=0 we find that the linear perturbations q 2 y=1 and mixed boundary conditions for q 2,3,4 y=1 . Of course, in this linearization procedure about the uniform brane, we insert the boundary conditions (2.23) and (2.28) of the leading solution; in particular, we impose Q 1 y=1 = 1 and Q 2 y=1 = α Λ .
Summarizing this second step, the above perturbation procedure at O( ) finds the critical zero mode of the Gregory-Laflamme (GL) instability of uniform branes with energy densitiesÊ B <Ê <Ê A . That is to say, it finds the dimensionless critical wavenumber k GL = k GL /Λ for the onset of the GL instability, and thus the minimum length L GL Λ = 2π/k GL above which the uniform brane is unstable. This critical valuek GL =k GL (T ) is only a function of the dimensionless temperatureT = T /Λ and is plotted in Fig. 7. We see that k GL = 0 at the endpoints A and B of the intermediate uniform branch whereT =T A andT =T B . These two branes are effectively stable sinceL GL → ∞ at these two temperatures. However, intermediate branes withT SoL GL is parametrized byT , and the energy density of uniform branes is also only a function of the temperature,Ê =Ê(T ). It follows that we can identify the onset GL curve 8 The superscript (n) here and henceforth always denotes the order n of the perturbation theory, not order of derivatives. 9 We have some freedom in the choice of the perturbation scheme. For example, an alternative perturbation scheme would be to keep the length L fixed by absorbing the L factor in the metric component gxx of the ansatz (2.18) into a new coordinatex. That is to say, we would change the x ∈ [0, 1] coordinate of (2 In this case, the U (1) dependence of the perturbation would be cos(kx) which would introduce the wavenumber k = 2π L in the problem. These two schemes are equivalent. This follows from the observation that the two sets of Fourier modes are equivalent: cos(η kx) = cos η 2π L L 2 x = cos(η π x). Further recall from the discussion above (2.8) that our solutions have Z2 symmetry: the solution inx ∈ [−L/2, 0[ can be obtained by simply flipping our solution over the x = 0 axis (computationally this is useful/efficient since we deploy a given number of grid points to study the range [0, L/2] instead of [−L/2, L/2]). This is why we have just a factor of π and not 2π in the arguments of our Fourier cosines. of uniform branes in a plotÊ vsL. This is done in Fig. 8. This plot is effectively a stability phase diagram for the uniform branes since the black dotted GL onset curve separates the region where the uniform branes are unstable -namely, the parabola-like shaped interior regionÊ B <Ê <Ê A withL >L GL -from its complementary region where branes are stable against the spinodal instability. In this figure note that the energy densityÊ =Ê A andÊ =Ê B corresponds to the energy densities of the uniform solutions A and B in Fig. 3 and note thatL GL → ∞ when the energy density of the black dashed GL onset curve approachesÊ A orÊ B . To summarize, Fig. 8 shows that intermediate uniform branes with a given energy densityÊ B <Ê <Ê A are unstable if their dimensionless length is higher that the GL critical length,L >L GL . Not less importantly, in a phase diagram of solutions, the GL onset curve also signals a bifurcation to a new family of solutions that describes nonuniform or lumpy branes. That is to say, the GL onset curve is a merger line between the uniform and lumpy nonconformal branes. Perturbation theory at order O( ) identifies this merger or intersection line (see Fig. 8) of two distinct surfaces in a 3D phase diagram but it cannot describe the properties of the lumpy brane surface as we move away from the merger line (roughly speaking, it cannot describe the "slope of the lumpy surface" in a 3D phase diagram). For that, we need to proceed to higher order O( n ) in the perturbation theory, as we do in in the next subsection. Figure 8. Stability diagram for uniform nonconformal branes. The interpretation of the two yellow square points Σ 1 and Σ 2 will be given when discussing

Lumpy branes: perturbative solution at O(n)
To find the solution at order O( n ) we expand the metric functions and wavenumber in powers of : In this expansion we have made the identification k (0) ≡ k GL and we have already found the n = 1 contribution in the previous section. Recall that this {k GL , q j } contribution was found by solving a homogeneous eigenvalue problem for k GL . The expansion (2.61) is such that at order O( n ) we solve the BVP to find the coefficients {k (n−1) , q (n) j }. Further note that, as explained above, our choice of perturbation scheme is such that the length L is corrected at each order n (see also footnote 9). That is, one has where the coefficients L (n−1) can be read straightforwardly from (2.61b). This also means that in our choice of scheme, the periodicity of the x circle allows us to introduce a separation ansatz for the perturbation coefficients q (n) j (x, y) whereby they are expressed as a sum of Fourier modes (with harmonic number η) in the x direction as (2.63) So here and onwards, 0 ≤ η ≤ n identifies a particular Fourier mode (harmonic) of our expansion at order O( n ). At order n ≥ 2, the perturbation EoM are no longer homogeneous. Instead, they describe an inhomogeneous boundary value problem with a source S (n,η) . Not surprisingly, this source is a function of the lower order solutions . This source can always be written as a sum of Fourier modes of the system. We find that at order n ≥ 2, the maximum Fourier mode harmonic that is excited in the source is η = n. This is due to the fact that at linear order we start with the single η = 1 Fourier mode and the n th polynomial power of this linear mode, after using trigonometric identities to eliminate powers of trigonometric functions, can be written as a sum of Fourier modes with the highest harmonic being η = n. This property of our source implies that the solution of the O(n) EoM can only excite harmonics up to η = n and this explains why we capped the sum in (2.63) at η = n.
To proceed, at each order O(n), we have to distinguish the Fourier modes η = 1 from the other, η = 1. This is because this particular Fourier mode η = 1 is the only one that is already excited at linear order n = 1.
Start with the generic case η = 1. Then the differential operator -call it L H -that describes the associated homogeneous system of equations, L H q (n,η) j = 0, is the same at each order n and for any Fourier mode η: it only depends on the uniform brane Q j (y) we expand about and k GL . The ODE system of 4 inhomogeneous equations is thus of the form It follows that the complementary functions of the homogeneous system are the same at each order n ≥ 2 and η. But, we also need to find the particular integral of the inhomogeneous system and this is different for each pair (n, η) since the sources S (n,η) differ. The general solution q , all of which follow from (2.27). Consider now the exceptional case η = 1. In this case, at order n ≥ 2, our BVP becomes a (non-conventional 10 ) eigenvalue problem in k (n−1) . That is to say, the ODE system of 4 inhomogeneous equations is now of the form where K jm is a diagonal matrix whose only non-vanishing components are K 11 = 1 = K 44 . Recall that L H is an operator that describes two second-order ODEs for q 1 , q 4 and two firstorder ODEs for q 2 , q 3 and this justifies the presence of this particular K jm in our eigenvalue term. We now have to solve (2.65) (subject to boundary conditions that are motivated as in the η = 1 case) to find the eigenvalue k (n−1) and the eigenfunctions q (n,1) j (y). To have a full understanding of the EoM of our perturbation problem one last observation is required. As pointed out above, the highest Fourier harmonic that is excited in our system at order O( n ) is η = n. This is because the n th polynomial power of the single Fourier mode that is present at linear order, after using trigonometric identities to eliminate powers of trigonometric functions, can be written as a sum of Fourier modes with the highest harmonic being η = n. But this trigonometric operation also indicates (as we explicitly confirmed) that not all Fourier modes with η ≤ n are excited. More concretely, for even n ≥ 2 we find that only even 0 ≤ η ≤ n modes are present in our system. And for any odd n ≥ 3, only odd 0 ≤ η ≤ n modes are excited. Therefore, up to order n = 5 we find that the modes that are excited in our system are: j (x, y) = q (5,1) j (y) cos(2 π x) + q This last property of our system, together with the previous observation -see the discussion of (2.65) -that Fourier modes with η = 1 are those that give the wavenumber correction k (n−1) at order n, immediately allows us to conclude that k (n−1) = 0 if n is even. At even n order the cos (π x) Fourier mode is not excited by the source and thus the only solution of (2.65) is the trivial solution.
Finally, note that the η = 0 harmonics are of particular special interest. Indeed note that modes with η = 0 do not contribute (since the integral of a cosine vanishes) to the total thermodynamic quantities of the solution such as the energy E, the entropy S, etc. It follows from the discussion of (2.66) that odd order n modes do not contribute to correct these thermodynamic quantities. 10 It is not a standard eigenvalue problem because the eigenvalue k (n−1) is not multiplying the unknown eigenfunction q (n,η) j . Instead, it multiplies an eigenfunction that was already determined at previous n = 1 order.
We can finally summarize the key aspects of the general flow of our perturbation theory as the order n grows: 1. even orders O( n ) introduce perturbative corrections to thermodynamic quantities like energy, entropy, pressure, etc., but they do not correct the wavenumber, k (n−1) = 0 (and thus do not correct L).
2. odd orders O( n ) give the wavenumber corrections k (n−1) but do not change the energy, entropy and pressure.
We complete this perturbation scheme up to order O( 5 ): this is the order required to find a deviation between the relevant thermodynamics of the lumpy branes and the uniform phase.
Once we have found all the Fourier coefficients q (n,η) j (y) and wavenumber corrections k (n−1) up to n = 5, we can reconstruct the four fields q j (x, y) using (2.61). We can then substitute these fields in the thermodynamic formulas of Sec. 2.2 to obtain all the thermodynamic quantities of the system up to O( 5 ). We find that all of them, as well as the wavenumber, have an even expansion in n , with the only exception of the temperature that is simply given by (2.42). Now that we have the thermodynamic description of lumpy branes up to O( 5 ), we can compare it against the thermodynamics of uniform branes and find which of these two families is the preferred phase. We are particularly interested in the microcanonical ensemble, so the dominant phase is the one that has the highestσ for a given pair (L,ρ). Let Q u and Q nu denote thermodynamic quantities Q for the uniform and nonuniform branes, respectively. When comparing these two solutions in the microcanonical ensemble, one must haveL nu =L u , andρ nu =ρ u . (2.67) Given a lumpy brane with (L nu ,ρ nu ) we must thus identify a uniform brane whose Killing densityρ u satisfies (2.67). Equivalently, we can impose that the energy density of the uniform brane obeys Both sides of this equation are known as a perturbative expansion in . This is because the energy density E u is a function of the dimensionless temperatureT u which is corrected at each order asT u =T 0 + 2T (2) + 4T (4) + O( 6 ) in our perturbation expansion. Similarly, the Killing energy densityρ nu ( ) and the lengthL nu ( ) of lumpy branes are also known as a Taylor expansion in . Therefore, in practice equation (2.68) becomes Taking the Taylor expansion of ρ u (T u ) we must impose Having thisT u we can now compute the entropy density of the uniform braneŝ u (T u ) and the Killing entropy densityσ u (T u ) =L uŝu (T u ). More concretely, a Taylor expansion in of this equality yieldŝ which allows us to find the entropy correction coefficientsσ (i) u and thus the Killing entropy densityσ u (T u ) up to order O( 6 ) of the uniform brane that has the same (L,ρ) as the particular lumpy brane we selected. This procedure (2.67)-(2.70) can now be repeated for all lumpy branes.
We are now ready to discuss our higher-order perturbative findings. First, in Fig. 9 we plot the wavenumber corrections k (2) (left panel) and k (4) (right panel), as defined in (2.61b). The fact that these higher order quantities grow large as one approachesT A and T B tells us that our perturbation theory breaks down in these regions. We will come back to this below.
Second, in order to determine the dominant phase, we are interested in the entropy difference between a nonuniform and a uniform brane when the two have the same length L and Killing energy densityρ. This is given by By construction ∆σ (0) ≡ 0 since the leading order of our perturbation theory describes the merger line of lumpy branes with uniform branes. Moreover, the first law for the Killing densities (2.59) can be rewritten, in the perturbative context, as ∂ ρ =T ∂ σ +p L ∂ L and has itself an expansion in that must be obeyed at each order. The leading-order term of this expansion implies that ∆σ (2) ≡ 0, a condition that we actually use to test our numerical results. Therefore the first non-trivial contribution to ∆σ( ) occurs at fourth order, namely This is the reason why we have to extend our perturbation analysis up to O( 5 ). We conclude that, for given (L,ρ), if ∆σ (4) > 0 then the lumpy branes are the preferred phase; otherwise the uniform branes are the dominant phase. We should thus plot the coefficient ∆σ (4) of (2.72) as a function ofL andρ. However, we find it clearer to plot instead ∆σ (4) as a function of the temperatureT u of the uniform brane that has the same (L,ρ) as the lumpy brane we compare it with. This is done in Fig. 10. Recall that uniform branes can be GL-unstable only in the rangeT A ≤T ≤T B andÊ B <Ê <Ê A , see Fig. 3. It follows that lumpy branes bifurcate from the uniform branch at the GL zero mode for temperatures in the rangeT A ≤T ≤T B . Fig. 10 plots this range of temperature and shows that forT Σ 1 <T <T Σ 2 , where the values ofT Σ 1 andT Σ 2 are identified in the caption, the lumpy branes are the preferred thermodynamic phase since ∆σ (4) > 0. However, for T A <T <T Σ 1 andT Σ 2 <T <T B , we have ∆σ (4) < 0 and thus uniform branes dominate over the lumpy phase when they have the same dimensionless lengthL and Killing energy densityρ. Going back to Fig. 8, for completeness we have also identified these points Σ 1 and Σ 2 in the associated GL merger curve.
Figs. 9 and 10 also illustrate the regime of validity of our perturbative expansion. For example, in Fig. 10 we see that ∆σ (4) grows arbitrarily negative as we approach the endpoints A and B of the intermediate branes with temperatureT A 0.387944 andT B 0.405724 (see also Fig. 3). But once the associated entropy correction becomes of the order of our expansion parameter, ∆σ (4) 4 ∼ , perturbation theory breaks down. So we should not trust our perturbative results close to the endpoints A and B.
Even away fromT A andT B , our perturbation theory is certainly valid only for 1. Therefore we expect it to describe accurately the properties of lumpy branes close to their GL merger line with the uniform branes (where = 0) but not far away from this merger. To learn what happens further away, we need to solve the full nonlinear BVP using numerical methods. This is what we do in the next subsection.

Full nonlinear solutions and phase diagram of nonconformal branes
To find accurately the lumpy branes and thus their thermodynamics in the full phase space where they exist, one needs to resort to numerical methods to solve nonlinearly the associated BVP, which was set up in Sec. 2.2. It consists of a coupled set of four quasilinear PDEs -two second-order PDEs for q 1,4 (x, y) and two second-order PDEs for q 2,3 (x, y) -that allow us to find the brane solutions (2.18) that obey the boundary conditions (2.20)-(2.27).
We solve our BVP using a Newton-Raphson algorithm. For the numerical grid discretization we use a pseudospectral collocation with a Chebyshev-Lobatto grid and the Newton-Raphson linear equations are solved by LU decomposition. These methods are reviewed and explained in detail in the review [20] and used in e.g. [21][22][23][24][25]. As explained in Sects. 2.1 and 2.2 (see in particular footnote 1 and the associated discussion) our gauge was judiciously chosen to guarantee that our solutions have analytical polynomial expansions at all the boundaries of the integration domain. In these conditions the pseudospectral collocation guarantees that our numerical results have exponential convergence with the number of grid points. We further use the first law and the Smarr relations (2.59)-(2.58) to check our numerics. In the worst cases, our solutions satisfy these relations with an error that is smaller than 1%. As a final check of our full nonlinear numerical results, we compare them against the perturbative expansion results of Sec. 2.4.
As usual, to initiate the Newton-Raphson algorithm one needs an educated seed. We use the perturbative solutions of Sec. 2.4 as seeds for the lumpy branes near the GL merger line with the uniform branes. The uniform branes are a 1-parameter family of solutions parametrized by the dimensionless temperature T /Λ. In contrast, the lumpy branes are a 2-parameter family of solutions that we can take to be T /Λ and the dimensionless length LΛ. This means that we need to scan a 2-dimensional parameter space. Our strategy to do so follows two routes. In one of them we follow lines of constant-temperature lumpy branes as their length LΛ changes. The temperature T is given by (2.42) where the constant α Λ and Λ (to build the dimensionless ratio T /Λ) are read from the uniform solution at the GL merger. The minimum length of these branes is the GL lengthL GL computed in Sec. 2.4.2, and constant-temperature branes exist for arbitrarily large LΛ. In a second route, we generate curves of lumpy branes that have fixed dimensionless length LΛ. In this path the temperature T /Λ of the branes changes but at the GL merger with the uniform branes, see e.g. Fig. 7, we know both the temperatureT and the associated GL lengthL GL (T ). Altogether these two solution-generating procedures allow us to construct a grid of two "orthogonal-like" lines of solutions that span the phase space of lumpy branes. Further, recall that once we have the numerical solutions q j (x, y), the thermodynamic quantities of the lumpy branes are read straightforwardly from the expressions discussed in Sec. 2.3.
After these preliminaries we are ready to discuss our numerical nonlinear findings. A first important plot is shown in Fig. 11, where we show the dimensionless average energy densityÊ = E/Λ 4 as a function of the dimensionless temperatureT = T /Λ. Recall that for uniform branesÊ coincides with the dimensionless energy density,Ê, which is constant across the entire system. This plot contains again the uniform-brane spinodal curve (blue circles) already shown in Fig. 3 but this time we also show some representative examples of lumpy brane solutions (all other lines/curves).
As illustrated in Fig. 11, a first non-trivial conclusion of our study is that lumpy branes exist only in the temperature windowT A ≤T ≤T B whereT A 0.387944 and T B 0.405724. That is, they exist only in the temperature range where the GL-unstable, intermediate branch of uniform solutions (the curve ACB) exists. Of course, we should have anticipated that lumpy branes merge with the uniform branes of the intermediate branch and thus in the windowT A ≤T ≤T B . However, it was a logical possibility that, away from this merger, lumpy branes might exist also for temperatures outside the rangê T ∈ [T A ,T B ]. We have generated considerably more solutions than those shown in Fig. 11 in order to test this possibility and, as stated above, we have found that it is not realised.
To continue interpreting Fig. 11, it is convenient to discuss separately the regionŝ T A ≤T <T c andT c <T ≤T B , i.e. the regions to the left and to the right, respectively, of the vertical dashed line DCE. Recall that this auxiliary line identifies the critical temperatureT =T c at which the first-order phase transition for uniform branes  Figure  place (see right panel of Fig. 3).
So consider first lumpy branes that exist in the windowT A ≤T <T c : 1. For a given temperatureT in this range, lumpy branes exist with a dimensionless length that satisfiesL GL ≤L ≤ ∞. In particular, the vertical lines of orange circles of Fig. 11 are lumpy branes at constantT that haveL =L GL (T ) when they bifurcate from the intermediate uniform-brane branch AC. Then they extend for arbitrarily largeL. More precisely, forT <T c , constant-T lumpy branes extend upwards (i.e. towards higherÊ) asL grows. However, we find that for a given step increase inL, the increase inÊ gets smaller and smaller asL grows, i.e. (∂Ê/∂L) T is a monotonically decreasing function ofL. This is explicitly observed in the vertical lines that we display: away from the merger each two consecutive orange circles are separated by the same step inL but the step increase inÊ is significantly decreasing as we move upwards. Due to the large hierarchy of scales that develops it is difficult to construct lumpy branes withL → ∞. But the above behaviour strongly suggests that lumpy branes withT <T c are precisely bounded by the heavy uniform branch segment AD whenL → ∞, i.e. we conjecture that Let us now follow these constant-L curves as they flow into the second relevant region, namelyT >T c . Fig. 11 shows that, if this was not already happening for smallerT , all these curves have a drop in theirÊ as they approachT c from the left. For very largeL this drop is dramatic with an almost vertical slope (see e.g. the magenta, ⊗ curve). Therefore, as best illustrated in the inset plot of Fig. 11 that zooms in the region around point C, all constant-L curves pile up around point C in a way such that: 1. AsT →T − c (approaching from the left) all curves haveÊ >Ê C . In particular, this means that these curves do not intersect the uniform branch AC near C.
2. Once atT >T c , all constant-L curves that bifurcated from the uniform branes in the trench AC cross the uniform brane branch curve between C and K. Recall that K describes the uniform brane solution that has the largest GL wavenumber k GL or, equivalently, that has the lowestL GL = 2π/k GL ; see Fig. 7. After this crossing, the constant-L lumpy branes keep extending to higherT with an energy density lower that the intermediate uniform brane with the sameT . This keeps happening until they merge again with the uniform brane in the trench KB at a critical temperature that is again the one predicted by the GL zero-mode analysis, i.e. at the highestT that satisfies the conditionL =L GL (T ), see again Fig. 7. Lumpy-brane curves with higher constantL merge with the uniform branch KB at a point that is closer to B. In the limit whereL → ∞ this merger occurs precisely at point B in Fig. 11, in agreement with the GL linear results of Fig. 7.
3. There are constant-L lumpy branes with very smallL that bifurcate from the uniform brane branch only in the trench CK (instead of AC). Then they extend to higher T , initially withÊ higher that the uniform branes with sameT before they cross the uniform branch CK at a temperatureT <T K and proceed to higherT below point K until they merge again with the uniform brane branch but this time in the trench KB (at a point very close to K). This happens for fixed-L branes whenever L GL (T K ) <L <L GL (T C ).
The three features of the lumpy branes just listed are compatible with the following interpretation that merges our nonlinear findings, summarized in Fig. 11, with the GL linear results of Sec. 2.4.2, summarized in Fig. 7. Indeed, let us go back to Fig. 7 and consider an auxiliary horizontal line at constantk GL , i.e. at constantL GL . This line intersects the curvek GL (T ) at two points. These are the two merger points of constantL lumpy branes with the uniform brane that we identify in Fig. 11. One of the mergers -let us denote it simply as the "left" merger -hasT A ≤T ≤T K and the other -the "right" merger -haŝ T K ≤T ≤T B . Since the maximum of the GL wavenumber occurs at a temperature that is higher than the one of the first-order phase transition of the uniform system,T K >T c , it follows that the "left" mergers of lumpy branes with constantL GL (T K ) <L <L GL (T C ) are in the trench CK of Fig. 11. But, forL >L GL (T C ), the "left" merger is located in the trench AC, with theL → ∞ "left" merger being at A. On the other hand, the "right" merger is always located in the trench KB of Fig. 11, with the "right" merger of theL → ∞ lumpy branes being at B. Our nonlinear results summarized in Fig. 11 further conclude that there are no lumpy branes withL <L GL (T K ). AsL approachesL GL (T K ) from above, lumpy branes exist only in a small neighbourhood around point K in Fig. 11, with the characteristics described in item 3 in the list above. 11 We stress again that lumpy branes exist only in the temperature rangeT A ≤T ≤T B . Our nonlinear results of Fig. 11 give strong evidence that constant-T <T c branes extend 11 Note that for other values of the (super)potential parameters φM and φQ in (2.2) (we have picked φM = 1 and φQ = 10), or in similar spinodal systems, it might well be the case thatTc >TK or, for a fine-tuned choice of potential, evenTc =TK . If that is the case our conclusions should still apply with the appropriate shift of K to the left of C in Figs. 7 and 11. Note that for this exercise we only need to find the uniform branes of the system and solve for static linear perturbations of these branes which determine the zero-mode GL wavenumber and thus the location of its maximum K with respect toTc. That is to say, we just need to complete the tasks described in sections 2.4.1 and 2.4.2.
to arbitrarily largeL with lim L→∞Ê constT →Ê AD u (T ) , (2.75) see the heavy uniform brane trench AD in Fig. 11. On the other hand, our results also strongly indicate that constant-T >T c branes extend to arbitrarily largeL with see the light uniform brane trench ED in Fig. 11. Moreover, asL grows arbitrarily large, the constant-L lumpy branes intersect (without merging with) the intermediate uniform brane branch AB at a point that is arbitrarily close (from the right) to point C in Fig. 11 and with a slope (∂Ê/∂T ) T c that grows unbounded, that is In the limitL → ∞ we thus conjecture that lumpy branes are limited by the curve ADCEB with two cusps connected by the vertical DCE line in Fig. 11. To argue further in favour of this conjecture, it is important to explore better the properties of the system in thiŝ L → ∞ limit and the associated limiting curve ADCEB. For that it is instructive to look at the energy density profileÊ(x) of the lumpy branes as a function of the inhomogeneous direction x.
In Fig. 12 we first consider lumpy branes with the sameT but different lengthsL. 12 In the left panel we have the profile of 3 lumpy branes withT 0.394948 <T c ; in the middle panel we have the profile of 3 lumpy branes withT 0.395894 T c (i.e. almost atT c 0.3958945); and, finally, in the right panel we show the profile of 3 lumpy branes withT 0.397308 >T c . In all panels, the blue diamond lines have a length only slightly aboveL GL (T ). Therefore, the profile of these lumpy branes is almost flat and very close to the horizontal dashed line that represents the intermediate uniform brane withÊ AC u (T ) (left/middle panels) orÊ CB u (T ) (right panel). Then, the green square curves have a length of roughlyL ∼ 1.25L GL (T ). We see that the profile starts becoming considerably deformed with one of the "halves" pulling well above (below) the uniform brane profile with the sameT . Finally, the red disk curves represent lumpy branes that have a lengthL(T ) that is considerably higher thanL GL (T ) (exact values in the caption). We see that the profile of lumpy branes withT 0.394948 <T c (left panel) is, in a wide range of x (x 0.7), very flat withÊ(x) ∼Ê AD u (T ), i.e. with an energy density that is the same as the one of the heavy uniform brane in the trench AD that has the sameT (upper horizontal dashed line). Then, for x 0.7,Ê(x) falls considerably towards the energy densityÊ light u (T ) of the light uniform brane that has the same temperature (lower horizontal dashed line). Still in Fig. 12, the middle panel shows that asT approachesT c and for largeL (red disks), the profileÊ(x) describes a domain-wall solution that interpolates betweenÊ AD u (T ) (for small 12 When interpreting these figures recall, from the discussion above (2.  TheL → ∞ limit of lumpy branes and its association with the limiting curve ADCEB is further revealed when we complement Fig. 12 with an analysis of the energy density profileÊ(x) of a constant-L family of branes for different values of the temperature. One such analysis is done in Fig. 13 where we fixL 11.501849: this picks the fourth constant-L curve (from bottom-left) in the plot of Fig. 11. For clarity we single out this curve and reproduce it -this time only the relevant zoomed in region of Fig. 11 -in the left panel of Fig. 13. We pinpoint a total of seven solutions with seven different temperatures (each one with its own distinctive plot marker shape and colour). The first ( ) and the last ( ) solutions are the two mergers with the intermediate uniform brane, the second ( ) and sixth ( ) solutions are the two extrema ofÊ(T ) L , and the third ( ), fourth (•) and fifth ( ) plot markers identify three solutions withT at or very close toT c . As in Fig. 12, we see that the profile of the two lumpy branes at the merger is flat:  12 The same shape/colour code is used. the uniform branes. As we move to the "extrema" solutions with plot markers and we see, like for similar solutions in Fig. 12, that the profile is considerably deformed. More important for our purposes are the solutions withT ∼T c , e.g. , •, . We see that for such cases the profile reaches its maximum deformation in the sense that the solution clearly interpolates between to regions that are fairly flat. Importantly, the small-x flat region is approaching the energy densityÊ D u (T c ) of the heavy uniform brane that hasT =T c (see the upper, horizontal, dashed, blue line labelled by D). Similarly, the large-x flat region is approaching the energy densityÊ E u (T c ) of the light uniform brane that hasT =T c (see the lower, horizontal, dashed, blue curve labelled by E). We further see that the closer we are toT − c (T + c ), the closer we get toÊ D u (T ) (Ê E u (T )). The plot of Fig. 13 is for a moderate value ofL. Combined with the findings of the discussion of Fig. 12 we conclude that asL grows large andT →T c , the flat regions get more extended in x and the domain wall that interpolates between them atÊ ∼Ê D u (T c ) andÊ ∼Ê E u (T c ) gets narrower. Altogether, the findings summarized in Fig. 12 and Fig. 13 lead to the following conclusion/conjecture. In the double limitL → ∞ andT →T c our results support the conjecture that That is, in this double limit we have a family of lumpy branes that fills up the segment DCE of Fig. 11. All this segment describes infinite-length lumpy branes that are sharp/narrow domain wall solutions interpolating (along 0 ≤x ≤ ∞) between two flat regions: one withÊ(x) =Ê D u (T c ) and the other withÊ(x) =Ê E u (T c ). These are the phase-separated configurations discussed above. As we move up from C to D, the region ofx withÊ(x) =Ê D u increases while as we move down from C to E, the region ofx withÊ(x) =Ê E u increases. We have infinite domain wall solutions that interpolate between the two uniform phases of the system at T = T c . Moreover, keeping the limitL → ∞, but relaxing the condition T →T c , the results summarized in Figs. 12 and 13 give evidence to conjecture that infinitelength lumpy branes exist only forT A ≤T ≤T B and are exactly at the line ADCEB of Fig. 11.
We will now discuss the thermal competition between lumpy and uniform nonconformal branes in the microcanonical ensemble. Recall that we keep the dimensionless lengthL and the Killing energy densityρ fixed and the relevant thermodynamic potential is the Killing entropy densityσ. Again, uniform and lumpy branes co-exist for temperatureŝ T A ≤T ≤T B . So in the microcanonical ensemble, given a lumpy brane with (L,ρ), our first task is to find the uniform brane (i.e. the temperatureT which parametrizes this family) that has the same (L,ρ) as the chosen lumpy brane. Once this is done, we can compare the Killing entropy densitiesσ of the two solutions at the same selected (L,ρ) pair. As for the perturbative analysis of Sec. 2.4.3 -see e.g. the discussion of (2.72) -we compute the entropy difference between the two phases when they have the sameL andρ: As before, the subscript "nu" stands for the nonuniform (lumpy) brane and "u" denotes the uniform brane. From our perturbative analysis recall that at the merger curve (ACB in Fig. 11 or the black dotted line in Fig. 8) between uniform and lumpy branes one must have ∆σ sameL,ρ = 0. Moreover, in the perturbative analysis leading to Fig. 10 we found that lumpy branes that bifurcate from uniform branes in the trench Σ 1 Σ 2 of Figs. 11 or 8 do so with a positive entropy difference slope. In other words, slightly away from the merger curve we have ∆σ sameL,ρ > 0. This means that lumpy branes emanating from the GL merger dominate over the uniform branes with the sameL,ρ, at least "initially". On the other hand, in the complement of Σ 1 Σ 2 , i.e. for lumpy branes bifurcating from AΣ 1 or Σ 2 B in Fig. 11, the perturbative analysis of Fig. 10 shows that ∆σ sameL,ρ < 0. That is, in this case for a given (L,ρ) lumpy branes have less Killing entropy density than the uniform solutions and thus the latter are the preferred phase in the microcanonical ensemble. Now that we have the full nonlinear solutions, our first task is to naturally compare these results with the perturbative results of Sec. 2.4.3 that led to Fig. 10. On the one hand this will check our numerical results. On the other hand it will identify the regime of validity of the perturbative analysis, i.e. how "far away" from the merger curve it holds. To illustrate this comparison, in Fig. 14 we show ∆σ sameL,ρ as a function ofρ for a family of lumpy branes that have a fixed temperatureT . Note that, asρ changes, so doesL in order to keepT fixed (this effect is better illustrated in Fig. 15, as we explain below). The plots in the top row of Fig. 14  T >T Σ 2 (right). In these cases the perturbative analysis summarized in Fig. 10 predicts that lumpy branes bifurcate from the GL merger with ∆σ sameL,ρ < 0, as indicated by the dashed red curves in Fig. 14. The numerical nonlinear results, shown as blue dots, indeed confirm this, and they are in excellent agreement with the perturbative results near the merger with the uniform brane. The numerical nonlinear results then show the regime where the perturbative analysis ceases to be valid and that ∆σ sameL,ρ decreases monotonically withρ (we have extended the computation to much higher values of ρ than those shown in the plot). The plot in the bottom row of Fig. 14 illustrates what happens for a constant-T lumpy brane family that bifurcates from an intermediate uniform brane In this case the perturbative analysis (see Fig. 10) tells us that the bifurcation occurs with ∆σ sameL,ρ > 0. Again the full nonlinear analysis confirms this is the case and is in excellent agreement with the perturbative results near the merger. However, in this case the nonlinear analysis provides new crucial information away from the merger: it shows that, although ∆σ sameL,ρ initially grows away from the GL merger, at a certain point it reaches a maximum and then it starts to decrease until it becomes negative. We have extended the computation to much larger values ofρ than those shown in the plot and we have found that, beyond this point, ∆σ sameL,ρ becomes more and more negative asρ becomes larger and larger. Since bothσ nu andσ u are non-negative, the reason for this is clearly thatσ u becomes arbitrarily large. In turn, this is due to the fact that, on a constant-T curve,L becomes larger and larger asρ increases (see Fig. 15), which causes the integral over the x-direction of the entropy densityŝ to diverge. At the value of (ρ,L) where ∆σ crosses zero, there is a phase transition between lumpy and uniform branes. This is a first-order phase transition since, for example, the temperature changes discontinuously. We emphasize that, at the qualitative level, this behaviour is the same for all constant-T families of lumpy branes that bifurcate from uniform branes in between points Σ 1 and Σ 2 in Fig. 11.
To further understand this phase transition, in Figs. 15 and 16 we reproduce again the stability diagram of Fig. 8, but this time we also plot a few constant-T or constant-L lumpy branes that depart from the GL merger curve. We use two plot marker codes: The solid blue markers, no matter their shape, represent the trench where ∆σ sameL,ρ > 0, while the empty orange markers, no matter their shape, describe the region where ∆σ sameL,ρ < 0. For reference, recall that For brevity, henceforth we will use the notation The main conclusions from Figs. 15 and 16 are as follows: 1. Recall that lumpy branes that bifurcate from the GL merger line at a temperaturê T <T Σ 1 (i.e. above Σ 1 in the figures) or at a temperatureT >T Σ 2 (i.e. below Σ 2 in the figures) have ∆σ sameL,ρ < 0 no matter how largeL is, as pointed out when discussing Fig. 14 (recall that Σ 1 and Σ 2 were introduced in Fig. 10). In Fig. 15 we display one family of lumpy branes in each of these classes. One has constant T × 0.390711 <T Σ 1 and is always described by empty orange colour markers, which means that the solutions indeed have ∆σ sameL,ρ < 0. All the curves that bifurcate from the GL merger above Σ 1 have this feature and they always bifurcate towards higherÊ and higherL with respect to the merger point (Ê GL ,L GL ). The other family has constantT ⊗ 0.405141 >T Σ 2 and is again always described by Figure 15. Same stability diagram as in Fig. 8 with the inclusion of some lumpy-brane curves that bifurcate from the GL merger curve. These curves have constantT given by {T × ,T ⊗ ,T ,T • ,T ,T } {0.390711, 0.405141, 0.395420, 0.395894, 0.396367, 0.397307}. Solid blue markers (empty orange markers), no matter their shape, indicate positive (negative) ∆σ sameL,ρ . empty orange markers, which means that the solutions indeed have ∆σ sameL,ρ < 0. All the curves that bifurcate from the GL merger below Σ 2 have this feature and they always bifurcate towards lowerÊ and higherL with respect to the merger point (Ê GL ,L GL ).
2. The situation is less monotonous for lumpy branes that bifurcate from a point on the GL merger curve that lies between Σ 1 and Σ 2 . Recall that these have, close to the merger, ∆σ sameL,ρ > 0. In Fig. 15 we display four curves in this class, namely: the family with constantT 0.395420 <T c ; the family with constant Figure 16. Same stability diagram as in Fig. 8 with the inclusion of some lumpy-brane lines that bifurcate from the GL merger curve. These lines have constantL given by {L ,L ,L ,L • } {5.299674, 6.004224, 6.900924, 11.501849}. Solid blue markers (empty orange markers), no matter their shape, indicate positive (negative) ∆σ sameL,ρ . Note that some orange circles are on top of some blue disks. This describes the region around the cusps of Fig. 17 and 18.
T • 0.395894 T c (so, very close toT c 0.3958945); the family with constant T 0.396367 >T c ; and the family with constantT 0.397307 >T c (slightly belowT K 0.397427). These curves withT Σ 2 <T <T Σ 2 bifurcate towardsL >L GL with ∆σ sameL,ρ > 0. Then, ifT Σ 1 <T <T c (e.g. the curve with diamond plot markers ) they typically move to higherÊ asL increases and ∆σ sameL,ρ changes from positive into negative when the plot markers change from solid blue into empty orange . On the other hand, ifT c <T <T Σ 2 (e.g. the curves initially with and ), the constantT -curves typically plunge into lowerÊ asL increases and ∆σ sameL,ρ changes from positive into negative when the plot markers change from solid blue into empty orange (i.e. → or → ).
3. WhenT ∼T c the properties described in the two previous points hold but the constant-T curves do not escape to largeÊ (ifT T ) or smallÊ (ifT T ) so quickly asL grows. A good example is given by the dotted (•) curve withT • 0.395894 T c .
The closer one is ofT c the longerL must be for the constant-T curve to cross the GL merger line again and then acquire ∆σ sameL,ρ < 0. Our results suggest that in the exact limitT →T c the curve extends toL → ∞ without ever leaving the window of energy densities [Ê B ,Ê A ].
To complete our understanding of the microcanonical phase diagram, in Fig. 17 we plot the Killing entropy density difference ∆σ sameL,ρ between lumpy and uniform branes with the same (L,ρ) as a function of the Killing energy densityρ for three families of lumpy branes that have constantL. Note thatL K <L Σ 1 <L Σ 2 (see e.g. Fig. 16). The three panels of Fig. 17 describe representative examples of the following three possible cases: (1) when the lumpy branes haveL K <L <L Σ 1 . We see that in this range ofL, lumpy branes (yellow inverted triangles) bifurcate from the uniform brane (blue line) at lowρ with ∆σ sameL,ρ > 0 and, asρ increases, the entropy difference grows until it reaches a maximum and then it decreases monotonically until the lumpy brane merges again with the uniform brane at higherρ. Since in this range of (L,ρ) one always has ∆σ sameL,ρ > 0, lumpy branes are the preferred phase in the microcanonical ensemble.
2. The top-right panel is for constantL 6.004224 solutions and illustrates what happens in the 3-dimensional microcanonical phase diagram ∆σ sameL,ρ versus (ρ,L) when the lumpy branes haveL Σ 1 <L <L Σ 2 . As in the previous case, in this range ofL, lumpy branes (brown diamonds) also bifurcate from the uniform brane (blue line) at lowρ with ∆σ sameL,ρ > 0 and, asρ increases, the entropy difference grows until it reaches a maximum. Then it again decreases monotonically but, this time, ∆σ sameL,ρ becomes negative at a certainρ. This first-order phase transition point is best seen in the inset plot that zooms into this region. The entropy difference keeps decreasing as ρ grows until it reaches a cusp. Then, asρ decreases, ∆σ sameL,ρ becomes less negative until the lumpy brane with constantL merges again with the uniform brane.
3. Finally, the bottom panel is for constantL 11.501849 solutions (whose energy density profile was discussed in Fig. 13). It illustrates how the 3-dimensional micro- Figure 17. Microcanonical phase diagram: Killing entropy density difference ∆σ sameL,ρ between lumpy and uniform branes with the same (L,ρ) as a function of the Killing energy densityρ for 3 families of lumpy branes that have constantL given byL 5.299674 <L Σ1 (top-left),L 6.004224 which is in the rangeL Σ1 <L <L Σ2 (top-right) andL 11.501849 >L Σ2 (bottom). For reference, L Σ1 5.618133 andL Σ2 6.592316. The three families are those with constant {L ,L ,L • } already displayed in Fig. 11; we use the same shape/colour coding for the markers execpt that here they are all solid. The horizontal blue line with ∆σ sameL,ρ = 0 describes the uniform brane family.
canonical phase diagram ∆σ sameL,ρ versus (ρ,L) looks like when the lumpy branes haveL >L Σ 2 . In this range ofL, at both GL mergers with the uniform brane (blue line), lumpy branes (orange disks) bifurcate with ∆σ sameL,ρ < 0. Then, as we move along the constant-L line away from the merger points, there are first two cusps (the left one is shown in more detail in the inset plot) and two first-order phase transition points where ∆σ sameL,ρ changes sign and becomes positive. Forρ in between these two transition points, one has a lumpy brane with ∆σ sameL,ρ > 0 and thus these lumpy branes are the preferred microcanonical phase. Otherwise, uniform branes dominate the microcanonical ensemble.
To complement this discussion, it is useful to plot the Killing entropy density difference ∆σ sameL,ρ between lumpy and uniform branes with the same (L,ρ) as a function of the average energy densityÊ for some families of lumpy branes that have constantL. Recall that the average energy density and the Killing energy density are related throughÊ =ρ/L. This means that comparing the entropy of uniform and nonuniform brane at the same (L,ρ) is the same as comparing them at the same (L,Ê). Fig. 18 shows this comparison for the four constant-L families {L ,L ,L ,L • } that were plotted in Fig. 16. Fig. 18, together with the projections to the (L,Ê)-plane shown in Fig. 16, is the key figure to understand the microcanonical phase diagram because it provides four representative slices of this plot at constantL. Gluing slices of this type together along theÊ-axis one obtains the threedimensional plot of ∆σ sameL,ρ versus (L,Ê). TheL → ∞ limit of the curves of Fig. 18 is the curve in Fig. 6(right).
So far we have discussed the lumpy branes only in the microcanonical ensemble. This is the most interesting ensemble because nonuniform branes can dominate this ensemble for certain windows of the parameter space and in a time evolution we typically fix the length and the average energy density of the solutions (i.e. the latter is conserved). But we may also ask about the role played by the lumpy branes in the canonical ensemble. In this case, we want to fix the length LΛ and the temperature T /Λ of the solutions and the dominant solution is the one that has the lowest Killing free energy density f /Λ 3 .
To address this question it is useful to first recall what happens when we consider only the uniform brane solutions. In the right panel of Fig. 3 we have already seen that the light uniform branch (the lower branch in the left panel of Fig. 3) is the preferred thermal phase forT <T c , while for fixedT >T c the heavy uniform branch (the upper branch in the left panel of Fig. 3) dominates the canonical ensemble. The intermediate uniform branch (between A and B in Fig. 3) is never a preferred thermal phase of the canonical ensemble. For this reason it is sometimes stated in textbooks that, atT =T c ∼ 0.3958945, there is a first-order phase transition at which the system jumps discontinuously between the light and heavy uniform branes. However, at infinite volume there is actually a degeneracy of states atT =T c because the average free energy density of any phase separated state is the same as that of the homogeneous states atT =T c . The reason for this is that the interface between the two phases in a phase-separated configuration gives a volumeindependent contribution. This contribution is therefore subleading in the infinite-volume limit with respect to those of the two coexisting phases, whose free energy densities are equal to each other and to those of homogeneous states atT =T c . Therefore, in the infinte-volume limit the system can transition between points D and E along a sequence of constant-temperature, constant free-energy, phase-separated states. The fact that all these states have the same free energy is the content of Maxwell's construction. Killing entropy density difference ∆σ sameL,ρ between lumpy and uniform branes with the same (L,ρ) or, equivalently, with the same (L,Ê), as a function of the average energy densityÊ. We show the same four families of lumpy branes with constant {L ,L ,L ,L • } {5.299674, 6.004224, 6.900924, 11.501849} already displayed in Fig. 11. We use the same shape/colour coding for the markers as in Fig. 11 except that here they are all solid. The families with , , • were also shown in Fig. 17, but the family was not. The horizontal blue line with ∆σ sameL,ρ = 0 describes the uniform brane family. The grey vertical lines indicate the turning points A and B in the phase diagram of Fig. 3. The labels "a" and "b" indicate the lumpy solutions withL • 11.501849 that lie away from the merger curve but have the same entropy density as the corresponding uniform branes. In other words, these are the points away from the merger curve at which ∆σ sameL,ρ crosses zero. The average energy densities at these points arê Figure 19. Canonical phase diagram: Dimensionless Killing free energy density difference ∆f sameL,T as a function of the dimensionless temperatureT for the six lumpy-brane families at constantL already shown in Fig. 11 (with the same colour/shape code). Namely, from the bot- In contrast, at finite volume the inhomogenous, phase-separated states are never thermodynamically favoured. To show this, in Fig. 19 we compare the free energy of lumpy branes with the light uniform branes ifT <T c , and with the heavy uniform branch ifT >T c . More concretely, we compute the difference between the Killing free energy of the nonuniform branef nu and the light (heavy) uniform Killing free energŷ f u whenT <T c (T >T c ) that has the same length LΛ and temperature T /Λ, i.e. ∆f sameL,T = f nu −f u sameL,T . Fig. 19 shows that for any temperatureT A ≤T ≤T B where nonuniform branes exist, one always has ∆f sameL,T > 0. That is to say, the Killing free energy density of the lumpy branes is always higher than the free energy of the relevant (light or heavy) uniform brane and thus lumpy branes never dominate the canonical ensemble. 13

Excited static lumpy branes: beyond the ground state solutions
So far we have discussed only the "ground state" lumpy branes of our spinodal system. The profile, for example that of the energy density E(x), of these fundamental branes has a single maximum and a single minimum, see Figs. 12 and 13. The phase diagram of the theory also contains infinitely many more lumpy brane phases whose profiles E(x) have η maxima and η minima for natural integer η. However, these are "excited states" of the theory in the sense that, as we will show below, for given (L,ρ) they always have lower Killing entropy densityσ than the ground state lumpy branes that we have constructed above. In other words, lumpy branes with η > 1 are subdominant phases of the microcanonical ensemble. In particular, this suggests that they should be dynamically unstable and evolve towards the fundamental lumpy brane if slightly perturbed. In the case of largeL this was explicitly verified in [2].
In principle, excited lumpy branes can be constructed using the perturbative method of Secs. 2.4.2 and 2.4.3. At linear order we would have to start with a Fourier mode that describes the η th harmonic of the system, namely with instead of the η = 1 case of (2.60). However, it is not necessary to perform this construction since the properties of these excited states can be obtained from those of the fundamental ones using extensivity. 14 Indeed, given a solution with η = 1 in a box of sizeL we can obtain a solution with η > 1 in a box of size ηL by taking η copies of the initial solution.
Once we know all solutions with η = 1 in boxes of any size, as we do, this procedure gives us all possible solutions with η > 1 maxima and minima in all possible boxes. Clearly, if the Killing energy and entropy densities of the initial solution areρ andσ, respectively, then those of the new solution are ηρ and ησ. In contrast, the average energy densitŷ E =ρ/L remains invariant. We must now compare the Killing entropy density of the solution with η maxima and minima with that of the corresponding η = 1 brane in a box of size ηL. Since the average energy density is invariant when taking copies of the initial solution, this comparison is most easily done by consideringσ as a function ofÊ andL. Therefore we must compare the entropy of the excited braneσ η (Ê, ηL) ≡ η ×σ(Ê,L) with that of the fundamental braneσ(Ê, ηL). It follows that if the entropy at fixedÊ grows withL faster than linearly then the fundamental brane always has higher entropy than the excited brane. This is indeed the case, as can be seen by taking constant-Ê slices of 13 For completeness, we have verified that the Killing free energy densityf of lumpy branes is always lower than the Killing free energy density of the intermediate branes AB, and that they become equal to one another precisely when the merger of these two branches occurs. In any case neither branch is ever preferred at finite volume in the canonical ensemble. 14 Similar arguments where used to find the thermodynamics of excited nonuniform black strings of the original GL system [26,27]). We can also start our linear order analysis with two (or more) harmonics with different amplitudes. This allows to construct lumpy branes with two (or more) maxima that have different amplitudes (in the spirit of [28]). Fig. 20 we do this forÊ = 0.85 and we compare ∆σ η=1 of the fundamental (η = 1) nonuniform brane (orange •) against ∆σ η=2 and ∆σ η=3 of the η = 2 (blue ) and η = 3 (green ) excited branes. For a given (L,ρ) or, equivalently, for a fixed (L,Ê), we see that the Killing entropy density decreases as η grows: in agreement with the most naive intuition, the fundamental lumpy brane has the highest Killing entropy density and therefore it dominates the microcanonical ensemble over any excited brane. The discussion above applies to excited states that can be obtained as copies of a single configuration. Therefore states of this type with η maxima and minima have a Z η discrete symmetry. There exist more general excited states with maxima and minima of different heights, but we expect these to be subdominant too. In the case of largeL this was explicitly verified in [2].

The spinodal (Gregory-Laflamme) timescale
In Sec. 2.4.2 we saw that intermediate uniform branes withÊ B <Ê <Ê A (see left panel of Fig. 3) can be GL-unstable. To find when the instability appears, we took the uniform branes Q j (y) of section 2.4.1 and considered static Fourier perturbations of the form (2.60) about this background, namely q j (x, y) = Q j (y) + q (1) j (y) cos(πx). This allowed us to find the minimum length L GL Λ = 2π/k GL (see Figs. 7 and 8) above which the uniform brane is unstable. This was enough for our purposes of Sec. 2, where we were just interested in finding the static lumpy branes. In particular, we found large regions of the microcanonical phase diagram where lumpy branes coexist with and are favoured over uniform branes. This suggests that, if we start with initial data that consists of a uniform brane that is GL unstable plus a perturbation, the system should evolve towards a lumpy brane with the same lengthL and Killing energy densityρ, and hence also the sameÊ =ρ/L. The initial stages of this time evolution should be well described by the linear GL frequencies. It is thus important to compute the GL timescales of the system.
Consider again the uniform branes constructed in Sec. 2.4.1 in the regimeÊ B <Ê <Ê A (see left panel of Fig. 3). Denote the collective fields byψ(y) = {ḡ µν (y),φ(y)}. We will now allow for time-dependent perturbations of this background. More concretely, we will use the fact that ∂ t and ∂x are Killing vector fields of the uniform brane background to Fourier decompose the time dependent perturbations as ψ(t, x, y) =ψ(y) + δψ (1) (y)e i kx e −i ω t . (2.83) This introduces the wavenumber k conjugate to the spatial directionx = x L 2 ∈ [0, L/2] and the frequency ω of the perturbation. Let δg µν ≡ h µν be the metric perturbations and δφ the scalar field perturbation. Perturbations δψ (1) (y) = {h µν (y), δφ(y)} that break the symmetries indicated in (2.83) excite a total of 8 fields, namely: δφ, h tt , h ty , h tx h yy , hx y , hxx and h x 2 We have not yet fixed the gauge freedom of the problem. Instead of doing so we construct two gauge invariant-quantities that encode the most general perturbations of the form (2.83) as described in [29]. The linearized Einstein equations then reduce to (and are closed by) a coupled system of two linear, second-order ODEs for these two gauge-invariant variables. The perturbations must be regular at the horizon in ingoing Eddington-Finkelstein coordinates and preserve the asymptotic AdS structure of the uniform background. This is a non-polynomial eigenvalue problem for the frequency ω where we give the uniform background and the wavenumber k and find ω. The GL modes of the uniform brane system have purely imaginary frequency.
In Fig. 21, as an illustrative example, we plot the dimensionless dispersion relationω(k) for a particular uniform brane with (τ ,Ê) (0.395894, 0.867966) that is very close to point C in Fig. 3, for which (τ ,Ê) C (0.3958945, 0.867956). 15 In Fig. 21, the red circle curve describes the dispersion relation of the fundamental harmonic η = 1. Not surprisingly, this curve starts at (k,ω) = (0, 0) and, ask increases, the dimensionless frequency Im ω/Λ first grows until it reaches a maximum and then starts decreasing. Precisely at the GL critical wavenumberk =k GL = 1.322499, as computed independently in Fig. 7, one has Im ω/Λ = 0 and fork >k GL the uniform brane is stable. For 0 <k <k GL the uniform brane is GL unstable and the maximum of the instability occurs at (k,ω)| max (0.626902, 0.120249 i).
Besides the fundamental GL mode, the uniform brane has an infinite tower of integer η spatial Fourier harmonics. Uniform branes are also unstable to these higher harmonics but the minimum unstable GL lengthL GL,η for the η th harmonic increases with η or, equivalently, the critical GL wavenumberk GL,η decreases with η. As examples, in Fig. 21  we also plot two other curves that describe the dispersion relation of the second (η = 2, blue ) and third (η = 3, green ) harmonics. Note that the dispersion relation of these higher harmonics can be obtained straightforwardly from that of the fundamental harmonic. Indeed, note that we can unwrap the S 1 and change the periodicity of its coordinatex from L to L η = η L, for integer η [26,27]. This also changes the wavenumber from k = 2π L into k η = k η . Altogether this leaves the phase of the Fourier mode e ikx invariant. But this means that the frequency ω η of the η th harmonic is related to the frequency of the fundamental harmonic simply by ω η (k) = ω(k/η) and that the critical GL zero mode of the η th harmonic isL GL,η = ηL GL ork GL,η =k GL /η. These properties, namely ω η (k) = ω(k/η) ,k GL,η =k GL /η , (2.84) are indeed observed in Fig. 21.
The linear results of Fig. 21 also provide a guide to the full nonlinear time evolution of nonconformal branes. In a microcanonical ensemble experiment, imagine that we start with a uniform brane in the regimeÊ B <Ê <Ê A where it can co-exist with lumpy branes, for example withÊ 0.867966. We want to perturb it to drive it towards a lumpy brane with the sameL andρ and thus sameÊ. What should we do? We certainly have to consider a Fourier perturbation withk <k GL as read from Fig. 21 or from Figs. 7 and 8. In these circumstances we still have different options that will result in substantially different time evolutions. Indeed, if we start with ak GL,2 <k <k GL where only the fundamental harmonic is unstable then the system will evolve "quickly" towards an η = 1 lumpy brane (the quickest evolution should occur ifk ∼k| max ). More generically this will still be the case also for ak <k GL,2 as long as it is higher than the criticalk ∼ 0.4205 where the curves for η = 1 (•) and η = 2 ( ) meet, see the right-most dotted vertical line in Fig. 21. (In this discussion we assume that the initial amplitudes of all modes are similar.) If instead 0.2513 k 0.4205, i.e. in between the two vertical dotted lines of Fig. 21, then the time evolution of the uniform brane should first approach an η = 2 lumpy brane before finally moving towards the fundamental η = 1 lumpy brane, which has a higher Killing entropy density. Finally, if the uniform brane is perturbed with ak 0.2513 mode then the system will first evolve towards an η ≥ 3 lumpy brane before being driven towards its fundamental lumpy brane endpoint. These expectations were explicitly verified in the case of large boxes in [2].

Real-time dynamics
Above we have constructed inhomogeneous static solutions using purely static methods to solve the Einstein equations. We will now examine several aspects of these solutions using real-time dynamical methods. We will first reproduce the static solutions obtaining excellent agreement. Then we will use the dynamical methods to address two novel aspects not studied above: the local dynamical stability of the inhomogeneous static solutions, and the full time evolution, including the end state, of the unstable solutions. The reader interested in the numerical methods that we use can consult e.g. [1,2,10,30,31].

Reproducing the static solutions from real-time dynamics
In Fig. 22 we compare the Killing entropy density of the static inhomogeneous solutions obtained with dynamical methods (black dots) and with static methods (orange dots) for a system withL 11.501849. In Fig. 23 we compare the average energy density-versustemperature relation. As is clear from the figures we find excellent agreement.
The process that we follow to reproduce the static inhomogeneous solutions from realtime dynamical evolution makes it natural to distinguish three cases: (I) Lumpy branes whose (L,ρ) lies inside the GL merger curve of Figs. 15 and 16. An example is given by the solution labelled as X a in Fig. 22.
(II) Lumpy branes that are outside the merger curve and have the largest entropy among lumpy branes with the same (L,ρ). An example is given by the solution labelled as X b in Fig. 22.
(III) Lumpy branes outside the merger curve with the smallest entropy for a given (L,ρ). An example is given by the solution labelled as X c in Fig. 22.
We follow different strategies to find each of these types of solutions. Solutions of type I are reproduced by following the full evolution of the spinodal instability, as in [2,1]. The initial state is a homogeneous brane with the same (L,ρ) of the lumpy brane that we want to obtain plus a small sinusoidal perturbation corresponding to the lowest Fourier mode that fits in the box. As this solution lies inside the GL merger, this perturbation is unstable and grows with time. 16 Upon dynamical evolution the system eventually enters the nonlinear regime and finally relaxes to the inhomogeneous solution. In Fig. 24 we show an example of one of these evolutions (top-left) and the comparison of the solution at asymptotically late times with the solution obtained via static methods (top-right), with excellent agreement. . Average energy density versus temperature for static inhomogeneous solutions obtained with dynamical methods (black dots) and with static methods (orange dots) for a system witĥ L 11.501849. Blue (red) curves indicate locally stable (unstable) solutions. Orange dots are exactly as in Fig. 13(left).
The previous procedure fails to produce solutions of type II because the homogeneous system is locally stable, so small perturbations decay in time and the system returns to the initial homogeneous state. Indeed, we consider a uniparametric family of perturbations, not necessarily sinusoidal, with the parameter given by the amplitude A of the perturbation, and find that if A is smaller than a certain critical value A * then the system evolves back to the homogeneous state. In order to obtain a lumpy brane as a final state we must start with a homogenous brane plus a perturbation that is so large that the system finds itself directly in the non-linear regime. This is indeed what happens if A > A * . In this case the system evolves in time towards the globally preferred state, namely towards a lumpy brane like the one labelled X b in Fig. 22. An example of this evolution is illustrated in Real-time evolution leading to a type I (top-left), a type II (middle-left) and a type III (bottom-left) lumpy brane withL 11.501849, labelled X a , X b and X c in Fig. 22, respectively. On the right panels we compare the energy density profiles at late times (continuous black lines) with those obtained by static methods (orange dots).
Finally, if the amplitude of the perturbation is tuned to be exactly A * then the system evolves in time towards a type III solution like the one labelled as X c in Fig. 22. An example of this evolution is illustrated in Fig. 24(bottom-left). The fact that A must be precisely tuned in order to reach the type III solution suggests that these solutions are locally dynamically unstable. We will verify this explicitly below. Since numerically it is impossible to tune A with infinite precision, this means that if we were to evolve the configuration in Fig. 24(bottom-left) for sufficiently long times we would see that either it falls back to the homogeneous state (if A is slightly smaller than A * ) or it evolves towards a type II configuration (if A is slightly larger than A * ). We will confirm this in Fig. 28.

Local stability
In this section we study the local stability of the static inhomogeneous solutions by using real-time dynamical methods. We consider an initial state given by the static inhomogeneous solution plus a small perturbation and study its time evolution. The system is said to be locally stable if all possible linear perturbations decay in time. If at least one of the perturbations grows in time, then the system is said to be locally unstable. In order to establish which is the case one must decouple, i.e. diagonalise, the full set of linearized equations around the inhomogeneous solution (note that all Fourier modes are indeed coupled to one another because the inhomogenous state breaks translational invariance). Each eigenmode then evolves in time as e −iωt , with ω the corresponding eigenfrequency. If the imaginary part of all the eigenfrequencies is negative the system is locally stable. If at least one of the eigenfrequencies has a positive imaginary part then it is locally unstable.
Rather than performing the exercise above, we will use our numerical code to obtain the time evolution of a generic small initial perturbation of the inhomogeneous state. Since the perturbation is generic we expect that it will be a linear combination of all the eigenmodes of the system. Thus, after some characteristic time, the eigenmode with the largest imaginary part of omega will dominate the evolution leading to a well defined exponential evolution. We have identified this region of exponential behaviour in all the time evolutions of the perturbed system that we have studied, and we have obtained the real and imaginary parts of omega for the dominant mode by performing fits. Note that this will not result in a mathematical proof of local stability. For example, our generic perturbations may accidentally have a very small projection on some unstable mode, or the positive imaginary part of the frequency of this mode may be exceedingly small and hence go unnoticed, etc. While these possibilities cannot be excluded with absolutely certainty, the detailed searches that we have performed, together with the consistent emergent physical picture, make us confident that they are highly unlikely.
Let us illustrate the procedure with the two examples in Fig. 25. The top panel corresponds to the relaxation to equilibrium at late times of the simulation presented in Fig. 24(top-left). Specifically, we take the spatial profile of the energy density at some late time, we subtract from it the profile of the inhomogeneous static solution (which we denoted as X a in Fig. 22), and we decompose this difference into Fourier modes. The time-dependent amplitude of the first few of these Fourier modes is shown in Fig. 25(top). We see that all of these modes oscillate and decay exponentially in time with the same frequency. This is as expected since the evolution is dominated by the single eigenmode with the slowest decay. The fact that there is no growing mode indicates that the type mode has leading amplitude. The dotted line corresponds to a fit to the envelope, from which we extract the imaginary part of omega. The fact that no mode grows in time indicates that X a is locally dynamically stable. (Bottom) Time evolution of some Fourier modes of a perturbation around the type III, inhomogeneous, static configuration X c of Fig. 24(bottomright) withÊ 0.651,L 11.501849. In this case the leading mode grows exponentially in time, indicating that X c is locally dynamically unstable. I, inhomogeneous, static solution to which this configuration asymptotes at late times (namely, X a in the current simulation) is locally dynamically stable.
The second example shown in Fig. 25(bottom) corresponds to the evolution presented in Fig. 24(bottom-left), but extended to longer times. Here we plot the Fourier modes corresponding to the difference between the spatial energy profile at a given time and the spatial energy profile of the inhomogeneous, static configuration that we denoted as X c in Fig. 22. We observe a first relaxation in which the stable modes decay but, this time, at later times the system is dominated by an exponential growth of an unstable mode. This confirms that the type III lumpy branes such as X c are locally dynamically unstable, as anticipated above. Recall that the initial state in this time evolution is a homogeneous brane plus a large perturbation of amplitude A that is tuned to be close to a critical value A * . This tuning is what suppresses the initial amplitude of the unstable mode, hence allowing the time evolution to drive the system close to X c for some time. Thus, intuitively, this solution behaves like a saddle point in configuration space with some stable and some unstable directions (i.e. a metastable configuration).
We have performed a scan to determine the real and imaginary parts of omega for the dominant mode of static inhomogeneous solutions withL 11.501849 and varying energy densities. The result is shown in Fig. 26. We find that type I and II solutions have negative imaginary parts of omega, and so they are locally stable, while the type III solutions have positive imaginary parts of omega, and so they are locally unstable. The imaginary part of omega crosses zero precisely at the "cusps" of Fig. 22, that is, at the static solutions lying precisely at the boundary between type II and III solutions. The real part of omega is non-vanishing in the locally stable cases, whereas it vanishes in the locally unstable cases, going to zero also at the cusps. In Figs. 22, 23 and 26 we show locally stable solutions in blue and locally unstable solutions in red.
In this section we have discussed the (in)stability of what we called "ground-state" or "fundamental" solutions in Sec. 2.6, namely of solutions whose spatial energy density profile has a single maximum and a single minimum. Here we have not explicitly investigated the case of "excited" solutions, namely those with multiple maxima and minima. However, some configurations of this type were studied in [2], and in all cases they were found to be locally dynamically unstable. As discussed in Sec. 2.6 and Sec. 2.7 the reason is that the entropy density can be continuously increased by moving two of these maxima or minima towards each other. Since this seems to be a generic feature, we expect all excited configurations to be locally dynamically unstable.

Full time evolution of the unstable solutions
In the previous section we studied the local stability properties of the inhomogeneous static solutions, finding some regions of local instability. A natural question is therefore what is the end state of the evolution if these locally unstable solutions are perturbed. In this section we perform the full time evolution of the system and determine the end state.
Given a locally unstable solution there are two natural possibilities for the end state of the evolution. In Fig. 27 we present a concrete example where we show the three static solutions with the same (L,Ê) (11.501849, 0.651): X c , X d and X e , where X c is the static solution presented in Fig. 24(bottom-right). For the locally unstable solution X c , the two possible candidates for the end state of the evolution are the homogeneous solution X d and the inhomogeneous solution X e , since both of these have larger entropy than X c . By performing full time evolution we confirm that both solutions X d and X e can be the end  Fig. 22, whereL 11.501849. X c , X d and X e are the three static solutions with the same average energy densityÊ 0.651. The nonlinear time evolution of Fig. 28(left) corresponds to an evolution from X c to X d , and Fig. 28(right) corresponds to an evolution from X c to X e . state of the evolution, and that which one is reached depends on the initial perturbation.
In order to illustrate this we essentially extend the range of the time evolution shown in Fig. 24(bottom-left). Recall that in that figure we dynamically generated a solution very close to X c by fine-tuning the amplitude of the initial perturbation to be close to the critical value A * . Since the amplitude we choose is close but not exactly equal to A * , the result of this time evolution at intermediate times is not exactly the solution X c but X c plus a small perturbation. Since the perturbation is small the system spends a sizeable amount of time in a very slowly evolving configuration close to X c , as can be seen from the intermediate-time behaviour in Fig. 28. However, if the exact amplitude is slightly smaller than the critical one then further time evolution eventually drives the system back to the homogeneous solution labelled as X d in Fig. 27. This is the case in Fig. 28(left). If instead the amplitude is slightly larger than the critical one then the system eventually evolves towards the stable, inhomogeneous solution labelled as X e in Fig. 27. This is the case in Fig. 28(right). Note that the evolution from the unstable to the stable solutions can be viewed as a long, approximately linear regime (when the unstable solution is perturbed), followed by a fast non-linear regime, further followed by another long, approximately linear regime (when the system relaxes to the corresponding stable solution). We have verified that, at the qualitative level, these results apply to all the -60 -E Λ 4 x tΛ E Λ 4 x tΛ Figure 28. Extension to longer times of the time evolution shown in Fig. 24(bottom-left), whose initial state is a homogeneous configuration plus a large perturbation of amplitude A. At intermediate times this generates the unstable solution X c plus a small perturbation. (Left) If the amplitude of the initial perturbation A is slightly smaller than the critical value A * then the system eventually evolves back to the homogeneous solution labelled as X d in Fig. 27. (Right) If the amplitude of the initial perturbation A is slightly larger than the critical value A * then the system eventually evolves towards the stable, inhomogeneous solution labelled as X e in Fig. 27.
unstable solutions withL 11.501849 and varying energy densities that we have studied.

Discussion
We have considered a bottom-up, five-dimensional gravity model that, at infinite volume, possesses a first-order, thermal phase transition in the canonical ensemble. As usual, we expect this model to be holographically dual to a strongly coupled, large-N gauge theory in four dimensions. We have placed the system in a box and, for simplicity, we have imposed translational invariance along two sides of the box. We have varied the volume by varying the sizeL = ΛL of the third side, where Λ is the microscopic scale of the gauge theory. We have then constructed what we believe is the complete set of all possible homogeneous or inhomogeneous equilibrium states at finiteL. On the gravity side these correspond to uniform or lumpy branes, respectively. Although we do not have a mathematical proof that this set is indeed complete, we have found no evidence to the contrary in our extensive investigations based both on static and dynamical methods.
The first effect of the finite volume is that some homogeneous states between points A and B in Fig. 1 become locally dynamically stable, as illustrated in Fig. 8. The reason for this is that the spinodal instability is a long-wavelength instability. IfL is below a certain energy-dependent value, then the potentially unstable mode does not fit in the box and the corresponding homogeneous state is actually stable. The unstable states are those in the region inside the parabola in Fig. 8. We see that asL → ∞ we recover the fact that all states with energy densities between A and B are unstable, but that at finiteL some of them are stable. In particular, there is a valueL =L K below which all homogeneous states are locally dynamically stable since none of them can accommodate an unstable mode.
The parabola in Fig. 8 is a curve of marginal stability. Therefore we expect a branch of static, inhomogeneous states to emanate from each point on this curve. One should think of the extra direction in which these branches emanate as the entropy relative to that of the homogeneous state, ∆σ. The union of all such inhomogeneous branches is therefore a surface in the three-dimensional space parametrized by the average energy density in the box E, the size of the boxL, and the entropy ∆σ. We will refer to this surface as the "entropy surface". The intersection of this surface with theÊ-L plane contains the parabola in Fig. 8 (as well as other points such as the points a and b of Fig. 18). The curves in Fig. 15 are the projections on this plane of constant-T slices of the entropy surface. Similarly, the vertical lines in Fig. 16 are the projections on this plane of constant-L slices. The same slices projected onto theÊ-T plane are shown in Fig. 11. This last figure makes it clear that inhomogeneous states only exist in the range of temperaturesT A ≤T ≤T B .
The structure of the entropy surface is most easily understood by thinking of it as the union of constant-L slices for allL >L K . The shape of each of these slices as a function of the energy density is shown in Figs. 17 and 18. We see that, at the qualitative level, there are three possibilities depending on the value ofL in relation to the following hierarchŷ (4.1) These three length scales are an intrinsic property of the theory at finite volume and their values are given in (2.80). IfL K <L <L Σ 1 <L Σ 2 , then ∆σ is always positive for all the values of the energy for which inhomogeneous states exist. This is the case illustrated by Fig. 17(top-left) and by the bottom curve with beige inverted triangles in Fig. 18. If insteadL K <L Σ 1 <L <L Σ 2 then the ∆σ curve becomes negative and develops a cusp near its endpoint on the right-hand side. This is the case illustrated by Fig. 17(top-right) and by the second-from-the-bottom curve with brown diamonds in Fig. 18. Finally, if L K <L Σ 1 <L Σ 2 <L then the ∆σ curve becomes negative and develops cusps near both of its endpoints. This is the case illustrated by Fig. 17(bottom) and by the two top curves with orange circles and green squares, respectively, in Fig. 18.
The shape of the entropy surface that we have just described determines the structure of phase transitions in the microcanonical ensemble. Recall that in the limitL → ∞ the set of globally preferred, maximum-entropy states are those indicated by the black curves (with arrows) in Fig. 1 (see also the discussion around Figs. 5 and 6). The direction of the arrows in Fig. 1 indicates what happens as the energy decreases from an arbitrarily high value. As the energy density decreases towards point D the preferred states are homogeneous branes of decreasing temperature. At D there is a phase transition into inhomogeneous states of constant temperature T c . SinceL → ∞ these are phase-separated configurations in which the homogeneous phases D and E coexist. At E there is another phase transition, in this case from inhomogeneous to homogeneous states. The fact that the fraction of the total volume occupied by each phase varies continuously between 0 and 1 as the energy density varies between D and E suggests that these transitions are continuous in the microcanonical ensemble. Continuity can also be seen more formally as follows. For fixed length and source, the first law (2.49) takes the form 1 In the microcanonical ensemble the total entropy S is the relevant thermodynamic potential, the total energy E is the control parameter and T is a derived quantity. At the points D and E the temperature is continuous, but its derivative dT /dE is not, because this is positive on the homogeneous branch but it vanishes on the inhomogeneous one. This picture is modified at finiteL. Note that in this case the system may still exhibit phase transitions since the planar limit that we work in, N → ∞, acts effectively as a thermodynamic limit. Consider first Fig. 29 which illustrates the structure of phase transitions for a lengthL 11.501849 such thatL K <L Σ 1 <L Σ 2 <L (this is the value corresponding to the orange circles in Fig. 11). In this case the preferred states lie on the homogeneous branch until the energy density reaches that of point a in Fig. 18. At this point a first-order phase transition takes place between the homogeneous state and the state a on the inhomogeneous branch, as indicated by the top horizontal arrow in Fig. 29. Note that this transition can take place before (as in the case of the orange circles in Fig. 11) or after (as in the case of the green squares Fig. 11) the turning point A is reached. The reason that this is a first-order transition is that the temperature changes discontinuously. As the energy decreases further, the preferred states are those on the inhomogeneous branch until the energy density reaches that of point b in Fig. 18. At this point another first-order phase transition takes place between the inhomogeneous state and a state on the homogeneous branch with the same average energy density, as indicated by the bottom horizontal arrow in Fig. 29. Note that this state is below the turning point B (i.e.Ê b <Ê B ). As the energy is further decreased the preferred state remains on the homogenous branch.
Thinking of the infinite-volume case as the limitL → ∞ of the situation described in Fig. 29 sheds light on the order of the phase transition at infinite volume. AsL increases the point a in Fig. 29 moves up and to the right. This means that the homogeneous and the inhomogeneous states between which the transition takes place become closer to one another. In the limitL → ∞ the point a tends to the point D (of Figs. 1 or 3) and the transition takes place between two states at the same temperatureT =T c . Since the discontinuity inT disappears the phase transition becomes second-order.
Consider now Fig. 30, which illustrates the structure of phase transitions for a lengtĥ L 5.299674 such thatL K <L <L Σ 1 <L Σ 2 (this is the value corresponding to the inverted beige triangles in Fig. 11). In this case the preferred states lie on the homogeneous branch until the merger point with the inhomogeneous branch is reached. At this point a transition between the homogeneous and the inhomogeneous branches takes place. Since the transition happens at the merger point, the temperature is continuous and the transition is second-order. As the energy is further decreased the preferred states remain on the inhomogeneous branch until this merges again with the homogeneous branch. At this point another second-order phase transition takes place. Below this point the preferred state lies on the homogeneous branch.
In the intermediate range of lengthsL K <L Σ 1 <L <L Σ 2 the structure of transitions is a hybrid between those described in Figs. 29 Figure 29. Phase diagram in the microcanonical ensemble for a lengthL 11.501849 such thatL K <L Σ1 <L Σ2 <L. The blue and the orange curves correspond to homogeneous and inhomogeneous states, respectively. The orange curve is the same as the curve of orange circles in Fig. 11. Solid segments indicate locally dynamically stable states; dashed segments indicate unstable ones. The black curves with arrows indicate the sequence of globally preferred, maximumentropy states as the average energy density decreases. The points a and b corrrespond to those in Fig. 18. The first and fourth (from top to bottom) dashed horizontal lines indicate the energy densities at these points, whereas the second and third lines indicate the energy densities at the turning points A and B. The phase transitions between the homogeneous and the inhomogeneous branches, indicated by the horizontal arrows, are first-order. first a first-order phase transition between the homogeneous branch and the inhomogeneous branch. In our model the point on the homogeneous branch lies between A and Σ 1 , and the point on the inhomogeneous branch is the analog of point a in Fig. 29. As the energy is further decreased the preferred state remains on the inhomogeneous branch until this merges with the homogeneous one. At this point a second-order phase transition occurs in which the preferred state becomes the one on the homogeneous branch. This second transition is analogous to that in Fig. 30 Figure 30. Phase diagram in the microcanonical ensemble for a lengthL 5.299674 such that L K <L <L Σ1 <L Σ2 . The blue and the beige curves correspond to homogeneous and inhomogeneous states, respectively. The beige curve is the same as the curve of beige inverted triangles in Fig. 11. Solid segments indicate locally dynamically stable states; dashed segments indicate unstable ones. The black curves with arrows indicate the sequence of globally preferred, maximum-entropy states as the average energy density decreases. The phase transitions between the homogeneous and the inhomogeneous branches take place at the points where these branches merge and they are second-order. this branch.
Finally, for lengths such thatL <L K <L Σ 1 <L Σ 2 , no inhomogeneous states exist and all the homogeneous ones are dynamically stable. In this case no phase transitions occur as the energy decreases from infinity to zero, as illustrated in Fig. 31.
In the figures above we have used continuous and dashed segments to distinguish between locally dynamically stable and locally dynamically unstable states. In the case of homogeneous states these properties can be established via a perturbative analysis. In the case of inhomogeneous states we used a numerical code for time evolution to study the behaviour of small perturbations. The results are shown in Fig. 22 and can be succinctly summarised as follows: all states on the upper part of the curve, shown in blue, are stable, whereas those on the lower part, shown in red, are unstable (see also a relevant zoom in Fig. 27). Once an unstable state is slightly perturbed, its full time evolution and its end state depend on the perturbation. For example, the locally dynamically unstable state X c in Fig. 22 and Fig. 27 can decay to either of the stable states X d or X e . The corresponding time evolutions are shown in Figs. 28(left) and 28(right), respectively.
For completeness, note that, unlike in the microcanonical ensemble, at finite volume lumpy branes are never the dominant thermodynamic phase in the canonical ensemble, as illustrated by Fig. 19.
In our analysis we have benefited from two simplifying assumptions. The first one is that we imposed translational invariance along two of the three spatial directions of the box. Lifting this restriction will generically allow for inhomogeneities to develop in all three directions. It would be interesting to study this more general setup, in particular the possible interplays between different length scales in different directions. Hopefully, the "one-dimensional building blocks" that we have investigated will be useful to understand the three-dimensional case.
The second simplifying assumption is the fact that we have worked in the N → ∞ limit, which ensures that the system is in the thermodynamic limit despite the finite volume. In particular, it guarantees that true phase transitions may occur. At finite N these transitions will turn into cross-overs. However, the latter can be made arbitrarily rapid by making N sufficiently large. This means that our results should be a good approximation to the physics at finite but large N .