Free Energy from Replica Wormholes

Euclidean wormholes -- geometries which connect disconnected boundaries -- present a challenge to a standard quantum mechanical interpretation of the theory. One potential resolution is that the gravitational path integral computes the ensemble average of many theories. The connected topologies contribute to the simplest possible observable: the free energy, which is computed using a replica trick. This is distinct from the replica trick used to compute entanglement entropies, and appears in the computation of any extensive quantity. We argue that both JT gravity and a simplified version of CGHS admit a regime where the contribution of connected replica wormholes to the free energy is larger than that of disconnected topologies. In both theories we find evidence of replica symmetry breaking, which is reminiscent of the behavior of certain spin glasses. We discuss possible insights about ensemble averaging in gravity from this perspective.


Introduction
The process by which information escapes from the black hole hole interior is a pivotal question in the study of quantum gravity. Recent work has brought to light two important points on this front. The first is that at least in some gravitational models the Euclidean gravitational path integral (GPI) exhibits traces of unitarity: a GPI calculation of the entropy of the Hawking radiation reproduces the unitary Page curve [1][2][3][4][5][6][7][8][9]. This hinges crucially on the contribution from Euclidean replica wormhole saddles that connect disconnected boundaries. The inclusion of such wormholes implies that absent some further UV effects, the GPI would not factorize across disconnected boundaries [10][11][12][13][14]. In this case, the GPI cannot be interpreted as computing the partition function of a standard quantum mechanical theory. One possible explanation is that the GPI should instead be interpreted as computing the ensemble average of many different quantum theories (see also [15][16][17][18] for related discussions). Another possibility is that additional contributions should be included, which would lead to the expected factorization of the Euclidean partition function on disconnected surfaces.
Our goal in this paper is to understand more systematically how Euclidean wormholes influence the physics of the GPI. We will put aside for the time being any further potential UV effects (such as certain doubly non-perturbative effects in JT gravity) which might be necessary to describe the dual of an individual quantum theory. We will investigate the contribution of Euclidean wormholes to a more general -and in a sense simpler -class of observables than the entropies described above. We find that, completely independently of any considerations of black hole physics, these wormholes make important (and apparently indispensable) contributions to the dynamics of the theory.
To understand how Euclidean wormholes contribute, let us imagine computing a Euclidean GPI where we sum over geometries with a particular choice of boundary B: (1.1) We will take B to be a connected surface, so this is usually interpreted as giving the gravitational computation of a partition function Z(B). One can also consider the integral over geometries with boundary B m = B ∪ · · · ∪ B: If Euclidean wormholes contribute, then P(B m ) = P(B) m and the resulting partition function (or, more generally, correlation functions) do not factorize. One potential interpretation is that the GPI computes ensemble averages: where the overline denotes the average over a family of unitary quantum theories, and Z(B) the partition function of a member of this family 1 . Here we will remain agnostic on whether this ensemble average is genuinely a feature of the gravitational theory, or whether it merely appears as an approximate contribution to the low-energy effective description of some UV-complete theory. Nevertheless, we will continue to interpret the GPI as an ensemble average, bearing in mind that this interpretation may only be valid in some effective description. We will revisit these issues in more detail in Section 6.
Our first observation is that if Euclidean wormholes contribute to the GPI, then they should contribute to even the most basic observable of the theory: the free energy F = −T ln Z evaluated at a particular temperature T . In particular, let us imagine computing the free energy via a GPI, where T enters through the choice of B (for example, in a two-dimensional theory of gravity, B is a circle of length β ≡ 1/T ). Naïvely, of course, one might try to compute it by simply taking F = F ann ≡ −T ln P(B) = −T ln Z(B). (1.4) This, however, is in tension with the ensemble interpretation: since Z(B) involves an integration over the random variables defining a particular instance of the ensemble, we may interpret Z(B) as the partition function of a theory in which the random variables themselves are permitted to fluctuate and come into equilibrium. In condensed matter systems, the free energy F ann defined above is therefore interpreted as an annealed free energy. Instead, what one is really interested in is the quenched free energy, in which the random variables defining a particular instance of the ensemble are not allowed to equilibrate. In other words, the free energy F = −T ln Z(B) is computed in a particular instance of the ensemble, and then the average is taken: (1. 5) In general the annealed and quenched free energies will be different. Indeed, from the gravitational point of view one might expect that ln Z(B) = ln Z(B) whenever Euclidean wormholes are present in the theory, for the same reason that Z m = Z m .
In order to understand exactly how Euclidean wormholes contribute to (1.5), one needs to compute F from the GPI using a replica trick that involves considering the GPI on m copies of the boundary B and then analytically continuing to m = 0. This replica trick is distinct from the one that is employed to compute the von Neumann entropy (which instead considers the GPI defined by an n-sheeted boundary manifold and then continues to near n = 1), and a completely consistent calculation of entanglement entropy must implement both replica tricks. This version of the replica trick will be reviewed in section 2, and is common in the condensed matter literature, especially in the study of spin glasses. In fact, although we have focused on the free energy, this new replica trick will apply to the computation of any extensive observable. For example, in the calculation of the Renyi entropy S n of a pure state from the GPI in [5], the result vanishes only to leading order if this additional replica trick is not implemented: the Renyi entropy vanishes identically only when the calculation correctly implements both replica tricks.
This additional replica trick means that F becomes sensitive to the contribution of wormholes connecting the replicas, and leads to the conclusion that it is not consistent to simultaneously interpret P(B) as computing an ensemble average and to compute the free energy (or more generally, any extensive obervable) without including contributions from Euclidean wormholes. If the free energy computation is dominated by the disconnected topology, then the ensemble averaging leaves no visible footprint, and the quenched free energy coincides with the annealed free energy: F ≈ −T ln P(B). However, if in some regime replica wormholes contribute nontrivially to F , then ensemble averaging is important for the computation of any observable in that regime. Failure to properly compute the free energy via the replica trick above will erase subtle signatures of the ensemble.
Of course, the skeptical reader may be concerned that replica wormholes might never actually make an appreciable contribution to F , at least in those regimes in which we have some control over the gravitational theory. Indeed, although it has now been verified that replica wormholes are important in the study of black hole entropy, it need not follow that such wormholes will be important in the computation of F .
To address this potential concern, in Sections 3 and 4 we compute the free energy in two different models of 2D gravity. We find that the naïve calculation of the annealed free energy F ann exhibits pathological behavior at sufficiently low temperature. Specifically, it is non-monotonic with temperature, implying a negative thermodynamic entropy S = −∂F/∂T . We then use the replica trick to investigate the contribution of replica wormholes to F , finding that this contribution becomes larger than that of the disconnected topology when the annealed free energy exhibits its unphysical behavior. The inclusion of wormholes ameliorates the pathological behavior of the free energy at low temperature, at least with a certain implementation of the replica trick.
The gravitational systems that we consider are CGHS [23,24] and JT gravity [25,26], and importantly we compute the free energy using the full GPI (computed for CGHS in [27] and JT gravity in [19]), rather than a saddle-point approximation. In both models, we find that replica wormholes substantially change the behavior of the free energy at sufficiently low temperature. Interestingly, in JT gravity, we find that the temperature at which the pathological behavior of the disconnected free energy manifests, and the temperature at which contributions from replica wormholes dominate, both scale like e −2S 0 /3 (where e −S 0 controls the JT gravity genus expansion). Since the gravitational theory is only under control for large S 0 , one might be concerned that the contribution of the replica wormholes happens in a regime of the theory in which we have no perturbative control. In fact, working at large S 0 but with T e 2S 0 /3 of order unity puts us in the so-called Airy limit, where the system is controlled by the universal behavior of the edge of the classical density of eigenvalues ρ 0 (E) 2 . In this limit, the genus expansion can be summed, providing a handle on doubly-nonperturbative corrections (in S 0 ). We find that these corrections are unimportant in part of the regime where replica wormholes dominate, so we can conclude that they genuinely do contribute even when doubly-nonperturbative corrections do not. This story is entirely analogous to the replica wormholes narrative in the context of black hole evaporation: some parameter k parametrizing the entropy of matter must become nonperturbatively large in S 0 in order for replica wormholes to dominate, and this transition happens right at the edge of validity of the semiclassical approximation. In our context, the parameter that must become large for wormholes to dominate is instead the inverse temperature β.
In an intriguing turn of events, while the replica wormholes do mitigate the pathologies in the free energy, we cannot show that they remove them entirely. We argue that this is due to the inherent ambiguity in the analytic continuation that defines F . To gain more insight into this ambiguity, in Section 5 we point out that an extremely similar phenomenon happens in spin glass systems, where a quenched disorder can allow for the spontaneous coupling of replicas used to calculate F . In that context, we review the Sherrington-Kirkpatrick (SK) model of spin glasses, and note that similar to our gravity calculations, at high temperature the free energy is dominated by a paramagnetic phase in which the replicas are uncorrelated, while at sufficiently low temperatures the system enters a spin glass phase in which the replicas correlate 3 . A replica-symmetric analysis of the spin glass phase exhibits the same sorts of pathologies that we see in the quenched free energy of CGHS and JT gravity; it turns out that in the SK model, replica symmetry breaking (RSB) is the key structure that "fixes" the analytic continuation in the replica trick and gives the correct free energy down to zero temperature. Motivated by the parallels between spin glasses and our gravitational results, we conjecture that the same sort of RSB is needed in the gravitational case to fully capture the correct behavior of F at low temperature. Importantly, the RSB that we discuss is notably different from the sort of RSB ordinarily discussed in the context of gravitational calculations of Renyi entropies. We make more exploratory comments about possible parallels between gravity and spin glasses in Section 6, but also note that our results should not necessarily be interpreted as indicative of a literal gravitational spin glass phase.
Relation to prior work: In the context of JT gravity, preludes of the transition in which we are interested can be found in analyses of the two-point correlator Z(β 1 )Z(β 2 ), which is relevant for studies of the spectral form factor. For instance, [28,29] find that at temperatures lower than O(e −2S 0 /3 ), the contribution of the cylinder topology to this correlator can become larger than that of the disk; see also [30,31] for the same behavior in nonperturbative completions of JT gravity, without needing to work at large S 0 . See also [32] for an analogous transition in a Gaussian matrix model. Our purpose here is specifically to investigate the contributions of connected topologies to the quenched free energy via the replica trick for ln Z.
While we emphasize that we do not claim a bona fide spin glass phase in JT gravity, the behavior is sufficiently similar that further comment is warranted given recent studies on SYK. These investigations show that SYK does not exhibit a spin glass phase; that is, a saddle-point analysis of the replica trick in the large-N limit (see e.g. [33,34]) indicates that no saddles correlating different replicas dominate the correlators Z m at any temperature [33,[35][36][37][38][39][40][41][42][43]. Here we point out that (i) we do not work in a saddle-point approximation, and in fact we expect that the behavior we study would be invisible in such a limit; and (ii) JT gravity is only dual to a low-energy regime of SYK, and as shown in [39] an appropriate IR limit of SYK can exhibit a different phase structure than the full SYK system. Hence there is no tension with our results.
More generally, attempts to model spin glasses holographically, such as e.g. [44,45], typically manually turn on a correlation between the different replica boundaries in order to induce a spin glass phase transition; this is analogous to the correlation between replica boundaries that occurs in computations of the entropy of Hawking radiation (due to tracing out a subsystem), or to the coupling of two boundaries in the traversable wormhole setup of [46,47]. Here we are specifically interested in the contribution of replica wormholes to the GPI P(B m ) defined by m completely uncoupled boundaries: the coupling happens entirely spontaneously and is an inevitable consequence of replica wormholes.
On a more tangential note, let us finally point out that there has been an ongoing discussion of the relevance of spin glasses to the physics of eternal inflation as well as to the landscape of string vacua. See e.g. [48] as well as [49] for an excellent review, and also [50] for more recent work. In a similar vein, [51] discussed these topics in the context of AdS 2 , and [52] and [53] studied a spin glass phase of black hole microstates (without external coupling). It would be interesting to explore connections to our present work.

The Replica Trick for ln Z
The purpose of this section is to discuss in more detail the replica trick necessary for the computation of the free energy F , and more generally the ensemble average of the generating functional ln Z considered as an arbitrary function of sources. Since such an average appears in the computation of Renyi entropies S n , and hence also of the von Neumann entropy, we will also discuss the relation to the replica trick used in the computation of von Neumann entropy.
The key point is that if the GPI is interpreted as the ensemble average of a partition function as per (1.3), then it cannot directly compute the ensemble average of any extensive quantity, such as ln Z. The replica trick relates such extensive observables to non-extensive objects via where the first term comes from summing over geometries that leave all the replica copies of B disconnected from one another, while the sum represents integrals over geometries that connect two or more copies of the boundary (i.e. replica wormholes) 4 . We therefore generically have P(B m ) = P(B) m . However, in certain cases one topological sector may dominate over others. If the dominant contribution is disconnected, then we have In this case, using (2.1) we see that ln Z ≈ ln Z, so the replica trick has no appreciable effect; in the condensed matter language used in Section 1, the quenched free energy and the annealed free energy approximately coincide. In particular, we may compute the gravitational free energy by just just taking F ≈ −T ln P(B), as usual. On the other hand, if a topology connecting multiple copies of B dominates, then we should expect that so the quenched and annealed free energies should not even approximately coincide, and a proper computation of the gravitational free energy will not coincide with the annealed free energy: F ≈ −T ln P(B).
Let us now exhibit how the replica trick (2.1) relates to the one used to compute the von Neumann entropy. This latter replica trick defines the von Neumann entropy of a subsystem (say a region R ⊂ B) as a limit of Renyi entropies: where the Renyi entropies S n are given by with B n an n-sheeted geometry consisting of n copies of B cut along the region R and then cyclically identified along this cut; see Figure 1. If R is empty, B n is just B n , consisting of n copies of B.
Suppose we now wish to evaluate the Renyi entropies via a gravitational path integral, under the interpretation that it computes an ensemble average of S n (and hence also of the von Neumann entropy). Such a computation requires the ensemble averages ln Z(B n ) and ln Z(B), which in turn requires use of the "extra" replica trick (2.1): where B m n consists of m separate copies of the n-sheeted geometry B n , as shown in Figure 1. A correct calculation of the von Neumann entropy therefore requires taking the double limit m → 0, n → 1.
A key distinction to note here is that the m replicated boundaries B m n are completely disconnected ; any geometric connection between them must come spontaneously from the GPI. On the other hand, when R is non-empty, the geometry B n is a single connected geometry, due to the identification of the n sheets along the cut R. In fact, m n B m n = Figure 1. A computation of the Renyi entropy (2.7) from the GPI requires an additional replica trick, involving computing the GPI with the boundary B m n shown here. Each of the columns is an n-sheeted geometry B n constructed by slicing n copies of B along the region R and then identifying these copies cyclically along the cut. B m n consists of m copies of this multi-sheeted geometry. The disorder-averaged von Neumann entropy is computed in the double limit m → 0, n → 1.
when R is the empty set, the Renyi entropies must vanish exactly (since we are computing the entropy of a pure state); it is precisely the auxiliary replica trick over m that guarantees this. To see this, note that if R is empty, B m n = B nm , and hence so the two terms in (2.7) cancel identically, giving S n = 0 for all n. Importantly, the vanishing of the Renyi entropy is independent of the dominant topology contributing to the path integral. This should be contrasted with, for example, the computation of Renyi entropy performed in [5], which (working in a semiclassical regime) claimed that the entropy a pure state vanishes because in that case the GPI is dominated only by disconnected topologies. The trouble with that interpretation is that even when the disconnected topology dominates, the path integral will still receive subdominant corrections from connected topologies which would lead to a nonvanishing (but small) Renyi entropy. The double replica trick makes clear that the Renyi entropy of a pure state vanishes exactly, and even when the dominant geometry is a replica wormhole. Of course, the claim of [5] that (at least in their JT gravity model) the disconnected topology dominates the gravitational path integral in a semiclassical limit when R is the empty set might lead to a concern: even if replica wormholes make subdominant contributions to the free energy, they might never be dominant in a regime in which the gravitational theory is under control. If so, the extra replica trick (2.1) will in practice never be necessary for computing leading-order effects. To address this concern, we will now explore explicit examples of gravitational models in which connected saddles Figure 2. The only topologies that can appear in the CGHS path integral are the disk and the cylinder.
do make dominant contributions when the theory is at least somewhat under control, focusing specifically on computations of the free energy F .

Free Energy in CGHS
We begin the investigation in gravity with a variant of standard CGHS dilaton gravity [23], introduced as the CGHS model in [54] (following [55]). This model is given by the Euclidean action where S ∂ is a boundary term. The equation of motion for A µ fixes Ψ to be constant, and it is this constant value that sets the temperature of black hole solutions, while the equation of motion for Φ sets R = 0. In fact, even in the path integral the integration over Φ means that only strictly flat geometries contribute. Hence the only contributions can come from the disk or the cylindrical topology, corresponding to one and two boundaries, respectively; see Figure 2. It is this simplification that will allow us to make definitive statements about the structure of the replicas and free energy in this model, without needing to worry about nonperturbative effects arising from highergenus contributions. This section is therefore a warmup for the JT gravity calculation in Section 4, which is complicated by contributions from all topologies.

Path integrals in CGHS
The path integrals of the disk and cylinder in CGHS were computed in [27]. For the disk with boundary length β, the result is 5 Already we can deduce the need for a phase transition. If the disk were to dominate the free energy, we would have F = −T ln P disk , which is clearly a non-monotonic function of temperature: it has a local maximum at T max = 1/ √ 2πe, corresponding to a negative thermodynamic entropy −∂F/∂T when T < T max (in fact, the entropy is logarithmically divergent at T = 0). We might hope that the contribution of the cylinder will rectify this low-temperature behavior.
To that end, the path integral on the cylinder (each of whose boundaries has length β) is Let us use P m (β) to denote the GPI defined by m boundaries of length β. This path integral receives competing contributions from the disk and the cylinder; the completely disconnected topology gives a contribution of while the topology that connects m/2 pairs of boundaries with cylinders (temporarily taking m to be even) gives a contribution At temperatures larger than T c ≡ 2 −1/3 , the contribution from the disk topology is larger, while for temperatures smaller than T c , the contributions from the cylinder topology is larger. So already at the level of this rough analysis we see a transition: the high-temperature behavior is controlled by the disconnected topology, while the lowtemperature behavior is controlled by a connected one 6 . Importantly, T c > T max , so the contribution from the cylinder modifies the free energy in the temperature regime in which the annealed free energy F ann ≡ −T ln P disk was pathological. Now let us be more thorough and compute P m (β) exactly, therefore attempting to obtain the free energy via the m → 0 limit (2.1). Defining r ≡ P cyl /P 2 disk , we have where the sum counts contributions from all aways of connecting an even number 2m of boundaries together via cylinders, the binomial coefficient counts the ways of choosing 2m boundaries from the full set of m, and the double factorial counts how many distinct ways there are of connecting those 2m boundaries pairwise with cylinder topologies. Expressing the double factorial as we find The sum can be evaluated using the identity 7 for any positive integers m and M , resulting in (3.10) To compute ln Z, we want to now continue to m → 0. 7 (3.9) can be shown by expanding the binomials on the right-hand side and then using the identity for sums of roots of unity: M, k ∈ Z and k = 0 (mod M ) .

Continuing to non-integer m
The result (3.10) can be naturally continued to non-integer m, but it exhibits a curious feature: because the second term 1 − √ 2tr will always become negative somewhere in the region of integration, for non-integer m this term need not be (and is not) real. Invoking the replica trick (2.1) at this stage would then yield a complex free energy, which is manifestly unphysical. Evidently, the obvious analytic continuation of (3.10) to non-integer m cannot be the correct one for the replica trick. A more well-behaved alternative can be obtained by noting the following. For any analytic function f (z) of a complex variable z, let f * (z) be the function obtained by complex-conjugating the Taylor series coefficients of f (z); then by construction the function f r (z) ≡ (f (z) + f * (z))/2 is also analytic, and is real whenever z is. If f (z) is real when z is a positive integer, then f r (z) = f (z) when z is a positive integer, and both f (z) and f r (z) therefore give admissible analytic continuations from the positive integers to general complex z. For this reason, for the purposes of computing F via the replica trick we are free to simply use the real part of (3.10) when m is real, which gives It may seem that we have pushed the replica trick to a breaking point. Of course there was always an infinite amount of freedom in how to continue the path integral P m (β) from positive integer m to non-integer m near zero, but the implied hope was that a "natural" analytic continuation should present itself, and that this continuation should be the correct one for getting the physically correct free energy. But the natural continuation of (3.10) gives a complex free energy, and we had to introduce a rather ad hoc procedure for modifying the continuation to obtain (3.11). What prevents us from, say, adding g(T ) sin(πm) to P m (β) with g(T ) an arbitrary function of temperature, and therefore getting whatever free energy we want?
This discomfort is well-justified, for there is an even more serious problem with the continuation of either (3.10) or (3.11) to general complex m. In order to consistently interpret P m (β) as giving the disorder average Z m of some power of the partition function, its behavior for purely imaginary m = iα must be bounded since where we have assumed that the disorder average is defined by a proper probability distribution (i.e. one that is positive and normalized). But while the terms on the first line of (3.11) are bounded when m is imaginary, the term on the second line is not, and indeed it grows arbitrarily large for large imaginary m. So (3.11) cannot be interpreted as the analytic continuation to complex m of an ensemble average Z m with respect to a positive and normalized probability distribution. In principle we should therefore look for a different analytic continuation that is well-behaved for imaginary m and hope that, say, Carlson's theorem is sufficient to ensure uniqueness of this continuation 8 . However, the growth of (3.11) at large real m excludes this possibility. To see why, note that (3.11) grows faster than exponentially in m at large real integer m, which can be seen easily by, say, keeping only the m = m/2 term in the sum (3.6). To try to prove that the analytic continuation to noninteger m must be unique (once we impose boundedness for imaginary m), suppose we had two different analytic continuations P (1) m and P (2) m , and let us try to show that their difference ∆P m must vanish. This difference of course vanishes on the positive integers, and must also be bounded on the imaginary axis if both P (1) m and P (2) m are. To invoke Carlson's theorem to conclude that ∆P m must vanish identically, we therefore only need to guarantee that ∆P m grows no faster than exponentially in the right half-plane; but this is not a condition we can enforce via any constraint on P (1) m and P (2) m due to their superexponential growth for integer m, and hence Carlson's theorem cannot be invoked.
The ambiguity in finding the "correct" analytic continuation is a substantial obstacle that we will address in much more detail in Section 5; it will be interpreted as a signature of replica symmetry breaking. For the time being, we will forge ahead by just using (3.11), assuming that the temperatures at which the quenched free energy is sensitive to contributions from the cylinder coincide with the temperatures at which the P m (β), and therefore the free energy obtained from (3.11), are. In proceeding in this way, we will be unable to determine what the correct form of the quenched free energy F actually should be, but we can still investigate when contributions from the cylinder cause the quenched and annealed free energies to differ.
With this important caveat in mind, the free energy obtained from (3.11) is At hight temperature T T c , r is small, so the second term is suppressed like O(r) and the free energy is controlled by the disconnected topology. On the other hand, at low temperature T T c , r is large and the integral can formally be expanded in powers of 1/r, with the leading contribution given by (1/2) ln r. Hence the behavior of the quenched free energy is where the ellipses denote subleading terms of order unity. At high temperatures, the free energy is the annealed free energy −T ln P disk sensitive only to the the disk topology, while at low temperature the leading-order behavior is modified thanks to the cylinders. Note that F is still not monotonic in temperature, even with the cylinder contribution. In particular, while the cylinder contribution decreases the severity of the logarithmic divergence (in reducing the prefactor of 2 to a 1/2), it does not eliminate it entirely. As discussed above, since the calculation of P m (β) for integer m was exact and involved no approximation, the culprit for this unphysical behavior is the analytic continuation away from integer m 9 . This should come as no surprise, as we have already established that the analytic continuation given by (3.11) does not behave correctly for imaginary m; clearly it needs to be modified to remove the pathological behavior entirely.
Nevertheless, the key point is that the replica trick is required to see that F receives large corrections from the cylinder topology right around the temperature where the annealed free energy is badly-behaved. Without properly understanding how the analytic continuation to non-integer m is to be perfored, we cannot know in precisely what way these additional corrections modify the free energy; the analytic continuation given in (3.11) is insufficient to remove the low-temperature pathology entirely, but we expect that the correct continuation should give a monotonic free energy that yields a vanishing entropy −∂F /∂T at zero temperature. We will revisit this issue in Section 5.
We will first do a preliminary analysis of the role of Euclidean wormholes in the replica computation of the free energy, beginning with a brief review of the salient features of the JT gravity path integral (using specifically the results of Saad, Shenker, and Stanford [19]). The (Euclidean) JT gravity action is where volume elements are left implied and K is the extrinsic curvature of ∂M . When ∂M consists of a single circle, the boundary conditions take the length of ∂M to be β/ and set the dilaton φ| ∂M = γ/ there; after the introduction of an appropriate counterterm, the limit → 0 is understood. For simplicity, we will work in units where γ = 1; this amounts to working with the dimensionless rescaled inverse temperature and free energy β/γ, γF respectively. When ∂M consists of several circles we may specify boundary conditions separately on each, but for our purposes it will suffice to take all boundary components to have the same length β/ . The path integral over the dilaton fixes the path integral over geometries to only include those with constant negative curvature; this space of topologies is of significantly richer structure than its flat counterpart and leads to the organization of the path integral in a genus expansion. For example, if P conn,2 (β) is the path integral over geometries that connect two boundary components (both of which have length β/ ), pictorially we have P conn,2 (β) = Explicitly, the path integral P conn,m (β) over geometries that connect m boundary components is given by where the objects Z g,m (β) are and V g,m (b 1 , . . . , b m ) are the volumes of the moduli spaces of Riemann surfaces with m geodesic boundaries of lengths b 1 , . . . , b m (we work in the convention where the normalization α of these volume forms is one, corresponding to V 0,3 = 1). The V g,m can be computed algorithmically using, for example, Mirzakhani's recursion relation [56]; a table summarizing the data for small g and m can be found in [57]. The genus expansion, as well as the contribution of topologies that connect arbitrarily many boundary components, makes the story for JT gravity substantially more involved than for CGHS. Nevertheless, even at this heuristic level we can now see that connected topologies must be included in, and will upon inclusion significantly affect, the low-temperature behavior of the free energy: for example, if we were to only consider the contributions from the disk topology Z 0,1 and the "double trumpet" Z 0,2 , the analysis would proceed just as in the CGHS case, and we would expect the double trumpet contribution to the free energy to compete with that of the disk whenever Z 0,2 /(e S 0 Z 0,1 ) 2 is order unity or larger. For large S 0 , this will occur at temperature T e −2S 0 /3 , so that at sufficiently small temperatures failure to include the connected topologies yields a result that is manifestly wrong, as those topologies contribute at least as much as the disconnected ones.
This observation raises a potential concern. The parameter e −S 0 is supposed to suppress the contributions from higher genus, as well as from topologies that connect more boundary components. But at low temperature β 1, the leading-order behavior of the Z g,m scales like β (3/2)(2g+m−2) , so contributions from higher genus and more-connected topologies are controlled by β 3/2 e −S 0 . The regime in which Euclidean wormholes contribute to the free energy therefore corresponds to the parametric regime in which we lose perturbative control of the genus expansion. What do we make of this?
From the perspective of the Euclidean wormholes, the story is completely analogous to that of quantum extremal islands in the computation of the entropy of Hawking radiation [4,5]. In that case, there is an auxiliary parameter k parametrizing the entropy of matter fields 10 , and replica wormholes lead to the presence of a quantum extremal island when k is nonperturbatively large: the Page transition happens at k ∼ e S 0 . In the present context, the inverse temperature β plays the role of k. On the other hand, from the perspective of the genus expansion we are justified in being concerned, because without control of the connected path integral P conn,m we cannot expect to make any substantive claim regarding the contribution of Euclidean wormholes. Fortunately, the regime we are discussing -that is, taking S 0 large but keeping β 3/2 e −S 0 of order unity -recovers the so-called Airy case of random matrix integrals, in which the partition function Z(β) is governed by the behavior at the edge of the spectral density ρ(E). This simplification makes it possible to resum the genus expansion to include doublynonperturbative (in S 0 ) effects, which we can use to assess how well-behaved the genus expansion is. Before proceeding, it will therefore be useful to discuss this regime in more detail.

The Airy limit
Before diving into the details of the Airy case 11 , let us first do a rough analysis of the behavior of the genus expansion in the regime β ∼ e 2S 0 /3 where we expect contributions from Euclidean wormholes to become important. Recall that the genus expansion (4.3) is asymptotic, meaning that it does not converge even when β 3/2 e −S 0 is small. Nevertheless, as with any asymptotic series, the partial sums in the genus expansion can be used to bound the free energy. When β 3/2 e −S 0 is not too small, the genus expansion can still be "under control" in the sense that the first few terms in the series (4.3) decrease, so that the partial sums provide a tight bound on the free energy. To that end, using (4.4) and the explicit forms of V g,m found in e.g. Appendix B of [57], in Figure 3a we plot the annealed free energy F ann ≡ −T ln P conn,1 (corresponding to the disconnected topology free energy −T ln Z) for S 0 = 7 where we include topologies only up to genus g = 5. The first few partial sums of the genus expansion do indeed provide accurate approximations to the free energy for T e 2S 0 /3 0.3, which crucially includes a local maximum. This is suggestive that this maximum should also be present in a full nonperturbative computation of F ann -but as discussed above, such a maximum is an unphysical feature of the free energy, which we expect to be resolved by the inclusion of connected topologies, indicating that inclusion of the latter is indeed necessary.
To proceed more carefully, we can in fact exchange the asymptotic genus expansion for an asymptotic low-temperature expansion with T e 2S 0 /3 fixed, verifying that it reproduces the behavior exhibited in Figure 3. To do so, note that the Weil-Petersson volume forms V g,m appearing in (4.4) are polynomials in the b i , and therefore the Z g,m are polynomials in β of order (3/2)(2g + m − 2), as mentioned above: where (up to various constants) the leading-order terms P 0,g,m are the intersection numbers of Chern classes (more generally, the P ,g,m are intersection numbers of the  first Miller-Morita-Mumford class with Chern classes [28,58]; more explicit expressions can be found in Appendix A). Inserting this expression into (4.3), for certain m the sum over genus can be performed as described in [28,29] to produce a low-temperature asymptotic expansion; for example, for m = 1 we have where the first fewz (h) are given explictly in [28]. For βe −2S 0 /3 of order unity, this asymptotic expansion is under control for large β. In Figure 3b we show the annealed free energy computed using (4.7) for S 0 = 7, and find that as expected, the lowtemperature expansion agrees with the first few partial sums of the genus expansion in the region T e 2S 0 /3 0.3. This allows us to conclude that the unphysical peak in the free energy at T e 2S 0 /3 ≈ 0.7 cannot be eliminated by either higher order terms in the genus expansion or by doubly non-perturbative effects.
In fact, there is more we can say in this low-temperature limit. Sincez 0 (h) = 1, the leading-order term in (4.7) is given by This is precisely the partition function in the Airy case of random matrix theory and topological gravity, where the Airy density of eigenvalues is given by [59,60] ρ (4.10) The leading-order behavior (in e −S 0 ) of ρ Airy (E) is just which is the universal behavior of the leading-order density of eigenvalues near the edge of the of the spectrum in the double-scaled matrix models of [19]. Hence the lowtemperature expansion (4.7) can be thought of as an expansion about the low-energy edge of the spectrum, with the subleading terms capturing deviations from the exact form (4.11). Concretely, it corresponds to taking S 0 → ∞ while keeping βe −2S 0 /3 fixed. The contribution to P conn,m from this leading-order low-temperature behavior can be summed over genus for any m using the results of [61]; we summarize the relevant results in Appendix A, and the relevant expression for P conn,m is given by (A.10).
The fact that the low-temperature limit in which we are interested is dominated by the universal behavior (4.11) means that we may gain some qualitative insights into the competition between connected and disconnected topologies by considering particularly simple matrix models. For example, the Gaussian matrix integral has a leading-order density of eigenvalues given by the Wigner semicircle which recovers (4.11) in the double-scaling limit E → E − a followed by a → ∞ [62]. The exchange of dominance between connected and disconnected topologies in the Gaussian matrix integral was studied in [32], where it was found that the connected correlator Z(β) 2 conn becomes larger than the disconnected correlator Z(β) 2 at temperatures lower than ∼ N −2/3 (or ∼ e −2S 0 /3 using JT terminology). So the behavior we are exploring is a general feature of random matrix models.
The upshot is that the low-temperature regime in which we are interested is quite well-understood; importantly, the contributions of higher genera (and their associated doubly-nonperturbative corrections) are insufficient to eliminate the pathological behavior of the annealed free energy. Therefore, we now turn to a computation of the quenched free energy via an analytic continuation to near m = 0.

The continuation in m
To compute the contribution of Euclidean wormholes to the quenched free energy via the replica trick, we need the JT gravitational path integral P m (β) defined by m disconnected boundary circles, each of length β/ . These are related to the connected path integrals (4.3) by the usual relation In order to continue to near m = 0, we need to express P m (β) in a form analytic in m; this is difficult because the Weil-Petersson volume forms V g,m , and consequently the coefficients Z g,m (β) in the genus expansion, are not known analytically in m. This is true also in the Airy limit discussed in Section 4.2 where although explicit formulas are known (see Appendix A for a review) they are not written as analytic functions of m. We will therefore proceed in an alternative fashion: we define a "truncated" path integral P m,M to be the JT gravity path integral including only topologies that connect up to M boundaries, with M some fixed integer (this amounts to truncating the sum on the right-hand side of (4.13) to m ≤ M ). We then analytically continue P m,M to non-integer m with M held fixed, defining a truncated free energy If as M → ∞ the analytic continuation of P m,M (β) to non-integer m converges to a function P m,∞ (β) which is also analytic in m, we may take P m,∞ (β) to define the analytic continuation of P m (β) to non-integer m. We can then express the free energy In practice, we will compute the truncated free energies F M for some relatively small values of M , which by the argument above we might expect to give us an approximation to the exact free energy F . In particular, F 1 is just the annealed free energy shown in Figure 3, so we are interested in modifications to the behavior of F M as M is increased, specifically in the regime T e 2S 0 /3 0.3.
To obtain the aforementioned continuation of P m,M (β) to non-integer m, we proceed inductively: noting that for M = 1 we have P m,1 (β) = P conn,1 (β) m , we will suppose that for arbitrary M we may write where C is any contour that encloses z = 0. Both of these equations are correct for integer m ; for m not an integer, the first expression is of course just the definition of the gamma function Γ(M m + 1) (for Re(M m ) > −1), but due to the branch cut of z −(m +1) along the negative real axis, the second only coincides with 1/Γ(m + 1) if C is chosen to be a Hankel countour 13 . But since (4.22) are only required to hold when m is a positive integer, there is no need to require C to be a Hankel contour, and in the freedom in choosing C we already see a foreshadowing of the freedom that will manifest in the analytic continuation to near m = 0. Using the identity (3.9), we may evaluate the sum over m to obtain Iterating these from the base case M = 1 (for which the sum (1) I is empty and A (1) = P conn,1 ), we therefore find The analytic continuation to near m = 0 is now straightforward; bearing in mind that as in the CGHS case we must take the real part, we find As already noted, this free energy depends on the choice of contours C k for the integrals over z k introduced in the analytic continuation (4.22). Specifically, the integrand of (4.27) exhibits branch cuts in the complex z k planes, and will therefore be sensitive to where the contour C intersects these cuts. This is not surprising: as discussed in Section 3.2, inferring the "correct" analytic continuation to near m = 0 is rather subtle.
We would now like to verify that the corrections from replica wormholes significantly alter and even dominate the behavior of the free energy in the regime e 2S 0 /3 T 0.3 with S 0 large in which we have shown we have perturbative control of the genus expansion. To that end, we again use (4.4) (along with the explicit forms of the V g,m ) to compute F M , incorporating contributions up to g = 2 and M = 5; the results are shown in Figure 4. Note that in Figure 4 we take the contour C in (4.22) to be the unit circle for simplicity. It is clear that in the regime e 2S 0 /3 T 0.3, the inclusion of replica wormholes can substantially modify the behavior of the free energy. The unphysical local maximum appears to be "softened" by the replica wormholes contribution, though  we should be careful not to draw any firm conclusions about the quantitative features of F M due to the ambiguity in the continuation to near m = 0 (including, for instance, whether the M → ∞ limit even exists). In short, we can ascribe meaning to the fact that the free energy changes when replica wormholes are included, but we cannot know its quantitative behavior until we know how to pick the "right" continuation. To highlight this point, in Figure 5 we compare the g = 0 free energies obtained from the analytic continuation (4.22) with C the unit circle to another analytic continuation in which we instead used the gamma function multiplication theorem to write (a) Using (4.22) with C the unit circle.  with the contour C taken to be the unit circle (this is the same as Figure 4a), while on the right we used (4.28). The qualitative features agree, but quantitative details do not. The blue, orange, green, red, and purple curves correspond to M = 1, 2, 3, 4, 5, respectively, and we take S 0 = 7.
and then expressed the gamma functions in the product in their integral form. The qualitative features of the free energy computed with these two different analytic continuations agree well, but of course they differ quantitatively. At this point we do not know how to specify the correct prescription, but for reasons that we will describe in the next section, we expect the answer will involve replica symmetry breaking in the m → 0 limit. As a final note, it is interesting to examine the behavior of F M using the the leadingorder low-temperature behavior of P conn,m discussed in Section 4.2 and Appendix A. Specifically, using equation (A.10) for the path integral in the Airy limit, we obtain the behavior of F M shown in Figure 6. While again we may not draw any definitive quantitative conclusions due to the ambiguity in the analytic continuation, we see that connected topologies affect the behavior of the free energy even when all all terms in the genus expansion are included.

Replica Symmetry Breaking and a Spin Glass Analogy
We have shown that in computing extensive quantities like the free energy in gravitational systems, the interpretation of the GPI as an ensemble average -requiring a Figure 6. The behavior of the quenched free energy F M for the Airy case, including contributions from all genera using the result (A.10). This amounts to taking S 0 → ∞ with T e 2S 0 /3 held fixed in the JT path integral. As in Figure 4, here we take the contour C in (4.22) to be the unit circle, and the blue, orange, green, red, and purple curves correspond to M = 1, 2, 3, 4, 5, respectively.
replica trick for the computation of the quenched free energy -can lead to a contribution from replica wormholes that exceeds that of disconnected topologies. The necessity of these corrections can already be inferred from the pathological properties of the lowtemperature behavior of the annealed free energy, computed just from disconnected topologies without resorting to a replica trick. We have also seen that the inclusion of replica wormholes remedies some of these pathologies but is not sufficient to remove them entirely; we interpret this as necessitating a clearer understanding of the correct analytic continuation to m = 0. Indeed, let us emphasize that in the simpler case of CGHS, the (nonperturbative) calculation that includes all of the allowed geometries still exhibits a pathological annealed free energy at low temperatures. This calculation had only one potential pitfall: the m → 0 analytic continuation. This immediately implies that it is the choice of the straightforward analytic continuation that is directly responsible for the incorrect result. All of these features -an annealed free energy with pathological low-temperature behavior, an improvement in this behavior under the inclusion of connected replicas in computing the quenched free energy, and the need for a careful analytic continuation to near m = 0 to eliminate the pathological behavior entirely -are exhibited in the well-studied context of spin glasses. In order to draw an analogy with these systems, we will now review one particularly well-known example: the Sherrington-Kirkpatrick (SK) model [63]. In this system, we will see that the non-uniqueness of the analytic continuation to m = 0 is due to a replica symmetry-breaking transition that occurs at m < 1, suggesting that a similar transition likely occurs in the gravitational sys-tems we have examined, and that it is unlike the usual Z n -replica symmetry breaking that is discussed in the context of gravitational calculations of the Renyi entropies. We will keep the review of the SK model limited to the bare essentials, but would recommend [64,65] and especially [66] for more comprehensive treatments.

Review of the SK model
The SK model is an infinite-ranged classical Ising model of N interacting spins σ i , with Hamiltonian where the sum runs over all distinct pairs of spins (ij). Each of the random couplings J ij is drawn from a Gaussian 14 distribution P (J ij ) with mean J 0 /N and and variance J 2 /N . As above, we will denote averages over the distribution P (J ij ) via an overline, so that, for instance, the ensemble average of the logarithm of the partition function is Note that ln Z is quite difficult to compute directly, but using the replica trick (2.1) requires us to simply compute the ensemble average of the m-replicated partition function where α is a replica index that labels m copies of the spins σ α , and the last trace is over all m replica systems. The last average is quite easy to express in terms of the moments J 0 and J of the distribution P (J ij ): The fact that the couplings J ij are correlated between the different replicas has led to the introduction of an effective coupling between replicas via the ensemble average. Moreover, by completing the squares in the sums over spin sites and introducing auxiliary variables s α , q (α,γ) with α = γ (sometimes called Hubbard-Stratonovich variables, 14 We could consider a more general distribution, but the important physics is captured by just the second moment of P (J ij ). collective fields, or mean fields), we may decouple the spin sites: where B is a prefactor that is sub-exponential in N (and therefore will be irrelevant in the thermodynamic limit N → ∞), the variables s α and q (α,γ) are all integrated over the real axis, and the notation (α, γ) denotes all distinct pairs of replicas. Here the effective Hamiltonian H eff is independent of N and given by where the trace is now over all m replicas of a single spin site and At this point (5.5) is still an exact equation, whose existence is made possible thanks to the all-to-all coupling of the SK model: the fact that the couplings between all pairs of sites are drawn from the same distribution allows for the factorization of different spin sites in (5.4) via the introduction of the variables s α and q (α,β) . We may now take the thermodynamic limit N → ∞, finding via a saddle point approximation that where now s α and q (α,γ) are solutions to the saddle point equations ∂H eff /∂s α = 0 = ∂H eff /∂q (α,γ) . It is easy to see that these conditions reduce to giving s α and q (α,γ) the interpretation of mean fields fixed by the self-consistency conditions (5.8). Importantly, the field q (α,γ) is interpreted as a coupling between replicas; a saddle with nonzero q (α,γ) indicates the spontaneous "turning on" of this coupling. Because this coupling is our main focus, from now on we will set J 0 = 0 so that Z m becomes independent of the mean field s α (this excludes the possibility of a ferromagnetic phase, in which we are not currently interested).
In order to now compute ln Z (and therefore F ) in the thermodynamic limit, we must analytically continue (5.7) to non-integer m near zero. Because the sums in H eff are only well-defined for integer m, this procedure requires positing some ansatz for the matrix q (α,γ) that is amenable to the analytic continuation to m = 0. Given the replica symmetry of the problem (corresponding to the permutation group S m ), it is natural to take the replica-symmetric ansatz q (α,γ) = q. (5.9) Indeed, for positive integer m, the dominant saddles do exhibit this symmetry [67].
The analytic continuation to near m = 0 is then straightforward, and the free energy becomes When βJ < 1 (i.e. at sufficiently high temperature), the only solution is q = 0, and hence the replicas are uncorrelated; this is the paramagnetic phase. The free energy obtained in this phase therefore satisfies ln Z = ln Z, i.e. we may average Z before taking the logarithm with no loss of information. Hence the replica trick does not introduce any novel behavior. As the temperature is lowered, however, a solution with nonzero q begins to exist once βJ > 1. This new solution dominates the free energy 15 , corresponding to the spin-glass phase in which the replicas spontaneously couple. While the field q was introduced in the context of the replica formalism, it has an interpretation in the m → 0 limit: it computes the so-called Edwards-Anderson order parameter q EA defined by the disorder-averaged square magnetization [68]: Here independence of the choice of lattice site i follows from translational invariance (after the disorder average), and the expectation value is a standard thermodynamic average taken with respect to a particular sampling of couplings: The non-vanishing of q in the spin glass phase therefore corresponds to magnetic order for any particular sampling of the couplings J ij . However, for J 0 = 0 the disorderaveraged magnetization vanishes: σ i = 0. Since this disorder-averaged magnetization measures the ferromagnetic order of the system, we see that the spin-glass phase corresponds to a cooperatively frozen magnetic state but with no ferromagnetic order.

Replica symmetry breaking in the SK model
As can be seen directly from (5.10a), the free energy of the paramagnetic phase q = 0 is pathological if we extend it to arbitrarily low temperature: at large temperatures it scales like −T , while at low temperatures is exhibits a −1/T divergence. These behaviors imply that it is non-monotonic, with the thermodynamic entropy becoming negative at sufficiently low temperatures (and in fact diverging at zero temperature). As shown in Figure 7, the turning on of the spin glass phase when T /J < 1 is necessary to alleviate these pathologies, rendering the free energy finite. However, it is still nonmonotonic: the zero-temperature entropy is S T =0 = −N/2π. Clearly the calculation remains incomplete; from our earlier discussion, we expect that this missing ingredient involves some nontrivial behavior of the analytic continuation from Z m at positive integer m to m = 0 16 . How do we understand what the correct analytic continuation is?
The answer can be gleaned by performing a stability analysis of the replica-symmetric ansatz (5.9). Indeed, though (5.9) does give the correct form of the saddles for computing Z m when m is a positive integer, it becomes unstable for sufficiently small m < 1: an eigenvalue of the Hessian ∂ 2 H eff /∂q (α,γ) ∂q (β,δ) evaluated on the ansatz q (α,γ) = q becomes positive in the limit m → 0 [69]. We must therefore invoke an alternative ansatz for q (α,γ) that avoids this instability as m → 0. The correct analytic continuation to m = 0 will then be determined by the behavior of the ansatz for q (α,γ) which remains stable down to m = 0; this behavior will undergo a phase transition at some critical m c (T ) < 1 [70] that was missed by just considering the replica-symmetric ansatz (5.9). The presence of this phase transition means that it is crucial to analytically continue the saddle-point equations ∂H eff /∂q (α,β) = 0 themselves down to m = 0, rather than first evaluating their on-shell value at integer m and then analytically continuing the results. Because the number of components of q (α,γ) is m(m − 1)/2 < 0 when m < 1, it is far from obvious how to construct a replica symmetry-breaking (RSB) ansatz that is amenable to analytic continuation. The answer is the well-established Parisi ansatz [71][72][73][74]. To get an idea of how this procedure works, consider splitting up the m replicas that define Z m into groups of m 1 , with m 1 an integer that divides m. We then write q (α,γ) in a block-diagonal form according to this grouping: where Q 1 and Q 2 are m 1 × m 1 matrices all of whose entries are q 1 and q 2 , respectively (in this example, we have m/m 1 = 4). This ansatz for q (α,β) can be analytically continued to m = 0 while leaving m 1 , q 1 , and q 2 free as variational parameters to be fixed by extremizing the free energy with respect to them (since 1 ≤ m 1 ≤ m, the analytic continuation of m also continues m 1 to be between zero and one). This procedure, called one-step RSB (or 1RSB), substantially improves the pathologies in the free energy shown in Figure 7, but the zero-temperature entropy is still negative (though substantially closer to zero) 17 .
To proceed further, we iterate this procedure: we introduce a new integer m 2 that divides m 1 and partition Q 2 into the same block-diagonal structure as (5.13), where Q 2 and Q 3 are m 2 × m 2 matrices all of whose entries are q 2 and q 3 , respectively (in this example m 1 /m 2 = 3). Repeating this process r times, we may then continue to m = 0, obtaining an expression for the free energy that depends on 2p+1 variational parameters, m i for i = 1, . . . , p and q i for i = 1, . . . , p + 1. After the continuation to m = 0 has been made, we may in fact take the limit p → ∞ which turns the (q i , m i ) into a continuous function q(x). The free energy is then a functional of q(x), and is obtained by a functional extremizaton with respect to q(x).
One way of understanding what the p → ∞ limit means is as follows. For positive integer m, the the ansatz (5.13) breaks the full replica symmetry group S m into the subgroup with the first factor corresponding to the permutation symmetry of each of the groups of m 1 rows and columns, and the second corresponding to the permutation symmetry of the m/m 1 groups amongst themselves. The iterative procedure outlined above amounts to breaking the subgroup further, into (with m p+1 ≡ 1), but of course we cannot take p arbitrarily large if the m i must all be divisors of m. However, if we analytically continue this group structure to m = 0, we obtain So we find that S 0 contains itself as a subgroup, which means we may continue to break the symmetry as much as desired by breaking the S 0 factor on the right-hand side. This is the feature that allows us to take p → ∞ in the Parisi ansatz after the continuation to m = 0 has been performed. The point is that RSB is contained in the structure of the Parisi function q(x): in the replica-symmetric ansatz (5.9) q(x) is just a constant q, so nontrivial structure e.g. the p-spin spherical model [75][76][77].
in q(x) is indicative of RSB. Because the Parisi ansatz changes the naïve analytic continuation to m = 0, we see that RSB is the mechanism reponsible for the phase transition at m < m c (T ), and it answers the question posed above: how do we correctly continue to m = 0?

RSB in Gravityà la Spin Glass
In Sections 3 and 4 we saw that in simple gravitational models, the introduction of replica wormholes alleviated some of the low-temperature pathologies of the disconnected free energy, but it did not remove them entirely; we interpreted this result as the statement that our anaytic continuation to m = 0 (which in the JT gravity case exhibited considerable freedom) was not correct. Having now reviewed spin glasses, there is quite an obvious analogy: since the paramagnetic and spin glass phases are characterized by correlated and uncorrelated replicas, respectively, we would like to interpret the "turning on" of replica wormholes in the gravitational free energy as the onset of spin glass-like behavior. It is important to note that the analogy will not be literal: perhaps the most important distinction is that a spin glass is a bona fide sharp phase transition that can be seen in the thermodynamic N → ∞, whereas we did not work in any saddle point approximation in our gravitational models (and in fact, the fact that the temperature at which connected topologies contributed was nonperturbatively small in S 0 suggests that the transition should be invisible to a semiclassical S 0 → ∞ analysis). The most relevant paralle we would like to highlight has to do with the allimportant analytic continuation: in the spin glass model, a replica symmetric ansatz remedies some low-temperature pathologies of the free energy, but it gives the incorrect analytic continuation, and RSB must be invoked due to a phase transition at small m. What does this analogy suggest for how to obtain the correct analytic continuation to m = 0 in the gravitational case?
One of the key lessons to draw from the spin glass example is that a naïve analytic continuation from the values of Z m for positive integer m to near m = 0 gives a wrong answer: we must first analytically continue the saddle point equations to near m = 0 with an appropriate ansatz, and only then do we solve them for the small-m behavior of Z m . In Sections 3 and 4, this is not what we did: we instead expressed the gravitational path integrals P m (β) for integer m, and then looked for an analytic continuation to m = 0. For the same reason as the spin glass, we might expect that in a gravitational theory we must look for RSB saddle points in order to perform the analytic continuation correctly.
Let us first be clear on what we mean by "replica symmetry breaking". There is a sense in which we could say that any replica wormhole breaks replica symmetry, since the symmetry group of m disconnected boundaries is S m , which is broken by any gravitational saddle that connects two or more of these boundaries. But the sort of RSB that appears in the spin glass example, and which we expect to determine the correct analytic continuation to near m = 0, is something more subtle: it is the breaking at m < 1 of a symmetry that is exhibited by the dominant saddles when m is a positive integer. For example, if the m-boundary gravitational path integral is dominated by disconnected saddles whenever m is a positive integer, the symmetry group is indeed S m , and we would say that RSB occurs if this group is broken for m < 1. But if the path integral for positive integer m is dominated by, say, a connected wormhole with Z m symmetry, we would not say that RSB occurs as m → 0 unless the Z m is broken for some m < 1. Now, since in Sections 3 and 4 we did not work in a saddle point approximation, no equations of motion were involved in our calculation. Hence it is not immediately clear what the analogue of the Parisi procedure might be in this models. It may instead be easier to consider working in the semiclassical limit of some more general gravitational theory, in which case probing the role of RSB, and computing the correct analytic continuation to near m = 0, requires us to look for a RSB ansatz for a gravitational solution that allows for the continuation of the gravitational equations of motion to near m = 0. This is still a difficult task, which is a natural starting point for future work. Instead, let us compare the approach we have in mind in this context with that of the Lewkowykz-Maldacena replica trick used to compute holographic von Neumann entropies [78]. In the latter case, we are required to compute the gravitational path integeral defined by an n-sheeted connected boundary manifold B n with Z n symmetry. Assuming the dominant bulk saddle also exhibits this symmetry, we may quotient the bulk geometry by Z n , after which the analytically-continued bulk equations of motion are just those on a manifold with boundary B 1 consisting of a single sheet, except with a conical defect proportional to (n−1) at the fixed point of the Z n isometry. For n near one, the bulk equations of motion can be expanded perturbatively around the smooth geometry with boundary B 1 and no conical defect, and the condition that the equations of motion hold near the (perturbative) conical defect reproduces the Ryu-Takayagani formula for holographic entanglement entropy [79]. In this context, the "usual" notion of RSB is the breaking of the Z n for n = 1 -but of course there is no breaking of replica symmetry for n = 1, since Z 1 is trivial. According to the alternative definition of RSB that occurs in spin glasses, RSB would require that the dominant saddles at positive integer n to exhibit Z n symmetry, but for the dominant saddles at small n, including n = 0, to break it.
Clearly the LM approach is along the lines we have in mind, as it continues the gravitational equations of motion to non-integer n. However, this continuation relies crucially on two properties. The first is the assumption of Z n symmetry, without which it would be unclear how to express the equations of motion on a manifold with a single boundary (just as in the SK model it was unclear how to generalize the replicasymmetric ansatz (5.9) until Parisi's breakthrough). The second is that there is a known n = 1 saddle around which the equations of motion can be perturbed to study the behavior near n = 1; there is no such saddle with n = 0. These are the two primary challenges that need to be overcome in order to properly understand the role of RSB in computing gravitational free energies, and more generally any extensive quantity.

Discussion
We have argued that the computation of extensive quantities via a gravitational path integral should be done using a replica trick which includes contributions from connected geometries. The inclusion of these connected saddle points dramatically changes the behavior of the theory at very low temperatures, and naturally accommodates the interpretation of semiclassical gravity as dual to an ensemble average rather than to a particular quantum theory. Let us now discuss open questions and natural directions for future work.
Ensemble Averaging in Higher Dimensions As alluded to in Section 1, UV corrections to the GPI may remedy the apparent lack of factorization that motivated the ensemble averaging interpretation in the first place, as discussed in the context of random matrix models and JT gravity in [19]. Such a picture becomes especially crisp in higher dimensions: for example, N = 4 SYM is a single theory, and AdS/CFT provides numerous other examples of unitary quantum theories of gravity without the need to ensemble average. If, however, one would like to apply the techniques of [4,5] to higher dimensions then we must include replica wormholes, whose most obvious interpretation is of an ensemble average. One possibility is that averaging is only genuinely necessary in certain low-dimensional theories (as was argued in e.g. [80]). For example, the low-temperature spectrum of higher dimensional gravity (and CFTs) is perfectly well-behaved, has a unique ground state, and does not resemble a spin glass. We do not expect to see replica wormholes or RSB dominating the free energy calculation at low temperature. Nevertheless, it is natural to speculate that replica wormholes will contribute to ln Z whenever we are in a regime where non-perturbative quantum gravitational corrections are important: for example, after the Page time [81] or at the Hawking-Page phase transition [82].
Another interesting possibility arises from the phenomenon of self-averaging: in a chaotic theory, the average over an ensemble of theories is often essentially harmless, as each individual instance of the ensemble is representative of the ensemble as a whole, at least for relatively coarse-grained observables. The ensemble average in this case is interpreted as a useful calculational trick to construct a universal effective theory which governs the dynamics at low energy, but of course the UV dynamics of each individual instance of the ensemble is that of a unitary quantum theory. Perhaps any gravitational theory which includes Euclidean wormholes should be understood as a low-energy effective theory in this sense; in this interpretation, the GPI plays the role of a convenient calculational trick for computing observables in a semiclassical limit. Such a possibility was discussed in various forms in [4,83,84].
Nonperturbative Completions At the end of Section 1, we briefly mentioned that although a large-N analysis of SYK does not exhibit a spin glass phase, [39] showed that in a large-coupling (or low-temperature) limit that reduces to an EFT of the low-energy dynamics of SYK, saddles that correlate replicas in the computation of Z m become dominant at both positive integer m as well as in the m → 0 limit, and therefore lead to a spin-glass like phase transition in this low-energy EFT. This observation may raise a concern: if a spin glass phase can only be obtained from the SYK model by excluding the UV, is the phase transition that we have found in JT gravity eliminated by a good UV completion? Our study of the Airy limit in Section 4.2 shows that a nonperturbative completion of JT gravity cannot eliminate the effect we have studied, since it is dominated by the universal behavior of the edge of the spectral density ρ (E). Indeed, the recent discussion of such completions in [31] explicitly finds that the twopoint correlator Z(β) 2 is controlled by the contribution of connected topologies at sufficiently low temperatures, even in a nonperturbative completion.
More generally, the results of [19] suggest that a good nonperturbative description of JT gravity should be available in the form of a matrix model (though this completion is not unique). Because the behavior we have studied in this paper is due to universal behavior at the spectral edge (at least at sufficiently low temperatures), we might investigate it more thoroughly by working in a toy matrix model like the Gaussian matrix integral investigated in [32]. To this end, it would be interesting to compute ln Z in such a model by expressing ln Z = ln Tr e −βH and then explicitly computing an average over the random matrix H, without resorting to a replica trick. We should expect to find a monotonic free energy all the way to zero temperature, with a free energy that agrees with the annealed free energy of the Airy case once the temperature becomes sufficiently (but not too) large.
The Emergence of Semiclassical Gravity A longstanding question in quantum gravity is how the (semi)classical metric g ab emerges from an underlying quantum theory. In the SK model, the partition functions Z m can be expressed exactly via the introduction of the mean fields s α and q (α,γ) in (5.5). In a large-N limit, the phase structure of the system is determined by the saddle point equations for these fields. Importantly, they appear purely as a consequence of the disorder average; they are not fundamental in the pre-disorder theory. (In the SYK case, the analogous fields are the auxiliary fields G αβ (τ 1 , τ 2 ) and Σ αβ (τ 1 , τ 2 ).) If we are to interpret the GPI as computing a disorder average (either genuinely or in an effective description for the purpose of probing appropriately coarse-grained observables), is there a sense in which the metric should then be thought of as a mean field, with the GPI analogous to the right-hand side of (5.5)? That is, rather than being a fundamental field of the underlying theory, is the metric a field whose existence relies fundamentally on the ensemble average? In such a case we would interpret the "turning on" of connected geometries between disconnected boundaries as analogous to the "turning on" of the matrix q (α,γ) in the SK model. This would give a clear meaning to the sum over topologies in the path integral, but even in the 2D models we have studied here it is unclear how this interpretation would incorporate a UV completion.

RSB and the Parisi Ansatz in Gravity
In the 2D models studied in this paper, the need for replica wormholes in the free energy (or more generally, any extensive quantity) is clear, and we have discovered hints of RSB. These suggest that gravity has some features analogous to a glassy phase just at the edge of semiclassicality. Since the gravitational path integral is in general -and in this regime in particular -of clear interest, clearly one important extension of our analysis would be the construction of a gravitational analogue of the Parisi ansatz for RSB. Of course, because we did not work in any saddle point approximation, we did not consider classical equations of motion. The resulting lack of any saddles to analyze for stability or to continue to m = 0 makes it difficult to explore the structure of RSB in any detail. In particular, the fact that (pure) JT gravity replica wormholes do not exist as solutions to any classical equations of motion suggests that there may be no way to study RSB in JT gravity in a way analogous to conventional spin glass systems (though admittedly the possibility of a phase transition at m < 1 means that the lack of on-shell wormholes for integer m does not necessarily exclude on-shell analytically continued wormholes for m near zero). A natural question, then, is whether there exist models of gravity that are sufficiently simple to allow for the continuation of classical equations of motion to m = 0, but sufficiently complex to still exhibit a phase transition. In other words, it would be valuable to find a gravitational model in which the effects of Euclidean wormholes can be disentangled from those of disconnected geometries with higher genus (analogous to the case of CGHS, in which higher genera don't appear at all). In such a model, we might imagine that the correct "gravitational" Parisi ansatz is a multi-branched wormhole connecting the various disconnected boundaries with wormholes of different sizes, with these sizes left as variational parameters with respect to which the free energy should be extremized. In the case of a near-extremal black hole (and consequently low temperature), the picture might be reminiscent of AdS fragmentation [85], in which the AdS 2 throat can fragment into many throats or disconnected universes. Understanding how this story works in gravity would be especially illuminating because the Parisi function q(x), which plays the role of an order parameter for the spin glass phase transition in the SK model, also probes the structure of microstates of the model. An analogous function in gravity could shed light onto the details of the underlying (that is, pre-disorder-average) theory. the ψ i are Chern classes, κ is the first Mumford-Morita-Miller class on M g,m , and we use the notation α = {α 1 , . . . , α m } and |α| = m i=1 α i ; see e.g. [57] for a review. The quantity in parenthesis in equation (A.1a) is the Weil-Peterson symplectic form on the moduli space of bordered Riemann surfaces. Because V g,m ({b i }) is a polynomial in the b i , when inserted into (4.4c) we may explicitly perform the integrations over the b i to obtain Z g,m (β) = α,p |α|+p=3g−3+m (2π 2 ) p p!(2π) m/2 β 3g−3+3m/2−p Mg,m ψ α 1 1 · · · ψ αm m κ p .
At low temperatures, the leading-order behavior of Z g,m comes from the terms in the sum with p = 0; keeping only these terms, (4.3) gives Mg,m ψ α 1 1 · · · ψ αm m + · · · , (A.3) where the ellipses denote terms that are subleading at low temperature.