Symmetry Classification and Universality in Non-Hermitian Many-Body Quantum Chaos by the Sachdev-Ye-Kitaev Model

Spectral correlations are a powerful tool to study the dynamics of quantum many-body systems. For Hermitian Hamiltonians, quantum chaotic motion is related to random matrix theory spectral correlations. Based on recent progress in the application of spectral analysis to non-Hermitian quantum systems, we show that local level statistics, which probe the dynamics around the Heisenberg time, of a non-Hermitian $q$-body Sachdev-Ye-Kitev (nHSYK) model with $N$ Majorana fermions, and its chiral and complex-fermion extensions, are also well described by random matrix theory for $q>2$, while for $q = 2$, they are given by the equivalent of Poisson statistics. For that comparison, we combine exact diagonalization numerical techniques with analytical results obtained for some of the random matrix spectral observables. Moreover, depending on $q$ and $N$, we identify $19$ out of the $38$ non-Hermitian universality classes in the nHSYK model, including those corresponding to the tenfold way. In particular, we realize explicitly $14$ out of the $15$ universality classes corresponding to non-pseudo-Hermitian Hamiltonians that involve universal bulk correlations of classes ${\rm AI}^\dagger$ and ${\rm AII}^\dagger$, beyond the Ginibre ensembles. These results provide strong evidence of striking universal features in non-unitary many-body quantum chaos, which in all cases can be captured by nHSYK models with $q>2$.


I. INTRODUCTION
Calculations in strongly interacting quantum systems are generically difficult. The spectrum of the Hamiltonian is arguably one of the least expensive quantities that can be computed numerically. Moreover, the spectrum is basis invariant and that, especially in many-body systems, is a substantial advantage. These facts help us understand the relevance of the Bohigas-Giannoni-Schmit (BGS) conjecture [1], which states that spectral correlations of quantum chaotic systems at the scale of the mean level spacing are given by the predictions of random matrix theory. It provides a simple, both conceptually and practically, but still very powerful link between the analysis of the spectrum and the quantum dynamics. More importantly, it points to a very robust universality of the late-time quantum dynamics under relatively mild conditions and with relatively well-understood exceptions due to phenomena like Anderson localization or integrability. Global anti-unitary symmetries, and not dynamical features, label different universality classes, which are rigorously classified using random matrix theory (RMT), the so-called tenfold way [2]. Although the Berry-Tabor conjecture [3] states a similar spectral characterization of integrable systems, its applicability is more restricted because, in some sense, each integrable system is integrable in its own way.
As an illustration of the aforementioned universality, the study of spectral correlations to characterize quantum motion encompasses multiple disciplines with no direct relation between them.
It started 70 years ago in the context of nuclear physics, where Wigner showed [4] that short-range spectral correlations of highly excited states of nuclei are well described by random matrix theory. In the 1970s and 1980s, the interest gradually shifted from nuclear physics to the dynamics of single-particle Hamiltonians in both deterministic [1] and random [5] potentials. Later, it was successfully applied [6,7] to describe the level statistics of the lattice QCD Dirac operator in the presence of gauge configurations. The chiral symmetry of the Dirac operator led to the proposal of universality classes in the context of random matrix theory, which were the direct precedent of the tenfold way [8]. In recent years, it has also found applications in quantum gravity after the proposal [9] that a certain universal bound on chaos is saturated in field theories with a gravity dual. More specifically, it is a bound on the growth of particular out-of-time-order correlation functions in the semiclassical limit. For quantum chaotic systems, the growth is exponential around the Ehrenfest time with a rate controlled by the classical Lyapunov exponent. In the early days of quantum chaos theory [10,11], this exponential growth was broadly employed to characterize quantum chaos in single-particle Hamiltonians. three chiral classes [61][62][63][64] that have found applications, e.g., in the context of Euclidean QCD at finite baryonic chemical potential. However, according to the full classification of symmetries in non-Hermitian random matrices [65,66], there are 38 non-Hermitian universality classes. It turns out that the Ginibre ensemble is only one of the three universality classes [67,68] for local level statistics in non-Hermitian random matrix theory. The other two, called AI † and AII † , are defined in terms of transposition symmetry instead of complex conjugation. For instance, they describe the spectral correlations of the QCD Dirac operator in two-color QCD coupled to a chiral U (1) gauge field [39,61].
The progress in the development of a full classification scheme has not yet fully translated into a systematic spectral characterization of non-Hermitian quantum chaotic systems, see Ref. [69] for a recent study focused on open fermionic quantum matter. First, there is no equivalent of the BGS conjecture, so the relation between dynamics and level statistics is unclear. There are also additional technical problems: Correlations of complex eigenvalues are weakened, and the necessary unfolding of eigenvalues may be problematic [47], in particular, when the eigenvalue distribution is not radially symmetric. However, these problems have been ameliorated in the last years with the introduction of spectral observables that do not require unfolding for short-range correlators, such as the ratio of spacings between nearest-neighbor eigenvalues [70], which have found applications in the study of collective-spin Liouvillians [71], non-Hermitian Anderson transitions [72][73][74], directed random graphs [75], nonunitary open quantum circuits [76,77], two-color QCD [39], and the classicalquantum transition [52]. The study of long-range correlators such as the number variance [78][79][80][81] or spectral form factor [82,83] (which requires unfolding) further suggests that some weakened form of spectral rigidity is still present in non-Hermitian systems and will be the subject of a separate publication [84].
To address the issues mentioned in the previous paragraph, one needs a many-body system where it is possible to test the results of random matrix theory, including the existence of many more universality classes; at the same time this system should be simple enough to be amenable to an analytical and numerical treatment so that it is possible to independently probe the nature of the non-Hermitian quantum dynamics. In this paper, we propose that this model is the non-Hermitian SYK (nHSYK) model. We focus our study on the local level statistics of this model, which probe the quantum dynamics for long timescales of the order of the Heisenberg time. For q = 2, we find that the spectrum is largely uncorrelated and well described by the equivalent of Poisson statistics.
For q > 2, there is excellent agreement with the predictions of non-Hermitian random matrix theory in short-range correlators like the complex spacing ratio [70] and microscopic correlators near the   [19,[85][86][87][88] in the literature. However, both the employed models and the focus of these studies are quite different from ours. In Ref. [19], the gravity dual of two complex conjugate copies of a q = 4 nHSYK model was identified as an Euclidean wormhole by an analysis of the free energy. The role of replica-symmetry-breaking solutions in the low-temperature limit of the free energy of two copies of an nHSYK has been investigated in Ref. [85]. The entanglement entropy and an emergent replica conformal symmetry were recently studied in chains of q = 2 nHSYK models [87,88]. The Page curve [89], describing the process of black hole evaporation [86], was computed by an effectively open SYK model, where the role of the environment is played by a q = 2

SYK.
We study a single non-Hermitian SYK model of N Majorana fermions in (0 + 1) dimensions with q-body interactions in Fock space, with complex instead of real couplings: where J i 1 ···i q and M i 1 ···i q are Gaussian-distributed random variables with zero average and variance [12,90] and {ψ i , ψ j } = 2δ ij . The Majorana fermions can be represented by 2 N/2 -dimensional Hermitian Dirac γ matrices, and below, we denote them by γ i . For our analytical considerations, we work in a basis where the odd-and even-numbered γ matrices are represented by real symmetric and purely imaginary antisymmetric matrices, respectively. We shall consider the cases q = 2, 3, 4, and 6, and the number of fermions, N , is always even.
The Hamiltonian has, in general, a complex spectrum. However, depending on q, it could have some additional symmetry. For instance, for q = 3, we see that it has a chiral symmetry where E i is, in general, complex.

III. SPECTRAL DENSITY
We start our analysis with the study of the spectral density of the nHSYK model. We compute its spectrum by exact diagonalization techniques for values of N ≤ 28. For each set of parameters, we obtain at least 10 6 eigenvalues. We recall that for real couplings and a fixed finite N 1, q > 2, or [91], the spectral density close to the ground state grows exponentially, a feature also typical of quantum black holes. For odd q, when the "Hamiltonian" is the supercharge of a supersymmetric theory, the system (i.e., the supercharge) has chiral symmetry. The spectral density close to zero energy, controlled by the chiral symmetry, is well described [22,27]   For q = 2 and q = 4-see Figs. 1 and 3, respectively-the spectral density has a maximum in the central region of small |E| with a monotonous decay towards a rather sharp edge for q = 4. In contrast, for q = 3, the spectral density depicted in Fig. 2 has a richer structure. Depending on N , we observe different degrees of suppression, or enhancement, of the spectral density, followed by a maximum at intermediate distances.
In some cases, we also observe universal oscillations close to zero on the scale of the level spacing due to a spectral inversion symmetry E → −E. For q = 6, see Fig. 4, we also observe the effect of inversion symmetry for N = 20 and N = 24, but otherwise, the spectral density is almost constant. In the next section, we study in more detail the relation between these special spectral features and additional global symmetries of the nHSYK Hamiltonian. This will eventually lead to a full match between different universality classes and  . As in the previous cases, the average spectral density is radially symmetric, but, unlike the q = 3 case, it is rather unstructured, with a broad maximum located in the central part, followed by a slow decay for larger |E| energies. Finally, as is the case for q = 3, a rather sharp edge is observed where the density vanishes abruptly. Because of the absence of inversion symmetry, the spectral density is not suppressed or enhanced in the |E| ∼ 0 region, as was the case for q = 3.
specific nHSYK Hamiltonians, depending on N and q, that implement them.
In the Hermitian case, the classification simplifies to the symmetries of Eqs. (3), (4), and (7). In Ref. [66], there is one more symmetry given by the substitution H → iH in Eq. (8) (pseudo-anti-Hermiticity), but this does not give any new classes for non-Hermitian Hamiltonians. In the case of pseudo-Hermiticity, Eq. (8), one would normally consider the modified Hamiltonian ηH, which is Hermitian, as (ηH) † = ηH. A well-known example of this is the Wilson Dirac operator, where the symmetry is known as Γ 5 -Hermiticity [97]. In the Bernard-LeClair (BL) classification scheme [65,94], Eqs. (3) and (4) involve complex conjugation of the Hamiltonian matrix H and are referred to as K symmetries; Eqs. (5) and (6) describe a transposition symmetry of H and are dubbed C symmetries; chiral symmetry, Eq. (7), is called P symmetry; and pseudo-Hermiticity, Eq. (8), is referred to as Q symmetry in Ref. [65]. As in the Hermitian case, chiral symmetry and pseudo-Hermiticity can be written as a composition of time-reversal and particle-hole symmetries, and they are only non-trivial in the absence of the latter. Furthermore, they can either commute or anti- For the nHSYK model, we consider the charge-conjugation operators [21,23,24,26] where K denotes the complex-conjugation operator, which square to The combination of these two operators yields the Hermitian operator, which squares to the identity. This operator is also known as Γ 5 or the chiral-symmetry operator.
The operators P and R act on Majorana fermions as the nHSYK Hamiltonian has the involutive symmetries Comparing with Eqs. (3)-(8), we see that P and R play the role of either C + or C − , while S either commutes or anti-commutes with the Hamiltonian. For the nHSYK model, the many-body matrix elements are manifestly complex, and no anti-unitary symmetries that map H back to itself exist. Then, only the involutive symmetries (5)-(7) can occur. [98] The nHSYK model thus belongs to one of the ten BL symmetry classes without reality conditions [65], which are in oneto-one correspondence with the Hermitian Altland-Zirnbauer (AZ) classes [8] and are summarized in Table I. The difference between the AZ classes and the BL classes without reality conditions is that complex conjugation is replaced by transposition and the Hermiticity constraint is lifted. For example, a Hermitian real Hamiltonian is to be replaced by a non-Hermitian complex symmetric Hamiltonian. Following the Kawabata-Shiozaki-Ueda-Sato nomenclature [66], they are dubbed Here, we have adopted a shorthand notation, where a subscript in the name of a class indicates that the chiral symmetry is commuting (subscript +) or anti-commuting (subscript −) with the time-reversal or particle-hole symmetries.
The nHSYK symmetry classification can be performed systematically by evaluating Eqs. (10) and (12)-(14) for different values of q mod 4 and N mod 8. Note that while the physical interpretation of the operators P, R, and S is different in the SYK and nHSYK models, the defining relations in Eqs. (10) and (12)- (14) are formally the same. It follows that the symmetry classification of the former [16,17,[21][22][23][24][25][26] also holds for the latter, provided that one replaces any reality condition by a transposition one. We now investigate in more detail the dependence of these symmetries of the odd or even nature of q in the nHSYK Hamiltonian Eq. (1).
Even q. According to Eq. (14), H commutes with S (which is proportional to the fermion parity operator), the Hilbert space is split into sectors of conserved even and odd parity, and the Hamiltonian is block-diagonal. There is no chiral symmetry. From Eqs. (12) and (13), we see that H transforms similarly under both P and R (when they act within the same block) and it suffices to consider the action of one, say P. We have the commutation relation • When N mod 8 = 2, 6, P is a fermionic operator that anti-commutes with S. Note that P is not an involutive symmetry of the Hamiltonian in a diagonal block representation, as it maps blocks of different parity into each other. The two blocks are the transpose of each other and have no further constraints (class A or complex Ginibre).
• When N mod 8 = 0, 4 and q mod 4 = 0, P is a bosonic operator that commutes with S. Each block of the Hamiltonian has the involutive symmetry, PH † P −1 = +H. If N mod 8 = 0, P 2 = Table I. Non-Hermitian symmetry classes without reality conditions, nine of which are realized in the non-Hermitian SYK model (i.e., all except class AIII † ). For each class, we list its anti-unitary and chiral symmetries, an explicit matrix realization [92], its name under the Kawabata-Shiozaki-Ueda-Sato classification [66], and the corresponding Hermitian ensemble. In the matrix realizations, A, B, C, and D are arbitrary non-Hermitian matrices unless specified otherwise, and empty entries correspond to zeros. In the last column, we list the AZ classes [8] that are in one-to-one correspondence with the non-Hermitian classes without reality conditions.
Matrix realization Class Hermitian corresp.
1, and we can find a basis in which the Hamiltonian is symmetric. This is the universality class of complex symmetric matrices, also known as AI † . If N mod 8 = 4, P 2 = −1, we can find a basis in which H = IHI −1 , with I the symplectic unit matrix. This class is AII † .
• When N mod 8 = 0, 4 and q mod 4 = 2, P is again a bosonic operator that commutes with P 2 = 1, we can find a basis in which the Hamiltonian becomes anti-symmetric and the universality class is given by that of complex anti-symmetric matrices (non-Hermitian class D); if N mod 8 = 4, P 2 = −1, and we can find a basis where H = −IHI −1 . Complex matrices satisfying this constraint belong to non-Hermitian class C.
Odd q. In this case, S is a chiral symmetry operator that anti-commutes with H, so that H acquires an off-diagonal block structure in a chiral basis. The operators P and R now act differently on H (one satisfies XH † X −1 = H and the other XH † X −1 = −H with X either P or R). Hence both must be considered if we use the square of P, R, and S to classify the matrices. However, for the derivation of the block structure, as given in Table I, we of course only need either P or R, and S.
• When N mod 8 = 0, both P and R are bosonic operators squaring to +1. Since they commute with S, and one of them satisfies Hamiltonian in a suitable basis is a complex matrix with vanishing diagonal blocks and offdiagonal blocks that are the transpose of each other (class AI † + ), irrespective of whether q mod 4 = 1 or 3.
• When N mod 8 = 4, both P and R are bosonic operators commuting with S and squaring to −1. Hence the off-diagonal blocks A and B of H are related by B = IAI −1 . This is the universality class AII † + , irrespective of whether q mod 4 = 1 or 3.
• When N mod 8 = 2, 6 we have that {S, P} = 0 and {S, R} = 0. So P and R have vanishing diagonal blocks, and depending on N , P 2 = 1 and R 2 = −1 (N mod 8 = 2), or P 2 = −1 and . We choose the operator that squares to 1. Since it anti-commutes with S, it has the block structure The complete symmetry classification of the nHSYK Hamiltonian for all q and even N in terms of BL classes is summarized in It is important to stress that, since the symmetry classification is algebraic, it is not necessarily related to specific features of spectral correlations that probe the quantum dynamics for long timescales of the order of the Heisenberg time. However, we see in the next two sections that this is the case. By studying local bulk and hard-edge spectral correlations for different values of q and N , we find not only excellent agreement with the predictions of non-Hermitian RMT for q > 2 but also that the different universality classes resulting from the symmetry classification can be characterized by the analysis of level statistics. Assuming that the relation between RMT level statistics and quantum chaos persists for non-Hermitian Hamiltonians, our results provide direct evidence that the nHSYK model is versatile enough to describe all possible quantum ergodic states to which quantum chaotic systems with a complex spectrum relax after a sufficiently long time.

V. LEVEL STATISTICS: COMPLEX SPACING RATIOS
We initiate our analysis of spectral correlations by studying spacing ratios, also called adjacent gap ratios,which are a spectral observable that has the advantage of not requiring unfolding. It was introduced to describe short-range spectral correlations of real spectra [102][103][104] but it has recently been generalized [70] to complex spectra. Because of their short-range nature, these ratios probe the quantum dynamics for late timescales of the order of the Heisenberg time. However, as was mentioned previously, a full dynamical characterization of level statistics, equivalent to the BGS conjecture [1], for non-Hermitian systems is still missing. Therefore, strictly speaking, we perform a around the origin and along the real positive semi-axis is observed for q > 2. This is also the random matrix theory prediction [70], which should indicate quantum chaotic dynamics.
comparison with the predictions of non-Hermitian random matrix theory, implicitly assuming that good agreement still implies quantum chaotic motion.
We define the complex spacing ratio as where E k with k = 1, 2, . . . , 2 N/2 is the complex spectrum for a given disorder realization, E NN k is the closest eigenvalue to E k using the standard distance in the complex plane, and E NNN k is the second closest eigenvalue. By construction, |λ k | ≤ 1, so it is restricted to the unit disk. This is the most natural generalization of the real spacing ratio to the complex case. We perform ensemble averaging, such that we have a minimum of 10 6 eigenvalues for each set of parameters N, q. The distribution of the resulting averaged complex spacing ratio λ k is depicted in Figs. 5 and 6. We observe qualitative differences between q = 2 and q > 2. For the former, it is rather unstructured (i.e., flat) with no clear signature of level repulsion for small spacing. In the real case, this is a signature of the absence of quantum chaos. By contrast, for q > 2, the complex spacing ratio is heavily suppressed for small spacings, especially at small angles, which is a signature of level repulsion. Indeed, a very similar pattern is observed for non-Hermitian random matrices [70].
We now show that the three universality classes of local bulk correlations-AI † , A (GinUE), and AII † [67]-with increasing level repulsion, can be clearly distinguished by the complex spacing ratio distribution.
In order to gain a more quantitative understanding of the spectral correlations, we compute the marginal angular, ρ(θ), and radial, ρ(r), complex spacing ratio distributions, where λ k = r k e iθ k .
The results, presented in Figs. 7 and 8, confirm the existence, depending on N and q, of the three universality classes mentioned above. However, only for q > 2, do we observe agreement with the random matrix prediction, which indicates that, as in the real case, this is a requirement for non-Hermitian many-body quantum chaos. We note that the full spectrum was employed in the evaluation of these marginal distributions. Therefore, ρ(r) and ρ(θ) cannot distinguish between, for instance, class A and classes C and D, as the last two only differ from class A in the region |E| ∼ 0 of small eigenvalues. In summary, the complex spacing ratios of the nHSYK Hamiltonian (1) distinguish universality classes A, AI † , and AII † . To get a more visual confirmation of the symmetry classification, we can characterize the complex spacing ratio distribution by a single number, say its first radial moment, r = dr rρ(r), as a function of N and q. The values of r for the three universal bulk statistics (A, AI † , and AII † ) are given in Table III. The results presented in Fig. 9 show that r computed numerically for the nHSYK model closely follows the predicted RMT pattern for q = 3, q = 4, and 6, while it goes to the Poisson value for q = 2.
These results confirm the predictions of Table II for the local bulk correlations (level repulsion).
In the next section, we study the distribution of the eigenvalue with the lowest absolute value.
The shape of this observable is expected to have a universal form for each universality class that should agree with the random matrix prediction provided the spectrum has one of the inversion symmetries studied previously. This enables us to identify additional universality classes depending on the type of inversion symmetry of the nHSYK Hamiltonian by the study of level statistics.  Table II in all cases.

VI. LEVEL STATISTICS: HARD-EDGE UNIVERSALITY
Through the use of complex spacing ratios, we can only distinguish three universality classes of bulk correlations. In addition, the classes with spectral inversion symmetry (chiral or particle-hole) show universal repulsion from the spectral origin, the so-called hard edge for real spectra. This universal behavior can be captured by zooming in on the eigenvalues closest to the origin, on a scale of up to a few level spacings, the so-called microscopic limit. In particular, the distribution of the eigenvalue with the smallest modulus, P 1 (|E 1 |), gives, when combined with the bulk complex spacing ratio distribution, a measure to uniquely distinguish the ten non-Hermitian symmetry classes without reality conditions. As an example, in Fig. 10, we show the distribution of |E 1 | for the nHSYK model with q = 6 and N = 20 and 24, and compare it with the prediction of non-Hermitian random matrix theory for classes C and D, respectively. In order to carry out a parameter-free comparison with the nHSYK model, we normalize the distribution to unity and rescale |E 1 | by its average. We thus see that, while the q = 6 nHSYK Hamiltonian has the same bulk statistics for all N , we can still resolve the Bott periodicity, which enables us to distinguish A convenient way to capture the hard-edge universality by a single number is to consider the ratio (normalized variance) following the proposal by Sun and Ye for Hermitian random matrices [26]. The values of R 1 for the seven non-Hermitian classes with inversion symmetry listed in Table I are tabulated in Table III.
To further confirm our symmetry classification, in Fig. 11, we show the value of R 1 as a function of N for the q = 3 and q = 6 nHSYK models. We again see excellent agreement with the random Table III. Universal single-number signatures of the non-Hermitian universality classes without reality conditions. The first radial moment r of the complex spacing distribution measures the bulk level repulsion of the three universal bulk classes A, AI † , and AII † , while the ratio R 1 in Eq. (18) gives the repulsion between the hard edge and the eigenvalue with smallest absolute value for the seven classes with spectral inversion symmetry. The values of r are obtained by numerical exact diagonalization of 2 15 ×2 15 random matrices of the corresponding class averaging over an ensemble of 2 8 realizations. Meanwhile, the value of r has been obtained analytically for class A [105] in the thermodynamical limit. Its value is equal to r = 0.73866, which is in agreement with our numerical value. In order to compute the ratio R 1 in Eq. (18), we numerically diagonalize 10 7 100×100 matrices of the corresponding universality class. Note that for the complex spacing ratio distribution (and its moments), it was shown in Ref. [70] that they have large finite-size corrections for Gaussian-distributed random matrices; hence, very large matrices have to be considered to converge to the universal result of the thermodynamical limit. In contrast, we have verified numerically that the smallest eigenvalue distribution (and thus R 1 ) converges to a universal distribution very rapidly with the system size; therefore, using relatively small matrices is justified. Hamiltonians in the nHSYK model. Indeed, there is only one universality class, AIII † (also known as chGinUE), of the original tenfold way still to be identified in nHSYK. In this section, we consider generalizations of the nHSYK that belong to AIII † and to several other universality classes.

A. Symmetry classification of chiral nHSYK models
We start by investigating models with a built-in chiral symmetry, where the Hamiltonian, termed non-Hermitian chiral or Wishart SYK (nHWSYK) [27], in an appropriate basis is represented by two off-diagonal blocks, where each block is either an Hermitian or non-Hermitian SYK. By tuning N and q in this model, we describe six new universality classes. More specifically, we consider the For all even q, if N mod 8 = 2, 6 and the couplings are complex, H and H have no symmetries, and the chiral symmetry is the only symmetry of H, which thus belongs to class AIII † . This is the only non-Hermitian class without reality conditions (see Table I Table IV, we list these extra classes, together with their defining symmetries and a matrix representation.  Table II. For q = 3, the Hamiltonian has a chiral symmetry for all even N , while for q = 6, there is a particle-hole symmetry only if N is a multiple of 4. We see excellent agreement between the nHSYK results and the RMT predictions for all available system sizes. If q mod 4 = 0, N mod 8 = 0, and the couplings are complex, then H and H belong to class AI † and have a transposition symmetry with anti-unitary P, PH † P −1 = +H, see Eq. (12), that squares to P 2 = +1. It then follows that H has a time-reversal symmetry C + H † C −1 + = +H with anti-unitary C + = antidiag(P, P) that squares to C 2 + = +1 and anti-commutes with Π. We thus see that H belongs to class AI † − (see also Table IV). Proceeding similarly, we find that if q mod 4 = 2 and N mod 8 = 4, H also belongs to class AI † − , while for q mod 4 = 0, N mod 8 = 4 and q mod 4 = 2, N mod 8 = 0, it belongs to class AII † − . These two classes were already implemented in the original nHSYK model.
If q mod 4 = 0, N mod 8 = 0, and the couplings are real, then H and H are Hermitian SYK models belonging to the AZ class AI. As in the preceding paragraph, H has transposition antiunitary symmetry C + H † C −1 + = +H with C 2 + = +1 that anti-commutes with the chiral symmetry operator Π, but it also has additional pseudo-Hermiticity implemented by η = antidiag(1, 1) that anti-commutes with Π and commutes with C + . The combined transposition symmetry and Table IV. Five non-Hermitian symmetry classes with reality conditions realized in the nHWSYK model. For each class, we list its anti-unitary and chiral symmetries and their commutation relations, an explicit matrix realization [92], and its name under the Kawabata-Shiozaki-Ueda-Sato classification [66]. The symbols ε XY indicate whether the two operators X and Y commute, e.g., T HT −1 = ε T H H, while the symbols η X denote the square of operator X, e.g., T 2 = η T . In the matrix realizations, A, B, C, and D are arbitrary non-Hermitian matrices unless specified otherwise, and empty entries correspond to zeros.
The symmetry classification of the chiral two-matrix model H with Hermitian and non-Hermitian SYK blocks is summarized in Tables. V and VI.
Level statistics in the chiral class AIII † Instead of an extensive numerical confirmation of these new universality classes in the SYK model, we focus on the AIII † universality class that completes the tenfold way. We note that only recently [27] have the correlations of the Hermitian version of AIII † , the chGUE universality class, been identified in the SYK model. We study level statistics for the chiral nHSYK Hamiltonian (19) with H, H blocks given by independent q = 4 nHSYK models. For q = 4, AIII † requires N mod8 = 2, 6 so we stick to N = 22, 26. We obtained more than 10 7 eigenvalues in each case by exact diagonalization techniques. In analogy with the previous cases, we start with the analysis of shortrange bulk spectral correlations represented by the radial ρ(r) and angular ρ(θ) distributions of the complex spacing ratio density. Results depicted in Fig. 12 confirm an excellent agreement with the AIII † random matrix prediction for both N . However, this is a bulk observable, so the result for AIII † is identical to that of the A universality class.
In order to differentiate these two universality classes, we study the distribution of the eigenvalue with the smallest modulus |E 1 |. Because of the chiral symmetry, we expect the distribution to be universal and very well approximated by the random matrix prediction for AIII † .
We first compute the normalized variance R 1 , Eq. (18). For N = 22, with around 20000 disorder realizations, we obtain R 1 = 1.131, while for N = 26, with 12000 disorder realizations, we get R 1 = 1.128. This is in excellent agreement with the random matrix result R 1 = 1.129. The agreement is not restricted to this low-order moment. The full distribution of |E 1 |, see Fig. 13, is very close to the numerical random matrix result. The comparison does not involve any fitting procedure, only a rescaling by the average in each case. This is further confirmation that this chiral nHSYK model belongs to the AIII † universality class.
For small |E|, the radial spectral density behaves as We conclude that the radial distribution of the smallest eigenvalue must also behave as |E 1 | 3 log |E 1 | for small |E 1 |, which is in agreement with numerical results. It is also possible to obtain an excellent analytical estimation of the random matrix result for the full distribution of |E 1 | (see also Refs. [108,109]). We consider a random matrix W from class AIII † with square D × D offdiagonal blocks (recall Table I). The 2D eigenvalues of W come in symmetric pairs and are denoted We introduce D new variables z j = E 2 j whose joint distribution is given by [110] Because the distribution (23) is invariant under permutations of eigenvalues z j , we can choose z 1 as the eigenvalue with the smallest absolute value. Then, by construction, the distribution of z 1 is obtained by integrating out all the remaining eigenvalues outside the disk centered at the origin and with radius |z 1 |: where we denote z j = |z j |e iϕ j in polar coordinates. We proceed by evaluating the integrals in , and changing variables back to the chiral eigenvalue |E 1 | = |z 1 |, we obtain where the normalization constant is equal to N = 32c. The arbitrary energy scale can be chosen as  In conclusion, the study of spectral correlations shows that the q = 4 chiral nHSYK model, Eq. (19), belongs to class AIII † , and its dynamics is quantum chaotic for sufficiently long timescales.
We expect that a similar agreement with the random matrix predictions will be obtained for the Table VII. Non-Hermitian symmetry classes with reality conditions but no pseudo-Hermiticity, four of which are realized in the complex-fermion nHSYK model (i.e., all except class AI − ). For each class, we list its antiunitary and chiral symmetries, an explicit matrix realization [92], its name under the Kawabata-Shiozaki-Ueda-Sato classification [66], and its name as a Ginibre ensemble (class AI − has no conventional Ginibre name). In the matrix realizations, A, B, C, and D are arbitrary non-Hermitian matrices unless specified otherwise, and empty entries correspond to zeros.
other chiral nHSYK models introduced in the previous subsection.

C. Classes with complex-conjugation symmetry in nHSYK models with complex fermions
In the previous subsections, we completed the tenfold BL classes without reality conditions.
There exist only five remaining classes without pseudo-Hermiticity (two non-chiral and three chiral classes), which are listed in Table VII. [111] These classes involve complex-conjugation symmetry but no transposition symmetry (because simultaneous complex-conjugation and transposition symmetries imply pseudo-Hermiticity). Four of these classes emerge naturally if we consider complexfermion nHSYK models, i.e., a Hamiltonian where w ij;k = −w ji;k = −w ij, k are independent random variables (real, complex, or quaternionic) and c † i and c i are the creation and annihilation operators of N c = N/2 complex fermions, respec-tively. The complex fermions can be expressed in terms of Majorana fermions as c i = ψ 2i−1 + iψ 2i , In a basis where the odd γ-matrices are real and the even ones are purely imaginary (which is our convention throughout this paper), the creation and annihilation operators are real, Since we have if we set w ij;k = w * k ;ij , H c is Hermitian and we recover the two-body random ensemble [30]. However, if we relax this condition and let w ij;k and w k ;ij become independent, H c is non-Hermitian and has complex eigenvalues. The independence of w ij;k and w k ;ij , together with Eq. (27), further precludes the existence of both transposition symmetries, Eqs. (5) and (6), and complex-conjugation particle-hole symmetry, Eq. • If w ij;k is a quaternion, i.e., with w ij;k and w (2) ij;k complex, there exists a basis where H c = IH * c I −1 , with I the symplectic unit matrix, because the matrix elements of c i and c † i are real (and a real number times a self-dual quaternion remains a self-dual quaternion). We see that H c belongs to the Ginibre Symplectic Ensemble (GinSE) or non-Hermitian class AII.
Proceeding as before, we can now consider a chiral or Wishart complex-fermion nHSYK model, i.e., a Hamiltonian given by Eq. (19) with off-diagonal blocks given by H c . These models again have a chiral symmetry implemented by the unitary operator Π = diag(1, −1). As before, the symmetry classification depends solely on the reality conditions of the couplings. • If the couplings are complex, the chiral symmetry is the only symmetry of the chiral Hamiltonian and it immediately follows that it belongs to class AIII † .
• If the couplings are real, each block has an anti-unitary symmetry realized by K (complex conjugation), which implies the chiral Hamiltonian satisfies T + HT −1 + = H, with T + = diag(K, K) squaring to T 2 + = +1. Because Π and T + commute, this defines the universality class AI + .
• If the couplings are quaternionic, the anti-unitary symmetry of the blocks, realized by I with which squares to T 2 + = −1 and again commutes with the chiral unitary Π. This is the universality class AII + .
The symmetry classification of the non-chiral and chiral nHSYK models with complex fermions is summarized in Table VIII. The final non-Hermitian symmetry class, AI − , consisting of chiral block matrices with zero diagonal blocks and off-diagonal blocks being the complex conjugates of each other (see last line in Table VII), is realized neither in the nHSYK nor in its chiral or complex-fermion extensions. It would be interesting to find a variant of the model that implements this symmetry class, and, in that way, complete the full fifteenfold classification of non-pseudo-Hermitian matrices by the nHSYK model.

VIII. OUTLOOK AND CONCLUSIONS
We have proposed a single-site SYK model with complex couplings as a toy model for manybody non-Hermitian quantum chaos. We have shown that 19 out of 38 universality classes in the symmetry classification of non-Hermitian quantum systems occur naturally in this SYK model and its chiral and complex-fermion extensions. In particular, we have reproduced all but one (AI − ) non-pseudo-Hermitian universality classes, namely, A, AIII † , AI † , AII † , D, C, AI † + , AII † + , AI † − , AII † − , AI, AII, AI + , and AII + , see Tables. II, V, VI, and VIII. A detailed spectral analysis of the original nHSYK model, involving short-range correlators such as the complex spacing ratios in the bulk and the distribution of the smallest eigenvalue near the hard edge, has revealed an excellent agreement with the RMT predictions for nine universality classes occurring for different choices of N and q in the Hamiltonian. Six additional universality classes were identified for the non-Hermitian version of the Wishart-SYK model. In particular, it realizes the class AIII † , the only non-Hermitian class without reality conditions not realized by the nHSYK model, thus completing our program of finding explicit realizations of the tenfold classification within the non-Hermitian SYK model.
Interestingly, and contrary to the Hermitian case, the nHWSYK model also realizes classes beyond the tenfold way, illustrating the richer classification of non-Hermitian random matrix models. More concretely, it also realizes the five pseudo-Hermitian classes AIII − , BDI † −+ , CII † −+ , BDI +− , and CII +− . Finally, the complex-fermion nHSYK model realizes non-pseudo-Hermitian classes with reality conditions, namely the Ginibre Orthogonal, Ginibre Symplectic, chiral Ginibre Orthogonal, and chiral Ginibre Symplectic Ensembles.
Additionally, we mention the possibility of realizing further symmetry classes in two-site nHSYK models. These models, relevant for wormhole physics [18][19][20]85] and also in the exploration of dominant off-diagonal replica configurations [85], provide an interesting playground to implement not only the universality classes already realized in the nHSYK but also going beyond by tuning the additional parameters appropriately. However, a potential issue in this case may be that the coupling needed to lift spectral degeneracies actually determines the RMT symmetry class and it may occur that no new classes are obtained without fine-tuning.
The study of local level statistics presented in this paper will be complemented by a companion work [84] on long-range correlations that explore shorter timescales. Deviations from random matrix universality allow us to use the nHSYK model as a toy model of both generic features of the universal quantum ergodic state reached around the Heisenberg time and non-universal, but still rather generic, properties of quantum interacting systems in its approach to ergodicity.
Other natural extensions of this work include the analytical calculation of the spectral density of the nHSYK model by combinatorial techniques and a short-time characterization of non-Hermitian many-body quantum chaos in this nHSYK model by the evaluation of the Lyapunov exponent resulting from an out-of-time-order correlation function [10,12,90] for times of the order of the Ehrenfest time. The latter could help dynamically characterize non-Hermitian quantum chaos. For real spectra, we have the BGS conjecture that relates dynamical and spectral correlations. However, despite the heavy use of terminology borrowed from the real-spectrum case, it is still unclear to what extent agreement with non-Hermitian random matrix predictions is related to quantum chaos in the original sense of quantum dynamics of classically chaotic systems. We plan to address some of these issues in the near future.