Symmetry Breaking in Coupled SYK or Tensor Models

We study a large $N$ tensor model with $O(N)^3$ symmetry containing two flavors of Majorana fermions, $\psi_1^{abc}$ and $\psi_2^{abc}$. We also study its random counterpart consisting of two coupled Sachdev-Ye-Kitaev models, each one containing $N_{\rm SYK}$ Majorana fermions. In these models we assume tetrahedral quartic Hamiltonians which depend on a real coupling parameter $\alpha$. We find a duality relation between two Hamiltonians with different values of $\alpha$, which allows us to restrict the model to the range of $-1\leq \alpha\leq 1/3$. The scaling dimension of the fermion number operator $Q=i\psi_1^{abc} \psi_2^{abc}$ is complex and of the form $1/2 +i f(\alpha)$ in the range $-1\leq \alpha<0$, indicating an instability of the conformal phase. Using Schwinger-Dyson equations to solve for the Green functions, we show that in the true low-temperature phase this operator acquires an expectation value. This demonstrates the breaking of an anti-unitary particle-hole symmetry and other discrete symmetries. We also calculate spectra of the coupled SYK models for values of $N_{\rm SYK}$ where exact diagonalizations are possible. For negative $\alpha$ we find a gap separating the two lowest energy states from the rest of the spectrum; this leads to exponential decay of the zero-temperature correlation functions. For $N_{\rm SYK}$ divisible by $4$, the two lowest states have a small splitting. They become degenerate in the large $N_{\rm SYK}$ limit, as expected from the spontaneous breaking of a $\mathbb{Z}_2$ symmetry.

physics, including the strange metals [3,[26][27][28][29][30][31]. With such applications in mind, it is interesting to study various dynamical phenomena in the SYK and tensor models. For example, phase transitions in such models have been studied in [32][33][34]. In this paper we identify a simple setting where spontaneous symmetry breaking can occur: two SYK or tensor models coupled via a quartic interaction. We take this interaction to be purely melonic (i.e. tetrahedral in the tensor model case), so that the symmetry breaking can be deduced from the large N Schwinger-Dyson equations.
In the random case, we will study two coupled SYK models with the Hamiltonian where, as usual, all repeated indices are summed over. The Majorana fermions are χ i 1 and χ i 2 with i = 1, . . . , N SYK , and J ijkl is a fully anti-symmetric real tensor with a Gaussian distribution. 1 We will show that the real parameter α may be restricted to the range −1 ≤ α ≤ 1/3 by a duality symmetry. This quartic Hamiltonian, which couples 2N SYK Majorana fermions, is invariant under an anti-unitary particle-hole symmetry [35][36][37][38][39][40] generated by P; see eq. (3.11). However, we will show that for −1 ≤ α < 0 this Z 2 symmetry is spontaneously broken when N SYK is divisible by 4 and taken to infinity. 2 In this limit the fermion number operator Q = iχ j 1 χ j 2 acquires an expectation value. This leads to a gapped phase in two coupled SYK models similar to that found by Maldacena and Qi [41] (for further results see [42]); however, instead of the quartic they assumed a quadratic coupling term µQ which breaks the Z 2 symmetry explicitly. This gapped phase was argued to be dual to a traversable wormhole in two-dimensional gravity [43,44], and our model (1.1) may have a similar interpretation for −1 ≤ α < 0.
As we show in section 2.4, a sign of instability of the conformal phase for −1 ≤ α < 0 is the presence of a complex scaling dimensions of the form 1/2 + if (α). Appearance of complex dimensions with real part equal to d/2 for some single-trace operators is a common phenomenon in large N models [45][46][47][48][49]. Via the AdS/CFT correspondence [50][51][52], such operators are related to scalar fields which violate the Breitenlohner-Freedman stability bound [53]. The fact that α = 0 is the lower edge of the conformal window is related to appearance of the marginal double-trace operator Q 2 there. For 0 < α ≤ 1/3 there are actually two fixed points connected by the flow of the coefficient of Q 2 , but at α = 0 they merge and annihilate, as explained for example in [54,55].
The complex scaling dimensions have been observed in bosonic tensor models [56,57], as well as in a complex fermionic model introduced in [6] following the work in [58]. This fermionic model is often called "bipartite" because of the two types of interaction vertices (black and white) arranged in an alternating fashion, since the propagator must connect different vertices. The bipartite model was further studied in [17] and shown to possess a complex scaling dimension of the operatorψ abc ψ abc . Here we generalize this tensor model to one with a continuous parameter α in such a way that the bipartite model corresponds to For α = 0 this describes two decoupled copies of the basic Majorana O(N ) 3 model with the tetrahedral interaction [6]. The coupling term proportional to α preserves its discrete symmetries and also has the tetrahedral structure, i.e. every two tensors have only one index contraction, so that the model (1.2) is melonic. It is the tensor counterpart of the coupled SYK model (1.1), and in the large N limit it is governed by the same Schwinger-Dyson equations for the two-point and four-point functions. 3 In section 2 we derive the Schwinger-Dyson equations and use them to study the scaling dimensions of various O(N ) 3 invariant fermion bilinears. We also exhibit a duality symmetry which allows us to restrict the model to the range −1 ≤ α ≤ 1/3. The nearly conformal phase of the theory is stable for 0 ≤ α ≤ 1/3, but it is unstable for −1 ≤ α < 0 as signaled by the complex scaling dimension of operator iψ abc 1 ψ abc 2 . The true behavior of the theory with negative α is the spontaneous breaking of the particle-hole Z 2 symmetry, as we demonstrate in section 3. In section 3.1 and 3.2 we numerically study the large N Schwinger-Dyson equations and exhibit the exponential decay of correlators at low temperature. We also ascertain the existence of second-order phase transitions by numerically computing the free 3 In [28,30] quartic interactions were added to SYK models, which have a "double-trace" structure and contain an additional random tensor C ij . These interactions do not have a tensor counterpart because ψ abc 1 ψ abc 1 is a c-number. energy. In section 3.3 we study the numerical spectrum of the coupled SYK model (1.1) via exact diagonalizations at finite N SYK . We observe that for −1 ≤ α < 0 there is a gap separating the two lowest energy states from the rest of the spectrum. For N SYK divisible by 4 there is also a small gap between the two lowest states, consistent with the fact that the ground state must be non-degenerate [35][36][37][38][39][40], but this gap decreases as N SYK is increased.
In the large N SYK limit, the two lowest states become degenerate and give rise to the two inequivalent vacua, which are present due to the spontaneous breaking of the Z 2 particle-hole symmetry.
This means that the low-temperature entropy is large for 0 ≤ α ≤ 1/3 but vanishes for −1 ≤ α < 0. It is tempting to suggest that the latter case is dual to a wormhole.
This senstivity to the sign of the interaction coupling two CFTs is like in [43], where the traversable wormhole appears only for one of the signs. 4

Schwinger-Dyson Equations and Scaling Dimensions
In this section we study the two-flavor tensor model with Hamiltonian (1.2). 5 It can be compactly written in the form where the capital letters are a shorthand notation for three tensor indices: I = a 1 b 1 c 1 , J = a 2 b 2 c 2 , etc, and the non-random tetrahedral tensor coupling consists of six terms The tensor J IJKL is antisymmetric under permutation of indices I, J, K, L and has a tetrahedron topology as shown in figure 1. In the form (2.1) the tensor model Hamiltonian is transparently similar to the SYK one (1.1). In terms of the complex tensors the Hamiltonian (2.1) assumes the form (2.4) Figure 1: Pictorial representation of the antisymmetric tensor J IJKL . The where A, B, and C are orthogonal matrices. In addition, it has a particle-hole Z 2 symmetry 6 generated by [35][36][37][38][39][40], where K is the anti-unitary operator which acts by The fermion number operator does not in general commute with H, but it is conserved mod 4. The particle-hole symmetry is not anomalous only if the total number of fermions 2N 3 is a multiple of 8, i.e. when N is even [35][36][37][38][39][40]. Even in this case, we will argue that in the large N limit the symmetry is spontaneously broken for −1 ≤ α < 0 because Q acquires an expectation value.

Duality in the Two-Flavor Models
In this section we show that the two-flavor models with different values of α can be equivalent.
We will demonstrate this explicitly in the tensor model case (2.1), but the SYK case (1.1) works analogously. Let us perform the following transformation on the Majorana fermions: It preserves the anticommutation relations, and turns the Hamiltonian (2.1) into 7 (2.10) Thus the energy levels are symmetric under the duality transformation (2.14) This means that the fundamental domain is −1 ≤α ≤ 1. Thus, we may restrict α to the The values of α outside of this domain are related to it by the duality. For α = −1 the transformation (2.11) maps the theory into itself, but with H → −H.

16)
7 Using antisymmetry of the tensor J IJKL one can operate with Majorana fermions as commuting variables but keeping order of I, J, K, L indices fixed. 8 For the original Hamiltonian (2.14 this transformation rescales the energy levels. Therefore, our results for dimensionful quantities, like energy levels and Green functions, will not respect the duality under (2.13).
where we introduced the complex tensor ψ I = 1 The theory with α = 1/3 is mapped into itself by (2.11). In this case the Hamiltonian is which has O(N ) 3 × U (1) symmetry. In the three-index notation This is different from the SU (N ) 2 × O(N ) × U (1) symmetric complex tensor model [6]; the latter involves taking only the first term in this Hamiltonian.

Feynman rules and two-point functions
At first we list the Feynman rules which follow from the Hamiltonian (1.2). In figures (2) and (3) we define propagators and interaction vertices for the given two-flavour tensor model. Since the interaction terms have a tetrahedral tensor structure the melonic Feynman  diagrams dominate in the large N limit. Let us define bare two-point functions where the sum over indices I is assumed. The leading melonic correction to the full two-point function G 11 is represented in figure 4. Using that we find

(2.21)
A similar expression can be derived for G 22 . Since there is a symmetry ψ 1 ↔ ψ 2 we can assume that G 11 = G 22 = G and obtain a Schwinger-Dyson equation for the full two-point function (see figure 5) where G(τ ) = 1 2 sgn(τ ) is the bare propagator. In writing this Schwinger-Dyson equation we implicitly made an important assumption that the two-point functions are zero G 12 (τ ) = G 21 (τ ) = 0. This follows from the Z 2 symmetry ψ 2 → −ψ 2 . As we will see below, the Z 2 symmetry can be spontaneously broken for some range of parameter α and dimensionless coupling βJ, where β = 1/T is the inverse temperature and J 2 = g 2 N 3 is effective coupling constant.
Let us first assume that Z 2 symmetry is not broken and analyze the SD equation (

Scaling dimensions of bilinear operators
We can use the large N Schwinger-Dyson equations for the three-point functions to deduce the scaling dimensions of four families of bilinear operators: where n = 0, 1, 2, . . . , and the sum over tensor indices is assumed. 9 Each of these operators is invariant under the O(N ) 3 symmetry, but they are distinguished by their transformations to discrete symmetry.
We take some operator O(τ ) and consider two three-point functions of the form where we assume summation over the index I. In the large N limit the functions (2.26) obey the melonic Bethe-Salpeter equations. They are schematically represented in figure 6.
In the conformal limit one can ignore the first diagram on the right and obtain  In the coupled SYK model (1.1) the same expressions for bilinear operators are applicable after replace- where assuming that G 11 = G 22 = G we find and K c is the kernel of the SYK model, defined in conformal limit as An arbitrary conformal three-point function of the form (2.26) with an operator of scaling and obviously must be antisymmetric under τ 1 ↔ τ 2 . This three-point function is an eigenvector of the kernel K c with the eigenvalue g(h): 10 To solve (2.27) one has to find eigenvalues of the matrix and equate them to unity. This gives an equation for possible scaling dimensions. It easy to see that this matrix acquires diagonal form in the basis of vectors v 11 + v 22 and v 11 − v 22 and we find two equations for the scaling dimensions The scaling dimensions of the operator O 2n+1 are independent of α. They are given by the well-known series h = 2.00, 3.77, 5.68, 7.63, 9.60, . . . which approaches 2n + 3 2 . These are the same scaling dimensions as in the basic O(N ) 3 tensor model [6] and the SYK model. On the other hand, the scaling dimensions of operators 1+3α 2 g A (h) = 1 and depend on α. As a check we note that for α = 0 the spectra of O 2 and O 1 are the same; this is as expected since the 10 To take the integrals one should use star-triangle identities twice [4]. two flavors are decoupled. Now consider the last possible three-point function The melonic Bethe-Salpeter equation for this three-point function is represented in figure 7 and in the conformal limit, neglecting the first diagram on the right we get In this case there are two general possibilities for conformal three-point function, namely anti-symmetric and symmetric cases (2.35) Therefore, we find equations which determine spectra of antisymmetric and symmetric op- The and indeed for α = −1 we get 6(α−α 2 ) 1+3α 2 g S (h) = g sym (h).
To summarize, we have found that scaling dimensions of the operators (2.25) can be obtained by solving equations g i (h) = 1, where The duality relation (2.11) is reflected in the behavior of functions g i (h), which define scaling dimensions of the operators O i . Using (2.38) and (2.11) one finds

Complex scaling dimensions
In this section, we examine if there exist any complex solutions of the equations g i (h) = 1 defined in (2.38). If such a complex root exists, then a conformal primary has a complex scaling dimension, which leads to a destabilization of the model. Indeed, a complex scaling dimension of the form 1 2 ± if corresponds to a scalar fields in AdS 2 whose m 2 is below the Breitenlohner-Freedman bound m 2 In such a case one may expect "tachyon condensation" in AdS space. In the dual CFT the operator dual to the tachyon acquires an expectation value, leading to symmetry breaking.
We will obtain some support to this picture.
First of all we notice that the functions g A (h) and g S (h) are real only if h is real or Using the fact that − 1 3 ≤ 1−α 2 1+3α 2 ≤ 1 (and the same for 2α(1+α) 1+3α 2 due to duality) we conclude that equations g i The plot of f (α) is shown in figure 8. For slightly negative α we find while f (−1) ≈ 1.5251 in agreement with the result for the bipartite model found in [17].
Thus, for −1 ≤ α < 0 there is an operator with the complex scaling dimension h = 1/2+if (α) or its complex conjugate: the fermion number operator Q = iO 0 4 . This makes the conformal large N limit unstable. operator Q 2 . In the large N limit, ∆ Q 2 = 2∆ Q . In one of the CFTs, ∆ Q = h + , so that ∆ Q 2 > 1. Since the operator Q 2 is irrelevant, this CFT is stable. There is RG flow leading to it, which originates from another large N fixed point where Q 2 is relevant [61,62]. At this UV fixed point, ∆ Q = h − . When α = 0 the two fixed points merge and annihilate, as in various other theories, for example [45,46,54,55]. For −1 ≤ α < 0 there are two different theories containing complex dimension ∆ Q = 1/2 + if (α) or its complex conjugate. They may be formally regarded as "complex CFTs" [49], but we will see in the next section that their true physics includes symmetry breaking, which leads to a gap in the energy spectrum.

Symmetry Breaking
In section 2.4, we showed that for the coupled tensor model (1.2) in the range −1 ≤ α < 0 the fermion number operator Q = iψ I 1 ψ I 2 has a complex scaling dimension, signaling an instability of the conformal phase of the model. In this section we show that this operator acquires a vacuum expectation value (VEV) in the true low-temperature phase of the large N model. Based on this, it is tempting to make the following conjecture.
Conjecture. If the assumption of conformal invariance in a large N theory leads to a singletrace operator with a complex scaling dimension of the form d/2 + if , then in the true lowtemperature phase this operator acquires a VEV.
In our case, the O(N ) 3 symmetry implies that where we used the short-handed notation I = abc, and A is of order 1 in the large N limit.
This leads to an exponential decay of correlation functions and signifies a gap in the energy spectrum. Furthermore, the VEV (3.1) implies that various discrete symmetries, including the particle-hole symmetry (2.6), the interchange symmetry between ψ 1 , ψ 2 , and the reflection symmetry ψ 2 → −ψ 2 , are spontaneously broken. Therefore, one should expect a second-order phase transition between the broken and unbroken symmetry phases. In addition, the spontaneously broken symmetry also implies a ground state degeneracy in the large N energy spectrum. 11 In this section we extensively analyze the phenomenon of symmetry breaking, sometimes Let us first demonstrate the connection between the tensor model and the SYK counterpart. For the one-flavor O(N ) 3 tensor model the analogous SYK model has the random tensor J ijkl which is fully antisymmetric. The mixed term A ijkl χ i 1 χ j 1 χ k 2 χ l 2 has only the symmetries which are the same as for the Riemann tensor. However, the full interaction term following Since A ijkl + A iljk + A iklj is fully antisymmetric due to (3.2), the mixed term has a fully antisymmetric random coupling. We will assume that it is proportional to the coupling in the diagonal term of (1.2), and are thus led to the random model (1.1). This model is the special M = 2 case of a periodic SYK chain model where the integer x labels the lattice site, and χ i M +1 ≡ χ i 1 . This can be obtained from the model of [26] by identifying the separate random couplings up to a factor of α.
Introducing the complex combination ψ j = 1 √ 2 (χ j 1 + iχ j 2 ), we may write the Hamiltonian (1.1) as As usual, we will assume that each variable J ijkl has a gaussian distribution with variance 6J 2 N −3 SYK . We will typically state energies in units of J, or equivalently set J = 1.
The fermion number operator does not in general commute with H, but it is conserved mod 4, just like in the Maldacena-Qi model [42]. For α = 1/3, however, we find the Hamiltonian Thus, we have enhanced U (1) symmetry ψ j → e iγ ψ j . 12 We note that, for α = 1/3 the symmetry generated by π 2 rotation R in χ 1 , χ 2 . any two broken symmetries that can be related by an unbroken symmetry do not produce extra ground state degeneracy, and therefore it is enough to focus on one of them.
Let us focus on the particle-hole symmetry [35][36][37][38][39][40] generated by It acts on the fermion number as For N SYK not divisible by 4, there is a two-fold degeneracy of the ground state in section 3.3, due to an anomaly in the particle-hole symmetry [35][36][37][38][39][40]. For N SYK divisible by 4 this symmetry is not anomalous, and we find a non-degenerate ground state, which is followed by a nearby state when −1 ≤ α < 0. The two lowest states become degenerate in the large N SYK limit, and they are separated by a gap from the remaining states. This leads to a spontaneous symmetry breaking through the formation of an expectation value of Q. We will demonstrate this effect by solving the large N SYK Schwinger-Dyson equations for the Green functions, and with diagonalizations at finite N SYK .

Schwinger-Dyson equations and the effective action
In this section we derive the large N SYK effective action of GΣ type, and the Schwinger-Dyson equations, for the coupled SYK model (1.1). Following [41], we introduce bi-local variables 13) and the corresponding Lagrange multipliers Σ ab (τ, τ ), where a, b = 1, 2. The effective action is given by

By translation invariance
(3.14) We also have the general properties The Schwinger Dyson (SD) equations for the two point functions assume the form 13

Solutions of Schwinger-Dyson equations and symmetry breaking
For 0 ≤ α ≤ 1/3 there are no operators with complex scaling dimensions, so it is consistent to assume that the discrete symmetries are unbroken and set G 12 = 0, and G 11 = G 22 , to obtain a nearly conformal solution in the low energy limit. However, the appearance of a complex scaling dimension for −1 ≤ α < 0 shows that such a conformal phase is unstable.
We will show that, in this range of α the true phase of the theory exhibits spontaneous symmetry breaking.
In order to exhibit it, we have to allow the possibility that G 12 (τ ) = 0. The underlying Z 2 symmetry of the Hamiltonian (1.1) implies that such solutions must come in pairs related by G 12 (τ ) → −G 12 (τ ) (in our numerical work we will typically exhibit only one of these two solutions). They correspond to working around the two inequivalent vacua, which we will call |0 + and |0 − . They are distinguished by the sign of the expectation value of operator 13 These equations are also valid in the two-flavor tensor model (1.2), where G ab (τ ) = 1 The unbroken symmetry R in (3.10) implies 18) and similarly for Σ ab . Using these constraints, we obtain for the effective action The Schwinger Dyson equations become and J −2 Σ 11 (τ ) = (1 + 3α 2 )G 3 11 (τ ) + 6α(1 − α)G 11 (τ )G 2 12 (τ ) , J −2 Σ 12 (τ ) = (1 + 3α 2 )G 3 12 (τ ) + 6α(1 − α)G 2 11 (τ )G 12 (τ ) . (3.20) may be written in momentum space as (3.22) These equations, together with (3.21), can be solved numerically using the method of weighted iterations used in [19]. 14 To trigger the spontaneous symmetry breaking, we start our iteration process with a tiny non-zero G 12 (τ ) which is purely imaginary. If we are in the unbroken phase, after the iterations G 12 becomes zero; whereas if we are in the broken phase we find a non-zero purely imaginary solution for G 12 .
The plots of G 11 and iG 12 for different values of α and βJ are shown in fig. 9. For each 14 In this case we find it more convenient to use a slow decay rate on the weight x.  are the ones with the lower free energy. As βJ decreases, |G 12 (τ )| decreases everywhere for the non-trivial solution (see figure 9, 10), and at the critical value becomes exactly zero. For βJ < (βJ) crit the only possible solution is G 12 (τ ) = 0. Thus, the Z 2 symmetry is restored, and this is a second-order phase transition. The plot of (βJ) crit vs. α is shown in figure 11; it blows up as α approaches zero from below. 15 Using the solutions of the Schwinger-Dyson equations we can numerically compute the large N free energy where the sum n log(−iω n ) is replaced by log(2). The energy can be computed with the and at low temperatures it should converge to the energy of the ground state E 0 divided by Symmetry broken phase Figure 11: Critical value of βJ as a function of α.
Now one can compare the free energy in the symmetry broken phase, F G 12 =0 with that of the symmetry unbroken phase, F G 12 =0 . In particular, the free energy of the latter is simply twice that of a single SYK with a rescaling J → √ 1 + 3α 2 J. As a check, in such cases one can consider a large q expansion [19,63], (3.25) where βJ (1 + 3α 2 )2 1−q q = πν cos πν

2
. The free energy of the symmetry unbroken phase F G 12 =0 is seen to agree well numerically with F 4 G 12 =0 .
In figure 12 we plot for α = −1 the free energy of the symmetry broken phase (3.23) as a function of βJ and compare it with that of the unbroken phase, obtained by setting G 12 = 0 in the SD equations (3.20) and (3.21). We also show the entropy as a function of βJ. The plot shows a clear second order phase transition at (βJ) crit ≈ 2.87, and the derivative of the entropy is discontinuous. We will systematically study the critical exponents in future work. The graph on the right shows the entropy; we can clearly see a second order phase transition, as there is a discontinuity in its derivative near critical temperature. Figure 13: Large N free energies at fixed β and J. We take β = 5, J = 1, and decrease α.
We observe also a second order phase transition.
We notice that at sufficiently large βJ, there is a range of τ where both iG 12 (τ ) and G 11 (τ ) decay exponentially and share the same decay rate. To explain this fact, let us study the T = 0 case and insert the complete set of states G 11 (τ ) = 0 + |e −Hτ χ 1 1 (0)e Hτ |n n|χ 1 1 (0)|0 + . (3.26) For large τ the sum is dominated by the lowest excited state, and we find Similarly, we find that the large τ behavior of G 12 is Thus the universal decay rate among correlators signifies a mass gap in the spectrum.
In the work of Maldacena and Qi [41] the functions G 11 and G 12 were also found to be exponentially decreasing for sufficiently large βJ. In fig. 14 we exhibit superimposed plots of the low temperature solutions to our system of equations and those from [41], with parameters chosen so that the solutions are close to one another for most of the range.
We observe a difference in the behavior of iG 12 (τ ) and iG LR (τ ) at small τ : in our case the function is smooth with a vanishing derivative at τ = 0, while in [41] its derivative is discontinuous at τ = 0; this is due to the fact that their Hamiltonian includes a quadratic term. Figure 14: Plot of solutions G LL and iG LR for the model in [41] superimposed with G 11 and iG 12 . Parameters chosen so that the solutions are close for most of the range of θ.
We can also study what happens at low temperatures (large βJ) as a function of α. In figure 15 we plot iG 12 (0), which is the expectation value of the order parameter Q for a large βJ. This quantity becomes small as α is increased towards zero. In figure 16 we plot the large N SYK limit of the energy gap E gap divided by J, calculated from the exponential decay of the Green functions. We also plot the ground state energy E 0 divided by JN SYK calculated using

Exact diagonalization for finite N SYK
In this section we present numerical results for the spectra of two coupled SYK models with Hamiltonian (1.1). We first check that the results from exact diagonalizations agree well with expectations: the spectrum for α = 0 and N SYK = 30, and the ground state energy of α = −1 for various N SYK concur well with analytical arguments, and with the results from 3.2. Then we present our results on the energy gap and broken symmetry.
The biggest number we are able to access via exact diagonalization of the coupled SYK models is N SYK = 16. In this case the discrete symmetry (3.11) is not anomalous, and the ground state is non-degenerate. However, for −1 ≤ α < 0 we observe a nearby excited state followed by a gap. We will interpret this as indication of approach to spontaneous symmetry breaking, which takes places in the large N SYK limit. We will also present spectra for N SYK = 15, where the discrete symmetry (A.5) is anomalous, so that the states are doubly degenerate. There is again a gap in the spectrum present for −1 ≤ α < 0. Furthermore, we will present numerical results on the VEV of operator iχ i 1 χ i 2 for N SYK = 14, which demonstrates that it is non-vanishing for −1 ≤ α < 0. First, let us consider α = 0, where we find the spectrum of two SYK model with the same random couplings. The density of states for this model is simply given by the convolution of that of the single SYK model: 16 ρ double (E) = deρ(e)ρ(E − e) (3.29) This in particular helps us determine the behavior of ρ double (E) near the ground state. Shifting the energy so that the ground state is at zero, we know that ρ(E) → A √ E for small E.
Therefore, for small E The numerical density of states, shown in figure 17 for N SYK = 30, is in good agreement with the E 2 dependence near the ground state. 16 We thank D. Stanford for a useful discussion about this.
Let us proceed to the spectra for non-vanishing values of α. In figs. 18, 19 we plot the spectra of energy divided by J for α = −1, −1/2, 1/3 and different values of N SYK . These energy distributions have interesting and unusual shapes. For the special values α = −1 and 1/3 we observe large numbers of states with E = 0; this creates the zero-energy peaks seen in the graphs. For α = −1 and odd N SYK we find that the E = 0 peak is separated by gaps from the remaining states, but for even N SYK it is not.
For N SYK = 15, due to the anomaly in particle-hole symmetry, there are two degenerate ground states, see fig. 18. 17 In fact, each energy level is doubly degenerate. Nevertheless, for −1 ≤ α < 0 we observe a gap between the lowest energy level and the next one, as expected.
On the other hand, for N SYK = 16 there is no exact degeneracy of the ground state, but the first gap is very small, indicating a tendency towards spontaneous symmetry breaking at large N SYK . We show the N SYK = 16 spectra for α = −1 and α = −0.5 in fig. 18. In both cases, for a typical sampling of the coupling constants J ijkl we observe two closely spaced states followed by a visible gap. For large N SYK the energy gap between the two lowest states is expected to decrease exponentially: For α ≥ 0 the low-lying spectrum is different -we observe many closely spaced low-lying states without large gaps, similarly to the standard SYK spectrum.
In fig. 20 figure 20 we also exhibit the energy gap between second and third states as a function of α. As α is increased from −1 to 0, the gap decreases as expected.
Exact diagonalizations also provide support for the statement that the fermion number Q    for a finite number of degrees of freedom [35][36][37]39,40]. In figure 21 the vacuum expectation value as a function of α is plotted for N SYK = 14. This is the finite N SYK analogue of fig.   15, where the large N SYK limit of the condensate is plotted. We also note the qualitative similarity of the plot 21 and that of the imaginary part of the scaling dimension of Q in fig.   8. A More on the discrete symmetries The model (1.1) has the anti-unitary particle-hole Z 2 symmetry generated by (3.11). The operator K is defined to take z →z, z ∈ C but acts as the identity on ψ orψ. It may be identified as a kind of time-reversal generator which satisfies K 2 = 1 [35][36][37]. It acts by and therefore, satisfies Note R 2 = (−1) F . There are also various reflection Z 2 symmetries that are spontaneously broken by the VEV of Q. In particular, we have the reflection symmetry: such that P P † = 1 , P χ i 1 P † = −χ i 1 , P χ i 2 P † = χ i 2 , P 2 = 1. (A. 6) In fact, R, P, and K are enough to generate all discrete symmetries of the model (1.1).
In particular, all the unitary discrete symmetries form D 4 , the dihedral group of order 8.
To see this, it's enough to check that the group presentation: R 4 = P 2 = (RP ) 2 = 1. The remaining reflections can be identified with RP, R 2 P and R 3 P. For a given unitary symmetry we can compose it with K to obtain an anti-unitary one.
In our case, when N SYK → ∞, although multiple Z 2 symmetries are spontaneously broken, we only expect a two-fold ground state degeneracy. In fact, any two broken symmetries that can be related by an unbroken symmetry do not produce any extra ground state degeneracy.
To see this, consider for example the reflection symmetry RP. Since R is unbroken, we may assume R |0 = |0 without losing of generality. Then RP |0 = RP R |0 = P |0 .
At finite N SYK , however, certain discrete symmetry can be anomalous and is responsible for an exact two fold degeneracy for certain N SYK . For example, the particle-hole symmetry P ∼ KP acts on the fermions as The fermion number operator (3.7) is odd under this symmetry: When N SYK is not divisble by 4, there are two degenerate ground states |0 ± , and the symmetry generator P maps them into each other [35][36][37][38][39][40]: In this case we can say that the particle-hole symmetry is anomalous.