We analyze a two-body non-Hermitian two-site Sachdev–Ye–Kitaev (SYK) model with the couplings of one site complex conjugated to the other site. This model, with no explicit coupling between the sites, shows an infinite number of second-order phase transitions, which is a consequence of the factorization of the partition function into a product over Matsubara frequencies. We calculate the quenched free energy in two different ways: first in terms of the single-particle energies and second by solving the Schwinger–Dyson equations of the two-site model. The first calculation can be done entirely in terms of a one-site model. The conjugate replica enters due to non-analyticities when Matsubara frequencies enter the spectral support of the coupling matrix. The second calculation is based on the replica trick of the two-site partition function. Both methods give the same result. The free-fermion partition function can also be rephrased as a matrix model for the coupling matrix. Up to minor details, this model is the random matrix model that describes the chiral phase transition of QCD, and the order parameter of the two-body model corresponds to the chiral condensate of QCD. Comparing to the corresponding four-body model, we are able to determine which features of the free energy are due to the chaotic nature of the four-body model. The high-temperature phase of both models is entropy dominated, and in both cases, the free energy is determined by the spectral density. The chaotic four-body SYK model has a low-temperature phase whose free energy is almost temperature-independent, signaling an effective gap of the theory even though the actual spectrum does not exhibit a gap. On the other hand, the low-temperature free energy of the two-body SYK model is not flat; in fact, it oscillates to arbitrarily low temperature. This indicates a less desirable feature that the entropy of the two-body model is not always positive in the low-temperature phase, which most likely is a consequence of the non-hermiticity.

## I. INTRODUCTION

In 1971, George Uhlenbeck asked Freeman Dyson “Which nucleus has levels distributed according to the semi-circle law?”^{1} As an answer to his criticism, Dyson published his last paper on Random Matrix Theory (RMT),^{2} in which he introduced a Brownian motion process to construct an ensemble of random matrices with an arbitrary level density, in particular the nuclear level density given by the Bethe formula,^{3}

Uhlenbeck could have asked a different question: “Is the nuclear force an all-to-all many-body interaction?,” which is the case for Wigner–Dyson ensembles. This question led French, Wong, Bohigas, Flores, and Mon^{4–8} around the same time to the introduction of the two-body [which in the Sachdev–Ye–Kitaev (SYK) literature is known as a four-body interaction] random ensemble, which reflects the two-body nature of the nuclear interaction. It took until 2015, through the seminal work of Kitaev,^{9,10} to realize that Uhlenbeck’s criticism could also have been addressed by this work. The two-body random ensemble, in particular the version of the model introduced by Mon and French,^{8} is now known as the complex Sachdev–Ye–Kitaev (SYK) model.^{11,12}

Since the pioneering work of Wigner, Dyson, Gaudin, and Mehta,^{2,13–20} random matrix theory has been applied to virtually all areas of physics and even outside of physics; see the comprehensive review by Guhr, Müller-Groeling, and Weidenmüller.^{21} In this paper, we study non-Hermitian random matrix theories, which were first introduced by Ginibre,^{22} but have also been applied to many areas of physics, for example, the distribution of poles of S-matrices,^{23–25} the Hatano–Nelson model,^{26–28} dissipative quantum systems,^{29–32} QCD at nonzero chemical potential,^{33–41} and Parity-Time (PT)-symmetric systems,^{42} to mention a few. Non-Hermitian random matrices were classified^{43–47} along the lines of the classification of Hermitian random matrices.^{2,20,48,49} A recent review of non-Hermitian physics was given by Ashida, Gong, and Ueda.^{50}

The possibility of a non-Hermitian version of the Sachdev–Ye–Kitaev (SYK) model was originally suggested by Maldacena and Qi^{51,52} as a two-site SYK model for Euclidean wormholes without explicit coupling between the two SYK models (the only coupling is through the randomness) and was studied in detail in subsequent papers.^{47,53–55} In this model, the Left (L) and Right (R) partition functions are complex conjugate to each other, and each of them has *N*/2 Majorana fermions so that energy levels of the *q*-body Hamiltonian

are given by ${Ek+El*}$ if the {*E*_{k}} are the eigenvalues of *H*_{L}. Therefore, the partition function factorizes as

and is necessarily positive definite.

One of the main conclusions from these studies is that the chaotic model has a first-order phase transition, which separates the low-temperature phase from the high-temperature phase. In the high-temperature phase, the average of the partition function factorizes

and the free energy follows from the eigenvalue density of the one-site Hamiltonian. Because of the complex phases of the eigenvalues, the average partition function is exponentially suppressed due to cancellations. In the case of maximum non-hermiticity when the eigenvalue density is isotropic in the complex plane, the partition function ⟨*Z*_{L}⟩ becomes temperature independent and the free energy is determined by the total number of states $F=\u2212TN2log\u20612$. It is clear that this result cannot be correct at low temperature when *F* → −|*E*_{0}|, with *E*_{0} being the ground state energy. The only possibility is that the low-temperature phase receives contributions from the correlations of the eigenvalues of the *L* and *R* Hamiltonians. Indeed, for the *q* = 4 SYK model, the free energy in the low-temperature phase is entirely determined by two-point correlations of the eigenvalues.^{54,55} The reason that this can happen is the exponential suppression of the single site partition function due to the complex phase of the eigenvalues. The dynamics of the *q* = 4 SYK model is chaotic with eigenvalue correlations in the universality class of the Ginibre model. The universal two-point correlations give rise to a temperature-independent free energy at low temperatures. This also explains that results obtained by solving Schwinger–Dyson equations are very close to the results for the Ginibre random matrix ensemble. We conclude that the quantum chaotic nature of the model is responsible for a nearly temperature-independent^{89} free energy in the low-temperature phase when an actual spectral gap is absent. As a consequence, the free energy of the low-temperature phase of the *q* = 2 SYK model, which is integrable and has no spectral gap, has to be different. The goal of this paper is to solve the *q* = 2 non-Hermitian SYK model to study the effects of chaos and integrability on the phase diagram.

As was the case for the *q* = 4 SYK model,^{54,55} in this paper, we evaluate the *quenched* free energy in two structurally different ways: First from the eigenvalues of the SYK Hamiltonian giving the quenched free energy and second from the solution of Schwinger–Dyson equations^{56,57} of the SYK model in the Σ*G* formulation giving the *annealed* free energy of the two-site SYK model. For *q* = 2, it is possible to perform the spectral calculation both analytically and numerically, while for *q* = 4, this could only be done numerically by an explicit diagonalization of the SYK Hamiltonian. In addition, the Schwinger–Dyson equation can be solved analytically for *q* = 2,^{58,59} while for *q* = 4, we had to rely on numerical techniques.

The Σ*G* formulation of the SYK model is based on the replica trick for the quenched free energy,^{60}

which is known to fail,^{61–63} in particular, for non-Hermitian theories.^{64} However, by now, it has been well understood how to refine the replica method so that its results can be trusted.^{33,65–71} For Hermitian models, the naive replica trick usually gives the correct result for a mean field analysis. This is also the case for the Σ*G* formulation of the SYK model.^{72,73} However, as should be clear from the arguments given above, the naive application of the replica trick to the one-site non-Hermitian SYK model gives an incorrect result for the quenched free energy of the low-temperature phase. Instead, the quenched free energy of the one-site Hamiltonian is given by the replica limit of the product of the partition function and its complex conjugate. In the case of the one-site partition function, the replica symmetry is broken in the sense that the conjugate replica emerges due to quenching and couples to the original replica. The emergence of conjugate replicas for quenched replicas is well-known in non-Hermitian RMT^{27,66,74} and QCD at nonzero chemical potential.^{33,34,70,75}

We start this paper with a short introduction of the SYK model and the replica trick. Since the *q* = 2 SYK model is a Fermi liquid, we can calculate the quenched free energy from a free-fermion formulation of the model; see Sec. III. In Sec. IV, we calculate annealed free energy for the Σ*G* formulation of the SYK model for one replica and one conjugate replica. The final result is in complete agreement with the result from the free-fermion calculation. The Σ*G* action of the *q* = 2 SYK model is quadratic in *G* so that it can be integrated out exactly. The result resembles a random matrix *σ*-model. In Sec. V, we show that this *σ*-model can also be directly obtained from the free-fermion description of the *q* = 2 SYK model. Concluding remarks are made in Sec. VI, and some technical details are worked out in Appendixes A and B. In Appendix A, we show that the occupation number representation also applies to the non-Hermitian SYK model, and in Appendix B, we work out the free-fermion calculation for the case that the non-hermiticity is not maximal. Some preliminary results of this paper were summarized in a recent letter.^{54}

## II. THE *q* = 2 NON-HERMITIAN SYK MODEL

The Hamiltonian of the *q* = 2 non-Hermitian SYK model is given by

where *ψ* are Majorana fermions and *J*, *M* are random couplings. A general *q*-body Hamiltonian would have products of *q* Majorana fermions. The tensor product structure on the second line of Eq. (6) follows from an explicit Dirac-matrix representation of *ψ*,

where *γ*_{i} are the Dirac matrices in *N*/2 dimensions and *γ*_{c} is the corresponding chirality matrix (assuming *N*/2 is even). The variances of the couplings are given by

where *v* is a dimensionful parameter that sets the physical scale. In a representation where the left gamma matrices are real and the right gamma matrices are purely imaginary, the Hamiltonian of a single SYK is anti-symmetric under transposition. The eigenvalues, which are representation independent, thus, occur in pairs ±*λ*_{k}. In this representation, we also have that $HL=HR*$ (with the minus sign from the sum included in *H*_{R}). The spectrum of Hamiltonian (6) is given by $\xb1\lambda k\xb1\lambda l*$ if *λ*_{k} are the eigenvalues of *H*_{L} with a positive real part. The partition function of this Hamiltonian is necessarily positive,

where *Z*_{L(R)} is the partition function of the *L*(*R*) Hamiltonian. The average partition function will be denoted as ⟨*Z*⟩ with the appropriate subscripts. Contrary to *q* > 2, the SYK Hamiltonian for *q* = 2 is a Fermi liquid with single-particle energies ±*ɛ*_{k} given by the eigenvalues of the coupling matrix, which is an anti-symmetric non-Hermitian random matrix. Therefore, the quenched partition function is given by

and

For the last equality, we have used that the average density of the single-particle energies satisfies $\u27e8\rho Lsp(z)\u27e9=\u27e8\rho Rsp(z*)\u27e9$. Therefore, the quenched free energy can be obtained from the one-site partition function, which is a direct consequence of the Fermi-liquid nature of the *q* = 2 SYK model.

When we evaluate the quenched free energy of the SYK partition function in the Σ*G* formulation, we will employ the replica trick,

For large *N*, the partition function can be evaluated by a saddle point approximation. As was argued in Refs. 54 and 55, since the partition function $Z=ZLZL*$ is positive definite, the replica trick is expected to give the correct mean field result with unbroken replica symmetry,

Therefore, the quenched free energy is equal to the annealed free energy,

and it can be calculated by evaluating the partition function for one replica of the two-site Hamiltonian and, in other words, one replica and one conjugate replica in terms of the one-site Hamiltonian. If the (left-right) replica symmetry is broken, we have that

so that it is not guaranteed that the quenched free energy can be obtained from the one-site annealed free energy. We will see in Sec. IV that is the case for the low-temperature phase of the Σ*G* formulation of the SYK model. In the literature, the coupling between replicas has been related to the formation of wormholes between black holes.^{76,77}

Before calculating the annealed free energy from the Σ*G* formulation of the SYK model, in Sec. III, we will evaluate the average partition function using the properties of the spectra of anti-symmetric non-Hermitian random matrices.

## III. THE FREE ENERGY OF THE TWO-SITE NON-HERMITIAN SYK FOR *q* = 2

Results for the *q* = 4 non-Hermitian SYK model^{54,55} suggest that quantum chaotic dynamics is responsible for a replica-symmetry-breaking (RSB) phase at low temperatures with a temperature-independent free energy. Indeed, the free energy of an integrable non-Hermitian model of random uncorrelated energies^{55} exhibits a distinct low-temperature phase with a temperature-dependent free energy. However, the random energy model lacks a natural interpretation as a many-body model. In this section, we test this hypothesis by evaluating quenched free energy of the two-site *q* = 2 non-Hermitian SYK model, which is a Fermi liquid with energies given by sums of single-particle energies. Therefore, the many-body eigenvalues obey Poisson statistics, which is consistent with a vanishing Lyapunov exponent^{78} (obtained by calculating the out of time order correlator). We will calculate the quenched free energy from the single-particle energy density, which is constant inside an ellipse in the complex plane. As argued in Sec. II, the quenched free energy can be obtained from the *one-site* partition function. In Sec. IV, we will see that the same result can be obtained from the SD equation of the *two-site* model.

After a change of basis (see Appendix A for a demonstration), the one-site Hamiltonian can be expressed as a free Fermi liquid with single-particle energies, which are the eigenvalues of the antisymmetric coupling matrix.^{90} To be concrete, the Hamiltonian *H*_{L} in (6) becomes

where *ɛ*_{k} are the eigenvalues of the coupling matrix *i*(*J*_{ij} + *ikM*_{ij})/2 with positive real parts and, hence, the sum runs up to *k* = *N*/4 instead of *N*/2. Just as in the Hermitian case,^{59} we have

with

Hence, we conclude that the many-body energies of *H*_{L} are given by filling *N*/4 free fermions into the single-particle states of (16). We note that generally *c*_{k} is not the Hermitian conjugate of $c\u0303k$, and hence, the energy eigenstates are not necessarily orthogonal to each other, which is consistent with the non-hermiticity of the Hamiltonian.

In this free-fermion representation, the quenched free energy of the one-site SYK model, which we shall see to be identical to half the two-site model, is simply

where ⟨*ρ*(*z*)⟩ is the averaged spectral density of the coupling matrix *i*(*J*_{ij} + *ikM*_{ij})/2. Note that since *ɛ*_{k} are only half of the levels of *i*(*J*_{ij} + *ikM*_{ij})/2 (those with positive real parts), the integral in the second equality should have only covered half of the support of *ρ*(*z*). However, since the integrand is invariant under *z* ↦ −*z*, we simply integrate over the whole support and compensate it by a pre-factor of 1/2.

At large *N*, the averaged spectral density ⟨*ρ*(*z*)⟩ is a constant inside the ellipse,^{29,79–81}

where

and *v* is the physical scale introduced in Eq. (8).

For *k* = 0, the eigenvalues are real with spectral density given by

The free energy was already calculated before^{10} and is given by

Using the Weierstrass product formula for cosh *x*, this can be expressed as

where *ω*_{n} are the Matsubara frequencies,

The integral over *x* is known analytically, resulting in

We recognize the expressions for Green’s function and the self-energy in the Σ*G* formulation,^{59} suggesting that this result can also be obtained from the solution of the SD equations; see Sec. IV.

Next, we discuss the case *k* = 1 where the ellipse becomes a circle with *ϵ*_{0} = *υ*_{0} (for the general elliptic case, see Appendix A). The free energy (19) can be expressed as

where $D\upsilon 0$ represents a disk of radius *υ*_{0} centered at the origin, and in the second equality, we have scaled the integral to the unit disk where

where *S*_{r} is a circle of radius *r*. We can directly evaluate *I*(*r*) by expressing the cosh as a product over Matsubara frequencies using the Weierstrass formula,

where

We note that the integrand of *I*_{n} has two cuts and no pole: the two cuts start at ±*iω*_{n}/2*υ*_{0} and extend horizontally to the negative infinity. If the circle *S*_{r} does not touch the cuts (*r* < *ω*_{n}/2*υ*_{0}), then *I*_{n} = 0; if *S*_{r} intersects with the cuts (*r* > *ω*_{n}/2*υ*_{0}), then we choose the contour to narrowly avoid the cuts (see Fig. 1) and then we have that the sum of *I*_{n} and the cut contributions vanishes so that the problem reduces to evaluating the contributions from the parts that surround the cuts. Hence, we conclude

Equation (29) finally truncates to

The free energy in (27) is, then, given by

In Sec. IV, we will show that this result can also be obtained from the solutions of Schwinger–Dyson equations.

We observe that a phase transition occurs at temperatures for which an additional Matsubara frequency enters in the sum of Eq. (33). This happens for $\upsilon 0\pi T\u221212=m$ (*m* being non-negative integers), resulting in a series of critical temperatures parameterized by *m*,

where $\upsilon 0=v/2$ according to Eq. (21). For *T* > *T*_{c,0}, there are no Matsubara frequency satisfying 0 < *ω*_{n} < 2*υ*_{0}, and no more phase transitions can occur.

We re-iterate that because of the complex-conjugation invariance of the *average* single-particle spectral density, we have

so that the two-site quenched free energy just doubles that of the one-site and the free energy density remains the same, that is, *F*/*N* = *F*_{L}/(*N*/2).

To determine the order of the phase transition, it is useful to study the first derivative of the free energy. Indeed, we observe, see Fig. 3, that it has kinks that point toward a family of second-order phase transitions, each time when a new pair of Matsubara frequencies enters in the sum. This can be shown analytically by expanding the free energy around *T*_{c,m}. We find that the contribution at this new critical temperature scales as $\u223c(T\u2212Tc,m)2$, so, as depicted in Fig. 2, the free energy (33) is smooth around this critical temperature, and therefore, the transition cannot be of first order. However, one can easily show that the derivative of the free energy is not smooth at *T*_{c,m} as is also clear from Fig. 3.

Using similar methods, the integrals can also be calculated for 0 < *k* < 1, see Appendix B,

Each term contributing to the sum in the first line is equal to the corresponding term in the free energy of the Hermitian SYK model (26) with $\u03f50\u2192\u03f502\u2212\upsilon 02$. In Figs. 2 and 3, we show the temperature dependence of the free energy and its derivative for various values of *k*. It is disturbing that the entropy becomes negative, which cannot be due to the failure of the replica trick since it is not used in the free-fermion method. Most likely, it is a consequence of the non-hermiticity, which will become more clear in Sec. V.

The critical temperatures are given by the same expression as for *k* = 1,

but with $\upsilon 0=k2v/1+k2$ according to Eq. (21). One can easily show that the derivative of the free energy is continuous at the critical points, while its second derivative is discontinuous. The appearance of an infinite number of critical points is a direct consequence of the factorization of the partition function into a product over positive Matsubara frequencies. For each Matsubara frequency, we have exactly one critical temperature.

## IV. FREE ENERGY FROM THE SCHWINGER–DYSON EQUATIONS

We now turn to confirm these results by an explicit large *N* calculation from the solutions of the Schwinger–Dyson (SD) equations of the SYK model. In this case, the calculation is also analytical though very different from the one carried out in Sec. III. However, we shall see that ultimately the expression for the free energy is the same. In the Schwinger–Dyson approach, replica symmetry breaking between conjugate replicas plays a crucial role. Indeed, we shall see that for *k* = 1, the free energy is determined by Green’s function *G*_{LR} (equivalently, its self-energy Σ_{LR}) related to the effective coupling of the two sites. However, it is assumed that the replica symmetry of a conjugate pair remains unbroken so that the quenched free energy can be obtained from just one replica and one conjugate replica.

The Euclidean Σ*G* action for the *q*-body SYK model takes the form^{51}

where the indices *a*, *b* can be equal to *R* or *L*. The integrations over *τ* variables are on the interval [0, *β*]. The factor *s*_{ab} is equal to 1 for *a* = *b* and equal to (−1)^{q/2} = −1 for *a* ≠ *b*. The couplings take the value $JLL=JRR=J$ when *a* = *b* and $JLR=JRL=J\u0303$ when *a* ≠ *b*. The term proportional to *ϵ* is included to break the symmetry that requires *G*_{LR} to vanish (more details are given in the analysis of the *q* = 4 model^{55}). In terms of the random coupling couplings of the *L* and *R* gamma matrices, $JijL$ and $JijR$, for *q* = 2, the constants $J$ and $J\u0303$ are defined by

where the left and right couplings are related to the couplings *J*_{ij} and *M*_{ij} by

For a discussion of the symmetries of *G*_{ab} ∼ Σ_{ab}, we refer to a study of the *q* = 4 model.^{55} As it stands, the integrals over Σ and *G* for the action (38) are not convergent, which is required to perform the integrations over *G*_{ab}. Convergence can be achieved by rotating

The rotations do not affect the saddle-point evaluation of the action integral, but as we shall see below, explicitly integrating out *G* for *q* = 2 simplifies the saddle-point analysis, and we prefer to use a convergent definition for this reason. The action that gives a convergent path integral is, then, given by

with *ξ*_{LL} = *ξ*_{RR} = 1 and *ξ*_{LR} = *ξ*_{RL} = *i*. Contrary to *q* = 4, the integrals over *G*_{ab}(*τ*_{1}, *τ*_{2}) are Gaussian and can be carried out exactly. This results in the effective action for Σ_{ab},

At the saddle point, Σ should be translation-invariant and, hence, only a function of *τ*_{1} − *τ*_{2}. For the purpose of saddle-point analysis, we can express the action in terms of Fourier modes of the Σ_{ab}, which is now a single-variable anti-periodic function,

where *ω*_{n} = (2*n* + 1)*π*/*β* are Matsubara frequencies defined already in Eq. (25). This results in

As already emphasized, contrary to the *q* = 4 case, in the *q* = 2 case, the SD equation simplifies to second-order equations (we took the limit *ϵ* → 0),

and other two equations with subscripts *L* and *R* interchanged. At the saddle point, we have that Σ_{RL}(*ω*_{n}) = −Σ_{LR}(*ω*_{n}) and Σ_{RR}(*ω*_{n}) = Σ_{LL}(*ω*_{n}). The saddle point equations couple positive and negative frequencies, but the solutions are simply related by

Using these relations, the saddle point equations are easily solved with a trivial solution given by (the symmetries of *G*_{ab} and Σ_{ab} are discussed in detail in the analysis of the *q* = 4 model^{55})

and a nontrivial solution that couples the left and right SYK models breaking the replica symmetry between them,

where we have introduced the critical frequency

This frequency will play an important role in the analysis of the partition function.

Pairing the positive and negative Matsubara frequencies allows us to write the free energy as a sum over only positive Matsubara frequencies,

For unbroken saddles, we have two solutions, and it turns out that one solution gives a larger −*S*_{E}, hence is the dominant saddle. This is the solution with the +sign(*ω*_{n}) term for Σ_{LL}. The free energy of this solution is equal to the free energy of two uncoupled SYK models^{59} [but with $J2=(1\u2212k2)v2$]. For each (positive) Matsubara frequency, this gives

while the free energy of the broken solution reduces to

The broken solution always gives the dominant action, but as we will see next, it does not always determine the free energy. The saddle point of Σ_{LL}(*ω*_{n}) is always purely imaginary, but the saddle point of Σ_{LR}(*ω*_{n}) switches from real to imaginary at *ω*_{n} = *ω*_{cr} where the free energy of the trivial and the nontrivial solution coincides. For *ω*_{n} > *ω*_{cr}, the imaginary part of the action is zero at the saddle point, but the action becomes complex along the integration manifold. In order to apply the steepest descent method, the integration manifold must be directed along the Picard–Lefschetz thimble. Otherwise, we will have large cancellations that may suppress the action of the saddle point that minimizes the free energy. It is a complicated problem to find the Lefschetz thimbles in a multidimensional space, but we can analyze the problem along the trajectory where Σ_{LL}(*ω*_{n}) and Σ_{RR}(*ω*_{n}) are at the saddle point, while for the off-diagonal Σ_{ab}(*ω*_{n}) variables, we restrict ourselves to the sub-manifold Σ_{RL}(*ω*_{n}) = −Σ_{LR}(*ω*)_{n} and Σ_{RL}(−*ω*_{n}) = −Σ_{LR}(*ω*_{n}), which intersects with the saddle-point. Combining positive and negative Matsubara frequencies, the action on this sub-manifold is given by

This action also arises in the study of a zero-dimensional Gross–Neveu-like model and its saddle point analysis,^{82,83} which we will apply here. For an illustration of the Lefschetz thimbles discussed below, we refer to Figs. 2 and 3 of the paper by Kanazawa and Tanizaki.^{82} The saddle points of this effective action are still given by Σ_{LR} of (48) and (49) of the full action, namely, Σ_{LR} = 0 and $\Sigma LR(\omega n)=\xb1J\u03031\u2212\omega n2/\omega cr2$. At all the saddle points, the action is real and the Lefschetz thimble of these saddle points is along the real axis if the saddle point solution for Σ_{LR} is real and along the imaginary axis when this saddle point is imaginary. In the latter case, i.e., for *ω*_{n} > *ω*_{cr}, the thimble ends are the zeros of the logarithm, and it is not possible to deform the real axis continuously into the thimble.^{91} Of course, we can deform the initial integration over the real axis to an integration path in the complex plane that goes over the saddle point on the imaginary axis. As long as we do not cross any singularities, by Cauchy’s theorem, the value of the integral along the deformed path will be the same despite that the integrand at the saddle point on the imaginary axis is much larger than that of the saddle points on the real axis. The phase of the integrand together with the Jacobian will assure that the contributions to the integral combine to the correct result. However, if we integrate only over the Gaussian fluctuations about the imaginary saddle point, we do *not* get the correct result. In order words, we cannot apply the saddle-point approximation to the imaginary saddle points. Instead, the integral can be evaluated at the trivial saddle point, which has its thimble on the real axis. For *ω*_{n} < *ω*_{cr}, the Σ_{LR}(*ω*_{n}) integral runs over both real saddle-points, but one of them is suppressed by the *ϵ* term in the action.

Strictly speaking, the Σ_{LR} = 0 saddle of the effective action (54) does not quite correspond to the Σ_{LR} = 0 saddle of the full action because in writing down the effective action, we already assumed Σ_{LL} of the form in solution (49), but the Σ_{LR} = 0 saddle of the full action belongs to solution (48) where Σ_{LL} takes a different form. Thus, the Σ_{LR} = 0 solution for the effective action should be viewed as a spurious saddle due to the sub-manifold constraint. Therefore, a more definitive analysis should be performed on the full action (the full action is quite similar, though not exactly the same, to the action of a zero-dimensional Nambu–Jona–Lasinio-like model^{82,83}). However, our analysis on the sub-manifold is indicative of the inaccessibility of nonzero imaginary Σ_{LR} solutions. We, thus, conclude that

The corresponding partition function is given by

The free energy of the replica diagonal solution is just the free energy of two decoupled SYK models.^{59}

Summing over all Matsubara frequencies, we obtain the total free energy

The first term can be evaluated using zeta function regularization,

This result gives the entropy of noninteracting Majorana particles. Our final expression for the free energy is given by

For *k* = 1, we have that $J=0$ and the last term vanishes. We note that the free-fermion expression (36) agrees with (59) after identifying the parameters by

The order parameter of the phase transition is given by Σ_{LR}. In Fig. 4, we show the analytical result (49) compared with a numerical calculation for *N*/2 = 256 that will be discussed in Sec. V.

The Σ model that is obtained after integration over the *G* variables resembles the usual random matrix theory *σ*-model. In Sec. V, we will derive basically the same *σ*-model directly starting from a partition function that is factorized into a product over Matsubara frequencies.

## V. NONLINEAR *σ*-MODEL FOR *q* = 2 PARTITION FUNCTION

From Eq. (16), we can express the partition function of the *q* = 2 Hamiltonian as

where *ɛ*_{k} are the eigenvalues of *i*(*J*_{ij} + *ikM*_{ij})/2 with positive real parts. Using the Weierstrass formula, this can rewritten in terms of a product over Matsubara frequencies $\omega n=2\pi (n+12)/\beta $,

where the product ∏_{n}*ω*_{n} has been evaluated to be $2$ by zeta function regularization [see Eq. (58)] and *h* is the *N*/2 × *N*/2 skew-symmetric matrix *J* + *ikM*. To leading order in 1/*N*, we have

Let us first consider a one-site SYK at frequency *ω*_{n}. The average partition function is given by

The determinant can be expressed^{84} as an integral over Grassmann variables *ϕ* and *ϕ**,

We recall that $\u27e8Jij2\u27e9=\u27e8Mij2\u27e9=J2/(1\u2212k2)$, so the Gaussian average over *h* can be performed by a cumulant expansion, resulting in

where we have used $\varphi i2=0$ due to their Grassmannian nature. Using the Hubbard–Stratonovich transformation,

we obtain

The saddle point equation is given by

The dominant solution is given by (we have *ω*_{n} > 0 to start with)

which leads to the one-site free energy,

This is exactly half of the Schwinger–Dyson result for the two-site model using only the replica-symmetric saddles. The total free energy is given by

This explicitly shows that the one-site annealed partition function gives the quenched result for the high-temperature phase of the non-Hermitian SYK model as is the case for the *q* = 4 SYK model. For *T* < *T*_{c}, the replica limit of the one-site partition function fails to give the quenched result. In order to get the correct result, we have to take the replica limit of the partition function and the conjugate partition function, which is well-known from the *σ*-model formulation of non-Hermitian random matrix theories^{27,66,74} and QCD at nonzero chemical potential.^{33,34,70,75}

Next, we consider the two-site non-Hermitian model and also assume that the ensemble average factorizes in the large *N* limit so that we can evaluate the partition function for a single frequency,

We can again average over *h* by a cumulant expansion using that

This results in the quartic action,

This action is invariant under

In addition to the Hubbard–Stratonovich transformation (67), we use the identity

to decouple the quartic terms. This results in

where $d\Sigma =d\Sigma LLd\Sigma RRd\Sigma LRd\Sigma LR*$ and $d\rho =d\rho LRd\rho LR*$ and we have used that the integral over the $\varphi ai$ variables factorizes into a product over *I*. Each of these factors is equal to

where *ϕ*_{a} no longer carry an index *i*. This integral can be evaluated by simply expanding the exponential in the integrand and collecting the terms that are proportional to $\varphi L\varphi L*\varphi R\varphi R*$. We get

The saddle points can be grouped in to the following three classes:

Now, the first two saddles are simply the ones we found in the SD equations in Sec. IV, and if the third saddle can be discarded, they would reproduce the same saddle-point analysis, which we recall here: although the second saddle (symmetry-breaking saddle) has a dominant action (larger −*S*_{E}) for all values of *ω*_{n}, thimble analysis requires us to pick the second saddle only for *ω*_{n} < *ω*_{cr}, and for *ω*_{n} > *ω*_{cr}, we must pick the first saddle (unbroken saddle). Let us now show why the third saddle can be discarded: the third saddle’s action (−*S*_{E}) is smaller than that of the second saddle for all *ω*_{n}, so we do not need to worry about it for *ω*_{n} < *ω*_{cr}; its action is smaller than that of the first saddle for *ω*_{n} > *ω*_{cr} (its action can be larger than that of the first saddle only for *ω*_{n} < *ω*_{cr}); hence, we do not need to worry about the third saddle for *ω*_{n} > *ω*_{cr} either. Thus, the third saddle drops out of our consideration for all values of *ω*_{n} and we reproduce the same free energy found in Sec. IV.

The order parameter of the phase transition of the coupled SYK model is given by the expectation value of Σ_{LR}. This is equal to

This is the chiral condensate corresponding to the spectral density of

It is given by the Banks–Casher formula,^{85}

The critical temperature is determined by the value of *ω*_{n} at which a gap opens and the spectrum becomes gapped for *T* > *T*_{c}. Another interpretation of the critical temperature is that *T*_{c} is the point at which *ω*_{n} enters the spectral support of *h*. The spectral density of $H$ can be obtained analytically and follows from the solution of a cubic equation.^{86,87} The phase transition is a typical Landau–Ginsberg phase transition with mean field critical exponents.

## VI. OUTLOOK AND CONCLUSIONS

In conclusion, the free energy of the integrable *q* = 2 SYK model is qualitatively different from the *q* > 2 case. Not only the order of the transition is different, but also there is an infinite series of transitions, while for *q* > 2, there is only one. This goes back to the factorization of the partition function into a product over Matsubara frequencies. Each of the factors undergoes a phase transition from a replica symmetric solution to a solution with broken replica symmetry with a critical temperature that depends on the Matsubara frequency. For the full partition function, this results in an infinite sequence of phase transitions.

We have calculated the quenched free energy in two structurally different ways: First, a quenched calculation based on the free-fermion description of the *q* = 2 SYK model; and second, an annealed calculation based on the solution of Schwinger–Dyson equations in the Σ*G* formulation of the SYK model using the replica trick. The two methods give the same result, which shows that despite the non-hermiticity of the model, the replica limit gives the correct result, provided that the starting point is the product of the one-site partition function and its complex conjugate before averaging. On the other hand, in the quenched free-fermion calculation, we did start from the one-site partition function without having to include the conjugate partition function. The reason that this gives the correct result is the factorization of the partition function into a product over single-particle energies.

The *q* = 4 non-Hermitian SYK model behaves quite differently. It is not a Fermi liquid, and the usual free-fermion description is invalid, which is most notable in the zero-temperature entropy, which is extensive.^{88} The non-Hermitian two-site *q* = 4 SYK model has a single first-order phase transition, which separates a low-temperature phase from a high-temperature phase. The high-temperature phase is entropy dominated, while the low-temperature phase is energy dominated. The free energy in the high-temperature phase follows from the many-body eigenvalue density and is entirely determined by the one-site partition function. This is also the case for the *q* = 2 SYK model. In the low-temperature phase, the replica limit of the one-site partition function breaks down and the quenched one-site partition function is given by the replica limit of the one-site partition function and its complex conjugate. In terms of the many-body spectral density, the two-point spectral correlation function determines the free-energy of the low-temperature phase. For this reason, there is a large difference between the *q* = 4 case and the *q* = 2 case. The dynamics of the *q* = 4 partition function is chaotic with universal eigenvalue correlations given by the Ginibre model, which as a consequence gives rise to a temperature-independent free energy in the low-temperature phase. On the other hand, the *q* = 2 SYK model is integrable with mostly but not entirely uncorrelated eigenvalues. We can distinguish two contributions from the two-point correlation function: one contribution is due to self-correlations and the second one is due to the many-body correlations resulting from the fact that the 2^{N/2} many-body eigenvalues are determined by *N*/2 single-particle energies. It is simple to evaluate the contribution from the self-correlations, but this only reproduces the free-energy at zero temperature and is smooth as a function of the temperature. This implies that the infinite series of second-order phase transitions are due to correlations of the many-body eigenvalues.

The phase transitions of the *q* = 2 non-Hermitian SYK model can also be understood in terms of the spectral properties of the coupling matrix. Using a random-matrix-theory-like *σ*-model calculation, we have related the order parameter of the phase transition Σ_{LR} to the formation of a gap of the hermiticized two-site Hamiltonian. In terms of the one-site Hamiltonian, this is the point where the Matsubara frequency enters the support of the spectrum of *H*_{L}. The starting point of the *σ*-model calculation is closely related to a random matrix model for the chiral phase transition in QCD where Σ_{LR} plays the role of the chiral condensate.

A natural question is whether the free energy of the *q* = 2 SYK model can also be understood in terms of the many-body spectral density and many-body spectral correlations. The cancellations that are responsible for the high-temperature phase of the *q* = 4 model are still at work for *q* = 2. For example, at high temperatures and for maximum non-hermiticity, the free energy of the *q* = 2 SYK model and the *q* = 4 SYK model is the same ($\u2212T2log\u20612$ per particle). From the solutions of the SD equations, it is clear that the two-point correlation function determines the low-temperature phase. In particular, the transitions observed in the low-temperature phase are due to the coupling between left and right sites but not by the dynamics within each of the sites.

In conclusion, the nature of the quantum dynamics plays a role for the replica symmetry breaking mechanism, which also induces phase transitions for the *q* = 2 non-Hermitian SYK model. However, while the replica dynamics of quantum chaotic systems is universal, there is a broad variety of dynamical behaviors associated with integrable systems, and we cannot conclude that the behavior we have observed for the non-Hermitian *q* = 2 SYK model is generic.

## ACKNOWLEDGMENTS

J.J.M.V. would like to thank Freeman Dyson for responding to him with a hand-written letter when he applied for a postdoc position in 1981. This response has encouraged his work on RMT throughout the years. Y.J. would like to thank Freeman Dyson for a conversation that happened in 2013 in Singapore, which has a lasting impact on Y.J.’s academic personality. Antonio García-García is thanked for collaboration in early stages of this project and providing eigenvalues of the *q* = 2 non-Hermitian SYK Hamiltonian. Gernot Akemann and Yuya Tanizaki are thanked for pointing out and explaining Refs. 79–83. Y.J. and J.J.M.V. acknowledge the partial support from the U.S. DOE under Grant No. DE-FAG-88FR40388. Y.J. was also partly funded by an Israel Science Foundation center for excellence grant (Grant No. 2289/18), by the United States–Israel Binational Science Foundation (BSF) under Grant No. 2018068, by the Minerva Foundation with funding from the Federal German Ministry for Education and Research, by the German Research Foundation through a German–Israeli Project Cooperation (DIP) grant “Holography and the Swampland,” by the Koshland postdoctoral fellowship, and by a research grant from Martin Eisenstein. D.R. acknowledges support from the Korea Institute of Basic Science (Grant No. IBS-R024-D1).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Yiyang Jia**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Dario Rosa**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). **Jacobus J. M. Verbaarschot**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: FREE-FERMION REPRESENTATION OF THE *q* = 2 SYK MODEL

We write the *q* = 2 one-site Hamiltonian as

where gamma matrices $\gamma \u20d7=(\gamma 1,\gamma 2,\u2026,\gamma M)$ with even *M* and *W* is a random antisymmetric complex matrix. In the main text, we have the convention that *M* = *N*/2. By matching with (6) and (7), we know

In the Hermitian SYK model (*k* = 0), *W* is an antisymmetric Hermitian matrix, and the Hamiltonian can be transformed into a fermion-filling form thanks to the fact that an *M*-dimensional antisymmetric Hermitian matrix has the normal form

where *O* is a real orthogonal matrix. Using a new basis for the *γ* matrices defined by $OT\gamma \u20d7$, one can easily write the Hermitian Hamiltonian in a fermion-filling form. Although a generic complex antisymmetric matrix does not have a normal form of Eq. (A3), there exists a parallel of it, which allows us to write the non-Hermitian Hamiltonian (A1) in a modified fermion-filling form. We will demonstrate this now.

We consider a diagonalizable complex antisymmetric matrix *W* with all eigenvalues being nonzero, which is almost always the case for our ensemble. The eigenvalues of *W* come in opposite pairs ±*ɛ*, so we can diagonalize *W* as

where

The column vectors of *S* are the eigenvectors of *W*, namely,

where

To completely fix the sign convention, we choose *ɛ*_{k} (*k* = 1, …, *M*/2) to have a positive real part. Because *W*^{T} = −*W*, we have

which implies that

Hence, if we scale the eigenvectors to redefine *S* as

we can easily see that

and hence, *S*^{−1} = Σ_{1}*S*^{T} (note, in general, *S*^{T}*S* ≠ *SS*^{T}). Substituting this into Eq. (A4), we obtain

We, thus, arrive at a normal form for complex antisymmetric matrices rather similar to that of the Hermitian antisymmetric ones (A3), with the difference that *S* is not orthogonal but satisfies *S*^{T}*S* = Σ_{1}. Now, we define a new set of operators ${c\u0303k,ck|k=1,2\u2026M/2}$ by

From the anti-commutation relation of *γ* matrices, we derive

This, in particular, implies that

Note that this is just the algebra for the ladder operators of *M*/2 spinless fermions, except that *c*_{k} and $c\u0303k$ are not related by a Hermitian conjugation. In terms of these ladder operators, the Hamiltonian (A1) becomes

Just as in the Hermitian case, we have

Hence, we conclude that the many-body energies of *H* are given by the filling of *M*/2 free fermions into the particle–hole symmetric levels of *W*: each fermion either occupies a “particle” level with energy *ɛ*_{k} or occupies a “hole” level with energy −*ɛ*_{k}. However, since the raising and lowering operators are not Hermitian conjugate to each other, the eigenstates are not necessarily orthogonal to each other (at least not with the original inner product ⟨*x*, *y*⟩ ≡ *x*^{†}*y*), just as one would expect for a non-Hermitian Hamiltonian.

### APPENDIX B: FREE ENERGY FROM THE FREE-FERMION REPRESENTATION FOR *k* < 1

Based on the free-fermion representation of the *q* = 2 SYK model of Appendix A, we proceed to the explicit analytical calculation of the free energy. The simpler spherical case *k* = 1 was already discussed in the main text. The free energy for *k* < 1 can be derived along the same lines, which is the purpose of this appendix.

For *k* < 1 (*υ*_{0} < *ϵ*_{0}), the large-*N* single-particle spectral density becomes a constant inside an elliptical disk as in Eq. (20). The elliptical disk can be parameterized by

We write

where

The subscript *E* denotes the “ellipse.” On its face, we cannot interpret *I*_{E}(*r*) as a complex contour integral as in Eq. (28) because *dz*/(*iz*) ≠ *dϕ* with the elliptical parameterization (B1). We can overcome this by considering the following conformal (Joukowski) transformation:

where

In terms *u*, the ellipse in Eq. (B1) at any given *r* becomes a unit circle,

We stress that *ϕ* here is the same *ϕ* as in the parameterization (B1). Now, we can write

where to obtain second equality, we applied Weierstrass factorization just as we did in the circular case, and the third equality simply defines *I*_{nE}. Note that *S*^{1} denotes the unit circle and the *r*-dependence of the integral comes from the *r*-dependence of *a* and *b*.

To analyze the cut and pole structures of the integral *I*_{nE}, we rewrite its integrand as

where *u*_{1±}, *u*_{2±} are the four roots of the equation,

namely,

Note that $4ab=(\u03f502\u2212\upsilon 02)r2>0$, and since *a* > *b*, *u*_{1±} are always inside the unit circle. We also note that for *u* = *e*^{iϕ},

Hence,

The integrand of *I*_{nE} has one pole at the origin with the residue

The integrand of *I*_{nE} has four branch cuts emanating from *u*_{1±}, *u*_{2±} horizontally to negative infinity. The *u*_{1±} cuts always intersect with the unit circle, whereas *u*_{2±} may or may not intersect with the unit circle depending on the values of *a* and *b*. To summarize, the *u*_{2±} cuts contribute to *I*_{nE} only if |*u*_{2±}| < 1 (which is to say *r* > *ω*_{n}/2*υ*_{0}), in much the same way as the ±*iω*_{n}/2*υ*_{0} cuts contribute to *I*_{n} in the circular case; what is new with the elliptical case are the *u*_{1±} cuts and the pole at the origin, which always contribute to *I*_{nE} regardless the value of *r*. Recycling the calculation done in the circular case, we obtain

if *r* < *ω*_{n}/2*υ*_{0}, which is the sum of one pole and two cut contributions. In addition,

if *r* > *ω*_{n}/2*υ*_{0}, which is the sum of one pole and four cut contributions. With these results, we arrive at

The result on the last line exactly matches with the SD calculation, provided that

## REFERENCES

^{†}and AII

^{†}

There are small deviations from −*E*_{0} when the temperature becomes closer to *T*_{c}. The nature of these deviations is not clear.

Since the coupling matrix is anti-symmetric, it eigenvalues occur in pairs as ±*ɛ*_{k}. This affects the level correlations close to *ɛ* = 0, but they converge rapidly to the level correlations of the Ginibre ensemble away from *ɛ* = 0. These correlations do not enter in the free energy discussed below.

To be precise, to get well-defined thimbles emanating from the zeros of the logarithm, the small symmetry-breaking term proportional to *ϵ* must be included, and each zero will give a separate thimble. But the basic conclusion remains the same: these two thimbles do not contribute to the path integral because the original contour of integration cannot be deformed into either of them.^{82,83}