Quantum Entanglement of the Sachdev-Ye-Kitaev Models

Quantum Entanglement of the Sachdev-Ye-Kitaev Models

Chunxiao Liu Department of Physics, University of California, Santa Barbara, CA, 93106-9530    Xiao Chen Kavli Institute for Theoretical Physics CA 93106-4030, USA    Leon Balents Kavli Institute for Theoretical Physics, University of California, Santa Barbara, CA 93106-4030, USA
July 17, 2019
Abstract

The Sachdev-Ye-Kitaev (SYK) model is a quantum mechanical model of fermions interacting with -body random couplings. For , it describes free particles, and is non-chaotic in the many-body sense, while for it is strongly interacting and exhibits many-body chaos. In this work we study the entanglement entropy (EE) of the SYK models, for a bipartition of real or complex fermions into subsystems containing real/ complex fermions and / fermions in the remainder. For the free model SYK, we obtain an analytic expression for the EE, derived from the -Jacobi random matrix ensemble. Furthermore, we use the replica trick and path integral formalism to show that the EE is maximal for when one subsystem is small, i.e. , for arbitrary . We also demonstrate that the EE for the SYK4 model is noticeably smaller than the Page value when the two subsystems are comparable in size, i.e. is . Finally, we explore the EE for a model with both SYK2 and SYK4 interaction and find a crossover from SYK2 (low temperature) to SYK4 (high temperature) behavior as we vary energy.

I Introduction

The Sachdev-Ye-Kitaev (SYK) models Sachdev and Ye (1993); Sachdev (2015); Kitaev (2015); Maldacena and Stanford (2016); Polchinski and Rosenhaus (2016) are large- solvable zero-dimensional systems of Majorana or complex fermions coupled through all-to-all -fermion interactions. They have drawn considerable interest in high energy physics for their emergent conformal symmetry and holographic duality, and excite the condensed matter community as solvable examples of strongly interacting quantum critical non-Fermi liquids, which lack any quasiparticle description Banerjee and Altman (2017). This is reflected in the solution of the Green’s function, the scaling dimension of Majorana fields, and the non-zero entropy at zero temperature Maldacena and Stanford (2016). At , the model is also maximally chaotic which can be seen from its quantum Lyapunov exponent defined from the out-of-time-order correlator Larkin and Ovchinnikov (1969); Shenker and Stanford (2015); Kitaev (2015); Maldacena and Stanford (2016); Polchinski and Rosenhaus (2016). Starting from the original SYK4 model, many extensions have been proposed and studied which retain the interacting and chaotic nature Banerjee and Altman (2017); Bi et al. (2017); You et al. (2017); Gu et al. (2017a); Jian et al. (2017); Gu et al. (2017b); Song et al. (2017); Chen et al. (2017); Jian and Yao (2017); García-García et al. (2017). On the other hand, the case is a free fermion model that can be completely solved by random matrix theory. It realizes a zero dimensional disordered fermi liquid in the large limit, and both many-body chaos and zero temperature entropy are wholly absent. Given these facts, it is natural to study the phase transition between these two regimes Banerjee and Altman (2017); Bi et al. (2017); Song et al. (2017).

In this paper we study the entanglement entropy (EE) of the SYK model and discuss the differences between the SYK2 and SYK models with from the entanglement perspective. EE, a concept borrowed from quantum information theory, can quantify the quantum entanglement of different subsystems. The von Neumann EE for subsystem is concisely defined as , where is the reduced density matrix of . For the SYK2 model, we show that the EE can be analytically computed in closed form by connecting to the concept of Jacobi ensemble in random matrix theory Forrester (2010). For SYK with , since the model is maximally chaotic in the large limit, one would naturally presume that it is also maximally entangled Fu and Sachdev (2016). We numerically calculate EE and find that in the limit when is much smaller than the total system size, this is correct: the ground state EE for the SYK4 is maximally entangled. However, this is also true for SYK2, which is not many-body chaotic at all. We then show that these observations can be explained analytically by using the standard replica trick in the path integral formalism. Deviations from maximal entanglement are seen when the two subsystems obtained from bipartition are comparable in size: we find that the EE for the SYK4 model is much larger than that for the SYK2 model, yet both are still smaller than the Page value, which can be considered as the EE for the Gaussian unitary random matrix ensemble (GUE). Further extrapolation indicates that the difference between the EE of the SYK4 and Page value persists to the thermodynamic limit, showing that the SYK4 model is not maximally entangled.

We further study a model which includes both SYK2 and SYK4 interactions. At low temperature, the physics is dominated by the SYK2 term, which may be understood from the relevance, in the scaling sense, of the SYK2 interaction at the SYK4 fixed point.Banerjee and Altman (2017); Bi et al. (2017); Song et al. (2017) This is consistent with our ground state EE result which shows that EE decreases as we increase and approaches the SYK2 value in the thermodynamical limit. For a highly excited state, we find that the EE is the same as that for SYK4 model and Page value, indicating there is crossover from low energy SYK2 physics to high energy SYK4 physics.

The rest of the paper is organized as follows: In Sec. II, we briefly review SYK model and discuss the methods to compute EE. In Sec. III, we analytically compute the EE of SYK2 model and analyze the scaling behavior of EE in various limits. In Sec. IV, we numerically study EE of SYK2+SYK4 model and explore the crossover behavior from SYK2 to SYK4 physics. We further discuss the replica trick method in Sec. IV.1. We conclude in Sec. V with some final remarks.

Ii Models and Methods

The Hamiltonians of the SYK models at and , both complex and Majorana fermion versions, are written as

 H2,c = N∑i,j=1Tijc†icj−μN∑i=1c†ici, (1a) H2,χ = N∑i,j=1iJijχiχj, (1b) H4,c = ∑{i,j,k,l}Tij;klc†ic†jckck−μN∑i=1c†ici, (1c) H4,χ = ∑1≤i

The matrix is hermitian and is antisymmetric, with all Gaussian entries with zero mean and variance , where is some constant; the tensor is complex with appropriate symmetry to ensure hermiticity and the tensor is real, both with Gaussian entries of zero mean and variance Maldacena and Stanford (2016); Fu and Sachdev (2016). For the complex SYK model, a chemical potential is included to tune the filling factor.

We study the EE for the models in Eq. (1). As SYK are zero-dimensional models, the bipartition is made by choosing subsystem to consist of random Majorana or random complex fermions, with the complement called . Since SYK includes all to all coupling between fermions, we expect that the disorder averaged EE for subsystem A has the form , describing a “volume law” for EE. There is no regime of “area law” due to the complete non-locality of interactions. The main purpose of this paper is to investigate the coefficient and the possible subleading correction . Below we will resort to both analytical methods, which involve random matrix theory and the path integral formalism, and numerical means, mainly exact diagonalization, to study the EE.

Iii Analytical results for SYK2

The complex fermion SYK2 model (Eq. (1a)) is a free Hamiltonian in which is a hermitian random matrix belonging to the Gaussian unitary ensemble (GUE) Forrester (2010). The eigenvalue distribution of GUE satisfies the “semi-circle” law, and the diagonalizing matrix , composed of column eigenvectors, is uniformly distributed in the space of unitary matrices according to the Haar measure.

For a free fermion model, the EE of a subsystem can be directly calculated from the reduced two point correlation function matrix with each entry for Vidal et al. (2003); Peschel (2003). A key observation is that belongs to the Jacobi random matrix ensemble at Forrester (2010) (see Appendix A), where is the -by- upper-left block of , and , are the filling number of energy bands and subsystem size of , respectively. The distribution of the eigenvalue of in the large limit can be derived from the joint probability distribution of the Jacobi ensemble and satisfies the Wachter law Wachter (1980); Forrester (2010),

 f(x,κ,λ) =12πλ√(λ+−x)(x−λ−)x(1−x)1[λ−,λ+] (2) +(1−κ/λ)Θ(λ−κ)δ(x),λ,κ∈(0,1/2],

where , and . The function exists when is singular. The continuum part of , shown in Fig. 1c), is nonzero only on the support . The Wachter law, Eq. (2), does not depend on the value of , and holds also for orthogonal () and quaternion matrices () – see Appendix A. The EE is therefore

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯SA(κ,λ)=m∫[−xlnx−(1−x)ln(1−x)]f(x,κ,λ)dx. (3)

Due to particle-hole symmetry, the EE of the Majorana SYK2 model has the same scaling behavior as that of the complex fermion SYK2 model at (Appendix B).

We verify the theoretical result for the EE in Eq. (3) by comparing it with a numerical calculation in a large system with [Fig. 1a)]. In fact the function has several limiting forms for which EE can be computed analytically. For instance, when the subsystem is much smaller than the total system, approaches the distribution [Fig. 1c) right] and the EE is

 ¯¯¯¯¯¯¯SAλ→0−−→m(−κlnκ−(1−κ)ln(1−κ)). (4)

This result has also been obtained in Ref. Magán, 2016. In particular, when , i.e. complex fermion at half-filling or Majorana fermion, the EE is simply equal to . This result agrees with simply counting the degrees of freedom of and exhibits maximal entanglement between and . However when is comparable to in size, disperses in the whole interval [Fig. 1c) left] and the entanglement deviates from the maximal value; in particular, when we have . This result is consistent with the bounds given for the quadratic fermionic Hamiltonian Vidmar et al. (2017).

Eq. (3) features a symmetry apart from the trivial symmetries and . Therefore in the zero filling factor limit, by analogy with Eq. (4). The quantity entanglement entropy density, defined as , measures the average EE for each fermion and is plotted in Fig. 1b). Clearly the maximal value is reached only when , and the EE density is a descending function of both and .

Iv Numerical and analytical results for SYK4+SYK2

The complex SYK4 model at half filling has been investigated in a small system in Ref. Fu and Sachdev, 2016, where the authors reported that this model exhibits maximal entanglement. Here via sparse matrix diagonalization we numerically calculate the ground state EE of both the complex fermion SYK4 for up to and the Majorana fermion versions for up to , where for the former we choose two filling factors. The results are shown in Fig. 2. We see that when , the EE density is equal to . This is the same as Eq. (4) for SYK2 and in fact is true for all SYK models (see below). As we increase , the EE for SYK4 is much larger than that for SYK2. However, from the results of , it is hard to predict the analytical form in the thermodynamical limit.

To better understand the scaling of EE around in the thermodynamical limit, we conduct finite size scaling analysis (see inset of Fig. 2) and compare EE of both Majorana and complex fermion SYK4 models at with Page value Page (1993)

 ¯¯¯¯¯¯¯¯¯¯SA,P=mln2−|HA|2|HB|, (5)

where and are the dimensions of the Hilbert space of and respectively. The Page value is the EE for a random pure state (eigenstate of GUE) and is equal to maximal EE up to a small constant, which is when . Fig. 2 shows that the Page value (red line) and SYK4 (green solid and dashed lines) are the same when but differ around . One naturally asks how much this lowering in SYK4 is due to finite size effects. To address this question, we further check the EE density of SYK4 at for different , denoted , and find excellent linearity between and (see the inset of Fig. 2). This suggests that at , with for Majorana (complex) fermions obtained from extrapolation. is smaller than , suggesting the deviation from Page state in the thermodynamical limit. On the other hand, and is larger than in the Page value. Both these results demonstrate that the SYK4 ground state is less entangled than Page state around . However, the difference should disappear as we increase since the Hamiltonian becomes less sparse and approaches a member of the GUE.

We also remark that the thermal entropy at zero temperature is smaller than the EE in the SYK4 model. It is known that Majorana SYK model (with ) has zero temperature residual entropy in the large N limit with this particular form , where is some function of Kitaev (2015); Maldacena and Stanford (2016). When , this value is around and is much smaller than the ground state EE density. This difference disappears as we increase .

We further numerically calculate the ground state EE of for different couplings , see Fig. 2. It is known that the SYK2 interaction is a relevant perturbation which induces a flow from SYK4 (a non-fermi liquid) behavior to a free fixed point Banerjee and Altman (2017); Bi et al. (2017); García-García et al. (2017). This physics is also reflected in the scaling of EE at finite , in two aspects: (i) for a fixed , the EE is lowered as is increased, which can be clearly observed in Fig. 2, (ii) for a fixed , as we increase , EE for the same subsystem size remains the same. This is consistent with the scaling dimension of fermion operator equal to .

On the other hand, as we increase the energy and move to a highly excited state, there is a “crossover” from SYK2 physics to SYK4 physics and we expect that the EE for the highly excited state should be the same as that for the pure SYK4 interaction.Song et al. (2017) To verify this statement, we first study the excited state EE for the pure SYK4 and the SYK2 models and show that they are different (see Fig. 3). We numerically find that the EE of a highly excited state for the SYK4 model is the same as Page value. In contrast, for the SYK2 model, both excited and ground state EE exhibit the same scaling since belongs to the Jacobi ensemble. As we introduce SYK4 interactions to the SYK2 model, the EE for the highly excited state drastically changes and becomes the Page value, the same as that for the pure SYK4 model.

iv.1 Replica trick method

The local observables, correlation functions, and thermodynamics of the SYK models are exactly obtained in the large limit using the path integral formalismKitaev (2015); Maldacena and Stanford (2016). It is natural to ask whether the EE of the SYK models may also be found exactly. Using the standard replica method and path integral formalism, one may write

 ¯¯¯¯¯¯¯SA=limn→1¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯lnTrA[ρnA]1−n=limn→1¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ln(Zn/Zn)1−n, (6)

where is the partition function involving copies of coupled with A as shown in Fig. 4b). The limit above is used to obtain the logarithm in the definition of the von Neumann entropy. Notice that for a disordered system, to compute we need in principle to perform a second replica trick. Luckily, the SYK in the large N limit is always in a replica symmetric phase and therefore we have .

After disorder average and a saddle point approximation valid at large , we have . The zeroth order physics is determined entirely by the saddle point equation (see Appendix C for the derivation):

 J2[λGab(τ,τ′)+(1−λ)gab(τ,τ′)]3 = Σab(τ,τ′), (7a) 1∂Aτδ(τ−τ′)δab−Σab(τ,τ′) = Gab(τ,τ′), (7b) 1∂Bτδ(τ−τ′)δab−Σab(τ,τ′) = gab(τ,τ′), (7c)

where are the replica indices, and are the Green’s functions, respectively, for the and subsystem. and have different boundary conditions

 Gab(β,τ′)=G(a+1)b(0−,τ′), Gnb(β,τ′)=−G0b(0−,τ′), (8) Gab(τ,β)=Ga(b+1)(τ,0−), Gan(τ,β)=−Ga0(τ,0−), gab(β,τ′)=−gab(0−,τ′), gab(τ,β)=−gab(τ,0−),

where , for the first line, , for the second line and for the third line. The differential operators are different for and which ensure the continuity of the equations (see Appendix C for explanation). These boundary conditions are illustrated in Fig. 4. Due to the diffferent boundary conditions for and , obtaining a general and explicit form for and turns out to be hard (Similar issue arises in entanglement dynamics in a SYK chain model in Ref. Gu et al., 2017c). However, the EE in the scaling limit does not require solving Eqs. (7) and can be readily obtained by counting degrees of freedom. In this limit, only the diagonal term in is nonzero and is the same as that in . Moreover, to first order of . As a result, all terms in cancel those in , except for the term that counts the degrees of freedom of subsystem . In , there is only one subsystem with time antiperiodicity , while there are copies in . Therefore we have

 (9)

This result is independent of and leads to Eq. (4). This shows the super-universality of maximal entanglement for small subsystems for all .

For around , solving Eqs. (7) is difficult. Notably, one must impose anti-periodicity in time on the subsystem, which couples to all -replicas. This coupling introduces an effective interaction between different -replicas, which leads to nonvanishing even for . In principle Eqs. (7) can be solved numerically but appropriate methods have to be developed and analytic continuation to needs to be understood. We leave this to future work.

V Discussion and Outlook

In this paper, we studied the entanglement aspects of the Sachdev-Ye-Kitaev (SYK) model at and . For the free model SYK, we showed that the reduced correlation matrix belongs to the -Jacobi ensemble at , and obtained from this an analytic expression for the EE of SYK. We further demonstrated that the ground state EE of the SYK model shows a derivation from the Page value which persists to the thermodynamic limit, indicating that SYK is not maximally entangled. Furthermore, we explored the EE for the SYK+SYK model, and found that there are two distinct regimes as we vary energy: for low energy states, the EE is dominated by the SYK2 term, while for highly excited states, the physics is dominated by the SYK4 term and the EE is the same as the Page value.

We were unable to obtain a full analytical expression for the EE of the SYK4 model for an arbitrary bipartition. This calculation requires coupling of different replicas and it would be very interesting to develop a field theoretic approach to this problem in the large limit, even if partly numerical.

Note added. Shortly after the submission of this draft to arXiv, another paper Huang and Gu (2017) gave the analytical estimation of the EE of the SYK4 model at subsystem/system=1/2, which is consistent with our numerical result at this ratio.

Acknowledgements.
We thank A. Ludwig, Z. Bi, L. Vidmarand, V. Rosenhaus and P. Lu for useful discussions. X.C. was supported by a postdoctoral fellowship from the the Gordon and Betty Moore Foundation, under the EPiQS initiative, Grant GBMF4304, at the Kavli Institute for Theoretical Physics. We acknowledge support from the Center for Scientific Computing from the CNSI, MRL: an NSF MRSEC (DMR-1121053). L.B. and C.L. were supported by the NSF Materials Theory program, grant number DMR1506199.

Appendix A β−Jacobi Ensemble and its eigenvalue scaling laws

Suppose matrices and are of size and , respectively, both with independently and identically distributed Gaussian entries over either the real, complex or the quaternion field ( respectively). Then by -Jacobi (or MANOVA) ensembles we mean the eigenvalue distribution of

 A†AA†A+B†B. (10)

The joint eigenvalue distribution is

 f(ε1,⋯,εm) =C(m1,m2,m)∏1≤i

In the limit while keeping the ratios

 λ=mm1+m2,κ=m1m1+m2

finite, the eigenvalue distribution satisfies the Wachter law [Eq. (2) in the main text].

The EE of the complex SYK2 model is calculated via the reduced correlation function . The eigenvalues of belong to the -Jacobi ensemble, which in the large limit satisfy the Wachter law. This is because the single-particle Hamiltonian of complex SYK2 is a GUE, so the eigenvector matrix is uniformly distributed in the space according to the Haar measure. The truncated correlation matrix therefore satisfies Eq. (10), provided we set , , and , where

 U=(Vn×mV′n×(N−m)W(N−n)×mW′(N−n)×(N−m))T.

The derivation from Eq. (10) to Eq. (11) is achieved by diagonalizing (10), which induces a basis transformation in the space of matrices; the “level repulsion” term comes from the Jacobian of the measure during transformation. Below we sketch the derivation from Eq. (11) to Eq. (2) in the main text, known as the Coulomb gas approach, and for simplicity we only discuss . For a complete proof, see Ref. Ramli et al. (2012).

The Coulomb gas approach expresses the joint distribution as a path integral. In the limit (keeping and finite) we have

 f(ε1,⋯,εm)→1Ze−NS[ε],

where , and

 S=−∬dxdyln|ε(x)−ε(y)|−∫dx(alnε+bln(1−ε)),

where , . Due to the presence of the large in front of the action, the leading order behavior is entirely determined by the saddle point equation

 ∫dy|ε(x)−ε(y)|+aε(x)−b1−ε(x)=0.

Under appropriate boundary conditions, this equation has unique physical solution, Eq. (2) in the main text.

The Wachter distribution for some generic values of are shown in Fig. 5.

Appendix B Calculating Entanglement Entropy of Majorana SYK2 Using Correlation Functions

The Majorana SYK2 model is

 H2,χ=iN∑i,j=1Jijχiχj=χTJχ, (12)

where are Majorana fermion operators with , and . Here is even and we set .

To ensure hermiticity, is antisymmetric with dimension . can then be “diagonalized” to the block form below by a orthogonal matrix :

 J=OΣOT,

where

 Σ=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝0⋱00λ1−λ10⋱0λr−λr0⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, (13)

where there are zeros on the diagonal block, and

 0<λ1≤λ2≤⋯≤λr.

Now let us denote , where (, ) denote the rotated Majorana operators which satisfy and . We have

 H=iχTJχ=iχTOΣOTχ=iγTΣγ=r∑i=12iλiγ1iγ2i. (14)

We further denote and , , where ’s are complex fermion operators with standard anticommutation relations, and

 γ1jγ2j=1i(fj+f†j)(fj−f†j)=1i(2f†jfj−1), (15)

then is diagonalized by ’s

 H2,χ=r∑i=12λi(2f†ifi−1), (16)

where states labeled by have finite excitation energy , and have zero excitation energy.

For the rest we will assume that , i.e. the Hamiltonian has no zero modes. This is justified for the SYK2 model since a random matrix has almost zero possibility in producing zero eigenvalues. The ground state of SYK2 is then , where is the state that is annihilated by all : for . The ground state expectation value is therefore , and

 ⟨γγT⟩=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝1i−i1⋱1i−i1⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠≡I+iK,

where we have defined

 K=−IM×M⊗(01−10).

Then, the correlation function matrix becomes

 C=⟨χχT⟩=⟨OγγTOT⟩=O⟨γγT⟩OT=I+iOKOT,

and the reduced correlation function matrix is (for subsystem with size )

 CA=I2m×2m+iO2m×2MK2M×2MOT2m×2M.

After diagonalizing we get , where

and we have Peschel (2003)

 SA=−m∑i=11−εi2ln1−εi2+1+εi2ln1+εi2. (17)

To prove Eq. (17), we assume that the reduced density matrix can be written in terms of an effective Hamiltonian , where is a matrix, by

 ρA=1Ze−HA,Z=TrA[e−HA].

Suppose is diagonalized in some Majorana basis and some complex basis (following the derivation for diagonalizing )

 HA=χTJAχ=m∑i=12ihiξ1iξ2i=m∑i=12hi(2ψ†iψ−1),

then , as a matrix that acts on the many-body Hilbert space of , has eigenvalues

 ˜εi=e−2hie2hi+e−2hi,i=1,2,...,m,

and finally, by definition, the EE is

 SA =−TrρAlnρA=−∑D∈D∏α∈D∏β∈A∖D˜εα(1−˜εβ)ln˜εα(1−˜εβ) (18) =−m∑i=1˜εiln˜εi−m∑i=1(1−˜εi)ln(1−˜εi).

To connect with the truncated correlation function matrix , note that is a block diagonal matrix, and

 ⟨ψ†jψj⟩ =∑{n}mi=1=0,1⟨{n}mi=1|e−m∑i=12hi(2ψ†iψi−1)ψ†jψ†j|{n}mi=1⟩∑{n}mi=1=0,1⟨{n}mi=1|e−m∑i=12hi(2ψ†iψi−1)|{n}mi=1⟩ (19) =∑nj=0,1⟨nj|e−2hj(2ψ†jψj−1)ψ†jψj|nj⟩e2hi+e−2hi =˜εi.

 CA=I−iY†AEAYA,

where

 EA=Diag(2˜ε1−1,⋯2˜εm−1)⊗(01−10),

thus the eigenvalues of are , . Plug these into Eq. (18) we get Eq. (17).

Appendix C Path integral formalism for the Entanglement Entropy of SYK4

The path integral formalism helps find the saddle point equations for the SYK models at leading order of . In this Appendix, we sketch the derivation for SYK4 but the generalization to all is straightforward. We introduce replicas to SYK, with index . Notice that

 Zn({Jijkl}) =TrAρnA (20) =TrA(TrBρ)n =∫dψA⟨ψA|(∫[dψaB]na=1⟨ψ1B|ρ|ψ1B⟩⋯⟨ψnB|ρ|ψnB⟩)|ψA⟩ =∫[Dχai]a,iexp⎡⎣∫β0dτn∑a=1⎛⎝−12χai∂aτχai−∑1≤i

where we use , , to denote a complete basis for the Hilbert space of subsystem , and a complete basis for the Hilbert space of subsystem . In the last line become real Grassmann fields with a continuous index and a discrete replica index . Because of the different trace rules in and , the boundary conditions for are

 χai(β)=χa+1i(0−) and χni(β)=−χ1i(0−)for i∈A,χbi(β)=−χbi(0−)for i∈B,

where and . The operators for and , accordingly, also differ, respecting the boundary condtions for and , respectively, which ensures they are continuous operators. To remind these differences, we also label with replica index , where () are the same for all ().

After averaging over we have

 ¯¯¯¯Zn =∫⎛⎜⎝∏1≤i

in which the fields for different sites are decoupled. Suppose has particles and has particles. We introduce Green’s functions and and self-energies and such that

 1 =∫[DGabDΣab]a,bexp[−12∫β0dτdτ′n∑a,b=1Σab(τ,τ′)(PGab(τ,τ′)−∑i∈Aχai(τ)χbi(τ′))], (22) 1 =∫[DgabDσab]a,bexp[−12∫β0dτdτ′n∑a,b=1σab(τ,τ′)((N−P)gab(τ,τ′)−∑i∈Bχai(τ)χbi(τ′))],

plug them into Eq. (21) and then conduct the integrals over the Majorana fields , we get

 ¯¯¯¯Zn =∫[DGabDΣabDgabDσab]a,bexp[PlnPf(∂Aτδ(τ−τ′)δab−Σab(τ,τ′))+(N−P)lnPf(∂Bτδ(τ−τ′)δab−σab(τ,τ′))] (23) ×exp[∫β0dτdτ′n∑a,b=1(−P2Σab(τ,τ′)Gab(τ,τ′)−N−P2σab(τ,τ′)gab(τ,τ′)+J28N3[PGab(τ,τ′)+(N−P)gab(τ,τ′)]4)],

which gives back Pfaffians of matrices and , with index pairs as rows and as columns. Also note that and respect the boundary conditions of the in and , which gives Eq. (8) in the main text. Now define , thus , and we have a large parameter that controls the behavior of . In the thermodynamic limit , we have

 ¯¯¯¯Zn=exp(−NScl), (24)

i.e. is dominated entirely by the classical action . is obtained by plugging in the solution of saddle point equations for , , and . The saddle point equations for , , and in are

 λΣab(τ,τ′)−λJ2[λGab(τ,τ′)+(1−λ)gab(τ,τ′)]3 = 0, (25a) (1−λ)σ