Phases of Holographic Hawking Radiation on spatially compact spacetimes

We study phases of equilibrium Hawking radiation in $d$-dimensional holographic CFTs on spatially compact spacetimes with two black holes. In the particular phases chosen the dual $(d+1)$-dimensional bulk solutions describe a variety of black funnels and droplets. In the former the CFT readily conducts heat between the two black holes, but it in the latter such conduction is highly suppressed. While the generic case can be understood in certain extreme limits of parameters on general grounds, we focus on CFTs on specific geometries conformally equivalent to a pair of $d \ge 4$ AdS${}_d$-Schwarzschild black holes of radius $R$. Such cases allow perturbative analyses of non-uniform funnels associated with Gregory-Laflamme zero-modes. For $d=4$ we construct a phase diagram for pure funnels and droplets by constructing the desired bulk solutions numerically. The fat non-uniform funnel is a particular interesting phase that dominates at small $R$ (due to having lowest free energy) despite being sub-dominant in the perturbative regime. The uniform funnel dominates at large $R$, and droplets and thin funnels dominate at certain intermediate values. The thin funnel phase provides a mystery as it dominates over our other phases all that way to a critical $R_{\mathrm{turn}}$ beyond which it fails to exist. The free energy of the system thus appears to be discontinuous at $R_{\mathrm{turn}}$, but such discontinuities are forbidden by the 2nd law. A new more-dominant phase is thus required near $R_{\mathrm{turn}}$ but the nature of this phase remains unclear.


Introduction
Understanding the thermodynamics of strongly interacting field theories remains a challenging task. This is no less the case when the theory is coupled to heat baths provided by black holes, in which context the results describe phases of the associated Hawking radiation. But for appropriate large N conformal field theories (CFTs), gauge/gravity duality [1] provides what one hopes may be a tractable description in terms of semi-classical gravity in asymptotically (locally) anti-de Sitter (AlAdS) spacetimes. Below we consider the bulk classical limit in cases where the bulk description may be truncated to Einstein-Hilbert gravity with a Λ < 0 cosmological constant. Early explorations [2,3,4,5,6,7,8,9,10,11] of both this setting and the related context of brane-world black holes found bulk solutions for which the dual CFT Hawking radiation behaved quite differently from that of a free theory. In particular, despite a large density of CFT states s CFT , the flux of energy to infinity remained quite small (order 1). However, inspired in part by [12], it was later recognized in [13] that such solutions represented only one possible phase of the Hawking radiation with other phases allowing heat transport of the same order as the CFT entropy density s CFT . Thus the CFT undergoes a conducting/insulating phase transition that is related by a conformal transformation to the more familiar confinement/deconfinement transition [14].
The properties of the above phases are readily seen from the dual d + 1-dimensional AlAdS bulk solutions. The induced conformal metrics on their conformal boundaries must all agree with that of the d-dimensional black hole spacetime on which the CFT is defined. Following [13], we thus refer to the CFT heat baths as boundary black holes. There will also be one or more black holes in the bulk whose horizons end on those of the boundary black holes. Bulk horizons conduct heat along themselves at the classical level [15] and thus at a level proportional to the CFT entropy of states s CFT . But two disconnected bulk horizons exchange heat only via bulk Hawking radiation, which remains an order 1 effect at large s CFT . The basic thermal conductivity properties of the CFT are thus determined by the pattern of bulk connections between the various boundary black holes. This argument, presented in [13], has been verified by direct studies of connected [16,17,18] and isolated [19] boundary black holes. See [20] for a simple solvable example describing heat transport along an analogous horizon.
Bulk horizons connecting two or more boundary black holes have become known as black funnels, while bulk horizons that connect to only one boundary black hole are called black droplets [13]. In this terminology it is standard to treat any asymptotic region of a globally hyperbolic CFT spacetime as an additional black hole, and thus to also apply the term funnel to bulk horizons that connect boundary black holes to the asymptotic regions of the CFT. Such asymptotic regions are in fact conformally equivalent to black holes in simple cases (see [21], based on [22,23,24]). More generally, the asymptotic regions of globally hyperbolic CFT spacetimes can be mapped to spacetimes with null singularities which one might consider to be singular black holes.
Our purpose below is to study phase transitions between AlAdS (d+1) funnels and droplets associated with a particular class of boundary black hole spacetimes conformally equivalent to a pair of (global) d ≥ 4 Schwarzschild-AdS d glued together along the AdS boundary. See figure 1 and section 3 for technical details. The resulting boundary spacetime is conformally equivalent to a pair of black holes in the Einstein Static Universe (ESU), generalizing the familiar way in which the ESU S d−1 × R is conformally equivalent to two copies of AdS d and gives the conformal boundary of global AdS d+1 . Thus we consider what one might call global droplets and global funnels. The analogous phase transitions were studied analytically for d = 3 with Bañados-Teitelboim-Zanelli (BTZ) boundaries in [25], where it was argued that all relevant phases for this case can be constructed by double Wick rotation of global AdS 4 and Schwarzschild-AdS 4 .
While a great many d ≥ 4 funnel and droplet solutions have by now been constructed [26,27,25,28,29,30,31,19,21,32,33,34,35,36], our work is the first to construct multiple such phases for the same boundary black holes. This allows a meaningful comparison of their free energies. In particular, although both droplet and funnel free energies F = E − T S receive divergent contributions from the infinite area of the non-compact bulk horizons, the fact that both solutions satisfy the same boundary conditions and that the AlAdS boundary is spatially compact implies that the difference ∆F must be finite and unambiguous as defined using any Fefferman-Graham regulator 1 . We consider equilibrium situations in which all bulk horizons are at the same temperature T , which we take to agree with the temperatures of all boundary black holes. It should be noted that more general 'detuned' equilibrium solutions should also exist, where the temperatures of the bulk and boundary black holes differ. Such solutions were predicted in [18] based on [22,23,24] and constructed with specific boundary metrics in [36]. As we discuss in section 3, the appearance of the Gregory-Laflamme instability in 5 or more bulk dimensions (i.e., for d ≥ 4) leads to additional funnel phases not seen in the d = 3 analysis of [25].
We begin with a brief overview in section 2 of the phases to be expected for general ESU boundary black holes. Here we simply characterize the boundary spacetime by the size R and temperature T of the boundary black hole and the ESU length scale d , describing the phases to be expected in various extreme limits of the dimensionless parameters R/ d and T d . This generalizes the general discussion of phases in [13]. Further insight into droplet phases is then obtained by recalling the d = 3 analysis of [25] (section 2.3), and insight into funnel phases is obtained by considering the Gregory-Laflamme instability (section 3) with some detailed expressions relegated to appendix A. Numerical methods and the framework for our calculations are described in section 4 while diagnostic machinery is described in section 5. Results for d = 4 are presented in section 6 and are supported by convergence tests described in appendix B. We close with further discussion in section 7. In particular, an argument based on the second law of thermodynamics predicts that a further new phase must exist with lower free energy must exist in certain regions of parameter space, but we find no natural candidates for this phase.

Dominant Funnels and Droplets for general ESU black hole boundaries
We begin by considering general compact static boundary spacetimes containing a pair of boundary black holes. We suppose here that the black holes lie on opposite sides of the boundary spacetime as shown in figure 1, though it is also interesting to consider phase transitions that occur as one varies the relative separation of the two black holes. We first discuss the relevant notion of the thermodynamic 2nd law in section 2.1 to establish that the most thermodynamically stable solutions are those with lowest ∆F . We then combine general arguments from [13] with observations from studies of the Gregory-Laflamme instability in section 2.2 to motivate the rough form of a general phase diagram for droplets and funnels. Finally, we review analytic results [25] for d = 3 with BTZ boundaries in section 2.3.

The second law for droplets and funnels
Before discussing the various possible phases, note that for general boundary metrics one may expect a variety of phases to exist, for which one will thermodynamically dominate over the others. Recall that we study a system in equilibrium with a heat bath provided by the black hole on the boundary. One thus expects the relevant notion of dominance to be determined by minimizing the free energy F = E − T S defined by the heat bath temperature T ; i.e., where T is fixed by the surface gravity of the boundary black hole. Indeed, for Lorentzian AlAdS spacetimes with such boundaries it was shown in [37] that finite processes cannot decrease the renormalized free energy F defined by using the entropy of the bulk event horizon. The same must hold for our ∆F since it is defined by subtracting off the free energy of a fixed reference solution. This gives a useful form of the 2nd law for our systems. This point is sufficiently important that it is useful to give an expanded discussion over the next few paragraphs. Many readers will find it most clear to consider the relevant thermodynamics from the viewpoint of a dual CFT living on the boundary spacetime. This boundary spacetime contains black holes. CFT entropy can disappear into the black hole, can also be emitted from the black holes, so the CFT outside the black hole is not a closed system and no law of thermodynamics can forbid ∆S from decreasing.
In contrast, if the black holes were dynamical, then we should find the total entropy S total = S CFT + S bndy BHs to be non-decreasing. That is to say, if we made gravity dynamical on the boundary, then there would be a boundary Generalized Second Law (GSL).
However, the boundary metric is not dynamical. We can think of this as taking the limit of a theory with dynamical boundary metric (with some boundary Newton constant G bndy ) and in particular taking G bndy → 0. In this limit, the Bekenstein-Hawking entropy of the boundary black holes diverges and, moreover, the boundary geometry does not change as energy flows in and out of the boundary black holes. The first fact makes the boundary GSL useless in the standard form (non-decrease of S CFT + S bndy BHs ), but the second fact comes to the rescue. Since changes in the boundary geometry are tiny, we can use the 1st law for boundary black holes to write dS bndy BHs = dE bndy BHs /T BHs . But energy conservation would The boundaries between phases are not shown as their locations are unknown. We cannot rule out the possibility that additional phases (such as thin funnels) also dominate in certain regimes, though comparison with known results for Kaluza-Klein black holes [38,39,40,41,42,43,44,45,46]. makes this seem unlikely. The absence of a Gregory-Laflamme instability for d = 3 boundary dimensions suggests that in that case the fat and uniform funnels are connected by a crossover, though in higher dimensions we expect a sharp phase transition.
also guarantee that dE bndy BHs = −dE CFT , so the boundary GSL in fact forbids decreases in the quantity That is to say, it states that the free energy F CFT must be non-increasing. This is just the usual way in which the 1st law allows us to rewrite the closed-system 2nd law as a useful 2nd law for open systems interacting with a heat bath. The above argument was the primary motivation for [37], which showed that the corresponding bulk spacetimes do indeed satisfy such an (open system) version of the 2nd law.

Droplets and Funnels for general boundary black holes
We are now ready to discuss natural phases of droplets and funnels for general boundary black holes. As argued in [13], one expects the radius R and the temperature T of the boundary black holes to be important in determining both allowed and dominant phases. The spatial scale d of the compact boundary will also play a role. Note that when the boundary spacetime is conformally equivalent to a pair of AdS d black holes, this d is also the boundary-AdS scale. The details of the associated transformation will be reviewed in section 3.
As we explain below, the arguments of [13] generalize readily to our current context and suggest the rough phase diagram shown in figure 2. For the moment we consider general independent R, d , T , though for Schwarzschild-AdS or BTZ boundary black holes the dimensionless combinations are all determined by a single parameter (e.g., R/ d ). Let us begin to understand figure 2 by analyzing the limit R/ d → 0 at fixed T d (i.e., along the far left side). Since the boundary black holes are then tiny compared to all other scales, they naturally support bulk horizons that extend only slightly into the bulk. This suggests a droplet phase similar to that of [30] near each droplet for small R/ d . As we will recall in section 2.3, the d = 3 analysis of [25] suggests that there can be more than one droplet solution for given boundary conditions but that droplets whose horizons stay closest to the boundary will dominate. Such solutions are called 'short' droplets in section 2.3 and we discuss only such droplets here.
The droplets in this regime make only minimal impact on the bulk physics. Thus the phase structure along the left side of figure 2 must be otherwise identical to that in the absence of boundary black holes and is dictated by the Hawking-Page transition [47]. At large T d , the dominant phase will contain two droplets as well as a large central AdS-Schwarzschild-like black hole. This phase is dynamically stable because the AdS gravitational potential inhibits attempts to merge the central black hole with either droplet. There will also be dynamically unstable phases involving more bulk black holes, carefully balanced between the central black hole and the droplets. But of greater interest are two dynamically stable sub-dominant phases: a pure droplet phase with no central black hole, and a phase having both droplets and a small central black hole. In the limit R/ d → 0, the transition occurs at the Hawking-Page temperature T HP = (d − 2)/(2π d ). At lower temperatures the pure droplet phase dominates, and at sufficiently low temperatures it is the only phase that exists (see figure 2).
We will now probe larger values of R/ d at high temperature, continuing clockwise around figure 2. As measured by a standard dimensionless Fefferman-Graham coordinate z (with z = 0 on the boundary), we expect the droplet to penetrate a distance z ∼ R/ d at small R/ d . This is a consequence of scale/radius duality (and thus of the bulk symmetries in empty global AdS). But the z-location of the large central black hole's horizon is set by the temperature, with z ∼ T d . So for RT 1 any large central black hole should merge with the droplets to form a funnel. The change in horizon topology requires some sort of phase transition, though the details remain to be investigated. If we are still at R/ d T d then, as measured by natural coordinates on the compactified spacetime, the droplets remain thin relative to the black hole and the resulting funnel will display the large bulge shown in left-most part of this region of the phase diagram. The solution will be well-approximated by the global bulk Schwarzschild-AdS black hole (with ESU boundary) over most of the spacetime, with significant departures only very close to the poles of the ESU. We call this the fat funnel phase, though at small R/ d this need not imply the existence of other funnels; it is possible that all sub-dominant phases are droplets. On the other hand, thinking of a fat funnel as in some sense composed of a more uniform funnel and a large black hole, at large R/ d (with still T d R/ d ) the Hawking-Page transition for pure ESU boundaries suggests the presence of (at least) two further funnel phases (subdominant at large T d ), which we call thin and uniform.
Increasing R/ d at fixed T d will gradually change the shape of the funnel to make the central bulge less pronounced. For large R/ d and small T d , one expects the funnel to reach down a distance ∆z ∼ d /R into the bulk. So the central bulge should disappear completely in the limit R/ d T d . As we will review in section 3 below, an analogy with the Gregory-Laflamme instability and the analysis of [48,49] suggests that for d ≥ 4 there is in fact a sharp phase transition from fat funnels to more uniform funnels, though the lack of a Gregory-Laflamme instability for d = 3 predicts a cross-over in that case. We will identify this phase transition below for d = 4 AdS-Schwarzschild boundary black holes. The uniform funnel phase should then persist as we decrease T d keeping R/ d sufficiently large. At small T d , decreasing R/ d will result in a transition back to the pure droplet phase.
This completes our discussion of the expected phase structure for general R, d , T. Although the dominant phases are clear in various asymptotic regimes, the details of the transitions and many features of possible sub-dominant phases remain to be explored. For example, both droplets and funnels may persist well into regimes where the other phase dominates.
For the rest of this work we therefore restrict attention to a particular 1-parameter family of boundary black hole metrics for each d. Section 2.3 will briefly review the d = 3 results of [25] using BTZ boundaries in order to provide further insight into the droplet phases. This structure seems likely to persist for general black holes with d ≥ 4. Section 3 then considers d ≥ 4, taking the d-dimensional boundary metric to be given by a pair of Schwarzschild-AdS d (SAdS) black holes. We focus in particular on insights from the Gregory-Laflamme instability, which introduce qualitative differences from d = 3. Since in both cases the boundary metrics are labeled by a single free parameter, not all of the above phases need arise.
The relevant curves through our phase diagram are drawn in figure 3. The curves lie largely in regions with intermediate values of parameters where the dominant phase is not yet clear. However, recalling that fat and uniform funnels should be connected by a crossover in d = 3, the left diagram for BTZ boundaries predicts a transition from droplets at small R, T to funnels at large R, T . While additional intermediate phase transitions are possible in principle, they would require a special symmetry that arises in this case to be spontaneously broken (see [25] or the review in section 2.3 below). And since the dropletplus-black-hole phase is beyond the scope of this work, the right diagram for d = 4 SAdS boundaries predicts that we will find fat funnels to dominate at small R and that either fat or uniform funnels will dominate at large R.

BTZ droplets and funnels: a brief review
It was argued in [13] that the most interesting droplet and funnel phases for BTZ boundary metrics could be found analytically due to an at-first-sight surprising SO(2,1) conformal symmetry of the boundary spacetimes. This constitutes a symmetry of the boundary conditions which corresponds to a bulk isometry of any phase in which it is not spontaneously broken. The phases preserving this symmetry can then be mapped via double Wick rotation The boundary metrics studied in [25] contain BTZ black holes while those studied below are Schwarzschild Anti-de Sitter (SAdS). Recalling that fat and uniform funnels should be connected by a cross-over in d = 3, the left diagram predicts a transition from droplets at small R, T to funnels at large R, T . Since the droplet-plus-blackhole phase is beyond the scope of this work, the right diagram predicts that we will find fat funnels to dominate at small R and that either fat or uniform funnels will dominate at large R. Other phases may of course be possible at intermediate values of parameters. To see this symmetry, we begin with the non-rotating 2 BTZ metric [50,51] where we have pulled out an overall factor of r 2 2 3 relative to the usual presentation. Introducing η = R 2 t/ 2 3 , τ = R 2 φ/ 2 3 , and sin θ = R/r this becomes Note that the factor in square brackets is just the θ ≥ 0 half of the static patch of dS 2 × S 1 . The full static patch is then obtained by gluing together two copies of the BTZ metric as shown in figure 1. Including the region behind the BTZ horizons leads to global dS 2 × S 1 , as is clear from the fact that this is the maximal analytic continuation preserving periodicity of τ. It is the SO(2,1) symmetry of global dS 2 that allows analytic control. Wick rotating dS 2 to S 2 , the boundary becomes S 2 × S 1 , which is just the thermal d = 3 Euclidean ESU. Applying the same operations to any bulk solution invariant under this SO(2,1) × U(1) symmetry gives a spherically-symmetric static solution asymptotic to empty (global) AdS 4 . So in the absence of spontaneous symmetry breaking the phases we seek can be obtained by double Wick rotation of empty (thermal) AdS 4 together with the large and small AdS-Schwarzschild black holes.
To understand the relation of these three geometries to funnels and droplets, consider a co-dimension 2 surface in each spacetime at some constant Killing time and constant azimuthal angle. As shown in figure 4, such surfaces terminate on horizons (solid lines) and at the azimuthal rotation axis (dashed lines). Double Wick rotation exchanges horizons with rotation axes, so one obtains spacetimes with BTZ boundaries where the solid line describes a rotation axis and the dotted lines describe horizons. The solution obtained from AdS 4 is a black funnel; it has no rotation axis and a horizon that runs from one side of the ESU 3 boundary to the other. Explicit calculation shows that it is in fact the BTZ black string of [52]. In contrast, double Wick rotation of the black hole spacetimes yields a rotation and two disconnected components to the horizon. These are the desired droplet solutions, which we call short and long based on the distance the horizon penetrates into the bulk.
Since double Wick rotation does not change the Euclidean action, the details of the phase transition are equivalent to those of Hawking-Page [47]. But the mapping of parameters is non-trivial. In particular, the high-temperature behavior of the Hawking-Page transition maps to physics of low-temperature BTZ boundary black holes. The funnel phase exists for all temperatures as a local minimum of free energy. It dominates the canonical ensmble below at T BT Z = 1 4π 3 . In contrast, the droplet phases exist only for . The short droplet phase locally minimizes the free energy and dominates below the transition temperature. The long droplet phase locally maximizes free energy and can be interpreted as mediating the transition between the other two phases. It is natural to expect a similar structure for droplets in higher dimensions, though funnels become more complicated as we discuss below.

Uniform and non-uniform phases of global black funnels
As mentioned above, further insight into the substructure of funnel phases can be obtained by considering the Gregory-Laflamme instability. After providing some general discussion below, section 3.1 proceeds to perturbative calculations that generalize and extend the work of [53,54,26]. Specializing to boundary spacetimes that describe global Schwarzschild-AdS d (or BTZ) black holes, we may use the observation of [25] that the global AdS-Schwarzschild (or BTZ) string solutions provide exact funnel geometries with the desired boundary conditions for all R/ d . The relevant bulk metrics are where d+1 , d are the bulk and boundary AdS length scales. Here z ∈ [0, π] and ds 2 d is a metric for (global) Schwarzschild-AdS d (BTZ): BTZ solution was obtained in [52] as a special case of the AdS C-metric while the d ≥ 4 solutions were studied in [53]. Such solutions may be constructed by first writing AdS d+1 in terms of AdS d slices and then replacing each slice with Schwarzschild AdS d (BTZ) of the same mass M , or equivalently with the same value of R/ d . Indeed (3.1) solves the vacuum bulk Einstein equations with cosmological constant so long as ds 2 d is an Einstein metric with Ricci scalar set in the usual way by 2 d . The construction thus extends to the rotating case, though for simplicity we fix all angular momenta to zero. The BTZ version was used in [18] as a basis for constructing funnels with heat flow.
It is convenient to use (3.1) to define boundary conditions. For each R/ d , we seek bulk metrics with the same leading-order behavior at z → 0. The ESU frame described in section 1 is defined by introducingz = d sinz √ r 2 + 2 d and rewriting (3.1) in terms ofz and ds where this ds 2 d is indeed the ESU metric for M → 0. The two boundary metrics at z = 0, π then combine to form a single boundary metric atz = 0. Since r diverges at the ESU equator, it is natural to replace r by ρ = sign (sinz) r so that the Z 2 symmetry z → −z acts on the boundary as ρ → −ρ. For odd d the metric is manifestly smooth, though for even d it contains M |ρ| d−1 and is only C d−2 .
While we will seek and find smooth bulk metrics in all cases, an interesting effect of the non-smooth boundaries for even d is that the energies of our solutions will diverge in the this ESU frame. As we will discuss in detail in section 5, this divergence is clearly associated with ρ = 0 where the two black hole solutions are patched together and not with the AdS d black hole horizons. In the particular case d = 4 studied below, the boundary metric is still C 2 . While it may at first seem surprising that a quantum field theory should have divergent energy simply because the spacetime on which it lives fails to be C 3 , for d = 4 this is in fact a natural consequence of both power counting and the conformal anomaly. Indeed, it is well known that the conformal anomaly contributes terms to the stress tensor involving two derivatives of the Ricci scalar which can manifestly diverge when the metric fails to be C 3 . Nevertheless, we will see that energy differences remain finite between distinct solutions sharing these boundary conditions. Returning to (3.1), inspection of the metric shows that ∂ z is a conformal Killing field. We therefore refer to this solution as the uniform black funnel. However, for R/ d 1 and d ≥ 4, the spacetime deep in the bulk near z = π/2 approximates that of the Λ = 0 Schwarzschild black string. As noted in [53,54,26], this should result in an instability like that found by Gregory and Laflamme [55]. We thus also expect to find non-uniform black funnels, analogous to the non-uniform black strings of [55,48,49,4], by following the zero-mode from the onset of the instability.
This zero mode fattens the funnel in some places and thins it in others. Two distinct nonuniform solutions may be found by following the zero-mode with either sign. For otherwise translationally invariant strings these solution are equivalent: thinning the string in one place is equivalent to fattening it in another. But our uniform funnels admit only a conformal Killing field. In particular, the Z 2 symmetry z → π − z gives a preferred middle (z = π/2) at which to compare the girth of the three funnels. We therefore refer to the non-uniform phases as fat and thin, based on their sizes relative to the unform funnel at this point.
As we discuss below, following the zero-mode any finite distance will require us to deform the value of R/ d away from the value R onset / d where the zero mode arises. Since our fat and thin funnels are distinct, following the zero mode in one direction should (at least initially) increase R/ d while following it in the other direction should (at least initially) decrease R/ d . One may thus wonder whether fat and thin funnels ever exist for the same value of R/ d . However, the Hawking-Page transition discussion in 1 suggests that both indeed exist at sufficiently small R/ d . This can be the case only if the branch of solutions (fat or thin) that moves from R onset toward larger R eventually turns around at some R turn and returns to small R. Such behavior is natural, as it agrees with that of the Hawking-Page transition where both big and small branches meet at the nucleation temperature , below which no black holes exist. Phases of Kaluza-Klein black holes are also well-known to display similar behavior [39,40,41,42,43,44,45,46].
It is natural to suppose that our funnel phase diagram near R turn resembles that of Hawking-Page [47] near T nucl . Thus the free energy F should decrease with increasing thickness of the funnel and the uniform funnels should dominate. We may also expect the behavior near R onset to resemble that found in [48,49] so that the non-uniform funnels near R onset have F > F uniform for small dimension d but F < F uniform when d is large. In [49] the transition 3 occurs between bulk dimensions d + 1 = 13 and d + 1 = 14. However, the perturbative analysis of uniform SAdS d funnels in section 3.1 below gives F < F uniform already in d = 4, the lowest dimension where non-uniform funnels arise! Finally, we expect the thin funnel phase to merge with some droplet phase in much the same way that thin non-uniform black string phases merge with phases describing Kaluza-Klein black holes in [38,39,40,41,42,43,44,45,46]. This expectation will be realized in section 6.

Perturbations of Uniform Black Funnels
Having established our general expectations above, we now proceed to calculations. Letḡ represent our background metric, and h its infinitesimal perturbation. Then, in the traceless transverse gauge where the metric perturbation satisfies the linearized Einstein equations take the form where Greek indices indicate boundary directions and Here, hatted quantities are computed with respect toĝ, the metric on ds 2 d . We are further interested in metric perturbationsĥ that preserve spherical symmetry and do not depend on time 4 , so thatĥ is given bŷ where the factors of f (r) inĥ tt andĥ rr and the factor of r 2 multiplying the d − 2 sphere, were introduced for later convenience. Remarkably, the gauge conditions (3.6) turn out to be algebraic in a(r) and c(r), which means they can be readily solved with respect to b(r) and b (r): where indicates derivative with respect to r.
We are now ready to determine the final equations resulting from tensor perturbations. First, we input the gauge conditions (3.8) and the ansatz (3.5) into the perturbed Einstein equations (3.4). This will give one second and two third order equations in r and z. This is not a surprise, since the gauge conditions are first order in r. The rr component of ∆ L h is second order in both z and r, and can be shown to solve the remaining two third order equations. Furthermore, as expected, we get two decoupled equations for Y and b, with a separation constant K. These are: where˙indicates derivative with respect to z. The equation for b has appeared, in a different but related context, in [56], where the negative mode of the Euclidean partition function of the Schwarzschild-AdS black hole was studied. There, K was identified as the Euclidean negative mode. The fact that we get the same equation is not a surprise, since we expect the local thermodynamic stability of the Schwarzschild-AdS black hole to be related to the dynamical stability of the corresponding uniform funnel. This plays a central role in the Gubser-Mitra conjecture [57,58]. It turns our that Eq. (3.9a) has a simple analytic solution in terms of Hypergeometric functions of the second kind. Here we choose the solution that is automatically regular at one of the two singular points, z = π: (3.10) The quantization of K comes from demanding a normalizable solution at z = 0, which yields We are thus left to solve Eq. (3.9b). Note that d+1 has completely decoupled from the problem, in particular, f only depends on r, M and d . For the sake of presentation, it is useful to parametrize our solution by the radius of the boundary black holes, normalized to d . This can be easily done by noting that on the horizon, located at r = R, f (R) = 0. This allow us to rewrite M as a function of R/ d . Furthermore, we introduce the following compact coordinate:x that takes values in the unit interval, i.e.x ∈ (0, 1), being 1 at asymptotic infinity and 0 at the horizon. Now we also note that: meaning that if we express f (r) in term of x, it depends on only. We are now ready to study the boundary conditions of b. In order to solve for b, we first need to investigate its boundary conditions. At the horizon, there are two possible and regularity demands C 2 = 0. On the other hand, close to x = 1, we find the following two possible behaviors ρ onset as a function of d: the blue disks represent our exact numerical data, and the dashed red line the analytic expression (3.19), valid at large d.
Normalizability at x → 1, requires A 1 = 0. Note that for p = 0 one might wonder whether A 1 = 0 is an allowable solution as well. However, one can check that, by using Eq. (3.8), this would correspond to a divergent a(r). Now that we have unravelled which boundary conditions we want to consider, we proceed to the actual numerical method we have used. First, we introduce a new function, q(x), that will actually be used in the numerics. This function is defined as (3.17) The equation for b, namely Eq. (3.9b), is secretly a quadratic Stürm-Liouville equation in q(x), with quadratic eigenvalue ρ. It takes the following schematic form where each of L's is a second order differential operator in x, that is ρ independent. Their explicit form can be found in appendix A. We can now use standard methods to solve this linear boundary value problem, see for instance [59]. Numerical results for R onset / d are plotted as a function of d in Fig. 5. As expected, R onset / d is a monotonic function of d. This is exactly what we expect, since for high d it should be easier to render the uniform funnel unstable. Also, we find that there is no instability for d = 3, in contrast to the expectations of [52].
We can repeat mutatis mutandis the 1/d expansion of [60] to find a closed form expression for √ ρ onset at large d

Fig. 5 also plots
√ ρ onset as given by Eq. (3.19) and finds excellent agreement at large d (for instance, for d = 55 we find that the relative error between the exact numerical value and our analytic expression is smaller than 6×10 −2 %). Note also that our onset curve shows that R onset is always smaller than R nucl (where large and small AdS black holes meet), regardless of the number of dimensions. We were also able to go take our perturbative approach beyond first order. The scheme is very similar to that used in [61], except that we cannot use the gauge invariant formalism of [62]. As such we need to start with a gauge choice, which we detail in the next section.
We start assuming that a general metric, sufficiently close toḡ ab , admits an expansion of the form (3.20) At each order in perturbation theory, Einstein equations always take the same form, namelỹ where k represents the quantum numbers of each of these perturbations. For instance, for tensor perturbations, k would be p in Eq. (3.11).The idea is to start at linear order with one of these modes, for instance the tensor perturbation discussed above, and compute T ab . We then decompose T (2) ab as a sum of tensors, vectors or scalars perturbations and compute h (2) ab . At second order one finds that there is a tensor metric perturbation, with the same quantum numbers as the one you started with, that cannot be made regular. It turns out, this singularity can be removed by promoting ρ to be ε dependent: where ρ (0) was determined from the eigenvalue problem detailed above. The singularity at second order can be readily removed by an appropriate choice of ρ (1) . Doing so and computing ∆F = F − F uniform yields the 2nd order perturbative result in terms of the parameter We use this parameter for convenience, but it also describes the rank of the gauge group of the dual super Yang-Mills theory in an AdS/CFT context [1]. As noted above, the d = 4 result (3.24) is a surprise when compared with the analogous study of the original Gregory-Laflamme instability in [49] which found ∆F < 0 only for bulk dimensions d + 2 ≥ 14.
We now turn to constructing and analyzing the corresponding non-perturbative solutions numerically. In section 6 we will present results that fit (3.24) well for T 5 ≈ 0.312617.

Numerical Procedure and Framework
Having set our expectations with analytic arguments, we now turn to numerics to explore details of the phase diagrams. We present new results only for d = 4, though for d = 3 our code also reproduces the analytic results of [25]. In addition, for d = 3 our code also allows for boundary black holes that are not BTZ. We will first introduce the general numerical technique that allowed us to determine the phase diagram of black funnels and droplets. This technique was first introduced in [63], discussed in great detail in [29] and reviewed in [64,59]. The idea is to solve a set of PDE's that are manifestly Elliptic, and whose solutions coincide with solutions of the Einstein equations in a certain gauge.
We deform Einstein's equations: by adding the following new term where ξ a = g cd [Γ a cd (g) −Γ a cd (ḡ)] andΓ(ḡ) is the Levi-Civita connection associated with a reference metricḡ. It is easy to show that any solution to G ab = 0 with ξ = 0 is a solution to G H ab = 0. However, the converse is not necessarily true. In certain circumstances one can show that solutions with ξ = 0, coined Ricci solitons, cannot exist [29]. The line elements discussed in this manuscript are exactly in that class. In particular, in order to ensure that ξ is everywhere zero, we need to ensure that all components of ξ are zero at any asymptotic end, and that the extrinsic curvature at fictitious boundaries, such as horizons, is zero. One can show that these conditions are only sufficient, and that ξ can be made zero for more general boundary conditions [29].
A choice of reference metric is equivalent to a gauge choice, which in the DeTurck formalism is given by a natural generalization of the Harmonic gauge ξ = 0 ⇔ x a = g cdΓa cd (ḡ). As a result, some reference metrics are more amenable to numerics than others. Our criteria for choosing the reference metric is simple: we demand it has the same axis and horizon locations as the metric we seek to find, and satisfies the same Dirichlet boundary conditions as g. We have used a standard pseudospectral collocation approximation on Chebyshev-Gauss-Lobatto points and solved the resulting non-linear algebraic equations using a Newton-Raphson method. This discretization is well know to have exponential convergence, so long as all functions are analytic in their integration domain. As we shall see below, this is not the case for the Einstein DeTurck equation. This issue is particularly relevant when reading asymptotic quantities such as the total energy of a given solution 5 .
Implementing this procedure in particular cases requires a choice of metric ansatz and boundary conditions. We describe these in turn for each of the situations we wish to study. Some readers may wish to skip directly to our diagnostics in section 5 or to section 6 which presents our numerical results.

Ansatz for black funnels in d = 4
The line element we use to describe thin and fat funnels takes the following form where we recall that ρ ≡ R 2 / 2 4 and A, B, F , S 1 and S 2 are function of x and y to be determined in our numerical procedure.
Here, (x, y) take values in the unit square, with x = 0 being the funnel horizon, y = 1 the conformal boundary x = 1 the point infinitely far away from the black funnel, and y = 0 the plane of symmetry that divides the funnel into two equal halves. Finally, for reference metric we take the line element above with A = B = S 1 = S 2 = 1 and F = 0.

Boundary conditions at the horizon x = 0
The metric (4.3) and the associated reference metric are regular at x = 0, if and only if A(0, y) = S 1 (0, y). The best way to understand this is to work in the Euclidean section, settingt = −iτ , and introducing the new coordinate which brings the line element 4.3 to the following form where we have expanded all functions around η = 0. We recognise the first two terms as being flat space, if and only if A(0, y) = S 1 (0, y) and if τ gets a periodicity of 4π √ ρ/g(0) and if F (x, y) = xF (x, y) for a smooth functionF (x, y). This in turn implies that our funnels have a temperature, measured in units of [t] −1 , parametrized by ρ, and given by Furthermore, the extrinsic curvature at η = 0 is zero, which is one of the conditions detailed in [29] for the nonexistence of DeTurck solitons. The remaining functions all have Neumanntype boundary conditions, which can be found by expanding the equations of motion close to x = 0. To wit, we find (4.8)

Boundary conditions at the asymptotic end x = 1
In [29] it was shown that asymptotic ends require ξ a = 0 on all components. If we impose A = B = S 1 = S 2 = 1 and F = 0 that turns out to be the case.

Boundary conditions at the reflection plane y = 0
The boundary conditions that we are interested at y = 0 are those of a reflection plane. Equivalently, we want the extrinsic curvature on the induced hyperslice defined by y = 0 to vanish. This can be easily achieve if we demand As alluded above, these boundary conditions also ensure that no DeTurck solitons exist.

Boundary conditions at the conformal boundary y = 1
Here we shall not only discuss the relevant boundary conditions, but also how to extract the corresponding stress energy tensor. At the conformal boundary we impose A(x, 1) = B(x, 1) = S 1 (x, 1) = S 2 (x, 1) = 1 and F (x, 1) = 0 . This automatically ensures that both ξ x and ξ y are zero at the boundary, and thus that no DeTurck solitons exist in our spacetime. One can also solve the equations off the conformal boundary up to any desired order. These take the following schematic form where all higher order terms depend on the five unknown functions {α(x), β(x), γ(x), δ(x), (x)} and their derivatives along x and Note that one can easily extract δ(x) and (x) by taking two derivatives of A and S 1 with respect to y, respectively, and evaluate them at the conformal boundary. All of these five functions are to be determined by requiring regularity in the bulk. It is not surprising that there are five such integration functions, since our PDE system is second order, and there are a total of five functions to solve for. We thus expect a generic expansion consistent with 10 free functions, half of which should be killed by our choice of boundary conditions, giving a total of five free functions. Note also that we expect ξ µ = 0 to emerge as the unique solution, but it is not a condition that one can see emerging locally by solving the Einstein DeTurck equation off the AdS boundary. One can, however, see what local conditions do come out by imposing ξ µ in the asymptotic expansion. As we shall see, these are related to the conservation of the holographic stress energy momentum tensor.
Requiring ξ µ to be zero order by order in a (1 − y) expansion further demands The first of these conditions ensures that the non-analytic piece that populates our expansion is pure gauge, since it can be reabsorbed via a redefinition of y. The second condition is related to the conservation of the holographic stress energy tensor, which we will extract next.
In order to read off the stress energy tensor, we closely follow the procedure first outlined in [65]. First we change to Fefferman-Graham coordinates and fix the conformal frame. Because we only know our functions numerically, we can only preform this coordinate change asymptotically. Up to O(z 5 ), we find that the relevant coordinate transformation is given byt which induces the following conformal boundary metric If we further define r = R/(1 − χ 2 ), we recover the line element (3.2) with d = 4. This confirms that our line element has as a boundary metric two copies of a Schwarzschild AdS 4 black hole.
The χ−dependent expectation value for the holographic stress energy tensor induced by the above coordinate transformation reads: where Ω i denotes any coordinate on the two sphere. The constant offsets in each of the components corresponds to the stress energy tensor of the uniform funnel, and the fact that is constant has been subject of intense study in the literature [26]. It results solely from the conformal anomaly, and is absent in uniform funnels that live in d + 1 even bulk spacetime dimensions. It is easy to check that the holographic stress energy tensor satisfies , and ∇ µ T µν = 0 (4.17) with the latter condition being enforced via the last equation in Eq. (4.13), and the first reproducing the standard four-dimensional conformal anomaly.

Ansatz for black droplets in d = 4
The metric ansatz for the black droplets takes the following form Alike in the funnels case A, B, F , S 1 and S 2 are function of x and y to be determined in our numerical procedure and (x, y) take values in the unit square. Here x = 1 is the droplets horizon, x = 0 is the plane of symmetry that divides the droplet into two halves, y = 0 is the axis of symmetry and y = 1 is the conformal boundary. For reference metric we take the line element above with A = B = S 1 = S 2 = 1 and F = 0.

Boundary conditions at the horizon x = 1
The boundary conditions at x = 1 are very similar to those of the black funnels at x = 0. In particular, the line element (4.18) and its associate reference metric are only regular at x = 1 if A(1, y) = S 1 (1, y) and F (1, y) = 0. This can be best understood if we again consider a coordinate transformation of the form and expand the metric (4.18) around η = 0. As before, we arrive at a temperature, measured in units of [t] −1 , given by We can also expand the line element (4.18) close to η = 0, to find that the extrinsic curvature at the horizon vanishes, which is one of the conditions detailed in [29] for the nonexistence of DeTurck solitons. For the remaining variables we find: A(1, y) = S 1 (1, y) , F (1, y) = 0 , and ∂A ∂x x=1

Boundary conditions at reflection plane x = 0
The boundary conditions at the reflection plane are perhaps a bit more involved. Suppose that we define your functions for x ≥ 0 but we wish to use the line element (4.18) to define a Z 2 -symmetric metric that extends to x < 0. Under what conditions is this extended metric smooth at x = 0? In general, the requirement can be stated as follows: consider Gaussian normal coordinates around the x = 0 hypersurface with λ the "normal" coordinate normalized to measure proper distance along geodesics that orthogonally intersect x = 0. We set λ = 0 when x = 0. Then the power series expansion of the spacetime metric in these coordinates should contain only even powers of λ, so long as we are away from the boundary.
One introduces Gaussian normal coordinates perturbatively in an expansion in λ. Schematically, they take the following form: (4.23) The {A i , B i } coefficients can be determined by demanding g λλ = 1 and g λY = 0. Note that there are no other cross terms. We can now look at the expansions in λ of the remaining metric components, and ask whether they are even in λ. If we impose the Einstein-DeTurck equations, together with the following boundary conditions at x = 0 that turns out to be the case. We went up to fifth order, and we are convinced this will happen to all orders. Note that that there are some rather complicated cancellations of the odd oder terms in λ which only occur if the equations of motion are used. This of course also implies zero extrinsic curvature at x = 0, and thus that no DeTurck solitons exist.

Boundary conditions at the axis y = 0
These boundary conditions are obtained demanding that the axis is smooth. In order for this to happen we need S 2 (x, 0) = B(x, 0) and F (x, 0) = 0, so that the last two terms in In this section we will not only detail the relevant boundary conditions at y = 1, but also how to extract the holographic stress energy tensor. We choose the following set of boundary conditions These automatically ensure that the leading terms of ξ a = O(1 − y), without the use of the equations of motion, and once more are within the class of boundary conditions studied in [29], for which DeTurck solutions can be ruled out.
One can solve the equations in an expansion off the conformal boundary, to any desired order, and in particular we find: are pure gauge, and the last enforces the conservation of the holographic stress energy tensor. Alike for the black funnels, we now change to Fefferman-Graham coordinates and fix the conformal frame. This can only be done perturbatively off the boundary, since we only know our functions numerically. We want to work in the same conformal frame we did before, since we want to explicitly compare the stress energy tensors of both phase. This is accomplished by the following coordinate transformatioñ was defined in Eq. (4.4). This coordinate transformation induces the following conformal boundary metric If we further define r = R/(1 − χ 2 ), we recover the line element (3.2) with d = 4. Alike the funnels, our droplet line element has as a boundary metric two copies of a Schwarzschild AdS 4 black hole.

Diagnostics
We will shortly compare the free energies of various droplets and funnels. However, since the bulk horizon extends to the boundary at infinity it clear that each of these quantities will be infinite. As such, we need to find a consistent way of regularizing them to make the comparison meaningful. In a nutshell, we will use the uniform funnel as a regulator.
Let us first explain the regularization of the energy. The energy is defined in the usual way via the boundary stress energy tensor T µν [66,67] (see [68,14] for reviews aimed at relativists): where K is the static Killing vector ∂ t , Σ t is a surface of constant t in the AlAdS boundary, with σ and T its induced metric and unit normal, respectively. The factor of two accounts for the fact that we have two copies of Schwarzschild-AdS 4 at the boundary. We can readily integrate in the angular coordinates, leaving only the integral over X.
For concreteness let us evaluate (5.1) on the d = 4 uniform funnel using eqs. (4.16) with δ(χ) = (χ) = 0. We find As foreshadowed in section 3, this is clearly divergent at the 'equator' χ = 1 where the boundary metric fails to be C 3 . This occurs even though we use the standard 'renormalized' boundary stress tensor T µν In order to deal with finite quantities, we will therefore systematically subtract the boundary stress-energy T uni µν ) of the uniform black funnel from both the non-uniform funnels and black droplets and compute We have explicitly checked that Eq. (5.3) gives finite results for both the droplet and funnels phase. However, we stress that this is a very difficult quantity to compute because we must (numerically) cancel this apparent divergent behaviour. For this reason we had to resort to extended precision in our calculations. To ease our numerical scheme, we used a standard Newton-Raphson algorithm to find a solution with double precision, and switch to a Broyden type method when performing extended precision calculations.
The astute reader will notice that in order for ∆E to be finite for the droplet phase, we must have We have checked this behaviour explicitly in our numerics. One can understand the origin of this simple pole by introducing funnel like coordinates near (x, y) ∼ (0, 1) and solving the corresponding equations off the conformal boundary. One might wonder whether this divergence can be removed by an appropriate choice of conformal frame, but one can show it is not possible. One might also wonder how this behaviour is consistent with reflection symmetry around x = 0. Indeed, near this special point it is natural to consider a "radial" variabler ≡ x 2 (2 − x) 2 + (1 − y 2 ) 2 , in terms of which all functions admit a regular Taylor series at least up to orderr 2 . For any y = 1, the reflection symmetry is present and that is why it is explicit in the bulk but not on the boundary. The computation of the entropy has very similar issues: since the horizons intersect the boundary, they are non-compact, and as such have infinite area. An additional complication is that we must find a gauge invariant way to subtract the divergent piece, again using the uniform funnel as a reference. The induced metric, for both funnels and droplets, on the intersection of the horizon with a partial Cauchy slice of constant t reads ds 2 Ht = 2 5 1 − y 2 Q(y) For the uniform funnel we have W = ρ and Q = 1. Clearly, the area diverges near y = 1 for all the phases, as anticipated. So we first introduce a cut-off y i , where i = {D, F, UF} labels droplets, funnels and uniform funnels, respectively. The question is then how to relate the cut offs when we look at the difference in the areas. We do this my demanding that the radius of the S 2 (a gauge invariant quantity) match as the cut off recedes to the boundary. For the black funnel phase, such a procedure demands For the funnel phase, one finds after some work that this procedure yields with the factor of 2 accounting for the two copies. In deducing this expression one has to ensure that the limit y F → 1 − exists and is finite, which one can do using the boundary expansions determined previously.
The droplet phase is more complicated, but the result is the same in spirit. To match the spheres we now have and we wish to look at where we have explicitly checked (by using our asymptotic expansion (4.27) around y = 1) that the remaining integral is finite and can be readily evaluated using our numerical data.
To compute the different in Helmoltz free energies we simply take ∆F = ∆E − T ∆S. This is the appropriate quantity to study when exploring the phase diagram at constant temperature, as we shall do below.
Finally, we will be using isometric embeddings when comparing funnels and droplets. These are specially useful if we want to visualize where the horizon bulges out. In a nutshell, we will embed the spatial cross section of the horizon into four dimensional hyperbolic space. We will foliate four-dimensional hyperbolic space using three-dimensional hyperbolic space: where˜ 3 and˜ 4 are the hyperbolic length scales of the embedding space. One then searches for an embedding of the form (R(y), Y (y)), which gives the following induced metric We can now compare this line element with the metric of a funnel or droplet induced on the intersection of the horizon with a partial Cauchy surface of constant t and read off a nonlinear first order equation forỸ (y). We fix the boundary conditions by demandingỸ (0) = 0 and choose the ratio˜ so thatR(1) = √ ρ 4 . The curve traced by (R(y),Ỹ (y)), as we vary y ∈ (0, 1), is the embedding diagram. Note that in this embedding, a uniform funnel is simply given bỹ R(y) = √ ρ 4 , as expected.
6 Droplets and funnels for d = 4 We now present the results of our numerics. We also encourage the reader to consult appendix B for evidence that our simulations exhibit the expected convergence with increasing numbers of grid points. We begin with the black funnels and then include droplets below. As described in section 3.1, the uniform black funnel becomes unstable for ρ > ρ onset . In this regime, at least perturbatively close to ρ onset , we can identify additional fat and thin branches of black funnels in our numerics. The former branch exists at ρ < ρ onset , while the latter exists at ρ > ρ onset . In fact, as is clear from the perturabtive discussion, the fat and thin funnels are really the same branch of (nonuniform) funnel solutions continued to opposite sides of ρ onset . We can then follow these branches numerically to larger values of |ρ − ρ onset |. As a measure of how fat or thin the funnels may be, we can monitor the minimum radius r min ≡ ρ S 2 (0, 0) of the S 2 along the horizon as a function of √ ρ. As shown in figure 6, the thin branch reaches a turning point at ρ turn ≈ 0.264553 (where r min ≈ 0.0534059) at which ρ then begins to decrease with continued decrease of r min , and beyond which funnels appear not to exist. This non-monotonic behaviour is very reminiscent of the behaviour of the transition between black strings and localised black holes in five and six dimensions, recently observed in [44,45,46]. We shall see below that the droplet phase also shows hints of merging with the thin funnel phase for values of ρ in the vicinity of ρ turn in direct parallel with known results [38,39,40,41,42,43,44,45,46]for Kaluza-Klein black holes. We can also see the predicted outward or inward bulging of the horizon by looking at the isometric embeddings found in our previous section. These are plotted in Fig. 7 for the smallest and largest values of ρ that we managed to probe.
For the droplets, we decided to monitor the proper distance P along the S 2 axis between the two horizons. The results are shown in (8). Our P diverges as ρ → 0 since the horizons recede to the boundary in that limit. The quantity P then generally decreases with ρ, though we again find a turning point √ ρ turn,D ≈ 0.5943 beyond which continued decrease of P requires ρ to decrease. Droplets appear not to exist for ρ > ρ turn,D . When two droplets exist, we may call them short and long in analogy with the d = 3 case [25] reviewed in section 2.3. As in that case, the short droplets have smaller ∆F as may be seen by comparing figures 8 and 15. The behavior is again akin to that of localized black holes in Kaluza-Klein theory [69,40,41,43,70,46,71,72,73]. Other properties of interest include profiles of the boundary stress-energy tensor along the AlAdS boundary in each phase. Due to spherical symmetry and time-reversal symmetry, this tensor has only three non-zero components. Furthermore, the trace anomaly fixes its trace in terms of the boundary curvature, and covariant conservation imposes a second constraint. As a result, there is only one independent component, with any two being determined by the third. We plot T t t in Fig. 9 as a function of the boundary coordinate χ for representative cases of interest again using the parameter N 2 from (3.25). For small enough boundary black holes, the fat funnel phase has a very negative T t t , while the droplet phase becomes increasingly positive away from χ = 1. At the largest value of ρ shown the quantity T t  Figure 8: The proper distance P between the two horizons in the droplet phase as a function of √ ρ. However, the most important quantities for each solution are the total entropy, energy and free energy, i.e. thermodynamic properties. As described in section 2.1, the dominant phase for each boundary metric is the one minimizing the free energy, and thus ∆F . For completeness, we nevertheless first show ∆E and ∆S separately in figure 10, 11, and 12 below. These quantities satisfy a first law in the usual sense: which follows directly from [74,75,76,77]. We have checked that our numerical data satisfies this very stringent relation, and we use it as a numerical check. All solutions presented in this manuscript satisfy this form of the first law, with a relative error of less than 1%. In all thermodynamics plots we represent droplets by orange squares and funnels by blue disks. We also comment that the geometry of a fat funnel near its central bulge remarkably like that a round spherical black hole. To make this quantitative, figure 11 compares the the entropy of a Schwarzschild-AdS 5 black hole with ∆S for a very fat funnel at the same tem- perature. Recalling that fat funnels exist only for small ρ, one sees that at high temperature subtracting the uniform funnel entropy in ∆S will have little effect on the finite part of the entropy so that this comparison will be meaningful up to some constant offset of order 1. Note that the plot for fat funnels is completely consistent with the SAdS 5 result ∆S ∝ T 3 at large T . As noted in the figure captions, the limit of very long droplets appears to merger with the limit of very thin funnels. In particular, the free energies of these families of solution appear to agree in that limit. Such a merger would again be in direct parallel with our understanding of Kaluza-Klein black holes. Indeed, as argued by Kol [38,42], the limiting solution can be expected to take the local form of a double cone. In the AdS 5 context one can find an exact analytic (Lorentz-signature) double-cone metric ds2 = dρ 2 + HdΩ 2 2 + HdΩ 2 2 (6.2) with H = L 2 3 sinh 2 ρ L , dΩ 2 2 the line element on the unit metric S 2 , and dΩ 2 2 the line element on a unit-radius 1+1 de Sitter space. While the boundary metric of this AlAdS spacetime does not match any of the ones we wish to consider, that does not prevent our solutions from approaching (6.2) locally in the bulk near the point where the horizons merge or pinch off. Comparisons of the horizon geometries in our numerical solutions with that of (6.2) are shown in figure 14. These plots are consistent with the idea that the funnels and droplets do indeed merge at a double-cone-like solution.
Finally, we discuss the free energies in detail. Our boundary conditions are specified by  Figure 15: Free energy of all all the competing phases with fixed boundary metric parametrised by ρ: the orange squares represent the droplets and the blue disks the funnels. The right panel shows a close-up of the region close to where non-uniform and uniform funnels merge. The finer scale used in the close up makes the merger less apparent, but given the results of figure 14 we expect that further numerics pushing closer to the presumed merger would fill in the gap between the droplet and funnel free energies. The black disk represents the transition between stable and unstable non-uniform funnels. The close-up also includes the region where droplets are expected to merge with funnels, and in particular the turning point ρ turn for the thin funnels. This turning point was easily visible in figures 10 and 12, showing finite differences S thin − S thinner < 0 and E thin − E thinner < 0 when comparing two branches of thin funnels at the same ρ. But the contributions nearly cancel when comparing free energies, so that the two parts of the thin funnel curve here are separated by much less than the width of the blue disks and the presence of two solutions near ρ turn cannot be seen; i.e., as the funnels become even thinner they are now moving toward smaller ρ and nearly retracing the curve defined by their somewhat-less-thin relatives. However, there is a tiny difference between the free energies of the two thin funnels at any given value of ρ. In particular, numerical computations find F thin − F thinner < 0 as expected. a single parameter ρ, which in turn determines the temperature T (ρ) of the boundary black holes. Since our bulk solutions are in equilibrium with the boundary black holes, T (ρ) is also the temperature of our bulk horizons. As a check on our computations, figure 13 compares numerical results for ∆F for fat and thin funnels near the transition with the 2nd order perturbative result (3.24) and finds strong agreement.
The free energies plotted in figure 15 indicate a rather intricate series of first order phase transitions as follows: For √ ρ 0.236037 the most dominant phase we study is a fat funnel. This is as expected from figure 3 (right) since droplet-plus-black-hole phases are beyond the scope of this work. At intermediate ρ (0.236037 < √ ρ < 0.477205), the most dominant phase is a droplet. Here figure 3 makes no definite prediction. At somewhat larger ρ (0.477205 < √ ρ < √ ρ turn ≈ 0.5143) there is a small region where our most dominant phase is a thin funnel. But as discussed above there appear to be no thin funnels for ρ > ρ turn . Instead, in this final regime the dominant phase we study is the uniform funnel. While figure 3 makes no definite prediction for the dominant phase beyond the small ρ regime, all of these results are certainly compatible with the general analysis of section 2.2.
A notable feature, however, is that the most dominant phase in figure 15 for ρ ρ turn (the thin funnel) does not exist beyond ρ turn . If there were no additional phases beyond those shown in the figure, increasing ρ slightly beyond ρ turn would thus cause ∆F to increase by a finite amount. As described in section 7, this would violate the 2nd law. As a result, there must be additional phases in this regime. See further discussion in section 7

Discussion
In the above work, we constructed and analyzed pure funnel and droplet AlAdS solutions numerically with d = 4 Schwarzschild-AdS boundaries. For each such boundary metric, the solution with lowest free energy F = E − T S should be considered most dominant, as the second law for such solutions [37] will prevent it from decaying to any solution with higher free energy. Here the analytic SAdS black string solution of e.g. [26] is interpreted as a funnel phase, which we call the uniform funnel due to its conformal symmetry. Nonuniform funnels branch off from this uniform funnel at the point where the uniform funnel becomes unstable. This process is closely akin to the Gregory-Laflamme instability of black strings [55] and our phase diagram 15 strongly resembles that of Kaluza-Klein black holes [69,40,41,43,70,46,71,72,73]. Here the funnels play the role of black strings and black droplets play the role of Kaluza-Klein black holes.
Taking R to be the radius of our SAdS boundary black holes, we find so-called 'fat' non-uniform funnels to dominate at small R, droplets to dominate for intermediate values of R, and uniform funnels to dominate at large R; see figure 15. Horizon shapes, energies, entropies, stress-tensor profiles were explored in section 6 and -with one exception described below -fit with general expectations reviewed in sections 2.2 and 3, and also with perturbative funnel computations performed in section 3.1. This supports the intuition developed in section 2 and suggests it will continue to hold for other boundary metrics. In particular, it supports the view of fat funnels as essentially large global SAdS black holes attached to the boundary by small pieces of uniform funnels; see e.g. figure 11. One interesting surprise captured by both numerics and perturbation theory is that even for d = 4 + 1 bulk dimensions, non-unform funnels have lower free energy than uniform funnels even close to the onset of the uniform funnel instability. In contrast, in the perturbative regime non-uniform black strings with Kaluza-Klein boundary conditions dominate over the corresponding uniform (translationally-invariant) string only for bulk dimension d + 1 = 14 or greater [49].
The above-mentioned exception concerns thin funnels, which we find to dominate in a regime between those dominated by droplets and that dominated by uniform funnels. In particular, they have lower free energy than our other phases near the point ρ turn = R 2 turn / 2 4 where the line of thin funnels turns around. However, the second law of thermodynamics requires that there be an additional more-dominant phase in at least the region very close to ρ turn . This may be argued by supposing that we begin with a thin funnel at ρ slightly less that ρ turn and dynamically change the boundary metric to increase ρ. The free energy ∆F should be a continuous function of time, so even at ρ slightly more than ρ turn our ∆F must remain lower than that of the uniform funnel. But as discussed in section 2.1, the 2nd law [37] for our systems prohibits ∆F from increasing when the boundary metric is static.
Our system is thus forbidden from evolving to the uniform funnel 6 . On the other hand, the system is clearly dissipative, so it should settle down to a stationary solution of lower free energy.
In considering possible candidates for this new phase, one should note that the funnel and droplet phases we study appear to merge in this regime at a ∆F greater than that of the thin funnels. Indeed, as shown in figure 14, our numerics is consistent with the idea that the merger is a double-cone-like transition directly analogous to that described in [38,42,46,71] for Kaluza-Klein black holes. Thus even if further turning points are found, the associated phases are expected to have higher free energy. Since the only remaining phase in figure 2 is the droplet-plus-black-hole phase, one may thus expect that this is the relevant new phase new ρ turn .
However, further study of the droplets near ρ turn suggests that this is also unlikely. To understand this point, consider adding a small black hole to a droplet solution. A sufficiently small black hole will act like a test particle and follow a geodesic, and by symmetry the static worldline at the origin is a geodesic in all droplet solutions. But the lowest free energy must be dynamically stable, and the central geodesic is unstable near ρ turn . This is shown by the results in figure 16, which plots the redshift γ(0) = −g tt (0) at the origin as well as the height γ(b) of any gravitational barrier that would stabilize the central geodesic with respect to perturbations along the axis of rotational symmetry. Here we note that since √ −g tt vanishes at horizons, it is bounded along this axis in a droplet phase. We may thus define b to be the point maximizing √ −g tt and compute the associated barrier height γ(b) = √ −g tt (b). The issue is then that near √ ρ turn ≈ 0.5143 figure 16 finds γ(0) = γ(b) (i.e., that the gravitational potential along the axis is maximized at the origin), so the central for further exploration of static phases or, even better, for dynamical real-time evolution of the system with time-dependent boundary conditions that increase ρ across the threshold at ρ turn . Furthermore, while our numerics show the expected convergence (see appendix A) and also the expected behavior near the funnel/droplet merger ( figure 14), it is difficult to exclude the suggestion that some subtle systematic effect might have led to small errors in our free energy computations that, when resolved will remove the need for this new phase. It would thus also be useful for our results to be reproduced independently, perhaps using other computational frameworks. Additional interesting directions to explore include the construction of droplet-plus-black-hole phases for both d = 4 with SAdS boundaries and d = 3 with BTZ boundaries, as in the latter case they would necessarily break the SO(2, 1) symmetry (see section 2.3) used to analytically construct and analyze pure droplets and funnels in [25]. A The differential operators L (i) , for i ∈ {0, 1, 2}: and In both cases we used √ ρ = 0.312.

B Convergence tests
Of all the quantities we computed in the main text, the energy and entropy variations with respect to the uniform funnel are the most difficult ones to extract. Note that once we know these two quantities, we can easily compute the corresponding free energy variation via standard thermodynamic relations. The entropy variation turns out to be very accurate to determine, and its numerical error is always below the 10 −8 %. So we focus on the energy extraction. We will also study convergence keeping the number of points in the x and y directions equal. That is to say, we set n x = n y = n. For each value of n we determine an energy variation ∆E n . To study its convergence we monitor the quantity Our results indicate a power law convergence for the energy, which is not surprising given the irrational decays we found in Eqs. (4.11) and Eqs. (4.27). For reference, we plot in Fig. 17 a convergence plot for the droplet (left panel) and funnel (right panel) for √ ρ = 0.312. For the droplets we find ∆E n ∝ n −4.7 (8) while for the funnels ∆E n ∝ n −4.9(0) . Other values of ρ exhibit similar convergence properties. In all plots in our manuscript, we used n = 250. We have also numerically checked the first law of thermodynamics shown in Eq. (6.1), and find that it is satisfied to 10 −3 % level.