Phase transitions in 3D gravity and fractal dimension

We show that for three dimensional gravity with higher genus boundary conditions, if the theory possesses a sufficiently light scalar, there is a second order phase transition where the scalar field condenses. This three dimensional version of the holographic superconducting phase transition occurs even though the pure gravity solutions are locally AdS$_3$. This is in addition to the first order Hawking-Page-like phase transitions between different locally AdS$_3$ handlebodies. This implies that the R\'enyi entropies of holographic CFTs will undergo phase transitions as the R\'enyi parameter is varied, as long as the theory possesses a scalar operator which is lighter than a certain critical dimension. We show that this critical dimension has an elegant mathematical interpretation as the Hausdorff dimension of the limit set of a quotient group of AdS$_3$, and use this to compute it, analytically near the boundary of moduli space and numerically in the interior of moduli space. We compare this to a CFT computation generalizing recent work of Belin, Keller and Zadeh, bounding the critical dimension using higher genus conformal blocks, and find a surprisingly good match.


Introduction
Three dimensional gravity has proven a remarkably rich testing ground for our ideas about classical and quantum gravity. Even though Einstein gravity possesses no local degrees of freedom, three dimensional theories of gravity nevertheless have many of the rich features of their higher dimensional cousins, including holography [1] and black hole solutions [2] whose Bekenstein-Hawking entropy can be computed microscopically [3]. Theories of gravity in AdS 3 are dual to two dimensional conformal field theories, allowing one to use CFT methods to gain insight into classical and quantum gravity in AdS. In this paper we will use CFT methods to motivate the existence of a new class of phase transitions in three dimensional gravity. We will then verify their existence directly in classical AdS gravity, and explore their features.
Our central result is simple, and is easiest to state for AdS 3 gravity in Euclidean signature. Such theories can be defined with a variety of boundary conditions: one can take the boundary of (Euclidean) space-time to be any smooth, two dimensional Riemann surface B. With appropriate boundary conditions [4] the theory will depend only on the conformal structure of B, so can be studied as a function of the conformal structure moduli of B. The bulk gravity path integral with these boundary conditions is, via AdS/CFT, equal to partition function of the dual CFT on the surface B. We will be interested in gravity theories in the semi-classical limit, where this bulk path integral is dominated by the classical geometry which minimizes the (appropriately regularized) gravitational action. For example, when the boundary is a sphere the dominant contribution comes from Euclidean AdS 3 , i.e. hyperbolic space H 3 , which is the unique constant negative curvature metric on the solid ball which "fills in" the boundary sphere. At higher genus, however, this path integral can have many saddle points, each of which correspond to a gravitational solution whose boundary is the surface B. For example, when the boundary is a torus the bulk saddles are constant negative curvature metrics on a solid donut which fills in the boundary torus. There are many such saddles, which are distinguished by which cycle in the boundary torus is contractible in the bulk (see e.g. [5][6][7]). For example, the geometry for which the Euclidean time coordinate is contractible is the (Euclidean) BTZ black hole, while the geometry where the angular coordinate is contractible is interpreted as the "thermal AdS" geometry used to compute finite temperature observables in a fixed AdS background. As one varies the moduli of the torus, these two saddles will interchange dominance in the bulk gravity path integral -this is the three dimensional version of the Hawking-Page phase transition [8] describing black hole formation in AdS.
We are interested in the case where the boundary B has genus g ≥ 2. Just as in the torus case, there are many different bulk solutions which give saddle point contributions to the partition function, and which can be characterized by a choice of cycles of the boundary surface B which become contractible in the bulk. The simplest of these solutions are handlebodies, where the bulk solution is the constant negative curvature metric on a solid genus g surface which fills in the boundary B. There will be phase transitions where these geometries interchange dominance: these are the higher genus versions of the Hawking-Page phase transition. As in the torus case, these handlebodies can be regarded as the Euclidean continuation of AdS 3 black holes; they are analytic continuations not of the BTZ black hole, but instead of multi-boundary black holes in AdS [9,10], as described in [11]. In the holographic context, these handlebodies describe contributions to the higher genus partition functions of holographic CFTs, which can be used to compute entanglement Rényi entropies [12], to constrain OPE coefficients [13][14][15], or as models of multi-party holographic entanglement [16]. 1 1 There are other "non-handlebody" solutions as well [17,18], which will not concern us in this paper. Near boundaries of moduli space -i.e. where cycles in the surface B become small -the handlebody solutions will always dominate [17,18]. Moreover, one can compute numerically the action of the non-handlebodies in the interior of moduli space, and -at least in the cases which have been studied -they are always subdominant compared to handlebodies [19]. We will therefore focus only on handlebodies in this paper.
The bulk solutions described above are all locally H 3 , so can be written as quotients of hyperbolic space of the form H 3 /Γ, where Γ is a discrete subgroup of the isometry group of H 3 . Indeed, Einstein gravity in three dimensions has no local degrees of freedom, so any solution of pure Einstein gravity must be locally H 3 . We are interested in more complicated theories of gravity, however, which have additional degrees of freedom. In this paper we will consider theories where we have an additional scalar field φ of mass m 2 . This means that the dual CFT has an operator O of dimension ∆, with m 2 = ∆(∆ − 2). 2 All of the solutions described above have φ = 0, and are dual to CFT configurations on the Riemann surface with O B = 0.
Our central result is the following: in some regions of moduli space, and for sufficiently light scalar fields, the handlebody solutions described above are unstable. This is because the kinetic operator (∇ 2 − m 2 ) will have a negative eigenvalue. Thus the solution with least action will not be a quotient of AdS, but rather a non-Einstein solutions with φ = 0. In the dual CFT, the expectation value of the scalar operator O B = 0 will be non-zero. This means that as the moduli are varied there will be phase transitions where these scalar fields condense. Although the general structure of these φ = 0 solutions is quite complicated -we expect the construction of these solutions to be a difficult numerical problem -we are able to prove rigorously the existence of the instability.
We will see that these phase transitions have several important features, including: Instabilities occur only when the dual operator is sufficiently light. In order for a given handlebody to be unstable, ∆ must be lighter than a certain critical value ∆ c which we will compute. The value of ∆ c will depend on the conformal structure moduli as well as on the choice of handlebody. Whenever the dimension ∆ < 2 (i.e. the bulk scalar φ has m 2 < 0) there is some region of moduli space where a given handlebody will be unstable. 3 Instabilities occur only in the interior of moduli space. At the boundary of moduli space one of the handlebody phases will always dominate, and will be stable against the condensation of any scalar field. In other words, ∆ c → 0 for the dominant handlebody as we approach the edge of moduli space. As one moves into the interior of moduli space ∆ c increases so the handlebody becomes more unstable to condensation of the scalar, until the first-order Hawking-Page transition is reached, and a topologically distinct handlebody becomes dominant.
For example, if B is a genus g = n − 1 surface constructed as an n-fold cover of the sphere branched over 4 points, parameterized by their cross-ratio x, the handlebody which dominates when x → 0 will become more unstable as x is increased. For this geometry, ∆ c is a monotonically increasing function of x. 4 Handlebodies become more unstable as the genus increases. Instabilities only occur when B has genus g ≥ 2. If B is the genus g = n − 1 surface constructed as an n-fold cover of the sphere branched over 4 points, then if we hold the cross-ratio x of the four points fixed, the corresponding handlebodies will become more unstable as n is increased. In other words, ∆ c is a monotonically increasing function of n. As we take n → ∞ with fixed x, ∆ c approaches a finite value which depends on x but is always greater than 1 2 .
This new phase transition in three dimensional gravity is quite similar to the holographic superconducting phase transition in higher dimensions [20,21]. There is, however, one crucial difference, which is that in the present case the solutions which become unstable are locally AdS 3 , and have no external potentials (aside from metric moduli) turned on. Our instability occurs because of global properties of the handlebody, not due to any local properties of the metric. 5 Although at first sight surprising, our results are a three dimensional version of a famous fact in two dimensions: the spectrum of the hyperbolic Laplacian on a Riemann surface depends not just on the local structure of the metric (which is always hyperbolic) but also on the moduli of the Riemann surface. Indeed, this spectrum is the central object of interest in the study of arithmetic and quantum chaos, and many of our results are borrowed from this literature.
Our results have important implications for entanglement entropies in two dimensional conformal field theories. For any state in a two dimensional CFT, one can consider the reduced density matrix associated to a particular spatial region. The Rényi entropies 1 1−n log Tr ρ n can then be used to characterize the spatial entanglement structure of this state. When the spatial region is collection of intervals, the Rényi entropy is -via the replica trick -equal to the partition function of a CFT on a higher genus Riemann surface whose genus depends on n (see e.g. [23]). For example, the Rényi entropy for a pair of intervals in the vacuum state is equal to the partition function on a genus g = n − 1 Riemann surface. The entanglement entropy is then computed by considering these Rényi entropies as an analytic function of n, and continuing to n → 1. In this procedure one assumes that the entropies are analytic functions of n. We have seen, however, that in holographic CFTs the Rényi entropies can undergo a phase transition as n is varied, at some finite value of n > 1. Thus the replica 4 It is important here that we are referring to the handlebody which dominates at small x. For the handlebody which dominates as x → 1, ∆c will be monotonically decreasing function of x between 0 and 1. 5 A rather similar phenomenon was observed in higher dimensions in [22], where hyperbolic black holes in AdS4 were observed to undergo similar phase transitions even though the solutions were locally AdS. As in the present case, the instability only arose because the hyperbolic black hole solutions differed globally from AdS4. Thus modes of the scalar field which are not normally present (since they are non-normalizable in global AdS) suddenly become normalizable and lead to a genuine instability of the locally AdS solution.
method for computing entanglement (von Neumann) entropies must be treated with care. 6 For example, if we consider the Rényi entropies for a pair of intervals, two handlebodies will interchange dominance precisely at cross-ratio x = 1/2 [12]; this is also exactly where the Ryu-Takayanagi formula for entanglement entropy [27] will undergo a phase transition. Our results imply, however, that if the CFT has a sufficiently light operator then the Rényi entropies will undergo a phase transition at cross-ratio x < 1/2. For example, the n = 2 Rényi entropy will undergo a phase transition if the theory has an operator with dimension ∆ < ∆ c = 0.189124 · · · . We note that this instability occurs for values of n which are strictly larger than one -we do not expect a non-analyticity in a neighborhood of n = 1. It would be interesting to revisit the arguments of [12,[28][29][30] in this context. Our results also make clear a sense in which higher genus CFT partition functions differ qualitatively from those on the sphere or torus. The torus partition function, for example, was shown by Hartman, Keller and Stoica to take a universal form at large central charge, provided one assumes that the spectrum of light states (i.e. those with dimension less than the central charge) does not grow too quickly [31]. This universal form is precisely that of a dual three dimensional theory of gravity which has a Hawking-Page transition between a thermal state and a BTZ black hole, and the sparseness condition is obeyed by any bulk local quantum field theory and even by string theories with string scale string AdS . At higher genus, however, we see that additional phase transitions are generic, and occur even for duals of local quantum field theories in the bulk. Thus at higher genus there is no analogous "universal partition function" at large central charge.
While our discussion will be entirely in the context of three-dimensional gravity, similar phenomena will also occur in higher dimensions. The most direct analogue is with locally AdS spacetimes, to which almost everything generalizes straightforwardly. In particular, the critical dimension ∆ c for an instability can be shown to be equal to the Hausdorff dimension of an appropriate limit set, just as we will see below for the three dimensional case. 7 This is slightly less natural than in three dimensions, because solutions to Einstein's equations need not be locally hyperbolic, so it is not clear when such geometries would dominate the path integral. More generally, the same mechanism of instability can apply, with global properties of the solution moving the critical mass above the naïve Breitenlohner-Freedman bound. Heuristically, scalars of negative mass squared can be stable because the reduction in action from the mass term in a finite region is compensated for by the positive contribution to the action from the gradient, required to match with the boundary conditions at infinity; it is important here that the volume of a region does not grow faster with size than its perimeter in negatively curved spaces. Without altering the local curvatures, nontrivial topology can upset this mechanism for stability by reducing the size of a region's boundary, and hence the gradient contribution to the action, for a given volume.
The discussion of the paper will be phrased in terms of the Euclidean solutions, but the 6 Similar phenomena were observed for spherical entangling surfaces in higher dimensional holographic CFTs in [24,25], and for the three dimensional O(N ) model in [26]. 7 Indeed, one can directly reinterpret the results of [22] in this context. Lorentzian signature. The relevant CFT states are  defined by a Euclidean path integral on a Riemann surface with one or more boundaries,  generalizing the familiar examples of the path integral on the disc preparing the vacuum  state, and on the cylinder preparing the thermofield double state on two entangled copies of the CFT Hilbert space. This defines the state at t = 0, which can be evolved in Lorentzian time.

results have interesting implications in
To find the semiclassical bulk dual of these states, we must first find the Euclidean solution that dominates the path integral on the 'Schottky double', the closed Riemann surface formed by gluing the surface to its mirror image along each of the boundaries, by construction producing a Z 2 -symmetric surface. The dominant solution is expected to respect this boundary time-reflection symmetry, so the bulk surface Σ fixed by the reflection acts as an initial data surface for Lorentzian evolution, and the quantum state of the bulk fields is the Hartle-Hawking wavefunction on Σ. For one possible solution, the t = 0 slice Σ is conformal to the original Riemann surface, describing a single-exterior black hole with topology hidden behind a horizon for a single-boundary case, or a multi-boundary black hole with an exterior region for each boundary, all joined by a non-traversable wormhole. Even in pure gravity, there are several phases of the dominant bulk solution, so depending on the moduli the bulk state can also be disconnected copies of pure AdS (but with fields in a state different from the vacuum), or something else. For a more detailed review of these states, see [19,[32][33][34]. Now, if there is a sufficiently relevant scalar operator in the CFT, there is an additional second-order phase transition to a dominant bulk solution with a nonzero classical value for the dual scalar field. This means that the initial data on Σ includes some scalar field configuration, which will evolve in time. The fact that these states are not stationary will then be visible even for a classical observer outside any horizon. When the phase includes a black hole, the scalar outside the horizon will rapidly decay away, falling into the black hole. A more interesting time evolution occurs when the dual state includes copies of pure AdS, which may now include some scalar configuration. When the amplitude is small, as will be the case close to the transition, linearised evolution will suffice, with the field bouncing around periodically, but eventually nonlinearities will likely become important, with resonances between different modes. Perhaps the most likely evolution thereafter is a turbulent cascade to excite higher and higher frequency modes, with the solution nonetheless being regular for all time, as evidenced by numerical studies of a massless scalar interacting only gravitationally [35]. This is different from the situation in higher dimensions, in which a black hole forms after finite time; this cannot occur in three dimensions, because there is a finite energy threshold between the vacuum and the lightest black hole 8 .
In section 2 we will review briefly a few salient features of three dimensional gravity, as well as the necessary aspects of CFT on Riemann surfaces. In section 3 we will give a CFT argument for the existence of an instability, inspired by recent results of Belin, Keller and Zadeh [36]. The main idea is that a free bulk scalar field is dual to a generalized free field in the boundary CFT, and we can compute the contribution of such a field to the higher genus CFT partition function. Using the asymptotic value of the OPE coefficients of multitrace operators built from a generalized free field, along with higher genus conformal blocks in the appropriate regime, one can show that these contributions diverge when the field is sufficiently light, which signals the phase transition. This argument allows us to bound the critical dimension of the scalar field; for example, for the genus two handlebody relevant for the computation of the third Rényi entropy of two intervals at cross-ratio x = 1 2 , we find ∆ c ≥ 0.189121 · · · .
In section 4 we will turn to the bulk instability. We will first review how the zero mode of the instability relates to various notions from the spectral theory of the Laplacian on a general bulk geometry M. We then specialize to the main case of interest, for which M is a quotient of hyperbolic space H 3 such as a handlebody, and find that the critical dimension has a rather beautiful mathematical interpretation. The quotient is by a group of Möbius maps, which has a limit set, a subset of the Riemann sphere (the boundary of Euclidean AdS 3 ). This limit set has a finite Hausdorff dimension δ > 0, which is sometimes referred to as the fractal dimension of the limit set. This Hausdorff dimension is precisely equal to the critical dimension of the scalar field, ∆ c = δ. In other words, a scalar is unstable if and only if its dimension is less than the Hausdorff dimension of the limit set. The calculation in section 3 can therefore be regarded as a CFT estimate of this Hausdorff dimension, which provides explicit lower bounds on δ.
In section 5 we turn to the explicit computation of the critical dimension, using an algorithm of McMullen for computing the Hausdorff dimension. We will use the algorithm to compute the critical dimension analytically, finding the asymptotic behaviour of δ as the boundary of moduli space is approached, for the handlebody which dominates the partition function. We also describe what happens to the instability at large genus. We will also use the algorithm to efficiently compute the critical dimension numerically. For example, for the genus two surface described above, the Hausdorff dimension is δ = 0.189124 · · · , close to our CFT bound. We use the numerical data to provide plots of this critical dimension as a function of moduli, and as a function of genus.

Review of higher genus partition functions in 3D gravity and 2D CFT
In this section we will review the description of higher genus Riemann surfaces, and the construction of solutions to three-dimensional gravity with such boundaries, which can be interpreted as saddle points for the higher genus partition function of a holographic CFT. In particular, we describe a class of symmetric surfaces that we will use as examples. We will also review the interpretation of these partition functions in terms of Rényi entropies.

Moduli spaces and handlebodies
We are interested in studying holographic two-dimensional conformal field theories, dual to three dimensional AdS gravity in Euclidean signature, in particular on a Riemann surface B of genus g ≥ 2. The partition function of the theory on such a surface, denoted Z g (τ ), will depend on the conformal structure of the surface B. Here τ is a collection of 3g − 3 complex coordinates which parameterize the moduli space M g of conformal structures on B. At genus one, τ can be identified with the usual torus modulus. At higher genus there are various different coordinates which can be used to describe the moduli τ , some of which we will now describe. 9 For many purposes in CFT and gravity, the most convenient way to realize a Riemann surface B is as a quotient of the Riemann sphere C * by a Schottky group Γ, B = Ω(Γ)/Γ. Here Γ is a discrete subgroup of P SL(2, C), which acts on C * in the usual way by Möbius transformations, and Ω(Γ) is the set of points on the Riemann sphere where this group 'acts nicely'. More precisely, Ω(Γ) is the set of points z ∈ C * which have some neighborhood U containing no other images of z under the group: γ · z ∈ U for γ ∈ Γ implies that γ is the identity. Equivalently, if we define the limit set Λ(Γ) to be the set of accumulation points of the action of Γ on C * (a set about which we will have much more to say later), then Ω(Γ) is just the Riemann sphere with those points removed: Ω(Γ) = C * − Λ(Γ). More specifically, a Schottky group Γ of genus g is a subgroup of P SL(2, C) that is freely generated by g loxodromic 10 elements of P SL(2, C), having as a fundamental domain the exterior of 2g closed curves (usually circles), such that each of the g generators of Γ maps one of these boundaries to another in pairs. Intuitively, to obtain a Schottky representation, we can cut the surface along g disjoint closed loops such that it stays in one piece and becomes a sphere with 2g holes, flatten it onto the complex plane, and build the Schottky group from the Möbius maps that glue the surface back together along its g seams. A given Riemann surface can be written as a Schottky group in many different ways, depending on the choice of g cycles to cut along. The presentation as a Schottky group is equivalent to the plumbing construction used in [37]. A more detailed review of Schottky uniformization can be found in [11,28,38]. A slightly different approach to calculations in the Schottky coordinates was used in [39].
A rather different presentation of the Riemann surface B is as an algebraic curve. In this 9 Because of the conformal anomaly the partition function will in addition depend on a choice of metric within a given conformal class. Thus Zg should not -strictly speaking -be regarded a function of τ alone. This dependence, however, involves only the central charge and not any of the other dynamical data of the CFT (such as operator dimensions or structure constants), so will not be important for us here. We will therefore suppress this dependence and simply indicate the dependence on the conformal structure moduli τ . 10 A loxodromic element γ of P SL(2, C) is one which is conjugate to q 0 0 q −1 for some q with 0 < |q| < 1.
case we represent B as the set of solutions to an equation such as Here, B is a genus g = (N − 1)(n − 1) surface, represented as n-fold branched cover over the Riemann sphere parameterized by the z-plane, with 2N branch points (u k , v k ). As the resulting Riemann surface automatically possesses a Z n symmetry (usually referred to as replica symmetry) where one permutes the n sheets, one cannot describe a general point in moduli space M g using this parameterization. Instead, this equation describes only a 2N − 3 dimensional slice of moduli space, a family of surfaces with an enhanced (Z n ) automorphism group. Except in special cases one cannot map out the full moduli space this way. 11 The advantage of this approach is that the moduli of this surface are easy to describe -they are the locations (u k , v k ) of the branch points. For physics, this presentation is natural for describing Rényi entropies (or certain correlation functions in orbifold theories), with the branch points u k and v k corresponding to the insertion points of twist and anti-twist operators. Except at genus one, or in special cases with very high symmetry, it is not possible to find an explicit map between the moduli of the algebraic curve and Schottky groups, or to find out whether two Schottky groups represent the same surface, sliced in a different way. However, the problem of finding a Schottky group associated with a particular algebraic curve, sometimes called 'Schottky uniformization', is equivalent to solving a monodromy problem, which we now briefly describe.
To do this, we begin by denoting the locations of branch points (u k , v k ) as z i (i = 1, . . . , 2N ). We would like to find the map w(z) from the algebraic curve coordinate z to the coordinate w of the complex plane on which the Schottky group acts. But w(z) is not single-valued, because there are many possible values of w related by elements of the Schottky group Γ. However, the Schwarzian derivative T c (z) = S(w)(z) = w w − 1 2 w w 2 is singlevalued, since the Schottky group consists of Möbius maps. If we take T c (z) as given, a simple calculation shows that solving T c (z) = S(w)(z) for w is equivalent to solving the ordinary differential equation with w(z) = ψ 1 (z) ψ 2 (z) being the ratio of two linearly independent solutions ψ 1,2 (z) to the ODE. This is not much use if we know nothing about T c (z). However, for the Riemann surface (2.1) T c (z) can be fixed up to a finite number of parameters, by using the fact that it is a meromorphic function of z which transforms like a stress-tensor: One important special case is the genus two moduli space. Every genus two curve is hyperelliptic, so can be represented as a 2-fold cover of the sphere branched over 6 points.
Here, we have assumed that the Schottky group respects the replica symmetry, so T c (z) is single-valued in z. The double poles are fixed by demanding smoothness in the y coordinate of eq. (2.1), and the γ i are free parameters, called 'accessory parameters'. It is also constrained by smoothness at infinity, which demands that T c (z) decays like 1 z 4 . This imposes three constraints on the γ i , leaving 2N − 3 free parameters.
It remains to fix these free parameters. To do this, note that if we go around a closed curve on the surface, a solution to the ODE will not usually come back to itself, but undergo monodromy, so the value of w will change by a Möbius map: The monodromies of the ODE form a representation of the fundamental group of the surface π 1 (B) by Möbius maps. But in the Schottky representation, not all the closed loops on the surface should take us to a different w, and a different copy of the fundamental domain for Γ: the g special loops that bound the fundamental domain should come back to the same value of w, and so correspond to trivial monodromy of the ODE eq. (2.2) (in fact, the monodromy matrix around these cycles is always minus the identity). Imposing these trivial monodromy conditions is precisely enough to fix the 2N − 3 free parameters. Once these parameters are fixed, we may solve eq. (2.2) to find the monodromy around g complementary cycles, which give the generators of the Schottky group Γ. In section 2.3 we go into more detail for a specific example, which we will subsequently use for analytic and numerical calculations.
This monodromy problem also appears in computations of the semiclassical limit of Virasoro conformal blocks [40], reviewed in [29,41], and described in generality for higher genus blocks in [13].

Schottky representations, handlebodies, and gravity
The Schottky representation has a very natural interpretation from the bulk point of view. To see this, note that the Möbius maps acting on the Riemann sphere can be extended into a bulk hyperbolic space H 3 , where they act as the orientation-preserving isometries. We can therefore extend the action of a Schottky group Γ into this bulk, obtaining a quotient of hyperbolic space with B as its boundary, M = H 3 /Γ. Representing hyperbolic space in the upper half-space model, this can be understood as taking the circles that bound the fundamental domain of the Schottky group and extending them as hemispheres into the bulk, giving a fundamental domain with the hemispheres identified by the generators of Γ.
A CFT partition function on the surface B can be computed holographically as the bulk gravity path integral over Euclidean geometries whose conformal boundary is B. Semiclassically, we just need to compute the action of a solution to the bulk equations of motion with boundary B, which will depend on the moduli of B. There are an infinite number of bulk solutions, which in general should include the contribution of matter fields, but a particularly simple class of solutions are those without matter fields turned on. Since pure gravity in three dimensions is locally trivial, Einstein's equations then imply that the bulk is locally H 3 , which means that it must be a quotient of hyperbolic space. The Schottky group quotients therefore provide a large class of solutions to the bulk problem, which are conjectured to dominate the path integral in pure gravity.
Topologically, the Schottky group quotients are handlebodies, obtained by 'filling in' the surface B along a choice of g cycles. These contractible cycles are precisely those we chose to cut the surface along to construct the Schottky group, or around which we imposed trivial monodromy. The Schottky group describes the remaining non-contractible cycles, in the sense that it is topologically interpreted as the fundamental group of the handlebody.
Some geometric properties of the bulk can be read off easily from the Schottky group, in particular the lengths of closed geodesics. A closed loop is represented topologically as a conjugacy class in the fundamental group, or equivalently in Γ, and since the eigenvalues of γ ∈ Γ are independent of the representative of the conjugacy class, the smaller eigenvalue q γ of γ ∈ SL(2, C) (0 < |q γ | < 1) is naturally associated with a closed curve. Writing q γ = e − 1 2 ( +iϑ) , is in fact the length of the closed geodesic, and θ is the amount the geodesic is twisted by (the angle a normal vector rotates by after parallel transport round the curve). Explicitly in terms of the trace, this length is For more detailed review and discussion, focussing on the Lorentzian versions of these geometries, see [32]. Because there are many possible Schottky groups corresponding to the same surface B, we must decide which geometry gives the correct semiclassical bulk dual for given moduli of the boundary surface (even before considering bulk matter fields). The naïve answer to this, and the one that reproduces CFT expectations, is that the handlebody with least action dominates the path integral, so is the dual bulk. It was shown how to compute this action in [11], from a particular higher-genus 'Liouville action' [42], depending crucially on the IR cutoff imposed on the bulk, and hence on the choice of metric on the boundary surface within the given conformal class. In a metric appropriate for Rényi entropies, flat away from conical singularities at branch points z i , and with a bulk preserving replica symmetry, the derivative of the action with respect to z i is proportional to the accessory parameter γ i [28]. In a constant curvature metric, for general surfaces, a numerical algorithm to compute the action was given in [19]. As a heuristic, to choose the dominant saddle point, the g shortest cycles of the surface should be filled in.
Given these tools, one can then attempt to construct the the higher genus partition function Z g (τ ) via a bulk path integral, as a sum over geometries. The handlebodies described above give semi-classical saddle point contributions to this bulk path integral, and the full partition function should be given be a sum over these semi-classical saddles along with a set of loop corrections. The loop corrections to these semi-classical contributions can be computed exactly at genus g = 1 and perturbatively at higher genus (see e.g. [7,43]). In pure gravity -i.e. in theories with no degrees of freedom aside from the metric -there is some hope that one could compute the higher genus partition function exactly [44]. We will be interested in more general theories, which contain scalar fields in addition to the metric. In this case the theory has local bulk degrees of freedom, and there is little hope of an exact computation. Nevertheless, the computations described above give contributions to the partition function of a holographic CFT which will be valid in the semi-classical (large c) limit.

A Z n symmetric family of genus n − 1 surfaces
We now illustrate this general discussion with an example, specifically a family of genus n − 1 surfaces with an enhanced Z n symmetry. As an algebraic curve, this family of Riemann surfaces is given by This is the N = 2 case of (2.1), where we have used P SL(2, C) transformations to put u 1 = 0, u 2 = 1 and v 2 = ∞. The remaining parameter is the cross-ratio x, which is the modulus of this family of Riemann surfaces. In general x can be any complex number, but for simplicity (and for the purposes of applications to Rényi entropies, described below) we will take it to be a real number between zero and one.
To find a Schottky group, or equivalently a bulk geometry, we can now solve the monodromy problem described above. Choosing to preserve the replica symmetry, the most general ansatz for the ODE eq. (2.2) is where we have imposed the constraints (which are slightly different, because there is a branch point at infinity), leaving the single accessory parameter γ. To fix this parameter, we must first choose the cycles around which we impose trivial monodromy. For our purposes, it suffices to take the cycles surrounding 0 and x; this gives a loop on the z plane enclosing one zero and one pole of y n , so remains on the same sheet of the branched cover, forming a closed loop on the surface. There are n of these, one on each sheet, but only n − 1 of them are in fact independent: the product (in the fundamental group) of the n loops, described by a loop enclosing z = 0 and z = 1, then moving to the next sheet, and repeating n times, is topologically trivial. Because we have imposed Z n symmetry already on our ansatz for T c (z), imposing trivial monodromy on any one of the sheets is sufficient.
Having chosen the accessory parameters to trivialize the monodromy around these cycles, we would like to read off the Schottky group. To do this, it is convenient to take full advantage of the symmetry of the situation, using the automorphisms of the surface (which are preserved by the handlebody). In the language of the quotient, an (orientation-preserving) isometry of the bulk is represented by an additional element γ ∈ P SL(2, C) (so it is an isometry on the covering space H 3 ) that commutes with the group Γ (γΓ = Γγ, so it has a well-defined action on the quotient). Including some such elements, we can form an extended groupΓ, of which Γ forms a normal subgroup. The largest possibleΓ, including all elements of γ ∈ P SL(2, C) such that γΓ = Γγ is the normalizer N (Γ) of Γ, and the isometry group of the bulk is then The most obvious extension providing an automorphism is the Z n replica symmetry R, represented as the monodromy around a loop containing 0. This comes back to a different sheet, so is not an element of Γ, but we can include it inΓ as an elliptic Möbius map of order n. From the point of view of the monodromy problem, these elements correspond to monodromy along curves that may not return to the same point on the surface, but go between some point and its image under the isometry. Near z = 0, the independent solutions to the ODE eq. (2.2) look like ψ ± (z) ∼ z n±1 2n , with corrections forming a power series in z and not affecting the monodromy around zero; choosing these as our basis ψ ± of solutions (w(z) = ψ + (z) ψ − (z) ), the loop around zero, enacting the replica symmetry, acts on the w coordinate as R : w → e 2πi n w. In fact, this family of surfaces automatically has more symmetry, containing an additional Z 2 extending the Z n to a dihedral group D 2n 12 . From the z coordinate, this can be understood as a map swapping 0 with 1 and x with ∞, z → x z , along with reversing the order of the sheets of the cover. It is straightforward to check that this leaves the ODE invariant, after transforming ψ as a weight − 1 2 field. In terms of the monodromy, this extra symmetry is enacted by taking the solutions ψ ± (z), following the solution from 0 to ∞, and reading off the coefficients of z n±1 2n in these solutions near ∞, giving some order two Möbius map S. In practice, except for the special case n = 2, finding S requires doing the calculation numerically, but we can deduce a lot about it, reducing the unknown parameters from the three numbers specifying a general Möbius map, to just one. Firstly, note that doing S twice corresponds to going round a loop with trivial topology, which implies that S is order two, S 2 = 1, which means it is specified by its fixed points. Secondly, without altering the form of R, we can change coordinates by rescaling and rotating in the w plane, and use this to remove one other parameter of freedom. We will use this freedom to set the product of the fixed points of S to be unity, which fixes S to act as S : w → w−ζ ζw−1 for some (in general complex) ζ. Now, the extended groupΓ is generated by just R and S (in fact, it is the free group generated by those elements with the only relations being those given by the orders of the elements:Γ = R, S|R n = S 2 = 1 ), so we are interested in the one-parameter family of groups generated by one Möbius map of order 2, and one of order n. The actual Schottky group Γ appears as a normal subgroup of this, generated by the loops that actually return to the same point on the surface (but are still non-contractible in the bulk), requiring an even number of S generators to appear, and also for the sheets to map back to themselves, rather than being permuted. The first of these is the element γ 1 = SRSR, taking a loop 12 Even further than this, these surfaces all have another additional Z2 commuting with this dihedral group, acting as z → x z−1 z−x , which is a hyperelliptic involution of B. We will not make use of this extra symmetry, but in the parameterization used below, it can be included as w → 1 w .
round zero by R, then going from 0 to ∞ by S, then a loop round to infinity by R again, and finally back to the starting point at zero, creating a closed loop surrounding zero and infinity, or equivalently x and 1. The remaining generators are similar, but starting on a different sheet, achieved by conjugating with R: γ k = R 1−k SRSR k . There are n of these, but they are not all independent, since γ 1 γ 2 · · · γ n = 1. Any n − 1 of these (a number equal to the genus) generate Γ. To relate this to the general discussion of symmetries above, the group of isometriesΓ/Γ described by this extension is the dihedral group of order 2n, since modding out by Γ is equivalent to imposing the additional relation This prescribes the family of Schottky groups we are interested in, parameterized by ζ = cos θ, though it is important to note that this only describes a Schottky group when ζ is sufficiently close to one (or θ close to zero). An alternative, more geometric parameterization is by q, defined as the smaller eigenvalue of γ 1 = SRSR ∈ SL(2, C) (noting that this is independent of the sign chosen for the matrix representatives of S and R), defined as above so that q = e − 1 2 ( +iϑ) gives the length and twist of a curve in the bulk geometry. In the case n = 2, q = tan 2 θ 2 is the usual elliptic nome of the boundary torus, lying in the punctured open unit disc, though for larger n it must be contained in a strictly smaller region.
In practice, to map from any given x to find the corresponding value of ζ (or θ or q), it is sufficient to compute the trace of any of the γ k (all are equal), from the trace of the monodromy of any loop containing x and 1. It is also possible to solve the monodromy problem perturbatively in small cross-ratio, as described in [38], and in our parametrisation, the result to leading order is θ = . Real x between 0 and 1 corresponds to real θ, or 0 < ζ < 1, or 0 < q < 1. Real negative x (or equivalently x > 1) also results in a Fuchsian group Γ, and corresponds to ζ > 1 (but bounded by sec π n so the group is Schottky). We can also find different handlebodies for the same surface by trivializing monodromy around some different cycle, but the only possibilities preserving the replica symmetry are much the same, the most obvious being to take the loop to surround x and 1 rather than 0 and x. There is a phase transition between the two corresponding handlebodies at Re(x) = 1 2 [12].

Relationship with Rényi entropies
The bulk computation of higher genus partition functions can be applied to the computation of Rényi entropies in holographic CFTs (see [45] for a review). For a density matrix ρ, the nth Rényi entropy is defined as In the n → 1 limit this becomes the von Neumann entropy S = − Tr(ρ log ρ). In order to probe the spatial entanglement structure of the theory, we can take ρ to be the reduced density matrix for a spatial region A in the vacuum state. Then ρ is defined by a Euclidean path integral on the sphere (with cuts introduced at A), and the Rényi entropy may be computed by gluing n copies of this sphere together along these cuts. Explicitly, we have Here Z n is the partition function on a manifold M n , which is the n-fold branched cover defined by gluing n copies of the original spacetime manifold along A, and the normalization constant Z 1 is the sphere partition function. If A consists of N disjoint intervals, this is precisely the n−fold cover of the sphere branched over 2N points (the endpoints) described above. We will focus on the case where A consists of two disjoint intervals [u 1 , v 1 ] and [u 2 , v 2 ]. Then the conformal structure of M n is completely determined by the cross ratio As this cross-ratio is varied, we sweep out a one (real) dimensional slice of the moduli space M n−1 of genus n − 1 dimensional Riemann surfaces. This is precisely the case described in the previous subsection. There are two handlebodies which compute the pure-gravity contribution to the Rényi entropies, which exchange dominance at the point x = 1/2. When n = 2 these saddles are precisely the thermal AdS and Euclidean BTZ black hole solutions, and the phase transition at x = 1/2 is the usual Hawking-Page phase transition. The scalar instabilities we describe in this paper will occur for n > 2 when the theory has an operator of dimension ∆ which is sufficiently light. In particular, we will find that there are two new phase transitions as x is varied, one at x = x c (∆) < 1/2 and one at x = 1 − x c (∆), where these two handlebodies will become unstable to the formation of scalar hair.

The phase transition from CFT
In this section we will make a CFT argument for the instability, by considering the contribution of a generalized free field -the boundary avatar of a free bulk scalar field -to the higher genus partition function of a CFT with large central charge. The result is that if the corresponding operator is sufficiently light, then the generalized free partition function will diverge somewhere in the interior of moduli space. This signals that the free approximation has broken down, so interactions become important, and the partition function will undergo a phase transition. We give an analytic lower bound on the critical dimension ∆ c in terms of the Schottky moduli of the surface. As the field becomes lighter, the phase transition will occur closer to the boundary of moduli space; in particular, for a sufficiently light field the corresponding bulk phase transition will occur before the usual "Hawking-Page" transition where (locally Einstein) bulk saddles are interchanged. In the bulk, this would be interpreted as the condensation of a bulk scalar field in a handlebody background. The discussion in this section is a refinement of the arguments presented in [36].

The partition function and conformal blocks
A higher genus partition function can, at least in principle, be computed from the basic dynamical data of the CFT, namely the spectrum of dimensions and spins (∆ i , s i ) of primary operators, along with their three-point coefficients C ijk . To do this, we can insert a complete set of states on a handle of the surface to reduce the computation to sum over two-point functions on a surface one genus lower, and repeat this (along with use of the OPE) until the computation has been reduced to three-point functions on spheres. A complete decomposition like this can be understood by cutting up the surface into pairs of pants: any genus g ≥ 2 surface can be decomposed (in many ways) into 2(g − 1) pairs of pants, joined along a total of 3(g − 1) cuffs. Along each of these cuffs, we can insert a complete set of states and, by the state-operator correspondence, the amplitude between three states defined by the path integral on the pair of pants is determined by a three-point coefficient.
Along with inserting complete sets of states in this way, we can use the fact that the states are arranged in multiplets of the Virasoro algebra, by summing up all contributions from a given multiplet appearing in the sums over states. The resulting object, collecting the contributions from a given primary on each of the 3(g − 1) cuffs along with all their descendants, is a higher genus conformal block F. This is determined by kinematics alone, depending only on the scaling dimensions and spins of the primaries chosen, the central charge, and the moduli of the surface 13 . In the partition function this conformal block will be multiplied by the product of 2(g − 1) OPE coefficients corresponding to the primaries.
Summing over all possible choices of primaries, the result is a general expression for the partition function of the form Here the sum is over all choices of the 3(g − 1) primary operators, and C 2(g−1) {i} denotes the product of 2(g−1) OPE coefficients corresponding to the primary operators propagating down the legs of each pair of pants. This expression may look very different for different pair-ofpants decompositions of the surface, but the result must be equal whichever decompositions is chosen; this is the statement of higher-genus crossing symmetry, which can be exploited to constrain CFT data [13][14][15].
In the case g = 2, there are two possible distinct types of decomposition into pairs of pants, depending on whether we choose to insert a complete set of states on a cycle dividing the surface into two pieces. Assuming we do not, the decomposition looks like two pairs of pants joined to one another along each of their three cuffs, as illustrated in fig. 1, and C 2(g−1) {i} is just C 2 ijk , where i, j, k denote the primaries chosen on each of these seams. In general, it is rather difficult to compute (3.1) explicitly, and even the conformal blocks cannot be calculated exactly in closed form. It is possible to calculate perturbatively in Figure 1: A genus 2 surface, which is cut into two pairs of pants glued together along the three black circles. Along each of the three circles we can insert a projection onto the descendants of a primary of dimension ∆ i (i = 1, 2, 3) to obtain the block F({∆ i }, c; τ ). The dual handlebody is found by 'filling in' the surface, as indicated by the shaded disks. The block F({∆ i }, c; τ ) can be computed in the bulk in a semi-classical approximation, valid in the limit 1 ∆ i c, by computing the action of the network of bulk geodesics indicated in red.
moduli of the surface, or using recursion relations exploiting the structure of degenerate representations [14,37], or in various semiclassical limits [13,46]. In particular, one needs to choose a conformal frame, and to account carefully for the way in which the surface is glued together from its constituents. Fortunately we will not need to work with this expression in generality, only requiring the blocks in a particular 'semiclassical global' limit.

Semiclassical global limit of higher genus blocks
We will require the blocks in a limit of large central charge, where the dimensions of exchanged operators are large also large (with ratios between different ∆ i fixed in the limit), but small compared to c: 1 ∆ i c. This limit has a dual holographic description in terms of semiclassical gravity coupled to particles in a probe limit, for which only the global conformal sl(2) subalgebra of the Virasoro algebra is important; hence the name 'semiclassical global' blocks.
In this limit the blocks simplify, becoming The functions S 0 and S 1 , depend on the moduli and (in the case of S 1 ) on the ratios of conformal dimensions, and have semi-classical gravity interpretations which we will describe below. The fact that the blocks exponentiate in this limit large c limit is most well known in the case of four-point functions [47], but has not been rigorously proven. It is, however, physically well-motivated, for example by considering a semiclassical limit of Liouville theory [41]. First, S 0 can be interpreted as the semiclassical vacuum block, i.e. the block for which all operators are taken to be the identity. It is equal to the on-shell action of the handlebody where each of the cycles in the pair-of-pants decomposition are chosen to be contractible. 14 This depends on a choice of conformal frame (the conformal anomaly precisely takes the form of a shift in S 0 ), and in general can only be computed exactly by numerics [19,28,29]. The frame-independent information contained in S 0 is the difference between its value in different channels, and this data is required to impose crossing symmetry. Luckily, for our purposes we will only need to express the partition function in one channel, so we can entirely disregard this piece.
For us, the more important contribution is S 1 , encoding the dependence on the dimensions. The important point is that in the limit c → ∞ with fixed ∆ i the contribution of the Virasoro descendants is unimportant after we factor out the vacuum block contribution S 0 . So the block reduces to a 'global' block, where only the L −1 descendants are kept [37] 15 . The gravitational interpretation is clear: as c → ∞, the backreaction from the matter and the loop corrections from the graviton can be ignored, and we need only the classical background action.
This global block is still tricky to compute at higher genus, but if we further assume that the internal dimensions ∆ i are large (but still much smaller than c), it simplifies to a 'semiclassical global' block. This can be determined by considering a network of geodesics in the handlebody spacetime, determined by the pair-of-pants decomposition. Specifically, for each pair of pants, assign a trivalent vertex in the bulk, and join these vertices by a geodesic for each seam joining the pairs of pants. This geodesic is interpreted as the worldline of a particle of mass m i ∼ ∆ i , determined by the dimension of the primary operator assigned to the corresponding cuff of the pants decomposition, and is assigned an action ∆ i i , where i is the length of that geodesic segment. We finally must specify the bulk locations of the vertices; these are chosen to minimize the total particle action i ∆ i i . The value of e − min i ∆ i i at this minimum obeys the semiclassical limit of the Casimir equations for the global conformal blocks, as shown in [46], so reproduces the correct moduli dependence of the blocks. 16 This 14 This only depends on the choice of g cycles, not the full decomposition, which is consistent because of the fusion rules of the identity: choosing the identity module along g cycles is enough to imply that the identity must be present in the other 2g − 3.
15 What constitutes a global descendant is a little ambiguous for higher-genus surfaces, since it is not invariant under general conformal transformations. The statement here requires a Schottky, or plumbing frame, for which all transition maps are Möbius maps. 16 For some values of the dimensions and moduli, the action can be minimized when one of more of the geodesics shrink down to zero size, in which case the block is given instead by some complexified saddle point. From the bulk point of view, this happens when the leading order amplitude in the large ∆ limit comes from double-trace contributions, rather than the case we would like to consider, where these are exponentially limit can be used as a starting point for a systematic perturbative expansion for the blocks, developed in terms of worldline quantum field theory coupled to gravity in [48]. The only thing that remains to fix is the overall normalization. The geodesic network prescription comes with an unambiguous normalization, but rather than being the canonical one, where we multiply by the appropriate OPE coefficients to find the contribution to the partition function, it comes with some nontrivial OPE coefficientsĈ(∆ i ) built in, depending on the dimensions of the operators meeting at each vertex. To compute these, consider taking the pinching limit in which all cuffs of the pairs of pants become small, suppressing the descendants so only the product of primary three-point functions remains. The function C(∆ i ) can therefore be computed by using a geodesic approximation to a three-point function, with three geodesics going from the boundary of AdS and meeting at a trivalent vertex [49]: Alternatively, the same result can be obtained from an appropriate limit of the DOZZ formula [50,51]. The function P is homogeneous of degree one in the dimensions, so gives a contribution scaling linearly with dimension in the exponential, as required. We will make particular use of the special case where all dimensions ∆ i are equal to ∆ p , for which P = − 3 2 log( 4 3 )∆ p . In the end, this gives the expressions for the blocks that we will use, applying in the limit 1 ∆ i c: We will apply this result in the specific case of the Z 3 symmetric, genus 2 handlebodies, with time-reflection symmetry, corresponding to the n = 2 version of the example in section 2.3, with x (or θ) real. The relevant geodesic network for the channel of interest is shown in red, in fig. 1. Furthermore we will take the dimensions of the three internal operators to be equal, ∆ 1 = ∆ 2 = ∆ 3 = ∆ p . In this case, finding the location of the vertices is straightforward, since they are fixed completely by symmetry, absolving us of the need to solve the minimization problem. It is now a simple exercise in hyperbolic geometry to work out the length of the geodesics connecting the two vertices, finding = log cot 2 θ 4 . This, along with P = − 3 2 log( 4 3 )∆ p as computed above, gives the result we will need for the block: The intuition behind the derivation of this expression relies on the operators in the internal channels being single trace operators, corresponding to single particle states, in a suppressed in ∆ relative to the single traces [48]. theory with semiclassical gravity dual. But because the blocks are kinematic objects, these restrictions are not required to apply the formula. We will use it in the case where the internal operators are highly composite multi-traces built from a primary of small dimension, for which the intuition behind the semiclassical blocks certainly does not hold.

Applying the blocks to generalized free fields
Consider a scalar O of dimension ∆, dual to a weakly interacting bulk field. As long as these interactions are unimportant, we can treat O as a generalized free field, which means that we can sensibly talk about composite 'multi-trace' operators built from products of O and derivatives, : ∂ # O · · · ∂ # O :. In the generalized free approximation, the dimensions of these products simply add, and they have vanishing connected correlation functions, so the correlators can be computed by Wick contractions. Now, let's try to compute the genus two partition function using the conformal block decomposition, accounting for such a free bulk scalar field. It is a slightly tricky prospect accounting for all the possible multi-trace exchanges, so we will make a slightly crude approximation, taking the contribution only of primary operators : O K : without derivatives, of dimension K∆, and also taking the same operator to propagate in all three legs. This gives us a lower bound for the partition function, since the OPE coefficients are real and the blocks are positive 17 : The OPE coefficients appearing can be computed from the combinatorics of the Wick contractions [36], and for three identical operators : O K :, in the limit of large K, the result is Putting this together with eq. (3.5) giving the blocks in the appropriate limit, the terms in the sum for large K look like But now, if ∆ is sufficiently small, these terms grow exponentially, causing the partition function to diverge! We can therefore put a bound on the critical dimension ∆ c at which Z GF (∆) diverges: We do not expect this to be exact, since we have dropped the contribution of so many operators, but we will see later that it becomes asymptotically equal to the correct value at small θ, corresponding to small x. For cross-ratio x = 1 2 , numerically solving the monodromy problem described around eq. (2.2), we find the corresponding value θ = .55128. This gives the bound ∆ c ≥ 0.189219 on the critical dimension (accurate to the number of quoted decimal places), in agreement with the analysis of [36].
A partition function should be well-defined for any surface, so it may seem puzzling to get a divergent answer. The resolution is that the partition function is not truly divergent, but our approximations on the spectrum and OPE coefficients do not apply when K is parametrically large. Even if we do not give a potential to the bulk field, it interacts through gravity, so the approximation of computing OPE coefficients of multi-trace operators by Wick contractions will cease to apply when K is of order √ c, though it could break down sooner if other interactions become important at a lower energy scale. When we pass the critical dimension, the sum over blocks will shift from being dominated by the vacuum, to being dominated by the multi-particle states at a scale set by the interactions. This signals a second order phase transition, which we will explain from the bulk as condensation of the scalar field.

The bulk instability
We have argued in Section 3 that the contribution of a free scalar to the genus two partition function will diverge for sufficiently small conformal dimension, ∆ < ∆ c . This divergence comes from the contribution of multi-trace states which are dual in the bulk to states with large particle number. It is therefore natural to expect that the divergence signals an instability where the bulk scalar field condenses to form a new solution with a nonzero classical value. This implies the existence of a second-order phase transition, below which the semiclassical bulk path integral is dominated by a new classical solution of the (nonlinear) bulk equations of motion: a 'hairy handlebody'.
The new classical solution will depend on the details of the theory, and in particular the interactions of the bulk field. These interactions give anomalous dimensions and couplings to the multi-trace operators in the theory, which become important above the scale of the interactions. 18 In particular, they will modify the asymptotic behaviour of the sum described in section 3 in such a way as to cure the divergence. The result is that the partition function will have some non-universal contribution at the interaction scale of the bulk field.
While the full nonlinear solution depends on details of the theory, the onset of the instability does not, and is sensitive only to the background geometry and the mass of the scalar field. It is characterized by the appearance at ∆ = ∆ c of a zero mode, a nonzero solution of the linearized bulk wave equation with source-free boundary conditions, which corresponds to a flat direction in the path integral. In this section, we will show that such a zero mode exists in quite general circumstances, and characterize the critical dimension ∆ c in terms of the bulk geometry.

The zero mode and spectral theory
In a d-dimensional holographic CFT, a single-trace scalar operator O of dimension ∆ is dual to a bulk scalar field φ of mass m 2 = ∆(∆ − d). The linearized bulk equation of motion (∇ 2 − m 2 )φ = 0 has two linearly independent solutions with different asymptotic behaviour near the boundary: Here J(x) is a source for the operator O in the CFT, which is fixed as a boundary condition, and the expectation value O(x) , in most circumstances determined uniquely by the boundary condition J and regularity, describes the 'response' of the scalar field in the presence of the linearized source J. We will be most interested in relevant operators ∆ < d, corresponding to masses which are naïvely tachyonic, but above the Breitenlohner-Freedman (BF) bound, − d 2 4 < m 2 < 0. In particular, we recall that for − d 2 4 < m 2 < − d 2 4 + 1, there are two possible choices of boundary condition for the scalar field φ with unitary duals [52]. These two different boundary conditions correspond simply to a choice of which of the two boundary behaviours in (4.1) one chooses to view as a source, and which as a response.
This linearized Klein-Gordon equation suffices to find leading order correlation functions (ignoring backreaction and other interactions), not just on pure AdS, but any asymptotically AdS geometry M obtained altering the boundary geometry or sourcing other fields, as long as φ = 0 on the background. 19 This corner of AdS/CFT therefore reduces to the theory of the Laplacian on the manifold M. Even for geometries that are locally AdS, this spectral theory can be rich and interesting, and we will import some ideas and results from the mathematics literature and explore the physical consequences.
Our key result is that, even in the absence of a source for the operator O, it is possible for φ to spontaneously acquire a nonzero classical expectation value. The second-order transition to this behaviour occurs when there is a nonzero solution of the bulk wave equation with vanishing source J = 0. In other words, as we vary the bulk solution M (or the dimension ∆), the solution will become unstable when there is an eigenfunction of the Laplacian with boundary condition J = 0 and eigenvalue ∆(∆ − d). For a given geometry M, we call the largest dimension for which this occurs the critical dimension ∆ c . Reducing ∆ further, this eigenfunction becomes a mode which decreases the action of the solution, so a given geometry is unstable to condensation of a scalar with ∆ < ∆ c . It is easy to see that the instability cannot happen if ∆ > d. In particular, for a scalar field with (∇ 2 − m 2 )φ = 0 we can use the standard argument for negativity of the Laplacian: In general, φ could have couplings to curvature or other nonzero fields which modify this linearized equation, but we will be interested primarily in geometries which are locally AdS. Thus all of these couplings may be incorporated into an effective bulk mass of the scalar field.
Here we have integrated by parts, and used the fact that the fast fall-off conditions (J = 0 and ∆ > d 2 ) imply that all integrals converge and boundary terms vanish. This instability is therefore ruled out for an irrelevant operator, but not immediately excluded for relevant operators, for which the boundary terms do not automatically vanish. We will see that such instabilities do occur, and are in fact quite generic.
Before giving our first characterization of the onset of the instability in terms of spectral theory, we should first clarify some mathematical terminology. We are seeking a solution φ of the equation ∇ 2 φ = ∆(∆ − d)φ, with boundary conditions J = 0. In the mathematics literature, this would be called an eigenfunction of the Laplacian only if ∆ > d 2 , since in this case the boundary condition J = 0 is equivalent to demanding square-integrability of the eigenfunction: φ ∈ L 2 (M). We will also be interested in the alternate quantization of the scalar field, corresponding to operators in the range 0 < ∆ < d 2 , where the boundary condition is imposed on the slowly-decaying mode. The dimension ∆ then corresponds to a resonance of the Laplacian, which is defined as a pole of the resolvent operator R ∆ = (∆(∆−d)−∇ 2 ) −1 , analytically continued in ∆. The resolvent is essentially the bulk Green's function (bulk-tobulk propagator) G ∆ on M; to compute the action of the resolvent on a function, integrate it against G ∆ , which satisfies the bulk wave equation with delta-function source: The critical dimension ∆ c will therefore show up as a pole in the bulk Green's function G ∆ . Another way to characterize the critical dimension ∆ c , is to note that the determinant of the resolvent det(∆(∆ − d) − ∇ 2 ) −1 is precisely the square of the one-loop partition function of φ. One can therefore find ∆ c by looking for a divergence in the one-loop contribution of a scalar field φ on the background M.
As a final characterisation of ∆ c , we can consider the linear response problem of turning on some small source J(x), solving the bulk wave equation with the corresponding boundary condition, and reading off the response O J . At generic values of ∆, this problem will have a unique solution, so defines a linear map S ∆ : J → O J between functions on the boundary B, known in the mathematics literature as the scattering matrix. If we tune ∆ to the critical dimension ∆ c , however, there is an ambiguity, as we can always add a multiple of the zero mode to the solution. The zero mode therefore also shows up as a pole of S ∆ , a scattering pole. 20 In the same way that the resolvent is related to the bulk Green's function, the scattering matrix is related to the CFT two-point function O(x 1 )O(x 2 ) B in the relevant background. This will also diverge as a function of ∆ as the critical dimension is approached (the familiar divergence in susceptibility at a second-order phase transition), with a pole at ∆ c , signalling the breakdown of the linearized bulk theory when the scalar becomes unstable. 20 Scattering poles do not coincide with the resonances for two reasons. The first is that there are also zeros of the scattering matrix, corresponding to solutions with a source but zero response, which may cancel a pole, giving a resonance without corresponding pole in S∆. Secondly, the scattering matrix has extra poles at half-integer values of ∆, even in pure hyperbolic space, related to the logarithms that appear in the boundary expansion 4.1 when the asymptotic powers differ by an integer, requiring additional counterterms. So far we have been quite general. We will now focus on the case of three dimensions, where we can make more concrete statements about ∆ c .

Locally hyperbolic spaces
Let us now consider the case where the bulk geometry M is a locally hyperbolic space of the form M = H 3 /Γ. We will be primarily interested in the case where M is handlebody, so we will take Γ to be a Schottky group of genus g > 1. In fact, the results of this section will apply in greater generality, to non-handlebodies, to some geometries containing conical defects, as well as to hyperbolic manifolds of general dimension. 21 We will consider a bulk scalar propagating on this geometry, and characterize the relevant spectral theory in terms of properties of the quotient group Γ. We will only motivate and explain the results here, referring to the appropriate mathematics literature for more details, precise statements, and proofs.
Consider first computation of the bulk two-point function of φ in the geometry M = H 3 /Γ. This can be computed using the method of images, by starting with the two-point function in H 3 and summing over all elements of the group, corresponding to sources at all image points. The result is where d(y, y ) is the geodesic distance between the points y and y , with respect to the H 3 metric. Formally, this gives a function invariant under the group Γ, hence well defined on M, and solves the Klein-Gordon equation with the appropriate source. However, this function will not be well-defined if the sum over images fails to converge. In particular, if the number of image points with d(y, γ · y ) less than some distance d grows rapidly enough as d → ∞ (for some fixed y, y ), then the sum will diverge. More specifically, if the number of image points with d(y, γ · y ) < d grows like e δd , then the sum will diverge for ∆ < δ. In fact, this is always the case for some δ > 0, as stated in the following result of Sullivan [53]: converges in the right half-plane Re ∆ > δ, where δ > 0, the exponent of convergence of Γ, is the location of the first resonance of H 3 /Γ. The Green's function G M ∆ (y, y ) (analytically continued in ∆) has a pole at ∆ = δ, and the residue of that pole is given by where φ 0 is the zero mode function, the solution of the free bulk wave equation with source-free boundary conditions.
As described in the previous section, this pole in the Green's function, the resonance, signals the onset of an instability. Thus the critical dimension ∆ c equals the exponent of convergence δ. We emphasize that δ is strictly positive given our assumptions on M, which implies that any handlebody of genus greater than one will be unstable if there is a sufficiently light operator in the spectrum.
We may also compute the CFT two-point function of O in this background, by taking the limit of the bulk Green's function as the points approach the boundary. The exact result for the two point function will depend on the conformal frame, which corresponds to a choice of regulator as we take the points to the boundary, but the convergence properties of the sum over images will be insensitive to this choice. We can write a general metric on the boundary as ds 2 = e 2σ(w) dwdw, where w is the complex coordinate on which Γ acts by Möbius maps. The conformal factor σ is defined on the regular set Ω of Γ, and defines a metric on the quotient manifold B = Ω/Γ under the condition σ(γ(w)) = σ(w) − log |γ (w)| for all Möbius maps γ ∈ Γ. The bulk computation of the two-point function on B gives a sum over images: In the summation on the right hand side, the denominator corresponds to the two-point function on the plane, and the numerator is the conformal factor appropriate for each image. Once again, this sum converges in the right half-plane Re ∆ > δ, and the divergence in the two-point function signals the onset of a second-order phase transition. In the mathematical literature, this CFT two-point function appears as the kernel of the scattering matrix [54]. This sum can be used to gain some intuition about the relationship between the exponent of convergence δ and the geometry of the group Γ. The first thing to note is that the tail of the sum, which controls the divergence, is closely related to the limit set Λ of the group Γ. The limit set is the set of points where the images γ(w) accumulate, for any starting point w. More precisely, a point is in Λ if every neighbourhood of that point contains infinitely many of the images γ(w). These are the places where the quotient by Γ acts 'badly', which we must remove to form the regular set Ω = C * − Λ, so we obtain the nice quotient space B = Ω/Γ. The tail of the sum is controlled by the limit set, since only a finite number of terms in the sum will lie outside of any arbitrarily small neighbourhood of Λ. In the simple case of BTZ, Γ is the cyclic group consisting of the maps γ n (w) = q 2n w for n ∈ Z, and Λ consists of the two points 0 and ∞. In most other cases, however, Λ is much more complicated.
For any limit point (that is, element of the limit set), there is a sequence of images of our starting point w that approach it, say γ n (w) for some γ n ∈ Γ (which are independent of w). As n increases, γ n will usually be a longer and longer word built out of the generators. For the images to tend to some limit, the γ n must eventually start with the same string of generators, because if they don't, they would map w to places that are separated by some (a) q = .8 + .3i, δ ≈ 1.09 (b) q = .8 + .44i, δ ≈ 1.29 Figure 2: The limit sets for two of the Z 3 symmetric genus two Schottky groups that arise when investigating n = 3 Rényi entropies. The parameter q defining the groups is an eigenvalue of one of the generators as specified in section 2.3. We give the value of the Hausdorff dimension δ for these two limit sets, computed using the methods of section 5.
finite distance: as the sequence γ n goes on, the words built out of the generators get longer and longer, and only change later and later on in the string. More precisely, the kth letter of the word γ n is constant after some sufficiently large n, for any k. For each limit point, we can in this way construct a unique semi-infinite word built from the generators of the group, a sort of decimal expansion, but using Möbius maps instead of digits. Such words are in one-to-one correspondence with elements of Λ. For g ≥ 2, this set is not only infinite, but uncountable. The 'rational numbers' in the analogy with decimals consist of strings of generators that eventually repeat, and are in one-to-one correspondence with the primitive elements P of the group Γ, that is, elements that cannot be written as γ n for any n > 1 (excepting the identity), corresponding to the attractive fixed point of that element.
The resulting set Λ, which controls the tail of the sum over the group, has a rich and beautiful fractal structure. For Fuchsian groups, generated by matrices with real entries, the limit set is a subset of the real line, and closely resembles (indeed, topologically, is homeomorphic to) a Cantor set. Allowing more general Schottky groups, the limit set moves into the complex plane, forming a twisting, intricate, self-similar pattern. Several examples arising when investigating the 3rd Rényi entropy of two intervals are illustrated in fig. 2. For many more images of limit sets, and a playful semi-popular account of the mathematics involved, we encourage a foray into [55].
Secondly, we note that the size of the terms in the sum (4.8) is controlled primarily by the factor |γ (w)| ∆ , which describes how things scale under the action of γ (in the flat or round metric on the Riemann sphere, not the metric pulled back from B). Given some small set near w, the characteristic length of its image under γ is scaled by |γ (w)|, and its area is scaled by |γ (w)| 2 , so it is natural generalise this, and say that |γ (w)| ∆ characterises the scaling in a ∆-dimensional notion of measure, where ∆ can be any positive real number. The convergence of the sum is therefore determined by the trade-off between the accumulation of many points at the limit set, and the shrinking of ∆-dimensional measure associated to images at those points. The critical dimension will occur when these two effects precisely balance, which is when the limit set itself can be assigned a ∆-dimensional measure invariant under Γ. Hopefully this discussion makes plausible the following theorem of Patterson [56], Sullivan [53,57] and Bishop-Jones [58], the precise statement of which uses the notion of Hausdorff dimension, a non-integer dimension defined for fractals in metric spaces.
Theorem 2 (Patterson-Sullivan). The exponent of convergence δ is equal to the Hausdorff dimension of the limit set Λ of Γ: This result connecting spectral theory and fractal geometry is certainly beautiful, which would be justification enough to include it in a mathematics paper, but amazingly enough it is also useful. Firstly, it gives us a new tool to intuit how the critical dimension ∆ c depends on the geometry, particularly in certain limits. But more importantly, it provides a method to accurately and efficiently compute δ for any group Γ, which is far better than naïvely solving the bulk Laplace equation numerically, or the method of extracting δ from the asymptotics of the terms in the sums introduced above. We will discuss an algorithm to compute δ in section 5.1, and use it to present both numerical and analytic results.

Divergence of the partition function
In this section we will offer one final perspective on the phase transition, to make a direct connection with the CFT argument discussed in section 3. In that section we summed up the contributions to the partition function from of a generalized free field, using the global limit of higher genus blocks with the spectrum. From the bulk point of view, this object is precisely the one-loop partition function of the bulk scalar field φ: This makes it apparent that the zero mode should again be visible as a zero eigenvalue of the operator ∇ 2 − m 2 (defined with suitable boundary conditions). In this way, the calculations of section 3 put a lower bound on ∆ c . The bulk computation of this one-loop partition proceeds much as the Green's function computation given above. In particular, one can compute this one-loop determinant using heat kernel methods and a sum over images [59]. For higher genus surfaces this one-loop determinant can be written as an infinite product Here, q γ is the smaller eigenvalue of γ, as previously introduced. The product is over primitive conjugacy classes γ ∈ P of the group Γ; these are conjugacy classes of elements which cannot be written as a power γ n of another element with n > 1, and q γ is the smaller of the eigenvalues of γ when written as an SL(2, C) matrix. Note that this definition counts γ −1 separately from γ, so that terms come in matching pairs. We have written the one-loop partition function terms of the Selberg zeta function ζ Γ associated to the group Γ, as defined in [54] 22 . We note that the product in (4.11) converges in the same right half-plane Re ∆ > δ as the image sums we have already introduced. In fact, the Selberg zeta function can be analytically continued to an entire function, with zeros precisely at the eigenvalues and resonances of the Laplacian on M, as expected [60]. In particular, the first resonance, corresponding to the phase transition of interest, leads to the one-loop partition function diverging as (∆ − ∆ c ) −1/2 . This is the divergence found from the CFT analyses of [36] and section 3.
The product (4.11) has a simple geometric interpretation in terms of the closed geodesics on the bulk manifold M = H 3 /Γ. Since Γ is the fundamental group of M, its conjugacy classes are in one-to one correspondence with homotopy classes of closed loops in the bulk, and in a hyperbolic manifold, there is a unique closed geodesic in each class. The primitive conjugacy classes P correspond to prime geodesics that do not trace over their image multiple times. The geometric parameters associated to a closed geodesic are its length γ , and its twist θ γ , the angle through which a normal vector gets rotated after being parallely transported around the curve, and are related to the associated conjugacy class of Γ by q 2 γ = e − γ +iθγ . The convergence of the product is therefore controlled by the asymptotics of the length spectrum of the bulk manifold. A precise statement of this is given by the prime geodesic theorem, so called because of its close analogy with the prime number theorem (provable using the analytic properties of the Selberg and Riemann zeta functions respectively): Theorem 3 (Prime geodesic theorem). The prime geodesic counting function π M ( ), defined as the number of prime geodesics of length at most , satisfies the asymptotic formula π M ( ) ∼ e δ δ as → ∞. (4.12) In this way, the instability of the scalar field is controlled by the asymptotic properties of the spectrum of very long geodesics. This relation between the spectrum of the Laplacian and the lengths of closed geodesics is a special case of the Selberg trace formula (or an appropriate generalization).
Consideration of the partition function leads to an alternative approach to the computation of δ, which we will not pursue further here, by numerically computing the Selberg zeta function, which can be done efficiently (though not directly from the product definition), and locating its zeros. 22 There are several closely related definitions of the Selberg zeta function. The definition we have given is appropriate for hyperbolic three-manifolds where Γ is a Kleinian group; another, more common definition is in the context of hyperbolic surfaces, where Γ is a Fuchsian group, so qγ is real, and the product overn is absent.

When Γ is Fuchsian
In the case when the group Γ is Fuchsian, i.e. when all elements are in SL(2, R) and so fix the real line (perhaps after conjugation with some Möbius map, for example fixing the unit circle instead), the discussion simplifies somewhat. Instead of requiring the full three dimensional geometry, all the main results discussed here can be reduced to the twodimensional slice Σ fixed by complex conjugation. In this section we briefly describe this reduction and its consequences.
Fuchsian groups are, in many circumstances, the most physically interesting cases, primarily because they correspond to geometries that have a real Lorentzian description. Interpreting the action of complex conjugation as a time-reversal symmetry, the slice Σ fixed by time-reversal has vanishing extrinsic curvature, and hence can be interpreted as an initial Cauchy surface for Lorentzian evolution. Very explicitly, the Euclidean bulk M can be written as where dΣ 2 is the hyperbolic metric on the χ = 0 slice Σ. The Lorentzian geometry (or, rather, a patch of it) is obtained by analytic continuation χ → it. This gives a locally AdS 3 solution to the equations of motion in an FRW-like coordinate system, where the spatial slices have constant negative curvature.
It is important to note that while all Fuchsian groups have a reflection symmetry, and corresponding Lorentzian interpretation, the converse is not true: a non-Fuchsian Schottky group may have a time-reflection symmetry and good Lorentzian continuation. To take one example, the pure entangled state on three copies of the CFT obtained by the path integral on a pair of pants is, for certain moduli, dual in the Lorentzian section to disconnected copies of pure AdS and BTZ [19], but the corresponding (connected) Euclidean geometry is not described by a Fuchsian group.
The bulk metric is not static, so to simplify the Laplacian it is not as straightforward as choosing a time-independent ansatz. But it is not much harder than that; instead, look for a separable eigenfunction F (σ, χ) = f (σ)g(χ), finding that if f is an eigenfunction of the Laplacian on the χ = 0 slice with eigenvalue ∆(∆ − 1), and obeys the correct AdS boundary conditions, then F (σ, χ) = (sech χ) ∆ f (σ) (4.14) is an eigenfunction of the full handlebody Laplacian with eigenvalue ∆(∆−2), with the correct boundary conditions. From this, the critical dimension of the handlebody is determined by the bottom of the spectrum of the slice Σ, and computing the actual profile of the zero mode is no longer a three-dimensional problem.

Results for the critical dimension ∆ c
We have seen that a scalar field on a handlebody H 3 /Γ will be unstable if the dimension is sufficiently small: ∆ < ∆ c . We now turn to an explicit computation of the critical dimension ∆ c , which will be a function of the moduli. A direct approach, where one studies the Laplacian directly on the geometry of interest, is a complicated numerical task. Our approach will be to instead use theorem 2 to calculate ∆ c from the Hausdorff dimension of the limit set of Γ.
Using this, we will obtain analytic results for ∆ c near the boundary of moduli space, as well as analytic bounds on ∆ c in the interior of moduli space. We will also obtain accurate numerical results. Our main tool will be an algorithm due to McMullen [61], which we now describe.

McMullen's algorithm
This section is somewhat technical and is not necessary to understand the results described in later sections. Readers who are not interested in the details of how the results are obtained can safely skip to section 5.2.
To begin, will we need to introduce an additional structure on the limit set Λ: a Γinvariant δ-dimensional measure µ. A measure µ (though we will not give a precise definition here) allows us to integrate functions on the limit set, in particular assigning a number µ(E) = E dµ ≥ 0 to subsets E ⊆ Γ, providing a measure of the 'content' of E. We require the additional property that it transforms as a δ-dimensional density under element of the group 23 : A nontrivial measure with this property exists when, and only when δ equals the Hausdorff dimension of Λ (in which case it is unique, up to normalization, for Schottky groups). The only feature of the right hand side that we require is that it is bounded by the measure of the set µ(E), times the extrema of the integrand |γ | δ : McMullen's algorithm works by splitting the limit set into a finite number of disjoint pieces E i , and attempting to approximate (or bound) the value of the measure µ on each of these pieces, µ i = µ(E i ). For a detailed explicitly worked example of this and the following, see fig. 3. We begin by imposing (5.1). Specifically, suppose we have some Möbius map γ i ∈ Γ, and one of the pieces E i , whose preimage under γ i is the union of some pieces E j 1 , E j 2 , . . . E jn (one of which may be E i itself): If we pick some points z j ∈ E j , then (5.1) implies that 23 Here, |γ | computes the local scaling of lengths under the map γ; we may use any metric on the boundary Riemann sphere for this purpose, for example a round metric, and the result for δ is insensitive to this choice. For practical computations the flat metric is often most convenient, in which case |γ | is the absolute value of the derivative of the Möbius map. It is then simplest to require that the point at infinity is not in Λ, which we will implicitly assume. Figure 3: An example of computing the transition matrix for McMullen's algorithm, in the case of a Kleinian group freely generated by two loxodromic elements g, h, so that C * /Γ is a genus two surface. In the figure, we have drawn a fundamental domain for Γ, the exterior of the four outermost circles (those corresponding to g −1 , h −1 are not shown in their entirety). Break the limit set into the four pieces E γ contained in each of these circles, labelled by γ = g, h, g −1 , h −1 corresponding to the element of the group that maps the fundamental domain to the interior of the circle, and choose points w γ ∈ E γ , for example the attractive fixed point of γ. The piece of the limit set E g can be broken up into three disjoint pieces, inside the circles labelled g 2 , gh and gh −1 , which are the images under g of E g , E h and E h −1 respectively. The scalings of these limit sets under the action of g go into the top row of the transition matrix: The other three rows repeat the same exercise for the other three regions, and finding δ such that the spectral radius of T δ is unity gives an approximation for the Hausdorff dimension. This can be refined by breaking the limit set up into the 3 × 4 n−1 regions E γ labelled by words of length n in g, h, g −1 , h −1 , and applying the δ-invariance imposed by considering the preimage of E γ under the first element (g, h, g −1 , or h −1 ) appearing in the word γ. Then T will be a sparse matrix, with three nonzero elements in each row and column, and the algorithm has error decreasing exponentially with n. The figure includes labels for words of length two, but also shows the images of circles under words of length three (unlabelled).
This is not an exact equality because the scale factor |γ i | is not constant on the limit set.
However, by taking the E i to be small |γ i | will be approximately constant, so the error will be small. To be more precise, we may replace the factors of |γ i (w j )| by upper or lower bounds on this scaling over the set E j , and use eq. (5.2) to replace the approximate equation by inequalities.
With an appropriately chosen partition {E i } and maps γ i , a similar argument can be repeated for every i. Our approximate formula (5.4) can then be written in terms of a square matrix T , the transition matrix, whose entries T ij are equal to |γ i (z j )| for each of the E j in the preimage γ −1 i (E i ), and zero otherwise. Invariance of the measure is then the statement that µ i is a unit eigenvector of T δ , where the power is taken element-wise. Since the matrix T δ has nonnegative entries it is guaranteed to have a unique eigenvector with positive components; furthermore this is the eigenvector with largest eigenvalue. 24 If µ is to satisfy the invariance criterion, this eigenvalue should be one. Thus to find the Hausdorff dimension, we find the value of δ such that the largest eigenvalue of T δ equals one. 25 An analogous result holds when we replace the approximate equations by inequalities, so by choosing upper or lower bounds on the transition matrix elements, we can obtain rigorous bounds on δ.
For a given partition {E i } of the limit set this gives an estimate for δ, and to obtain a more accurate estimate we can refine the partition into a larger number of pieces. With an appropriate refinement, the result converges rapidly to the Hausdorff dimension, and in practice it is sufficient to use a rather coarse partition of Λ. Although this discussion is rather abstract, the explicit implementation of this algorithm is quite straightforward; see figure 3 for a simple example.

Analytic results
Our first analytic results are for Fuchsian groups; this includes the surfaces described in section 2.3 with real cross-ratio 0 < x < 1. In this case we note that the limit points must all lie on the real axis of the w-plane, corresponding to the slice fixed by time-reflection symmetry. Since Λ is a subset of the one-dimensional line, it must have dimension δ ≤ 1. Thus ∆ c ≤ 1. So the only potentially unstable fields are those with ∆ < 1, which correspond to bulk scalars which are quantized using alternate boundary conditions.
In the rest of this section we will focus on the case of the genus g = n − 1 surfaces with Z n symmetry, described in section 2.3, relevant for computing the nth Rényi entropy of a pair of intervals. We begin by applying the above algorithm in this case at the coarsest level of approximation, to obtain analytic bounds. These are especially useful at the edge of moduli space where x → 0, because in this limit the different pieces of the limit set become well-separated and small, so the scale factor does not vary much over it. These bounds thus become tighter and tighter as x → 0.
As described in section 2.3, we can extend the group Γ toΓ, generated by R : w → e 2πi/n w and S : w → w−ζ ζw−1 , by including a dihedral group of holomorphic automorphisms. This extension of the group does not alter the limit set, and a Γ-invariant measure constructed on it will also be invariant underΓ. Using this, we will divide the limit set into n pieces, all related by the Z n symmetry R, and hence having equal measure, and use the mapping under S to constrain the dimension of this measure. This is somewhat simpler than using the original presentation of Γ such as in fig. 3.
There are n pieces of the limit set E k , each centred at a root of unity e 2πik/n , with size of order θ 2 for small θ, and related to each other by R. A simple way to show this is by constructing a fundamental domain forΓ, bounded by the radial lines from the origin at angles ±πi/n, related by action of R, and a circle C n mapped to itself by S, centred at sec θ with radius tan θ (recall ζ = cos θ). Then E n is the part of the limit set inside this circle, and the remaining E k are inside corresponding circles C k = R k (C n ) obtained by rotations by angle 2πk/n. This immediately bounds the size of the limit set by the size of the circles C k , of order θ, which is a sufficiently strong result for our immediate purposes 26 .
By the Z n symmetry, and the fact that the action of R doesn't scale (|R | = 1), the sets E k must all have equal measure µ k = µ. We can then apply the action of S, which maps E 1 , E 2 , · · · , E n−1 onto E n . The amount by which E k scales under S can be computed from This requires δ to tend to zero as θ goes to zero, and solving to leading order for small δ, we get the result that with the first corrections coming from solving the equation to higher order in δ, rather than the order θ 2 variation of the scaling relation giving the correction to eq. (5.6). With only minor modifications, this derivation continues to apply if we allow θ to be complex, which is why we have included the modulus in the result. The case n = 3, with real θ, was treated in [61], though instead of using S, that paper uses a reflection, which requires θ to be real 27 .
To facilitate comparison between different values of n, we write this in terms of the cross-ratio x. With the monodromy methods outlined in section 2.3, the map from x to the 26 To obtain an improved bound, note that En must be contained in a smaller set, namely the union of the interiors of S(C k ), which is concentrated in a region with size of order θ 2 .
27 For comparison, the parameter θ used in that paper is half of the θ used here.
Schottky parameter θ can be computed as a series expansion, with the leading order result that θ = √ x n (1 + O(x)): Note that this is the asymptotic behaviour for fixed n as x → 0, but clearly must break down if n is parametrically large. In the first instance, it is not self-consistent if nx is of order one, since in that case the leading order term in the expansion would be of order one. But we can take a different order of limits to understand what happens at large n.
Starting with the scaling relation eq. (5.6) at small fixed x, and naïvely taking the large n limit term by term, we arrive at where the factor of two is to count contributions both from fixed k and fixed n − k. The first thing to notice is that the tail of the series decays as k −2δ , so convergence of the sum immediately requires δ > 1 2 . Summing the series, we arrive at from which we see that the left hand side is large for small x, since δ cannot be small. The zeta function must therefore be close to the pole at δ = 1 2 , and we can find a perturbative solution: This result should be interpreted as the limit of δ as n → ∞, for fixed but small x. It turns out this naïve argument is, in essence, correct, and can be made completely precise by repeating the argument for the group generated by one parabolic element and one elliptic element of order 2, equivalent to Theorem 3.6 of [61] (for real cross-ratio). In fact, the result that the Hausdorff dimension does not go to zero as x → 0 is a consequence of a general result, that δ > 1 2 whenever the group in question contains a parabolic element (Corollary 2.2 in [62]). This Schottky group with R parabolic, instead of elliptic order n, corresponds to the 'n = ∞' version of the geometry described in section 2.3. This complex one-dimensional family of Kleinian groups is known in the mathematics literature as the 'Riley slice' of Schottky space. It is a little tricky to think about the n → ∞ limit of the handlebody, bounded by a Riemann surface of infinite genus, but it is rather simpler to understand the geometry after taking a quotient by the Z n replica symmetry, as suggested in [30]. The boundary of this geometry is just the original Riemann sphere, and the bulk has conical defects, of opening angle 2π/n, going from 0 to x and from 1 and ∞. Taking a formal analytic continuation of the geometry to n = 1, the conical defects become the Ryu-Takayanagi surface, lying on geodesics [30], but we are taking the opposite limit, in which the defects become cusps, in particular receding to infinite proper distance.
At this point, let us pause briefly to understand the physical consequences. This result means that the x → 0 and n → ∞ limits of the critical dimension do not commute, so that while for any fixed n, any dimension of scalar will be stable for sufficiently small x, if there is a scalar of dimension less than 1/2, it will be subject to the phase transition for any cross ratio, if n is taken sufficiently large. Note that this is all in a limit where we have taken c to infinity first, and new behaviour dominated by quantum corrections may take over when n is parametrically large in c. In particular, the large n limit of the Rényi entropy is controlled by the largest eigenvalue of the reduced density matrix, or the 'ground state energy' of the modular Hamiltonian H A = − log ρ A , but it is unclear whether the semiclassical description is sensitive to a single lowest eigenvalue, or a dense collection of parametrically many low-lying eigenvalues of H A .
To conclude the discussion of analytic results, let us briefly describe the other limit of the geometry, when x → 1, corresponding to the horizon sizes of the multi-boundary wormhole becoming small. This limit is of less direct physical interest, since it is well past the first-order phase transition of the partition function, so is not the dominant saddle-point geometry, but is nonetheless useful to hone intuition. In the case that x is real, so Γ is a Fuchsian group, it is sufficient to describe the bulk geometry in terms of the two-dimensional hyperbolic surface making up the t = 0 slice. In this limit, it is helpful to separate the geometry into the exterior pieces, lying between each boundary and a horizon, and the 'convex core' or 'causal shadow' region linking them together, bounded by the n horizons. In the x → 1 limit, the horizons become very small, and the centre of the geometry recedes down a long, narrowing neck. The core then approximates the geometry of the negatively curved metric on some compact surface with punctures, though the punctures do not quite pinch off, rather reaching a minimum radius at the narrow horizons where they join to the exterior funnels. On such a surface, the critical dimension approaches one, the maximal eigenvalue of the Laplacian (on the t = 0 slice, not the Laplacian in the three-dimensional bulk: see section 4.4) being small and positive. The corresponding eigenfunction is roughly constant on the convex core, and small in the exterior funnels, with an interpolation over the long, narrow necks connecting them. The asymptotic behaviour of ∆ c in this limit can be computed by directly approximating this Laplace eigenfunction (the zero mode of the instability) [63].
Taking x → 1 through complex values is much more complicated, so we will not be able to say much about it. Schottky space, perhaps parameterized in this case by the values of ζ corresponding to some complex x, itself has a complicated fractal boundary, and the features of the handlebody depend sensitively on how this boundary is approached. This is a deep and beautiful subject, but goes far beyond the scope of this article. In any case, the numerical computations we describe next show that it is possible to obtain dimensions δ > 1 in this limit, for example the limit sets illustrated in fig. 2.  Figure 4: The critical dimension ∆ c = δ as a function of cross-ratio x for the handlebodies corresponding to the Rényi entropies of a pair of intervals. From top to bottom, the curves correspond to genus 2, 3, 4, 5, 6, and finally the n → ∞ result in black. The shading visible on the right side of the plot indicates the bounds achieved by applying McMullen's algorithm at the crudest level of approximation.

Numerical results
These analytic results are very useful to understand the behaviour of the critical dimension at the edges of moduli space, and as genus is varied, but McMullen's algorithm is also useful to quickly compute the Hausdorff dimension numerically, to many digits of precision. We conclude by presenting the results of these computations.
Firstly, fig. 4 plots the Hausdorff dimension as a function of cross-ratio x for various values of n, including also the limit as n → ∞. With this parameterization, for generic values of x, the convergence of the algorithm is remarkably rapid. Indeed, the plot includes shaded regions to indicate the rigorous bounds obtained by applying the algorithm at the crudest level. These are computed by numerically solving the equation eq. (5.6), but replacing the terms on the left hand side with upper or lower bounds for |S (w)| δ over w ∈ E k , rather than the estimates used there. The allowed regions for δ are in many cases not even visible until x is rather close to 1. Refining further, the algorithm gives results with ten or more digits of precision in a fraction of a second on a laptop of modest specifications. In fact, by far the larger source of computing time and error comes from the conversion between the cross-ratio and Schottky variables, rather than the algorithm to compute δ from the group generators.  Figure 5: The critical dimension for the x = 1 2 Rényi surface as a function of replica number n. The asymptote is the computed limit as n → ∞. On the right is a log-log plot showing convergence to this value.
A physically motivated value to consider is at the boundary with the Hawking-Page phase transition x = 1 2 , which will give the maximum value of ∆ c within this class of geometries, while in the dominant phase. At genus 2 (n = 3), this value is ∆ c = 0.189124003, which is rather close to (and the correct side of) the bound ∆ c ≥ 0.18912109 obtained in section 3 from refining the CFT methods of [36]. Staying at x = 1 2 and increasing the genus, we find that the critical dimension increases rapidly at first, before slowly approaching the limiting value ∆ c → 0.599 as n → ∞, as shown in fig. 5.
In our last plot, fig. 6, we indicate how the Hausdorff dimension behaves for complex values of the cross-ratio, for genus two. Note that this is invariant under inversion in the circle of unit radius centred at one. This is because the extended groupsΓ corresponding to these geometries are the same, though Γ consists of different subgroups in each case. From the geometric point of view, taking the Z n quotient of the handlebody gives the same geometry, with conical defects at the fixed points, though the original geometries are distinct (being branched around the defects in different ways). This relates a cross-ratio 0 < x < 1 with a negative cross-ratio − x 1−x , which corresponds to swapping the location of twist and antitwist operators, relevant for computing Rényi negativity of two disjoint intervals [64]. The correspondence between the geometries implies a correspondence between the classical limits of Rényi entropy and Rényi negativity for two intervals, and taking an analytic continuation to n = 1, the logarithmic negativity of two intervals must vanish (to leading order in c) in the regime x < 1 2 where this geometry dominates the path integral. Finally, it is interesting to ask what the largest possible value of ∆ c could be for a geometry that dominates the path integral. A lower bound (precluding surprising new symmetrybreaking phases that dominate the path integral) comes from our numerics for the infinite genus limit, giving examples where ∆ c as large as .599 can be achieved. An interesting result that may bound this in the other direction comes from [65], showing in particular that every Riemann surface admits a uniformization by a Schottky group of Hausdorff dimension less than one. As a heuristic, matching our expectations in limits of moduli space, the dominant Figure 6: The Hausdorff dimension of the Z 3 symmetric genus two handlebody, as a function of the (complex) cross-ratio x. Note that the Hausdorff dimension goes to zero at the origin (x = 0) and approaches one as x → 1 − along the real axis.
saddle-point seems to be that with minimal Hausdorff dimension, so this is suggestive, though not conclusive, that there may never be a dominant geometry with ∆ c > 1.