Universality and its limits in non-Hermitian many-body quantum chaos using the Sachdev-Ye-Kitaev model

Spectral rigidity in Hermitian quantum chaotic systems signals the presence of dynamical universal features at timescales that can be much shorter than the Heisenberg time. We study the analog of this timescale in many-body non-Hermitian quantum chaos by a detailed analysis of long-range spectral correlators. For that purpose, we investigate the number variance and the spectral form factor of a non-Hermitian $q$-body Sachdev-Ye-Kitaev (nHSYK) model, which describes $N$ fermions in zero spatial dimensions. After an analytical and numerical analysis of these spectral observables for non-Hermitian random matrices, and a careful unfolding, we find good agreement with the nHSYK model for $q>2$ starting at a timescale that decreases sharply with $q$. The source of deviation from universality, identified analytically, is ensemble fluctuations not related to the quantum dynamics. For fixed $q$ and large enough $N$, these fluctuations become dominant up until after the Heisenberg time, so that the spectral form factor is no longer useful for the study of quantum chaos. In all cases, our results point to a weakened or vanishing spectral rigidity that effectively delays the observation of full quantum ergodicity. We also show that the number variance displays nonstationary spectral correlations for both the nHSYK model and random matrices. This nonstationarity, also not related to the quantum dynamics, points to intrinsic limitations of these observables to describe the quantum chaotic motion. On the other hand, we introduce the local spectral form factor, which is shown to be stationary and not affected by collective fluctuations, and propose it as an effective diagnostic of non-Hermitian quantum chaos. For $q = 2$, we find saturation to Poisson statistics at a timescale of $\log D$, compared to a scale of $\sqrt D$ for $ q>2$, with $D $ the total number of states.

ical universal features at timescales that can be much shorter than the Heisenberg time.
We study the analog of this timescale in many-body non-Hermitian quantum chaos by a detailed analysis of long-range spectral correlators. For that purpose, we investigate the number variance and the spectral form factor of a non-Hermitian q-body Sachdev-Ye-Kitaev (nHSYK) model, which describes N fermions in zero spatial dimensions. After an analytical and numerical analysis of these spectral observables for non-Hermitian random matrices, and a careful unfolding, we find good agreement with the nHSYK model for q > 2 starting at a timescale that decreases sharply with q. The source of deviation from universality, identified analytically, is ensemble fluctuations not related to the quantum dynamics. For fixed q and large enough N , these fluctuations become dominant up until after the Heisenberg time, so that the spectral form factor is no longer useful for the study of quantum chaos. In all cases, our results point to a weakened or vanishing spectral rigidity that effectively delays the observation of full quantum ergodicity. We also show that the number variance displays nonstationary spectral correlations for both the nHSYK model and random matrices. This nonstationarity, also not related to the quantum dynamics, points to intrinsic limitations of these observables to describe the quantum chaotic motion. On the other hand, we introduce the local spectral form factor, which is shown to be stationary and not affected by collective fluctuations, and propose it as an effective diagnostic of non-Hermitian quantum chaos. For q = 2, we find saturation to Poisson statistics at a timescale of log D, compared to a scale of √ D for q > 2, with D the total number of states.

I. INTRODUCTION
The study of quantum chaotic dynamics attracts a great deal of interest in different fields because of its robust universal features. For sufficiently long timescales, the evolution of very different quantum systems is qualitatively similar provided that the dynamics is quantum chaotic.
By contrast, the quantum dynamics of integrable systems is very sensitive to the details of the Hamiltonian. A central result in the theory of quantum chaos is the Bohigas-Giannoni-Schmit (BGS) conjecture [1] that states that spectral correlations of a quantum chaotic system are given by random matrix theory (RMT). In the context of single-body quantum mechanics, the conjecture has received strong analytic support [2-4] by using periodic orbit theory techniques. The BGS conjecture has been very influential because the spectrum of the Hamiltonian (or another relevant operator) is one of the least expensive quantities to obtain numerically even for quantum manybody systems. Therefore, a relatively straightforward spectral analysis is sufficient to determine the quantum chaotic nature of the motion for sufficiently long times of the order of the Heisenberg time-the timescale related to the (inverse) mean level spacing.
However, the agreement with RMT extends in many cases to substantially shorter times due to so-called spectral rigidity-directly related to the power-law tails of the two-level correlation function-which is responsible for the slow logarithmic growth of the number variance or the ramp of the spectral form factor. For sufficiently short times, level statistics of realistic quantum chaotic Hamiltonians deviate from the random matrix prediction. The timescale that marks the onset of these deviations and delimits the region of universal quantum chaotic dynamics-related to the so-called dip or correlation hole [5][6][7][8] of the connected spectral form factor, or to power-law deviations of the logarithmic growth of the number variance-depends on details of the dynamics.
For disordered systems, where it is called the Thouless time [9], it is related to the typical diffusion time needed for a single particle to cross the sample. However, also in the context of disordered systems, this timescale is sometimes determined by ensemble fluctuations not directly related to the type of motion [10,11].
So far, the discussion has been restricted to Hermitian quantum systems. A natural question to ask is to what extent these ideas and results are applicable to non-Hermitian quantum chaotic many-body Hamiltonians. The main goal of this paper is to address this question. We investigate long-range spectral correlations such as the number variance and the spectral form factor, which probe shorter timescales of the dynamical evolution of non-Hermitian systems.
The theory of non-Hermitian random matrices is well developed for some universality classes corresponding to the so-called Ginibre ensembles [12]. However, no equivalent of the BGS conjecture is known, and so the relation between dynamics and level statistics is less clear than in the Hermitian case. There are also technical problems: correlations of complex eigenvalues are weakened, and the necessary unfolding of eigenvalues may be problematic [13,14] when the eigenvalue distribution is not radially symmetric. Unfolding problems have been ameliorated in the last years for shortrange observables with the introduction of spectral observables such as the adjacent gap ratios [15] for complex spectra [16] that do not require unfolding. They have already been applied in a variety of non-Hermitian systems: phase transitions in many-body Liouvillians [17][18][19][20], non-Hermitian Anderson localization [21][22][23][24], nonunitary open quantum circuits [25,26], two-color QCD at imaginary chiral chemical potential [27], and, more recently, the Sachdev-Ye-Kitaev (SYK) model [28][29][30].
Long-range spectral correlators such as the number variance [31][32][33][34][35][36] or spectral form factor (SFF) [32,37] have already been investigated in the context of non-Hermitian systems but there are still problems with the unfolding procedure, which is necessary, especially for the SFF, for a correct dynamical interpretation of the results. Moreover, the role of spectral rigidity, if any, and the determination of the timescale that signals deviations from universality are still poorly understood in quantum chaotic non-Hermitian systems. We aim to shed light on this problem by computing these spectral observables for the non-Hermitian SYK model [29,[38][39][40]. This model is a natural building block of Euclidean [29,38] and Keldysh [41] wormholes, models for Lindbladian dissipation [30,42], and entanglement dynamics [43].
The SYK model, describing N fermions with infinite-range interactions in zero spatial dimensions, was introduced more than 50 years ago [44][45][46][47][48] in the context of nuclear physics as a toy model for nuclei. Later, it played an important role in the development of so-called many-body quantum chaos [49][50][51][52] and also in the description of certain aspects of spin liquids [53]. The revival of interest in the SYK model is motivated by its role in quantum gravity as a toy model for holography [54][55][56] and also due to the use of Majorana fermions, proposed by Kitaev [54], that simplifies the model allowing analytical calculations for some region of the parameters. For instance, it was possible [54,55] to demonstrate that the SYK model saturates a universal bound [57] on the exponential growth in time of certain out-of-time correlation functions that probe quantum chaos at short timescales of the order of the Ehrenfest time-the time for which quantum effects start to become relevant. It has also been shown that the SYK model, both Hermitian and non-Hermitian, is quantum chaotic with spectral correlations well described by RMT [28,[58][59][60]. By tuning q and N , the SYK model can also reproduce several of the different universality classes, controlled by the global symmetries of the system, in which a many-body quantum chaotic system can relax to ergodicity [28,[61][62][63][64][65][66]. It is therefore a natural choice for the problems we will be addressing.
We initiate our analysis in Sec. II with a description of long-range spectral correlations of non-Hermitian random matrices belonging to the Ginibre ensemble [12]. We will focus on the number variance and the spectral form factor of the real parts of the eigenvalues, which has recently been proposed as a measure of quantum chaos in non-Hermitian systems [37]. Some of the analytic results of this section were already derived in Refs. [32,37,67,68], but we present them here from a unified viewpoint. Then, in Sec. III we apply the same tools to investigate the emergence of random-matrix universality, and its limits, in the non-Hermitian SYK model.

RANDOM MATRIX MODELS
For real spectra, some widely used long-range correlators are the number variance, the ∆ 3 statistic, and the spectral form factor [69]. In this section, we study analog statistics for non-Hermitian spectra of the Ginibre unitary ensemble (GinUE) of random matrices [12].
A. The Ginibre ensemble Throughout this paper, we frequently refer to the Ginibre ensemble. Here, we collect some wellknown results that we use below. The Ginibre ensemble is the ensemble of D × D complex random matrices H with probability distribution given by [12] With this normalization, in the limit of large D, the eigenvalues z i of H are distributed uniformly inside the complex unit disk. The real parts E i = Rez i of the eigenvalues, therefore, follow the semicircular distributionρ normalized as At finite D, the connected two-point correlation function of the eigenvalues is given by [12,69] with the kernel given by: In the large-D limit, this simplifies to [12,69]

B. Unfolding
Depending on the long-range correlator, we may have to unfold the eigenvalues (i.e., reparametrize them such that the spectral density is constant in the new variables) in order to obtain universal results. Since the number variance is defined as the variance of the number of eigenvalues in a fixed interval, there is no need for unfolding, although it is convenient to do so if it is calculated by spectral averaging or a combination of spectral averaging and ensemble averaging. The spectral form factor is an observable that involves the entire spectrum so in this case unfolding is essential for making quantitative comparisons between different systems. Since the contribution of eigenvalue pairs with large spacings is suppressed by large phase oscillations, the main contribution to the spectral form factor is due to eigenvalues that are close. Therefore, it is possible to give a local definition of the spectral form factor that only includes the eigenvalues on a scale where the average spectral density is approximately constant, see Sec. II E below. Because our goal is to identify universal features of the quantum dynamics, we only consider connected two-point correlators in the analysis of the spectral correlations.
For the Ginibre ensemble, there is no need to unfold the spectrum because, well away from the spectral edge, the spectral density is constant and a rescaling is enough. However, for the nHSYK model, to be studied later, the average spectral density is not constant and nonuniversal as it is determined by the details of the phase space and the dynamics. Therefore, we first provide a detailed explanation of the unfolding of complex spectra. For simplicity, we focus on the case of radially symmetric spectra.
Because of unfolding ambiguities, we only unfold the eigenvalues in this case 1 and we only have to reparametrize the absolute value of the eigenvalues. The unfolding is performed using the average radial spectral density, which is a smooth function and therefore does not affect local statistics. If the average radial density is given byρ(r), the unfolded eigenvalues, z unf k in terms of the original eigenvalues are given by As a check of this transformation, we can take the flat densityρ(r) = 1/π which results in z unf The unfolded eigenvalues have constant density inside the unit disk. Below, we will see that the spectral density of the nHSYK model for q ≥ 6 is almost constant, and unfolding can be approximated by just rescaling the eigenvalues to the unit disk by When the eigenvalue density is not strictly isotropic, generically the spectrum is still locally isotropic. In that case, unfolding ambiguities can be resolved by requiring that local spectral isotropy is preserved.

Complex eigenvalues
The best-known spectral correlator to probe the existence of spectral rigidity is the number variance, defined as the variance of the number of eigenvalues in a spectral window containing a fixed number of eigenvalues on average. Let that window be a compact domain D inside the unit disk, containing n eigenvalues. In terms of the connected two-point correlation function, the 1 There is still an ambiguity. For example, we could have unfolded the eigenvalues z k so that the density of Re(z k e iθ ) (i.e., the projection of the eigenvalues along some axis with angle θ) becomes constant. However, since analytical results are only available for constant density inside the unit disk we unfold the eigenvalues this way.
number variance reads For spectra of Hermitian quantum chaotic systems and Hermitian random matrices, the growth of the number variance increases only logarithmically with the number of eigenvalues, while it grows linearly in the case of a nondegenerate uncorrelated spectrum typical of an integrable system. For non-Hermitian systems, the spectral rigidity takes a different form because the spectrum is now two-dimensional. One can imagine the D eigenvalues in the unit disk as small disks with area in Poisson statistics with Σ 2 (n) = n. It is clear that the coefficient of √ n depends on the length of the perimeter of D for n eigenvalue inside this domain. Below, we will see that for a disk [35] we have that Σ 2 (n) = n/π, while for a square we find Σ 2 (n) = 2 √ n/π, which behaves as the ratio of the perimeters for the same area.
Next, we calculate the number variance of the Ginibre ensemble for two simple geometries, the rectangle and the disk. Using Eq. (6), the number variance (10) for the rectangle with sides 2a and 2b and centered at the origin is given by with The second term in Eq. (11) factorizes into F (a)F (b) with where erf stands for the error function. Using Eq. (12), this result can be expressed in terms of n and the aspect ratio of the rectangle, α ≡ b/a, which has a well-defined large-D limit. F (b) is obtained by replacing α → 1/α. If the rectangle contains more than a few eigenvalues, this result is well approximated by This gives the number variance This approximation is already accurate to four digits for n = 2 and becomes rapidly more accurate for larger values of n.
A second interesting case is a rectangle with b independent of n and a ∼ b. We again have that 4abD/π = n, but now 4a 2 D = n 2 π 2 /(4Db 2 ) instead of nπ/α. In this case, If b is fixed we have in the large-D limit that The second term is subleading and will be neglected. This gives the number variance Note that this result for the number variance is nonuniversal: it depends on D and b. A stable double scaling limit is obtained by taking the large-D limit at fixed b √ D.
We now turn to the evaluation of the number variance (10) for a disk of radius R centered at zero. This can be done most simply by using the exact finite-D expression for the connected twopoint correlation function of the Ginibre ensemble written in terms of the kernel, Eq. (4). Inserting it into Eq. (10) we obtain with the kernel given by Eq. (5). Using polar coordinates, the integrals can be simplified to which was first obtained in Ref. [35]. In the large-D limit, we have which allows us to take the large-D limit of Eq. (21). Using the asymptotic expansion of the incomplete Gamma function, the number variance for finite but large n is approximated by

Real parts of complex eigenvalues
In this paper, we also consider level correlations of the projections of the eigenvalues on a line through the origin. If the spectrum is radially symmetric, these are just the level correlations of where we expect a rounding to occur when n ≈ 2 √ D.
The number variance of the real parts of the eigenvalues (or the eigenvalues projected on a ray through the origin) is not universal. It depends sensitively on the geometry of the eigenvalue support. For example, for the elliptic Ginibre ensemble supported on an ellipse with a x and a y as long and short semiaxis (a x a y = 1), order 2a y √ N eigenvalues can diffuse in and out of a segment. More generally, the number variance of the real parts depends on the vertical dimension of the eigenvalue support. The geometry of the eigenvalue support thus determines the fraction of eigenvalue pairs with projections that are much closer than their spacing in the complex plane and are, therefore, essentially uncorrelated. We refer to this as the Poisson admixture. For the same reason, we expect that the number variance of the projected eigenvalues is not stationary. Closer to the edge, there is less Poisson admixture resulting in a smaller value of the number variance.
The number variance of the real parts of the eigenvalues at energy E can be obtained from the number variance in a rectangle with fixed side b, Eq. (19). In this case, b is equal to the range of the imaginary parts, i.e., b = √ 1 − E 2 = πρ(E)/2D withρ(E) the density of the real parts of the This results in the number variance showing that the number variance is determined by a universal function g.
We shall employ these expressions in the comparison with the number variance of the nHSYK model in Sec. III.

D. Spectral form factor
For real spectra, the Fourier transform of the two-point correlation function, called the spectral form factor, although non-self-averaging [70], is also a popular probe of quantum chaos. The logarithmic growth in energy of the number variance, which signals spectral rigidity, translates into linear growth in time of the spectral form factor for intermediate times much larger than the Ehrenfest time but smaller than the Heisenberg time.
The connected spectral form factor of D (real) eigenvalues x k is defined by The normalization is such that The second term in Eq. (27) is proportional to the nonuniversal disconnected spectral form factor: Only K c (t) provides direct information on the quantum dynamics.
A spectral form factor for complex eigenvalues was first introduced in Ref. [32], with z 1 = x 1 + iy 1 and z 2 = x 2 + iy 2 . For s 1 = s 2 = 0, this becomes the spectral form factor of the real parts of the eigenvalues. In Ref. [32], the spectral form factor was calculated for the weak non-Hermiticity limit of the elliptic Ginibre ensemble but, as stated by the authors, their calculations are also valid in the case of interest here, termed the strong non-Hermiticity limit. They found: in the normalization where the support of the D eigenvalues is the unit disk. A derivation of this result from the finite-D Ginibre kernel is given in Appendix A, see also Ref. [37,71]. The spectral density of the real parts of the eigenvalues of the Ginibre ensemble is given by the semicircle distribution, Eq.
(2). We note that unfolding the eigenvalues from semicircular to constant density can almost entirely be absorbed by rescaling the argument of the exponential in Eq. (31) by 1.18.
More generally, setting t = τ cos θ, s 1 = −iτ sin θ, and s 2 = iτ sin θ in Eq. (30), we can define the spectral form factor of the eigenvalues projected onto the direction defined by the angle θ: If the spectral properties are axially symmetric, this spectral form factor does not depend on the angle θ onto which the eigenvalues are projected. By averaging over θ it is possible to increase the statistics of the form factor [37]. The spectral form factor of the projected eigenvalues was recently proposed [37] as a measure of quantum chaos in dissipative systems, see also Refs. [71,72].
It was dubbed the dissipative spectral form factor [37] not to be confused with other closely related quantities: the dissipative form factor introduced in Ref. [73] and the open-system spectral form factor put forward in Refs. [74,75].
The spectral form factor (32) measures the same long-range spectral correlations as the number variance of the real parts and therefore depends sensitively on the geometry of the spectrum. This can be seen explicitly from the relation between the number variance and the spectral form factor [5] which is valid if the spectral form factor is calculated for unfolded eigenvalues. Equation (33) shows that the number variance and the spectral form factor are complementary observables. In the normalization with D eigenvalues in the complex unit circle, the spectral form factor for τ > T determines the number variance for n < D/(2T ). Likewise, the spectral form factor for τ < T mostly contributes to the number variance for n > D/(2T ). This does not imply, however, that the asymptotic behavior of the number variance is given by the small-τ behavior of the spectral form factor. In particular, the spectral form factor for τ √ D increases quadratically This region contributes to the number variance for n > √ N /2 but it does not determine the large-n saturation value of the number variance.
Finally, we note that the limits D → ∞ and τ → 0 do not commute (see Appendix A for details). Indeed, the coefficient of τ 2 in the expansion of obtained by taking the large-D limit first, is not equal to the coefficient obtained by taking the small-τ limit first. The latter is given by K c,exact (0)/2, where K c,exact (τ ) is the exact finite-D result for the Ginibre ensemble, Eq. (A8). Its Taylor expansion to second order in τ is given by . It is actually a factor of 2 larger than the result obtained by taking the large-D limit first and is equal to the perturbative coefficient which can be easily checked numerically, for example, for an ensemble of 1000 of 100 × 100 Ginibre matrices. For a real spectrum belonging to the Gaussian Unitary Ensemble (GUE), we also find that for sufficiently small τ , the spectral form factor K(τ ) ∼ τ 2 with the perturbative prefactor given by the analogous expression.
E. Local spectral form factor

Real eigenvalues
Let us again start by considering the case of real spectra. Except for times much shorter than the inverse mean level spacing, the main contribution to the spectral form factor (27) comes from eigenvalue pairs that are sufficiently close. That is, we can define a local connected spectral form where ∆(x) is the average spacing of the eigenvalues at x, the level densityρ(x) = 1/∆(x), and |π(x)| is the length of an interval located at x that satisfies ∆(x) |π(x)| D∆(x). We have, for now, chosen a hard cutoff to enforce locality. Other choices are possible, however, and below we will find a Gaussian cutoff more useful.
Locally, unfolding is just rescaling the eigenvalues by the local level spacing ∆(x). We thus have Integrating over the entire spectrum gives We can also invert the relation (37): Using Eq. (33), the spectral form factor is related to the number variance by [5] where Σ 2 (x, n) is the local number variance. Using Eq. (37), it evaluates to The integral over t can be evaluated analytically resulting in the number variance This shows that the local number variance can be obtained by integrating the local spectral form factor. The spectral form factor of the entire spectrum is thus related to the spectral average of the number variance. This is not an issue if the number variance is stationary (i.e., independent of the point x) but, as we will see below, the number variance of the real parts of the eigenvalues of both the Ginibre ensemble and the nHSYK model are not stationary. On the other hand, the local spectral form factor turns out to be stationary.
The local spectral form factor (36) shows strong oscillations resulting from the Fourier modes of the hard cutoff. These oscillations can be eliminated by introducing a smooth cutoff. For a Gaussian cutoff, we obtain the local spectral form factor [76] K loc where N is a normalization factor chosen such that K loc c (x, t) asymptotes to 1 for large t and w is the width of the cutoff. The large-t behavior determined by the contribution of the self-correlations is given by This results in the normalization factor

Real parts of complex eigenvalues
So far, we have only considered the local spectral form factor of real eigenvalues. We now turn to the local spectral form factor of the real parts of complex eigenvalues, i.e., the local counterpart of Eq. (30). A straightforward generalization of Eq. (43) yields where z j = x j + iy j .
We now show that this local spectral form factor for eigenvalues unfolded to constant density inside the unit disk is stationary. The non-Hermitian two-point correlator is given by where R univ is the universal two-point correlator and z c = (z 1 + z 2 )/2. Changing to variables the universal contribution to the spectral form factor in Eq. (46) can be written as For eigenvalues unfolded to the complex unit disk, this becomes For the Ginibre universality class with universal two-point correlator (6), the second integral can be worked out: Collecting all terms and using that the normalization factor is equal to [see Eq. (45)] We find that, for eigenvalues unfolded to a constant density inside the unit disk and in the limit w 2 1/D, the local spectral form factor of the real parts of the eigenvalues is given by: This shows that the spectral form factor is stationary. In contrast, the global spectral form factor, given by Eq. (38), is the integral of a nonstationary quantity due to the semicircular distribution of the projected eigenvalues multiplying the local spectral form factor.

Relation to the local number variance
It is possible to calculate the number variance of the real parts of the eigenvalues from the spectral form factor using the relation (33) which is valid for unfolded eigenvalues [32]. If the local eigenvalue density isρ(E) (normalized to ρ(E)dE = D), the spectral form factor unfolded to the unit density is given by [see Eq. (37)] The resulting integral can be performed analytically [32] leading to This is in agreement with the previous result (25) obtained from the number variance of a rectangular geometry with vertical side b = √ 1 − E 2 . If we take the limit D → ∞ at fixed n, we get Poisson statistics: However, the expression for the number variance also has a nontrivial double scaling limit It is tempting to interpret the existence of this scaling limit as a signature of quantum chaos in non-Hermitian systems. However, a similar scaling behavior has been observed for the number variance of integrable systems at a finite distance above the ground state [2, [77][78][79].
As a check of the analytical result (55), we compare in Fig. 1 this expression to the numerically calculated number variance at the center of the spectrum for an ensemble of 8192 realizations of 1024×1024 complex Ginibre matrices. The normalization is such that the support of the eigenvalues is the complex unit disk but this does not affect the number variance. The discrepancy between the analytical and numerical results for large D is due to finite-size corrections. There is also a correction due to the semicircular shape of the spectral density, but because n D, this correction is much smaller and can be neglected. The spectral form factor for the other two universal bulk statistics, AI † and AII † , is much less understood and no analytical results are available. Some numerical results were presented in Ref. [72] for class AI † , while the spectral form factor for class AII † has not been investigated before.
We obtained them numerically and plot them in Fig. 2. As for class A, for classes AI † and AII † there is also an early quadratic growth, albeit with a different prefactor-for class AI † it is twice the prefactor of A, while for class AII † it is half. Note that the prefactor for class A, 1/(2D), is in agreement with our previous considerations, see Eq. (34). As can be seen from Fig. 2 (right), in all three cases, the prefactor of the τ 2 dependence decreases by a factor of 2 going from the perturbative to the nonperturbative domain [see the discussion below Eq. (34) for an explanation of this anomaly for class A]. The approach to the late-time plateau is slower for AI † than for A, while it is faster for AII † . Contrary to the Hermitian case, there is no nonanalyticity in the spectral form factor around the Heisenberg time.
The spectral form factor is a bulk observable, namely, it is defined as a sum over all eigenvalues, so in principle it can identify only three different classes A, AI † , and AII † . However, we will show in the next section, that for ensembles with chiral symmetry, the number variance calculated for a symmetric interval around zero is a factor of 2 different with respect to the Ginibre ensembles, so that it can be employed to identify non-Hermitian systems with chiral symmetry.

G. Spectral form factor for Poisson statistics
In this section, we calculate the spectral form factor for spectra with Poisson statistics (i.e., 2d uncorrelated points) unfolded to constant density inside the complex unit disk. The connected two-point correlation function is given by The connected spectral form factor is given by In Fig. 3, we show a numerical verification of the result of Eq. (59), finding perfect agreement.
We will see in the next section that although the q = 2 nHSYK model is integrable, its spectral form factor has more structure than plain Poisson statistics of completely uncorrelated random Table I. Classification of the nHSYK Hamiltonian into non-Hermitian bulk universality classes for all q and even N . Note that this does not correspond to the full symmetry classification, which is richer and goes beyond bulk level statistics [28].

MODEL WITH COMPLEX COUPLINGS
We now probe the dynamics of the non-Hermitian SYK (nHSYK) model for timescales shorter than the Heisenberg time by a detailed comparison of the unfolded spectral form factor and the number variance with the random matrix predictions worked out in the previous section. The nHSYK Hamiltonian is defined as [54,55] where N and q are integers (N is taken to be even), J i 1 ···iq and M i 1 ···iq are real Gaussian random variables with zero mean and variance and ψ i are Majorana fermions satisfying {ψ i , ψ j } = 2δ ij . To be precise, we note that for odd q, H corresponds to a supercharge (not Hamiltonian) operator.
The symmetry classification of the nHSYK model was put forward in Ref. [28]. It was found that, depending on q mod 4 and N mod 8, it belongs to nine out of the 38 non-Hermitian symmetry classes. However, the bulk correlators we are employing here only capture the local level repulsion, i.e., only distinguish the bulk universality classes A, AI † , and AII † . The bulk universality classes for different q and N [27,28] are tabulated in Table I.
We obtain the spectrum by exact diagonalization techniques. We carry out an ensemble average to suppress statistical fluctuations, reaching at least 10 5 eigenvalues for a given q and N . Since the spectrum is radially symmetric, the necessary unfolding is carried out as explained in Sec. II B.
Depending on q, we shall employ polynomials of different degrees to approximateρ(r). For q = 4, the radial spectral density is well approximated by a fourth-order even polynomial, while for q = 2 it is close to a Gaussian, and for q = 6 the radial spectral density is almost constant so that unfolding is basically a rescaling of the eigenvalues. For q = 3, the radial spectral density can be unfolded by an eighth-order even polynomial.
In what follows, we will compare the spectral form factor and number variance of the nHSYK model with the random matrix prediction in the corresponding universality class (Table I). We start our analysis with the spectral form factor.

A. Spectral form factor of the nHSYK model
We carry out the numerical evaluation of the ensemble-averaged spectral form factor corresponding to the real parts of the eigenvalues, K c (τ, θ), for various values of N and q. As before, the spectral form factor is normalized such that it asymptotes to 1 for large times. Since the spectrum is rotationally invariant, we additionally average the spectral form factor over 10 values of θ j = πj/5, j = 1, . . . , 10.
1. Ginibre universality for q = 4 In Fig. 4, we depict the results for the q = 4 nHSYK Hamiltonian for N = 22 and N = 26 (black curves) and compare them to the analytical result (31) for the Ginibre ensemble (red curves). We find agreement with the random matrix prediction for τ > √ D but the two results differ for smaller times. This is fully consistent with the results of short-range correlations [28] which are insensitive to these deviations. In Fig. 5, we show the results for N = 24 and N = 28. The agreement with the corresponding random matrix universal result is excellent for τ > √ D. In these two cases, which are in the universality class of AI † (N = 24) and AII † (N = 28), no analytical formula is available and the random matrix result was obtained numerically for D = 2048. The spectral form factor for other values of D can be obtained by using the scaling relation

The limits of universality: Collective scale fluctuations
To better understand the short-time behavior of the spectral form factor we have enlarged the region close to the origin in the plots of the right column of Figs. 4 and 5. The local minimum of K c (τ ) for τ > 0, usually termed correlation hole [5][6][7][8]81], defines, for a real spectrum, the maximum timescale for which the dynamics did not fully relax to the universal prediction of RMT.
In the Hermitian SYK model [62], it is determined by the collective fluctuations of the spectrum that arise because the number of independent matrix elements (∼ N q ) is much smaller than the number of matrix elements of the Hamiltonian (2 N/2 ). The same mechanism is at work in the non-Hermitian case, where we have the same mismatch in the number of matrix elements.
The oscillations for small times are mostly due to collective scale fluctuations [82]. They correspond to fluctuations in the overall scale of eigenvalues from one realization to the next: x n → x n (1 + ξ), for all n, where x n are the real parts of the eigenvalues. ξ is a random variable  Figure 5. The connected part of the spectral form factor of the real parts of the unfolded eigenvalues, Eq. (32), normalized by the number of eigenvalues D ∼ 2 N/2−1 , of the q = 4 nHSYK model for N = 24 (upper) and N = 28 (lower). The red curves are the random matrix prediction which was obtained numerically for the AI † and AII † universality classes. The right plots are a magnification of the left ones. We find excellent quantitative agreement without any fitting up to relatively small τ ∼ √ D, of the order of the correlation hole, which signals the timescale for which the quantum chaotic dynamics is universal.
with zero mean and gives rise to the scale fluctuation of the spectral density: whereρ is the ensemble-averaged spectral density. The connected two-point correlator for these scale fluctuations is given by where the prime denotes the derivative, we have used ξ = 0, and we have dropped all terms of order ξ 4 and above. Recall thatρ(x) is normalized as ρ(x)dx = D. This contributes to the spectral form factor as [83] δK The spectral density of the real parts of the eigenvalues is given bȳ where J 2 is a Bessel function, resulting in the contribution of the scale fluctuations to the spectral form factor which decreases as 1/τ for large τ . The analytical result for total spectral form factor in class A including the scale fluctuations factor is then given by The variance ξ 2 can be computed as with the moments defined as Here, E k are either the real eigenvalues or the real parts of the complex eigenvalues and the brackets · denote ensemble averaging. For the Hermitian SYK model, the moments can be evaluated exactly and we find that [82] For the real parts of the eigenvalues of a non-Hermitian matrix with spectral density unfolded to constant density inside the unit disk, ξ 2 cannot be obtained from traces of moments of the Hamiltonian, but its exact numerical value can be obtained from the real parts of the eigenvalues using the definition (70). For our non-Hermitian SYK model, we find that it is approximately equal In Fig. 6, we show the difference between the spectral form factor of the real parts of the eigen-  Fig. 7, we compare the full analytical spectral form factor, Eq. (69), with numerical results. Given that higher multipole collective fluctuations also contribute to the difference, the agreement with the analytical result is better than expected, in particular for small times. We thus conclude that most of the oscillatory behavior is due to the lowest-order multipole, i.e., the scale fluctuations. We emphasize that the results for δK c (τ ) are obtained without using fitting parameters. Note that the period of the oscillations does not depend on N and is close to the period of the oscillations of J 2 2 (τ ). The amplitude increases with D and also varies as the amplitude of J 2 2 (τ ). The location of the correlation hole can be obtained by equating the two contributions to the spectral form factor: the universal Ginibre contribution, Eq. (31), and the collective fluctuations contribution, Eq. (68). By replacing the oscillatory part of J 2 2 by its asymptotic average, the location of the correlation hole is thus given by the minimum of This condition cannot be solved analytically, but it gives the rough position of the correlation hole and can be studied numerically for small values of N . For the cases studied in this paper belonging to class A, we compare in Table II the position of the correlation hole obtained from the figures for the spectral form factor with the result given by the minimum of (74) and find good agreement between the two.
The condition (74) can be recast as This result is to be contrasted with the Hermitian case, for which the location of the correlation hole is roughly determined by the condition, Although ξ 2 is approximately the same as before, the Heisenberg time is now of order D (instead Finally, we note that it is possible to eliminate the collective spectral fluctuations by unfolding the spectrum realization by realization [76,82]. Then these oscillations do not show up in the spectral form factor.

Dependence of nonuniversal features on q
Results for the q = 3 and q = 6 nHSYK Hamiltonian, see Figs. 8 and 9, confirm the picture obtained for q = 4. Agreement with the random matrix predictions corresponding to the expected universality class is observed for τ > √ D. The area below the small-time peak decreases markedly for increasing values of q. This is expected since a larger q > 2 brings the nHSYK Hamiltonian closer to a random matrix, as more entries of the Hamiltonian are nonzero. The area below the peak is proportional to 2 N/2 / N q . For N = 24, it is given by 14.84, 2.02, 0.39, and 0.03 for q = 2, 3, 4, and 6, respectively.
The oscillatory behavior in the small-τ region for q = 3, although not qualitatively different from q = 4, has a much larger amplitude than in the q = 4 case (see the right panel of Fig. 8).
This results in a correlation hole that is shifted to a larger value of τ . On the other hand, for q = 6, the amplitude of the oscillations is very small, and we barely observe any deviation from the For the nHSYK model, this corresponds to a constant level density inside the eigenvalue disk so that the real parts of the eigenvalues are distributed according to a semicircle. For q = 6, we are likely in this asymptotic region. Indeed, this is confirmed by a comparison of the disconnected part of the spectral form factor, see Fig. 10, with the random matrix prediction K dis (τ ) = 4D 2 J 2 1 (τ )/τ 2 , the square of the Fourier transform of the semicircle law, where J 1 is a Bessel function. They are almost indistinguishable which explains why, for q = 6, the spectral density is very close to that of the Ginibre ensemble. This also suggests that the spectral correlations are very close to that of the Ginibre ensemble. In contrast, for q = 3, we observe larger deviations with respect to the semicircle law in the disconnected part. This is consistent with the fact that, by reducing q, the Hamiltonian is much sparser and, therefore, deviations from the RMT predictions should be more visible. Another issue is that there is a systematic difference between even q and odd q related to  Figure 9. The connected part of the spectral form factor of the real parts of the unfolded eigenvalues, Eq. (32), normalized by the number of eigenvalues D ∼ 2 N/2−1 , of the q = 6 nHSYK model for N = 22 (upper) and N = 24 (lower). We find excellent quantitative agreement with the random matrix prediction for the GinUE. We note that this is a global observable and therefore the RMT prediction for class D is indeed identical to class A because the two only differ for eigenvalues around E = 0. A remarkable feature of the q = 6 results compared to smaller values of q is that the correlation hole has almost disappeared.
cancellations that occur in the calculation of moments of eigenvalues of the supercharge [85] which we expect to persist in the non-Hermitian case.

Integrable behavior for q = 2
The q = 2 SYK and nHSYK models are both integrable with all energy levels determined by N/2 single-particle energies. In this case, we expect Poisson level statistics for sufficiently long times, but deviations from Poisson statistics may be observed for shorter times. Indeed, as illustrated in Fig. 11, the spectral form factor saturates to the Poisson limit, K c (τ ) = 1, at a scale of order log D, which is much shorter than for q > 2, where the scale is determined by √ D (for N = 24 the two scales are of the same order of magnitude and our data cannot really distinguish between the two). The analytical result for uncorrelated eigenvalues unfolded to constant density inside the complex unit disk, given by Eq. (59), saturates to Poisson statistics at τ = O(1) and does not match the numerical result (see the solid red curve in Fig. 11, left). A reasonable fit is obtained Although this is a nonuniversal observable related to the Fourier transform of the spectral density, for q = 6, we find excellent agreement with the random matrix prediction, K dis (τ ) = 4D 2 J 1 (τ ) 2 /τ 2 after an overall rescaling. However, for q = 3, the agreement is only qualitative. This is not surprising, as the spectral density for q = 3 is not constant while for q = 6 it is already almost constant resulting in a semicircular spectral density of the real parts. We also expect that, for larger N , some deviations will be observed in the latter case because the spectral density is less uniform.  Figure 11. The connected part of the spectral form factor of the real parts of the unfolded eigenvalues, Eq. (32), normalized by the number of eigenvalues D ∼ 2 N/2−1 , of the q = 2 nHSYK model for N = 24. The solid red curves show the result for Poisson statistics, Eq. (59), in the left panel, and the same expression with the horizontal axis rescaled by a factor of log D in the right panel. As was expected, the q = 2 SYK model shows correlations that are in between Poisson statistics and RMT statistics. The reason is that the integrable many-body spectrum is determined by a small number (N/2) of chaotic single-particle energies. The oscillations observed for small, but nonzero, τ are due to collective fluctuations of the spectral density.
by replacing τ → τ / log D (solid red curve in Fig. 11, right), but we have no rigorous argument for this substitution.
Physically, the saturation scale of the q = 2 spectral form factor is related to the fact that the model can be mapped onto free fermions with single-particle energies correlated according to RMT [60]. The short-time dynamics, controlled by the single-particle excitations, will be very different from that expected for a generic integrable system. However, for longer times of the order log D, multiparticle excitations will reveal the generic integrable nature of the quantum dynamics.
As is the case for q = 3 and q = 4, we find oscillations for small values of τ with the same period but with a larger amplitude. These oscillations, which dominate the quadratic τ -dependence, are due to scale fluctuations of the average spectral density.  (55), while in the case of class AII † we have to rely on a numerical calculation of the spectrum of the corresponding random matrix ensemble. We assume that the number variance of AII † still has the scaling behavior (26) obtained analytically for the Ginibre ensemble, whereρ(E) is the eigenvalue density on the real axis, for some universal function f .

q = 4 and nonstationarity
We first discuss the q = 4 case. In Fig. 12, we plot the number variance of the real parts, Σ 2 (n), versus the average number of levels, n, in an interval that is chosen to be symmetric around zero. No unfolding is necessary this way-if we would have unfolded the real part of the eigenvalues (times a phase factor) we would have obtained the same result. In order to suppress statistical fluctuations, we also average {e iθ |z k |} over ten values of θ as we did for the calculation of the spectral form factor.
The results are compared to the analytical result (55) for the Ginibre ensemble (red curves) and the result obtained by integrating the numerical spectral form factor for τ < √ D using Eq. (33) (blue curves) instead of using the analytical result of the spectral form factor all the way to τ = 0. For N = 22, this correction explains the difference between the number variance for the nHSYK model and the Ginibre ensemble, but for N = 26 a discrepancy remains. One issue, as was discussed above, is that by integrating the spectral form factor, we obtain the spectral average of the number variance, while in Fig. 12 we show the number variance for intervals centered about zero energy. The conclusion is that the number variance is not stationary, but this is also the case for N = 22 Remarkably, for E c = 0.6 we find agreement between the nHSYK model and the Ginibre ensemble for the entire range of n we looked at. Since the number variance of the nHSYK model ensemble overshoots the Ginibre number variance for E c > 0.6 and undershoots it for E c < 0.6, its spectral average will agree with the random matrix prediction up to a larger value of n, i.e., to n ≈ 50. This is the spectral average related to the spectral form factor of the real parts of the eigenvalues. We thus expect it will agree with RMT for τ > 4096/100 ≈ 40, which is consistent with Fig. 6. This calculation also shows that the good agreement we find for the spectral form factor for q = 4 is due to the fact that the stationarity behavior of the nHSYK model and the  Although a good agreement has also been observed in other systems [37,71,72], we do not expect that the spectral form factor of the real parts of the eigenvalues is a good observable for detecting universal random matrix behavior in generic quantum dissipative systems. The reason for that is that the statistical properties of the real parts of the eigenvalues depend sensitively on the shape of the spectral density which is not universal. In addition, as seen in the previous section, the contribution due to collective fluctuations may dominate the spectral form factor all the way up to the Heisenberg time.

q = 3 and the effect of chiral symmetry
We turn now to the q = 3 case. In Fig. 14 is due to the chiral symmetry of the AII † − symmetry class which doubles the number variance only for E c = 0. The agreement with the random matrix result decreases rapidly for increasing E c and the number variance is close to Poisson for E c = 0.9.
The same behavior is also visible in the nearest neighbor spacing distribution, P (S) where S is the absolute value of the distance between neighboring eigenvalues in units of the mean level spacing. It is illustrated by the log plot in Fig. 15, where we compare the spacing distribution of the unfolded nHSYK eigenvalues for E k ∈ [0.85, 0.9] (blue curve), E k ∈ [0.9, 0.925] (black curve), and the corresponding random matrix result for AI † (red curve). We find excellent agreement for S < 1.75, but for larger spacings the distribution becomes exponential (characteristic of Poisson statistics) instead of Gaussian (characteristic of random matrix statistics) as we move close to the edge. Note that the average spacing is 0.024 so that there are no edge effects for the interval

C. Local spectral form factor and stationarity
In this and the next section, we study the stationarity of the local spectral form factor and the relation with the number variance in the context of the nHSYK model. In Fig. 17, we show the local spectral form factor of the real parts of the eigenvalues, Eq. (43), for N = 24 and q = 6 (left) and N = 26 and q = 4 (right). The local spectral form factor is calculated atx = 0,x = 0.4 and x = 0.8. The width of the Gaussian cutoff is w = 0.05. As expected, the oscillations for small times are absent. Contrary to the number variance though, we observe only a weak dependence of the local spectral form factor onx. This is in agreement with the stationarity of the local spectral form factor shown in Sec. II E. We conclude that the local spectral form factor overcomes the shortcomings of the global spectral form factor (nonstationarity and a diverging nonuniversal collective-excitation peak) and is, therefore, a good diagnostic of non-Hermitian quantum chaos.
More surprising is the fact that, for q = 6, for which the spectrum has a reflection symmetry across the origin, the spectral form factor in the center of the spectrum is the same as away from the center. The explanation is as follows. The spectral form factor can be written as: Away from the center of the spectrum, the contribution of both terms is equal. However, forx = 0, the term containing the sine functions vanishes because the eigenvalues occur in pairs ±x k . In the cosine term, we can restrict the sum to the positive real parts only at the expense of an overall factor of four. For a finite correlation length, we then only sum over half the total number of eigenvalues so that the cosine term in the reflection-symmetric case is twice as large as without this symmetry.
Since the sine term does not contribute in the reflection-symmetric case, we find that the spectral form factor is the same in both cases. We thus conclude that chiral symmetry does not affect the local spectral form factor.

D. Local number variance
Contrary to the spectral form factor, the number variance for an interval centered about zero is quite different from the number variance in the bulk of the spectrum in the case when the spectrum is reflection symmetric (see Fig. 18, left), but they are the same when there is no reflection symmetry (see Fig. 18, right). The analytical result in the right panel of this figure is just the result for class A (GinUE) evaluated at the center of the spectrum (red curve). In the case of chiral symmetry in the left panel, the analytical result is obtained by the following approximation. For the interval [0, n] the number variance is half the number variance of class A, Σ 2 A , because the eigenvalues can only fluctuate in and out of the interval at n, but Σ A has to be evaluated at 2n, Σ 2 (0, n) = 1 2 Σ 2 A (2n).
Finally, to show consistency between the spectral form factor and the number variance we The analytical results are given by Σ 2 (0, n) = 1 2 Σ 2 A (n) and Σ 2 (−n/2, n/2) = 2Σ 2 A (n) in the left panel and by Σ 2 A (n) in the right panel. calculate the number variance from the spectral form factor using Eq. (33). This relation assumes translational invariance and is not applicable in the center of the spectrum when the spectrum is reflection symmetric, and we can only give results for the q = 4 case, see Fig. 19. In this figure, we compare the direct evaluation of the number variance to the result obtained from the local spectral form factor atx = 0 andx = 0.8 (dashed curves). The agreement between the two shows that the nonstationarity of the number variance is "kinematical" and can be eliminated by a proper rescaling, see Eq. (26).

IV. CONCLUSIONS AND OUTLOOK
In this paper, we have studied long-range correlations of the non-Hermitian SYK model by means of the number variance and the spectral form factor of the real parts of the eigenvalues with results for the Ginibre or Ginibre-like ensembles as a benchmark. To eliminate unfolding ambiguities we have only considered non-Hermitian SYK models with a radially symmetric spectrum. A feature of spectral correlations of the real parts of the eigenvalues is that eigenvalues that are many level spacings apart, and are essentially uncorrelated, can have real parts that are close. This results in Poisson statistics already after a timescale of ∼ √ D. The early onset of Poisson statistics has the consequence that the collective spectral fluctuations can no longer be separated from universal eigenvalue fluctuations.
For small times, the spectral form factor deviates from the Ginibre result and shows an oscillatory behavior with a period conjugate to the overall width of the spectrum and an amplitude decreasing with time that is very sensitive to the number q of interacting Majoranas, a parameter that controls the fraction of independent matrix elements of the Hamiltonian in Fock space. The area below the peak is proportional to 2 N/2 / N q and decreases rapidly going from q = 2 to q = 6. Averaging over the oscillations, the small-time behavior of the spectral form factor is similar to that of Hermitian systems with a correlation hole. Although this is a typical feature of strongly interacting quantum chaotic systems, in this case, it is caused by collective ensemble fluctuations rather than by the nonuniversal dynamics at that timescale. Therefore, the nHSYK model describes both generic features of the universal quantum ergodic state reached around the Heisenberg time and nonuniversal, but still rather generic, properties of quantum interacting systems in its approach to ergodicity.
Having said that, we note that the spectral correlations of the real parts of the eigenvalues are not stationary and show deviations from the Ginibre ensemble that depend on the region of the spectrum that is considered. Since the corresponding spectral form factor is an average over the complete spectrum, we expect that, in general, it is not universal with a result that depends on the nonstationarity of the spectral correlations. Remarkably, in the nHSYK model, deviations from stationary compared to those of the Ginibre ensemble seem to average out, resulting in a much better agreement than could be expected. At this point, we do not have a good understanding of this remarkable coincidence, but we hope to further explore this in future work. Most likely it is related to the stationarity of the local spectral form factor introduced in this paper. In any case, our results also point to intrinsic limitations of the global observables that we have investigated to describe dynamical features. Their local counterparts, on the other hand, overcome these shortcomings and could prove an effective diagnostic of non-Hermitian quantum chaos.
Other related problems that are worthwhile to pursue, and seem within reach, are to find analytical results for the spectral density of the nHSYK model and also to compute analytically the spectral form factor, and number variance, for other universality classes, such as AI † and AII † .
Our final expression for the spectral form factor is given by This result differs from the expression quoted in Ref. [37] by the factor 2 −|p−s| . Asymptotically, for large D, it simplifies to which is in agreement with Ref. [37]. We note that the large-D limit of the spectral form factor (A8) does not commute with the τ → 0 limit. The Taylor expansion of Eq. (A8) gives while from the large D result (A9) we obtain (A11)