Replica symmetry breaking in random matrix model: a toy model of wormhole networks

We study the replica symmetry breaking (RSB) in the Gaussian Unitary Ensemble (GUE) random matrix model. We find that the RSB occurs at the transition temperature $T_{\text{RSB}}\sim N^{-2/3}$ in the large $N$ limit. We argue that this transition originates from a landscape of almost degenerate local minima of free energy labeled by Young diagrams. We also discuss a possible implication of our findings for the multiply-coupled Sachdev-Ye-Kitaev models and their holographic dual Jackiw-Teitelboim gravity from the viewpoint of ER=EPR conjecture.


Introduction
The replica symmetry breaking (RSB) in disordered system is an important concept in the study of glass transition. In particular, the spin-glass model of Sherrington and Kirkpatrick [1] was solved by Parisi [2,3] based on the RSB ansatz (see [4] for a review).
In the disordered system a natural quantity to compute is the quenched average of the free energy log Z(β) , where Z(β) = Tr e −βH is the partition function and β is the inverse of temperature T β = 1 T . (1.1) The average · · · in the disordered system is defined by integrating over the distribution of random couplings in the Hamiltonian H. The replica method is commonly used to compute the quenched average log Z(β) by preparing n copies of the system (replicas) and taking the limit n → 0 at the end of calculation log Z(β) = lim In this paper, we consider the RSB in the Gaussian Unitary Ensemble (GUE) random matrix model, where the random couplings in the Hamiltonian are modeled by treating the Hamiltonian H itself as a random N × N hermitian matrix with the Gaussian distribution (1.4) Note that the RSB in the GUE random matrix model has been considered in [5], where the correlator of resolvent Tr(E − H) −1 was studied by the replica method. In this paper we consider the correlator of Z(β) = Tr e −βH and study the RSB characterized by (1.3). 1 We will not consider the limit n → 0 in (1.2); we study Z(β) n with positive integer n in its own right. We find that in the large N limit Z(β) n in the GUE random matrix model indeed exhibits the RSB at the transition temperature T RSB Z(β) n ≈ Z(β) n , (T > T RSB ), Z(β) n c , (T < T RSB ), (1.5) where the subscript c in Z(β) n c denotes the connected part of the correlator. A similar behavior was observed in the so-called random-energy model, where each energy levels are assumed to be Gaussian distributed but the correlations between different energy levels are neglected [6] 2 . It turns out that T RSB for the GUE matrix model scales as We will argue that the transition at T = T RSB originates from the existence of a landscape of almost degenerate local minima of free energy. This landscape of local minima arises as follows: When we expand Z(β) n in terms of the connected correlators, there appear many terms labeled by Young diagrams. Each term of this expansion can be thought of as a local minimum of free energy. We will also discuss a possible implication of our findings for the Sachdev-Ye-Kitaev (SYK) model [8,9] and its holographic dual Jackiw-Teitelboim (JT) gravity [10][11][12]. We examine the idea that the connected part of correlator Z(β) n c corresponds to a dual spacetime geometry where the n boundaries are connected by an n-pronged Euclidean wormhole. We speculate that around the RSB transition temperature the dual gravity side does not corresponds to a classical geometry, but a superposition of random wormhole networks.
This paper is organized as follows. In section 2, we first review the exact result of Z(β) n at finite N in the GUE random matrix model. In section 3, we study the RSB of Z(β) n in the GUE random matrix model using the exact result at finite N . In section 4, we discuss possible implications of our findings in the GUE matrix model for the multiplycoupled SYK models and their dual JT gravity.

Review of the exact results of GUE random matrix model
In this section, we review the known results of the correlation function Z(β) n in the GUE random matrix model (1.4). As pointed out in [13,14], this quantity is formally equivalent 1 The correlator Z(β) n in random matrix model has been studied for a long time in the context of 2d quantum gravity under the name of "loop operators" [7]. 2 In the high temperature regime T > TRSB, taking the limit in (1.2) we find log Z(β) ≈ log Z(β) .
On the other hand, the n → 0 limit in the low temperature regime T < TRSB is rather subtle and the analytic continuation to n = 0 is not straightforward [6]. As emphasized in [6], the correct free energy is not reproduced by simply taking the n → 0 limit of Z(β) n c in (1.5).
to the expectation value of the 1/2 BPS Wilson loops in 4d N = 4 U (N ) super Yang-Mills (SYM) theory under the identification √ λ ↔ 2β, (2.1) where λ is the 't Hooft coupling of N = 4 SYM. Thus we can immediately find the exact finite N expression of Z(β) n by borrowing the known result of N = 4 SYM. We can naturally decompose Z(β) n into the connected components Z(β) j c (j = 1, · · · , n). The relation between Z(β) n and the connected part Z(β) j c is compactly expressed in terms of the generating function, as usual where α is a formal expansion parameter. For instance, Z(β) n with n = 2, 3, 4 are expanded as In general, the expansion of Z(β) n is characterized by the partition of n labeled by a Young diagram Y For example, [1 n ] and [n] correspond to the totally disconnected and the totally connected part, respectively As shown in [15], we can write down the exact finite N result of Z(β) n c in terms of the N × N matrix A(β) where L α n (x) denotes the associated Laguerre polynomial. The one-point function Z(β) is simply given by the trace of A(β) Similarly, the connected part Z(β) n c is given by some combination of A(β) [15] (2.10) For instance, Z(β) 2 c is given by (2.11) In a similar manner, one can write down Z(β) 3 c and Z(β) 4 c as (2.12) Using these exact results (2.11) and (2.12), in the next section we will numerically study the behavior of Z(β) n (n = 2, 3, 4) in (2.3).
Next, let us briefly summarize the known results of Z(β) n in the large N limit. In the large N limit, the eigenvalues of H with the Gaussian measure (1.4) are distributed along the cut [E 0 , −E 0 ] with the eigenvalue density ρ(E) obeying the Wigner semi-circle law In our normalization of the measure (1.4), the "ground state energy" E 0 is given by (2.14) The eigenvalue density (2.13) is a property of a single eigenvalue. One can consider the correlation of n eigenvalues characterized by the correlation function ρ (n) (E 1 , · · · , E n ). It turns out that ρ (n) is written in terms of the kernel K (see e.g. [16]) where ψ (E) is the wavefunction associated with the measure (1.4), i.e. the wavefunction of harmonic oscillator in this case Here H (x) denotes the Hermite polynomial of order . Note that the eigenvalue density is given by the diagonal part of the kernel K (2.18) One can show that the decomposition (2.5) corresponds to the sum over the conjugacy class of σ ∈ S n in (2.18). The totally connected part Z [n] (β) = Z(β) n c comes from the cyclic permutation with length n The leading order behavior of Z(β) n c in the large N expansion can be easily found by borrowing the known results of N = 4 SYM in [17][18][19] via the dictionary (2.1) where I ν (z) denotes the modified Bessel function of the first kind. Note that in the large N limit Z(β) n c scales as where χ is the Euler number of a sphere with n punctures

RSB in the GUE random matrix model
In this section, we will study the RSB in the GUE random matrix model using the exact result reviewed in the previous section 2.

Numerical analysis of RSB
To see the RSB in (1.5), let us consider the ratio between the annealed average Z(β) n and the quenched average Z(β) n We can numerically study the behavior of R n (β) (n = 2, 3, 4) as a function of β using the exact result of Z(β) in (2.9) and Z(β) n c (n = 2, 3, 4) in (2.11) and (2.12), together with the relation (2.3). It turns out that R n (β) is conveniently described by the scaling variable β defined by In Figure 1, we plot R n (β) (n = 2, 3, 4) as a function of β for N = 60 (red curves) and N = 100 (blue dots). As we can see from Figure 1, the plots of R n (β) for N = 60 and N = 100 are almost identical in the region β 2. This justifies a posteriori our definition of the scaling variable β in (3.2). One can see from Figure 1 that R n (β) decays to zero about β ∼ 1. This indicates that the replica symmetry is broken in the low temperature regime β 1, and the critical value of β is roughly given by which corresponds to the temperature T RSB in (1.6). In the next subsection, we consider an analytic explanation of the N -dependence of the transition temperature in (3.3). For n = 2, the decay of R 2 (β) means that in the low temperature region T < T RSB the connected part Z(β) 2 c becomes dominant, and around the transition temperature T ∼ T RSB it takes place the exchange of dominance between the connected part Z(β) 2 c and the disconnected part Z(β) 2 . However, for n ≥ 3 what is happening around T ∼ T RSB is more involved since Z(β) n has many terms (2.5) when written in terms of the combination of connected correlators. We will study the detailed structure for n ≥ 3 in section 3.3.

Analytic estimate of the transition temperature
In this subsection, we will estimate the transition temperature using the large N result of Z(β) n c in (2.20). From the asymptotic behavior of the modified Bessel function As discussed in [20], the combination β 3/2 /N appearing in (3.5) has a natural interpretation as the effective genus counting parameter Indeed, the large N expansion of the exact finite N result (2.9), (2.11) and (2.12) beyond the leading order (2.20) are organized as the genus expansion with the effective string coupling in (3.6) [20]. In terms of this effective string coupling g eff , (3.5) is written as where χ is defined in (2.22) and E 0 = −2 (2.14). (3.7) should be compared with the disconnected part The connected part (3.7) and the disconnected part (3.8) become comparable when the effective string coupling becomes of order 1 From (3.6), this condition explains the N -dependence of the transition temperature β RSB in (3.3). The scaling (3.5) for n = 1 is easily understood from the behavior of the eigenvalue density near the ground state E ∼ E 0 . In the low temperature regime β 1 the dominant contribution to the integral (3.10) comes from the edge of the spectrum near the ground state E = E 0 . Introducing the excitation energy E from the ground state, the eigenvalue density (2.13) near the ground state E ≈ 0 behaves as Then (3.10) is approximated as which reproduces the scaling behavior in (3.5) for n = 1. For n ≥ 2, we could not find a simple scaling argument to derive (3.5) from the integral representation of Z(β) n c in (2.19). It would be interesting to understand the behavior (3.5) directly from (2.19).

Comparison with the spectral form factor
The spectral form factor (SFF), defined by Z(it)Z(−it) , is widely studied as a useful diagnostic of quantum chaos [21]. In SFF, the exchange of dominance between the connected and the disconnected parts also occurs at the so-called dip time t dip , which can be estimated as follows [23]. The disconnected part of SFF behaves as which describes the early time decay of SFF, the so-called "slope." The origin of the factor 1/t 3/2 in (3.14) is essentially the same as the calculation in (3.13). On the other hand, the connected part grows linearly in time, which is known as the "ramp": Equating (3.14) and (3.15) we find the time scale of "dip" The N -dependence of t dip is different from that of the transition temperature β RSB in (3.3). This is not a contradiction since Z(β) 2 and the SFF are sensitive to the different aspects of the two-point correlation ρ (2) (E 1 , E 2 ) of energy eigenvalues: (3.17) Namely, Z(β) 2 is the Laplace transform of ρ (2) with respect to the total energy, while SFF is the Fourier transform of ρ (2) with respect to the energy difference.

Landscape of free energy
The behavior of Z(β) n c in (3.5) can be easily generalized to arbitrary term Z Y (β) in the expansion (2.5) (3.18) When g eff ∼ 1, not only the totally connected and disconnected terms with Y = [n] and Y = [1 n ], but all terms in the expansion (2.5) become comparable. As we did in section 3.1, we can numerically study the behavior of Z Y (β) using the exact result in section 2. Let us define the free energy of Z Y (β) relative to the disconnected part In Figure 2, we plot F Y (β) for various Y 's with |Y | = 3 and |Y | = 4 for N = 60. We have checked numerically that the plot for N = 100 is almost identical to the case of N = 60 in Figure 2 in the regime β 1.5. As one can see from Figure 2, all F Y (β) become almost degenerate around the transition point β ∼ 1. For |Y | = 4 in Figure 2b, the contribution of Y = [2 2 ] is relatively small compared to the other partitions, but F [2 2 ] is still of order O(N 0 ) and it is not suppressed by a negative power of N . Near the transition around T = T RSB , all terms Z Y (β) in the expansion (2.5) contribute with almost equal weight. This transition is not a true phase transition in the sense of Ehrenfest classification: it does not involve any discontinuity in the derivative of free energy. This is in contrast to the confinement/deconfinement transition in 4d N = 4 SYM and its dual Hawking-Page transition [24,25], which is very likely a first order phase transition [26,27]. In that case the free energy changes from O(N 0 ) to O(N 2 ) at the transition, but in our case there is no such large change of free energy: the free energy stays the same order below and above the transition. In fact, the exponential factor e −nβE 0 of Z Y in (3.18) is common for all Y . Moreover, the number of terms (or the number of partitions p(n)) in the expansion (2.5) grows rapidly as n increases When the temperature is below the RSB transition T < T RSB , the totally connected contribution Y = [n] becomes dominant. One can see this explicitly for n = 3, 4 in Figure 2: the blue curve labeled by Y = [n] becomes the largest contribution in the low temperature regime β 1. Thus we confirmed numerically the behavior claimed in (1.5).

Low temperature regime
In the previous subsection 3. 3 we have seen that the totally connected part Z(β) n c becomes dominant in the low temperature regime T < T RSB . We find that in this regime Z(β) n c is further approximated by the first term in (2.11) and (2.12) Z(β) n c ≈ Tr A(nβ) = Z(nβ) , (T < T RSB ). (3.21) To see this, let us consider the ratio r n (β) between the both sides of (3.21) In Figure 3, we plot r n (β) for n = 2, 3, 4 as a function of β = β/N 2/3 for N = 60 (red curve) and N = 100 (blue dots). As we can see from Figure 3, r n (β) approaches 1 as β increases, which confirms the claim in (3.21). Finally, combining (3.21) and (1.5), in the low temperature regime we find

Possible implications for multiply-coupled SYK models
In this section, we will discuss possible implications of our findings in the GUE random matrix model for multiply-coupled SYK models and their dual JT gravity. The SYK model exhibits a chaotic behavior which is diagnosed by the SFF [23] as well as the out-of-time ordered correlators (OTOCs) [8,9,28]. The Lyapunov exponent extracted from the OTOC of SYK model saturates the chaos bound [29] which is consistent with the holographically dual black hole interpretation.
The level statistics of energy eigenvalues of the SYK model is well approximated by the random matrix models [22,23], which is also a characteristic feature of the quantum chaos. At low energy, the SYK model is described by the Schwarzian theory whose density of states is given by [30,31] ρ Sch (E) = sinh(2π Although this is different from the simple semi-circle law in (2.13), ρ Sch (E) exhibits the same square-root behavior ρ(E) ∼ √ E (3.12) near the edge of the spectrum E = 0. The above discussion suggests that the RSB transition observed in the GUE random matrix model in the previous section 3 also appears in the multiply-coupled SYK models, where the expectation value of Z(β) n is taken by the disorder average over the random coupling in the Hamiltonian of SYK model. Recall that the scaling of the effective string coupling g eff in (3.6) is sensitive only to the behavior of eigenvalue density near the edge of the spectrum. Thus the behavior of g eff is expected to be the same for both the GUE random matrix model and the Schwarzian theory, and the transition occurs when g eff ∼ 1.
However, the RSB and the role of replica non-diagonal saddle point in the SYK model are still controversial [32][33][34][35]; It is argued in [32] that RSB does not occur in the SYK model. In [33][34][35] the replica non-diagonal saddle points are found by numerical analysis but it is argued that such saddle points are always sub-dominant compared to the diagonal one. 4 In a recent paper [37], it is argued that a non-perturbative completion of Schwarzian theory and its dual JT gravity is given by a certain double-scaled matrix model, and it is suggested that the replica non-diagonal saddle point does play an essential role.
Interestingly, it is shown in [38] that by turning on a mass term connecting the fermions living on different replicas, the partition function of coupled SYM models exhibits a first order phase transition, which is very similar to the transition we observed for Z(β) 2 in the GUE random matrix model (see also [39] for the symmetry breaking in coupled SYK models).
In the rest of this section, we will have in mind a certain soft breaking of replica symmetry as in [38], and see what kind of dual gravity picture is expected based on the observation in the previous section 3. Our argument is merely a heuristic one and we have not performed an actual calculation of the multiply-coupled SYK models with n ≥ 3 replicas generalizing [38]. It would be very interesting to study such multiply-coupled SYK models.
Based on the fact that the level statistics of the SYK model is described by a random matrix model, we expect that the phase structure of the GUE random matrix model we observed in the previous section captures the essential part of the low energy dynamics of multiply-coupled SYK models.

Wormhole networks
The bulk gravity picture of n = 2 coupled SYK models in [38] is that the two terms in the expansion of Z(β) 2 in (2.3) correspond to the two geometries with different topology: the first term Z(β) 2 corresponds to the two disconnected disks associated with the two Euclidean black holes, while the second term Z(β) 2 c corresponds to the Euclidean wormhole connecting the two boundaries with the topology of annulus. This is depicted schematically as: This picture for n = 2 replicas suggests that the holographic dual geometry of the connected correlator is literally connected by an Euclidean wormhole, in the same spirit as the ER=EPR conjecture [40]. Assuming this interpretation, we can draw a cartoon of spacetime geometry for all terms Z Y (β) in the expansion (2.5). For n = 3 we have three terms Y = [1 3 ], [1 · 2], [3] in the expansion (2.3), where the corresponding spacetime geometries look like and Z [3] (β) = Z(β) 3 c = . (4.5) Note that the coefficient 3 in Z [1·2] (β) = 3 Z(β) Z(β) 2 c has a natural interpretation as the number of ways to connect two boundaries by a wormhole out of three boundaries, as shown in (4.4). One can also check that the coefficient appearing in the expansion of Z(β) 4 in (2.3) has a similar interpretation.
For general n, from the result of random matrix model in (1.5), we can draw the following dual spacetime picture: at high temperature T > T RSB the bulk geometry is n disconnected disks representing n disconnected Euclidean black holes, while at low temperature T < T RSB the corresponding bulk geometry is an n-pronged wormhole connecting the n boundaries (see (4.5) for the example of three-pronged wormhole). Around the RSB transition temperature T ∼ T RSB , we speculate that the dual gravity side does not correspond to a classical geometry, but it is described by a superposition of random wormhole networks. It would be very interesting to understand the implication of this RSB transition in the bulk gravity theory.

Thermofield n-tuple and ground state
In the previous section 3.4, we observed that Z(β) n is reduced to Z(β) n c ≈ Z(nβ) in the low temperature regime (3.23). Let us consider a possible interpretation of this behavior. We first notice that our starting point Z(β) n is written as a trace over the n-th power of the Hilbert space H of a single system Z(β) n = Tr H ⊗n e −βH ⊗ · · · ⊗ e −βH n . (4.6) Since Z(β) n is defined as a trace in H ⊗n , it is natural to express its low temperature limit Z(nβ) in (3.23) as a quantity involving H ⊗n . However, this seems impossible since Z(nβ) = Tr H e −nβH is a trace on the single Hilbert space H. This tension is resolved by assuming that the low temperature state in H ⊗n is highly entangled. We propose that Z(nβ) should be written as Z(nβ) = TF n |TF n , (4.7) where |TF n denotes the "thermofield n-tuple state" |TF n = i e − nβ 2 E i |i ⊗ · · · ⊗ |i n ∈ H ⊗n . (4.8) This is a generalization of the thermofield double state |TFD for n = 2 |TFD = |TF 2 . (4.9) In the low temperature limit, the partition function is approximated by the expectation value in the ground state. From (3.23) and (4.7), we expect that the ground state |Ψ gnd of multiply-coupled system is approximated by the thermofield n-tuple state |Ψ gnd ≈ |TF n . (4.10) This property is indeed observed for n = 2 in [38]. Our relations (3.23) and (4.7) are also consistent with the picture of spacetime geometry built by entanglement [41] and the ER=EPR conjecture [40]. In fact, the relation between the thermofield double state |TFD and the wormhole picture of Z(β) 2 c depicted in (4.2) is the motivation for the proposal in [40,41]. Then it is natural to generalize this picture to the case of multiple boundaries. For instance, the thermofield triple state |TF 3 corresponds to the three-pronged wormhole depicted in (4.5). 5 We emphasize that only the totally connected part is dominant in the low temperature limit, and hence the ground state |Ψ gnd ≈ |TF n would be dual to an n-pronged wormhole. In the spirit of [40,41], this n-pronged wormhole can be thought of as a geometric manifestation of the entanglement in the state (4.8).