To go or not to go with the flow: Hawking radiation at strong coupling

We construct the gravitational dual of a one-parameter class of states of strongly coupled $SU(N)$ $\mathcal{N}=4$ SYM at infinite $N$ and asymptotic temperature $T_{\infty}$, on a fixed Schwarzschild black hole background with temperature $T_{\mathrm{BH}}$. The resulting bulk geometry is of the flowing type and allow us to measure Hawking radiation at strong coupling. The outgoing Hawking flux is a function of the dimensionless ratio $\tau\equiv T_{\infty}/T_{\mathrm{BH}}$ and appears to be non-monotonic in $\tau$. At present, we have no field theory understanding for this behaviour.

extending from the boundary black hole into the bulk and out to an asymptotic region. In contrast, the horizon in a Black Droplet solution has two disconnected parts. One extends from and surrounds the boundary black hole (the "droplet"), and the other is a deformed planar horizon which is not connected to the boundary at all.
Given the fixed black hole geometry on the boundary with temperature T BH , the bulk solutions are characterized by two free parameters: T ∞ , the temperature of the bulk horizon in the asymptotic region, and T H , the temperature of the bulk horizon where it meets the boundary black hole. If T H = T BH then the Euclidean boundary geometry exhibits a conical singularity and the stress tensor diverges at the horizon. These solutions are important (the Boulware vacuum state is described by such a solution [14]) but we will not consider them in this paper. From this point on we shall always take T H = T BH .
Given our remaining freedom in fixing the parameters, there are two special cases of particular interest. The first is the choice of parameters T ∞ = T BH . This gives an equilibrium solution and corresponds to the Hartle-Hawking state of the field theory. Physically, a thermal state of the field theory is in equilibrium with a Schwarzschild black hole which emits Hawking radiation at the same temperature. The second case of particular interest is when T ∞ = 0. This is an out of equilibrium solution corresponding to the so called Unruh state of the field theory. Physically, the field theory approaches its natural vacuum state in the asymptotically flat region (where a natural choice of vacuum state exists) but we now expect to see an outgoing energy flux coming from the Hawking radiation emitted by the higher temperature black hole. This is a good approximation to the state of the field theory that would be obtained after gravitational collapse, but before the resulting black hole has had time to evaporate.
Bulk duals corresponding to the Hartle-Hawking state and to the Unruh state have been constructed previously. For the Hartle-Hawking T ∞ = T BH choice of parameters only black funnel solutions exist [15,16]. These solutions can be said to describe Hawking radiation in that they contain a black hole in equilibrium with a plasma, but there is no net energy flux. For the Unruh T ∞ = 0 choice of parameters, the only bulk solutions constructed so far have been of the droplet type [17]. Surprisingly, this means that they also exhibit no energy flux, at least at the leading order in N in which a classical description of the bulk spacetime is valid. It is expected that Hawking radiation will be present at next leading order in N, but the goal of observing something which can be interpreted as Hawking radiation in a classical solution to the Einstein equation has not yet been realised in the case of the Unruh state.
In this paper we construct for the first time numerical solutions of the black funnel type, containing a Schwarzschild black hole at the boundary, in which T ∞ = T BH . Such states have been considered previously in the weak coupling context by Frolov and Page in [18]. We will study the properties of the system as we vary the ratio T ∞ /T BH . We will start with T ∞ /T BH > 1 and decrease it passing through T ∞ /T BH = 1 into the region T ∞ /T BH < 1. While we are unable to take T ∞ all the way to zero, we can lower it significantly, and observe an energy flux corresponding to outgoing Hawking radiation from the black hole. This also means that our bulk horizon is not Killing, where we evade the zeroth law of black hole mechanics [19][20][21] due to the fact that the horizon is non-compact. Other "flowing funnel" solutions containing two boundary black holes have been obtained previously in the context of AdS boundary geometries [22,23] (see also [24][25][26][27][28][29][30][31] for closely related geometries).
The outline of the paper is as follows. Section 2 explains how to formulate the construction of flowing funnels with arbitrary T ∞ /T BH . In section 3 we detail how to compute some of the horizon properties of flowing geometries, such as its expansion and shear. Section 4 shows how one can extract the holographic stress energy tensor from the numerical solutions, and in section 5 we present our numerical results. We close with some final discussions in section 6.

Constructing the holographic dual
We will work in five bulk spacetime dimensions, and thus have a four-dimensional boundary spacetime. We expect that appropriate generalisations exist in higher dimensions. The solution we seek to construct necessarily exhibits temperature gradients across the horizon, the horizon generators cannot be associated with the integral curves of a Killing field [23,24], i.e. the horizon is not Killing. So we expect our solution to follow under the class of flowing geometries constructed in [24][25][26][27][28][29][30][31].
We begin by reviewing a coordinate system that is well adapted to the construction of black funnels [15,23]. In essence, a black funnel contains a single component horizon which connects the boundary black hole to an asymptotic region that is located infinitely far away from the boundary black hole (see left panel of Fig. 1). In this asymptotic region, the geometry should approach that of a standard five-dimensional planar black hole, which we denote on the left panel of Fig. 1 by P. Thus, black funnels admit a natural triangular integration domain with three boundaries: a horizon, the planar black hole metric, and the conformal boundary. Working with triangular integration domains can be tricky, so in [15,23] a new coordinate system was introduced, such that the point where the bulk future horizon H + meets the boundary, was blown up into a line. By a careful inspection of the region where the bulk horizon meets the boundary horizon, one can show that the metric approaches that of a hyperbolic black hole. In fact this has to be the case for any horizon that meets the boundary, the reason being close to the conformal boundary the geometry is hyperbolic, up to subleading terms that are related to the holographic stress energy tensor. Moreover, the most general cohomogeneity one line element that is static and manifestly exhibits hyperbolic symmetry for each constant time hyperslice is that of a hyperbolic black hole. Each of these hyperbolic black holes is then determined by a choice of boundary black hole temperature, and bulk horizon temperature at the hyperbolic black hole horizon. We are interested in the situation where the boundary black hole has the same temperature as the bulk horizon temperature at the hyperbolic black hole horizon. The only such black hole is the zero energy hyperbolic black hole (which is isometric to pure AdS), and is given by where here henceforth dΩ 2 2 will denote the unit radius round metric on S 2 . We want to find coordinates (x, y) where at x = 0 we approach the zero energy hyperbolic black hole with the same temperature as the boundary black hole T BH , at x = 1 we want to impose that we approach a Schwarzschild black brane with some temperature T ∞ and finally we want y = 0 to denote the conformal boundary. Ideally the horizon would be an hypersurface with y = 1, which was achieved in [23] by a careful choice of gauge. In this work we decided to follow more closely the approach devised in [24], but with some important differences which we will detail below.
Finally, at the conformal boundary we want to consider a geometry which is conformal to Schwarzschild, since this is the spacetime we want our field theory to live on. This means that close to y = 0 the solution has to approach a Schwarschild black string with temperature where we defined the Schwarzschild radius r 0 = 2M . Note that the temperature of the black string matches that of the Hawking temperature of a Schwarzschild black hole. The general idea in both [23] and [24] was to use the so called DeTurck method, which was first introduced in [32], and reviewed in detail in [33,34]. Recall we seek to construct solutions of the five-dimensional Einstein equation where L is the AdS length scale. We will take latin indices to run over bulk spacetime dimensions, and greek indices to run over boundary spacetime dimensions. In order to apply the DeTurck method, we modify Eq. (2.4) and consider instead where ξ a ≡ g bc [Γ a bc (g) − Γ a bc (ḡ)], withḡ being a reference metric which should obey to the same Dirichlet boundary conditions as the physical metric g we wish to find, Γ a bc (g) is the Christoffel connection associated with a metric g. Of course, solutions of Eq. (2.5) will only coincide with solutions of Eq. (2.4) if ∇ (a ξ b) = 0. Under certain special circumstances one can show that solutions with ξ a = 0 cannot exist (see for instance [17,35]), and as such solutions of Eq. (2.5) will necessarily coincide with solutions of Eq. (2.4). However, the solutions we seek do not satisfy the conditions of such theorems, so we need to check a posteriori if ξ approaches 0 in the continuum limit.
There are many other circumstances where one cannot show that ξ is necessarily zero on solutions of (2.5) (see for instance [36][37][38][39][40]). However, in many of these cases, one can show that the resulting system of equations is Elliptic, and as such one can trust local uniqueness of solutions to determine whether ξ vanishes or not. The situation here is, however, more delicate. One can show that the system of equations we wish to solve does not appear Elliptic, instead the system appears to be of the mixed Elliptic-Hyperbolic type. This in turn means that we cannot use local uniqueness to distinguish between solutions with ξ = 0 and solutions with ξ = 0. We note however that this problem has an alternative formulation [41] which does appear to have an Elliptic character.
We will now introduce our line element, and show that it reproduces all of the Dirichlet conditions detailed above. Our line element reads For the reference metric, we will use the line element above with Finally, we will also choose H( We shall see below that will control the temperature of the black brane infinitely far way from the hole, and thus fix the temperature of the field theory reservoir. We now discuss the thorny issue of boundary conditions. At y = 0 we demand q 1 = q 3 = q 4 = 1, q 2 = H(x) and q i = 0 for i ∈ {5, 6, 7}. The metric then reads (2.7) We now take the following coordinate transformations The line element above is that of a Schwarzschild black string in AdS (2.2). That is to say, the boundary metric is conformal to that of a Schwarzschild black hole. This is precisely what we want if the field theory is to live on a fixed Schwarzschild black hole background. At x = 0, we demand q 1 = 1 − y 2 , q 2 = q 3 = q 4 = 1 and q i = 0 for i ∈ {5, 6, 7}. This transforms (2.6) to We now define the following coordinate transformations which brings Eq. (2.10) to coincide with the large η limit of Eq. (2.1). We also recall that such black hole has temperature measure by the time defined in Eq. (2.8a). Finally, we come to the more delicate boundary located at x = 1. At this boundary we want the impose the metric of a planar Schwarzschild black brane with a different temperature than the line element (2.10). At x = 1, we demand q 1 = 1 − y 4 , q 2 = H(1) and q 3 = q 4 = 1 and q i = 0 for i ∈ {5, 6, 7}, which brings the line element (2.6) to Now take which transforms the metric given in Eq. (2.12) to (2.14) This metric describes a five-dimensional planar black hole with temperature  15) we see that there is a gradient of temperatures between the two horizons unless, H(1) = 2, i.e. = 1. We will be interested in situations where = 1, and in particular One might wonder whether we are still missing a boundary condition in the interior. However, that is not the case. One can understand this in the following manner: the solution we are seeking is regular in ingoing coordinates, but badly singular in outgoing coordinates. This is enough to show that there are two possible solutions in the interior: one that blows up and one that does not. By working with a Chebyshev-Lobatto grid we automatically assume enough smoothness to guarantee that we only capture the smooth solution. The price we pay of not imposing a boundary condition in the interior is that we do not know a priori where the horizon is. All we know is that it must join y = 1 at both x = 0 and x = 1, due to the Dirichlet boundary conditions detailed above.

The horizon and its properties
Horizons in general relativity are null-hypersurfaces, and can thus be written as the zero level sets of a functionĥ, which we choose to take the simple formĥ(x, y) ≡ y − P (x) [24] g ab ∇ aĥ ∇ bĥ = 0 , (3.1) i.e. dĥ is null. Eq. (3.1) yields and ODE for P , which one can readily solve. Furthermore, we know that P (0) = P (1) = 1, either of which we can use as boundary conditions. Since the horizons we seek to construct are not Killing horizons, they will have unusual properties. For example, the expansion Θ and shear σ of the horizon generators will be non-vanishing. Note, however, that since the horizon is a null-hypersurface, the rotation of the horizon generators can always be chosen to vanish even for non-Killing horizons. In principle, computing Θ and σ can be a relatively daunting task. However, in [23] a series of tricks where used to determine Θ and σ, as well as a choice of affine parameter λ for the horizon generators. Since the properties of such quantities will play an important test of our numerics, we will review (and generalise slightly) the construction detailed in [23].
The class of spacetimes we seek to find have a time translation symmetry and an S 2 symmetry. We shall use coordinates {v, x, y, θ, φ}, where ∂/∂v is the Killing vector field associated with time translations, φ and θ parametrise the S 2 , and y = P (x) on H + , the future event horizon. For the sake of notation, let us define the shorthand notation ∂ I , so that ∂ v ≡ ∂/∂v, ∂ θ ≡ ∂/∂θ and ∂ φ ≡ ∂/∂φ, and thus I = v, θ, φ.
H + is 4-dimensional, with a 3−dimensional space of generators. We want to fix an x coordinate on H + , and then uniquely label all horizon generators using the coordinates (v, θ, φ) of their intersection with this surface. This should be possible provided that H + is not Killing. We choose one horizon generator and pick an affine parameter λ. λ can then be extended to a scalar function on H + by requiring it to be independent of v, θ and φ, so that λ = λ(x). By symmetry, λ will serve as an affine parameter for each geodesic.
Let k a be the tangent vector to the horizon generators with affine parameter λ. Since ∂ I are all parallel to H + we have k ⊥ ∂ I , for all I. ∂ I are deviations vectors for the geodesic congruence, since k and ∂ I commute. Therefore we have In the coordinates {v, x, y, θ, φ}, we can identify h IJ with the IJth component of the metric tensor through The following is then true where B IJ denotes the IJth component of B.
To obtain the expansion and shear of the congruence, it is necessary to work in the 3 = 5 − 2 dimensional vector spaceV obtained by first restricting to vectors orthogonal to k, and then quotienting by an equivalence relation where two vectors are equivalent if they differ by a multiple of k. Tensors in the spacetime give rise to natural tensors inV if they obey the property that contracting any one index with k a or k a and the remainder with vectors or dual vectors having natural realisations inV , gives zero. B, g and ∂ I all have this property. Furthermore, since k and all ∂ I are linearly independent, we have that all∂ I are linearly independent onV .
We can thus pick {∂ I } as a basis forV . For any tensor T ab naturally giving rise to a tensor inV , i.e.
we can identifyT inV by reading off the components {v, θ, φ} of T . This is the case for the tensors g and B.
The expansion and shear can now be determined using the usual formulae in terms of the quantitiesĝ It remains to find λ as a function of x, but this can be achieved by solving Raychaudhuris equation written using the above expressions forĝ andB [23]. This gives an equation for λ as a function of x on H + : It is then possible to obtainB asB To compute the expansion and shear, we simply look at the trace and trace-free parts of B, namely

Extracting the Holographic-Stress energy tensor
One of the most interesting quantities to extract from these solutions is the associated holographic stress energy tensor FP|T µ ν |FP , where FP stands for Frolov-Page states [18]. We follow closely [42,43] whose starting point is to cast our solutions into Fefferman-Graham coordinates [44][45][46][47]. Fortunately, we only need to perform this coordinate transformation asymptotically.
We start by determining the behaviour of all qĨ , withĨ ∈ {1, . . . , 7}, close to y = 0 by solving Eq. (2.5) asymptotically. A careful analysis reveals the following asymptotic expansion 2 (x), q 3 (x), q 5 (x), q 5 (x), q 7 (x)} not being fixed by any local analysis of the Einstein-DeTurck equation, i.e. these coefficients correspond to data that can only be determined once regularity deep in the bulk is imposed and the corresponding equations of motion are solved. The remaining coefficients are all determined as functions of {q 5 (x), q 7 (x), H(x)} and their derivatives along x. Imposing ξ = 0 asymptotically imposes further constraints, which in turn imply the tracelessness and transversality of the holographic stress energy tensor. In addition, it imposes a linear relation between q 1 (x) and q 7 (x), which can be used to show that the terms proportional to y 2+2 √ 3 are pure gauge. Note, however, that in the DeTurck gauge they are present, and change the expected convergence of the spectral collocation methods we used to solve the Einstein-DeTurck equation from exponential to power law.
To determine the holographic stress energy tensor, we define the following asymptotic change of coordinates where V is the boundary Eddington-Finkelstein coordinate and r is the standard Schwarzschild tortoise coordinate defined as with w 1 ∈ (0, 1) a constant. The coefficients α i , β i , γ i , α, β, γ, α, β and γ are then determined by changing to Fefferman-Graham coordinates, where the five-dimensional line element takes a particularly simple form The metric components g µν (x ρ , z) themselves admit an expansion around z = 0 of the form We choose our conformal boundary metric, i.e.ḡ(x ρ ), to be Schwarzchild Since our choice of boundary metric is Ricci flat, one can show that C = 0 [43]. Indeed, all the log terms in Eq. (4.1) and Eq. (4.2) conspire to give a net C = 0. Furthermore, as we have commented above, the terms proportional to z 2+2 where A ≡ḡ µν A µν . Since the boundary metric is Ricci flat, there are no ambiguities in defining FP|T ρ |FP . Finally, following the standard AdS/CFT dictionary [5] we identify One quantity of interest to compute, once the holographic stress energy tensor is determined, is the heat flux Φ. This quantity is interpreted here as the energy flux integrated over a two-sphere of constant r, that is to say where ξ µ = (∂/∂V ) µ , √ −h is the volume element induced on a constant r slice and n ν an outward unit normal to S 2 r . By virtue of the conservation of the stress energy tensor and since ξ µ is a Killing vector, Φ is independent of r. With our conventions, an outgoing flux of radiation corresponds to Φ > 0 and an ingoing flux to Φ < 0. We expect Φ to be positive if T ∞ < T BH , i.e. > 1, to vanish for the Hartle-Hawking state = 1, and to be negative for < 1. In the subsequent section we will confirm this behaviour numerically.

Results
We start by studying the geometry of the future horizon H + . To help us gain intuition, we construct isometric embeddings of spatial cross sections of H + into hyperbolic space. These are specially useful if we want to visualize where the horizon bulges out. Since spatial cross-sections of H + are three-dimensional, we seek to construct embeddings into H 4 . We foliate four-dimensional hyperbolic space using three-dimensional flat space (i.e. Poincaré slicing of Euclidean AdS) One then searches for an embedding of the form (Z(x), R(x)), that is to say the pull-back of Eq. (5.1) to a parametrised surface (Z(x), R(x)), which gives the following induced metric We can now compare this line element with the pull back of the metric for the flowing funnel, i.e. Eq. (2.6), induced on the intersection of the the future event horizon H + , identified via solving Eq. (3.1), with a partial Cauchy surface of constant v and read off a nonlinear first order equation for Z(x). We fix the boundary conditions by demanding Z(1) = 1. In making the identification between line elements we also setL = L. The advantage of this embedding is that a black string, see Eq. (2.2), appears as a line of constant R, whereas a five-dimensional planar black hole, see Eq. (2.14), appears as a line of constant Z. Thus, a black funnel should naturally interpolate between these two curves.
In addition, we shall also be interested in constructing isometric embeddings of the ergosurfaces of our solutions. These surfaces are defined as surfaces for which the norm ∂/∂v 2 becomes spacelike. Inside the ergosurface, that is to say in the ergoregion, the character of Eq. (2.5) changes from elliptic to hyperbolic. This is the reason why the existence of these surfaces is at the core of the difficulties in trying to prove that our numerical method ensures the absence of DeTurck solitons. We shall follow the same strategy and also embed the ergosurfaces into four-dimensional hyperbolic space.
In Fig. 2 we can observe the embedding of a partial Cauchy surface of constant v with H + (blue disks) and of the ergosurface (orange squares). The ergoregion remains very narrow even when T ∞ /T BH = 2.5. The same is true in the opposite limit, i.e. T ∞ /T BH ∼ 0.6. As expected, the ergoregion shrinks down to zero size at the boundary (y = 0) and at the planar black hole asymptotic region (x = 1). The existence of the bulk ergoregion also raises questions about the stability of the solutions we have just found, specially in light of [48]. However, we note that in [48] it was essential that the spatial cross section of the boundary metric was (metrically) a round sphere. A comprehensive analysis of the linear-mode stability of our solutions is outside the scope of this manuscript, though we plan to return to it in the near future. Perhaps the most important quantity to extract is the stress energy tensor and its associated Hawking flux Φ defined in Eq. (4.9). Within our symmetry class, the holographic stress energy tensor has four non-zero independent components. In addition, it should be traceless and covariantly conserved, which gives three constraints amongst these four components. Thus, the full stress energy tensor is determined by, say, Using conservation of the holographic stress energy tensor, it is also simple to show that if FP|T V w |FP = 0 on H + , and the holographic stress energy tensor is smooth on H + , then T V V | H + = 0. This is what we expect to happen when T ∞ /T BH = 1.
In Fig. 3 we plot T V V for several values of T ∞ /T BH . On the left panel of Fig. 3 blue disks, orange squares and green diamonds correspond to T ∞ /T BH = 0.557103 , 0.844773 , 1, respectively. The right panel of Fig. 3 shows T V V | H + and, as anticipated, it is non-zero except when T ∞ /T BH = 1 (the Hartle-Hawking state). For the Hartle-Hawking state, we recover the results of [14,15]. We now turn our attention to the Hawking flux Φ defined in Eq. (4.9), which we plot in Fig. 4. The right panel shows a zoom of the left panel in the region T ∞ /T BH ≤ 1. As expected Φ < 0 for T ∞ /T BH > 1, meaning that there is an influx of radiation into the black hole. The temperature of the heat bath provides an energy reservoir that sources the geometry. For the Hartle-Hawking state, T ∞ /T BH = 1, there is no net flux of Hawking particles, so that Φ = 0. This is marked as a black disk in Fig. 4. Finally, we see a somehow surprising result when T ∞ /T BH < 1. This is indeed when we expect the black hole to radiate Hawking quanta, and as such Φ > 0. However, Φ does not seem to be monotonic in T ∞ /T BH . Instead, it initially grows with decreasing T ∞ /T BH , but it reaches a maximum at T ∞ /T BH ≡ τ max ≈ 0.698282. This value is marked in Fig. 4 as a vertical dashed black line. At the moment, we have no field theory understanding for this non-monotonic behaviour of Φ. Using the results of [49], we have some preliminary results indicating that the flowing geometry is Gregory-Laflamme unstable [50] in the region T ∞ /T BH < τ max . It is natural to expect that the endpoint of such instability is one of the black droplets found in [16]. Such evolution will necessarily involve a violation of the weak cosmic censorship conjecture [51], alike the non-linear evolution of the Gregory-Laflamme instability [52]. This is not the first time that violations of weak cosmic censorship conjecture play a role in the AdS context, see for instance [53][54][55][56][57], although it seems one of the easiest scenarios to actually sort out the evolution numerically using the methods outlined in [58,59].
We now discuss the properties of this novel class of horizons. For simplicity, we restrict our discussion to the region of moduli space T ∞ /T BH < 1. We begin by studying the behaviour of h IJ itself, and in particular its determinant h ≡ deth IJ . We plot this quantity on the left panel of Fig. 5 for T ∞ /T BH ≈ 0.557103. h is clearly a monotonically increasing function of x, showing that x increases towards the future along the future horizon, and thus that the past horizon lies at x = 0. We note that in this region the temperature of the black hole is higher than that of the heat reservoir. To understand the direction of the flow, we can also determine the horizon coordinate velocity Ω(x) along H + . To do so, we can look at the pull-back of our line element (2.6) to the future horizon H + (see Eq. (3.1)) where the last equality follows from the fact that H + is null,qĨ = qĨ (x, P (x)) and denotes differentiation with respect to x. The last equality in Eq. (5.3) implicitly defines the coordinate velocity Ω(x). We have plotted Ω(x) for T ∞ /T BH ≈ 0.557103 on the right panel of Fig. (5), and we find Ω(x) ≥ 0. Indeed, we find that Ω(x) ≥ 0 so long as T ∞ /T BH ≤ 1, and negative otherwise. Both the sign of Ω and the monotonicity properties of h, indicate that the past horizon is located at the black hole horizon (the hotter reservoir). This is to contrast with the coordinate choice made in [22,23], in which the cooler horizon appears to be closer to the past horizon H − . Next, we look at the behaviour of the horizon expansion Θ as a function of λ. According to Raychaudhuri's equation and the definition of Θ, it better be that Θ > 0, dΘ/dλ < 0 and that Θ approaches zero for large λ. On the left panel of Fig. 6 we plot the affine parameter λ as a function of x on the future event horizon. Recall that x is one of the coordinates introduced in Eq. (2.6). As expected, λ is a monotonically increasing function of x, with x = 0 locating the past horizon. On the right panel, we show the expansion Θ as a function of the affine parameter λ, and we confirm that dΘ/dλ < 0. In order to show that this is the case using Raychaudhuri's equation, one needs to use the equations of motion. We thus see this as the most solid test of our numerical procedures. On the right panel of Fig. 6 we also show, as a dashed red line, a linear fit in the (log 10 λ, log 10 Θ) plane, which yields log 10 Θ = 0.00195019 − log 10 λ. The fit seems to work for all ranges of x, so there is a clear indication that Θ is diverging as λ −1 at the past horizon (located at λ = 0), signalling the presence of a caustic. This caustic, in turn, gives rise to a curvature singularity. As noted in [23], the easiest way to see this is to note that for any Killing field, the Ricci identify for vectors implies that In particular, this will be true for any Killing vector on the two-sphere, say K = ∂/∂φ. It is easy to see that for any such Killing vector K 2 diverges at x = 0, i.e. when λ = 0. But since K obeys the second order partial differential equation (5.4), it can only diverge at finite affine parameter λ = 0 if R abcd diverges at λ = 0 in all orthonormal frames. Finally, we end with a comment regarding the behaviour of our solutions close to the minimal temperature ratio T ∞ /T BH = 0.557103 that we have reached. In Fig. 7 Figure 6. Left: Logarithmic plot of the affine parameter λ as a function of x on H + . Right: log 10 − log 10 plot of the expansion Θ as a function of λ. The dashed red line indicates a best fit consistent with log 10 Θ = 0.00195019 − log 10 λ. Both panels were generated with T ∞ /T BH ≈ 0.729262.
the behaviour of C 2 ≡ max C C abcd C abcd , as a function of T ∞ /T BH , where C stands for domain of outer communications and C is the five-dimensional Weyl tensor. Throughout moduli space, we find no evidence for a divergent behaviour. We suspect the reason why we cannot reach lower temperatures is purely numerical, and is related to the existence of large gradients inside the future event horizon.

Discussion
We constructed the holographic dual of the Frolov-Page states [18] for N = 4 SYM with gauge group SU (N ) and asymptotic temperature T ∞ , at large 't Hooft coupling and infinite N on a Schwarzschild black hole background with temperature T BH . When T BH = T ∞ , the Frolov-Page states reduce to the Hartle-Hawking state, and when T ∞ → 0 one recovers the Unruh state. The resulting bulk geometries are of the flowing type, and thus similar in nature to the ones first uncovered in [23,24]. They possess spherical symmetry and a stationary Killing field ∂/∂v as well as an horizon H + whose spatial cross sections are non-compact. However, the horizon is not a Killing horizon and in particular it is not generated by any linear combination of Killing fields. This is not in contradiction with the standard rigidity theorems [19][20][21], since the latter assume horizons with compactly generated spatial cross sections. When T BH = T ∞ we find a net flux of Hawking radition, being outgoing whenever T BH > T ∞ , and ingoing otherwise. For the Hartle-Hawking state, corresponding to T BH = T ∞ , the flux vanishes identically, as expected, and the results of [15] are recovered. These novel flowing horizons have unfamiliar properties, such as a non-vanishing expansion Θ. We have studied how Θ varies as a function of T BH /T ∞ , and find it is positive and extends to the future along each null generator. The horizon generators extend to infinite affine parameter λ in the far future, but reach a caustic (at finite affine parameter λ. We have computed C abcd C abcd and it seems to remain bounded in the domain of outer communications, indicating that the caustic located at λ = 0 is likely to be tidal singularity. We can now merge the results of [16] to the ones found in this manuscript, to infer a complete phase diagram for N = 4 SYM with gauge group SU (N ) and asymptotic temperature T ∞ , at large 't Hooft coupling and infinite N on a Schwarzschild black hole background with temperature T BH . We attempt to draw such diagram in Fig. 8 as a function of τ ≡ T ∞ /T BH . The black droplets of [16] were found to exist in the window τ ∈ [0, τ D ], with τ D ≈ 0.93. The solution with τ = 0, corresponding to the Unruh vacuum, was found previously in [17]. In the droplet phase, there is no O(N 2 ) Hawking radiation, and the corresponding state can be best understood in terms of a 'jammed phase', though no underlying understanding of this phenomena exists to date on the field theory side. The funnel solutions we found in this manuscript exist at least in the range τ > τ F ≈ 0.557103, though we expect such solutions to exist for even lower values of τ , which we are unable to probe with current numerical techniques. Most notably, the Hawking flux Φ in the flowing funnel phase does not appear to be monotonic with decreasing τ , reaching a maximum at τ = τ max ≈ 0.698282. Borrowing intuition from the Stefan-Boltzman law, this would suggest that the effective number of degrees of freedom in the field theory acquires a nontrivial temperature dependence. We conjectured that flowing funnels with τ < τ max are linearly unstable to the Gregory-Laflamme instability, and that the concomitant endpoint is one of the droplets found in [16]. Using results from [49], one can give numerical evidence in favour of this conjecture, which will be presented elsewhere. It is clear that our results will largely be stable to corrections in the 't Hooft coupling, which manifest themselves as higher derivative corrections in the bulk. However, this is not the case for finite N effects. At finite N , we expect Hawking radiation in the bulk, and thus a non-trivial flux of Hawking radiation even in the droplet phase. We stress, however, that this effect is subleading in N , perhaps being O(1), instead of the O(N 2 ) effect of the flowing funnel phase.
The most resounding mystery of this work still pertains the existence of the droplet phase [16,17], in which Hawking radiation seems trapped by the local geometry and does not escape to infinity, perhaps due to some form of local confinement yet to be explored on the quantum field theory side. convergence depending slightly on which component of ξ one looks at. This convergence properties is to be expected, since non-analytic terms were identified in an expansion off the conformal boundary (see Eq. (4.2)). In Fig. 9 we show a logarithmic plot of max C |ξ t | (left), max C |ξ x | (middle) and max C |ξ y | (right) as a function of N, where C stands for domain of outer communications. For τ = 0.665557, we find max C |ξ t | ∝ N −12.22 , max C |ξ x | ∝ N −9.89 and max C |ξ y | ∝ N −9.83 . . Convergence properties of the non-vanishing components of DeTurck vector ξ a , as a function of the number of grid points N. From left to right, we plot in a logarithmic scale max C |ξ t |, max C |ξ x |, max C |ξ y |, with C being the domain of outer communications. All plots were generated using τ = 0.665557.
Finally, we also investigated the convergence properties of Φ, by defining where Φ N stands for Φ computed on a mesh with N x = N y = N grid points. The results can be seen in Fig. 10 for τ = 0.7293 where we now find δ N Φ ∝ N −5.85 , which is very close to the non-analytic power reported in the expansion (4.2).