Complex spacing ratios: a signature of dissipative quantum chaos

We introduce a complex-plane generalization of the consecutive level-spacing distribution, used to distinguish regular from chaotic quantum spectra. Our approach features the distribution of complex-valued ratios between nearest- and next-to-nearest neighbor spacings. We show that this quantity can successfully detect the chaotic or regular nature of complex-valued spectra. This is done in two steps. First, we show that, if eigenvalues are uncorrelated, the distribution of complex spacing ratios is flat within the unit circle, whereas random matrices show a strong angular dependence in addition to the usual level repulsion. The universal fluctuations of Gaussian Unitary and Ginibre Unitary universality classes in the large-matrix-size limit are shown to be well described by Wigner-like surmises for small-size matrices with eigenvalues on the circle and on the two-torus, respectively. To study the latter case, we introduce the Toric Unitary Ensemble, characterized by a flat joint eigenvalue distribution on the two-torus. Second, we study different physical situations where nonhermitian matrices arise: dissipative quantum systems described by a Lindbladian, non-unitary quantum dynamics described by nonhermitian Hamiltonians, and classical stochastic processes. We show that known integrable models have a flat distribution of complex spacing ratios whereas generic cases, expected to be chaotic, conform to Random Matrix Theory predictions. Specifically, we were able to clearly distinguish chaotic from integrable dynamics in boundary-driven dissipative spin-chain Liouvillians and in the classical asymmetric simple exclusion process and to differentiate localized from delocalized phases in a nonhermitian disordered many-body system.

For hermitian systems, the by now universally accepted conjectures of Berry and Tabor [43] and of Bohigas, Giannoni, and Schmit [44] (see also Ref. [45]) assert, respectively, that classically integrable systems follow Poisson statistics of uncorrelated random variables, while systems with a chaotic semiclassical limit have statistics well described by Random Matrix Theory (RMT). Most astonishingly, many-body systems with no classical counterpart follow a similar rule. Poisson level statistics is found for integrable or (many-body) localized systems whereas RMT distributions are observed in generic thermalizing * lucas.seara.sa@tecnico.ulisboa.pt † ribeiro.pedro@gmail.com ‡ tomaz.prosen@fmf.uni-lj.si phases [46,47]. The power of the RMT approach relies on the fact that spectral fluctuations (measuring correlations of levels) are highly universal, depending solely on the symmetries of the system, and not on the details of particular models. For instance, the three classical Gaussian ensembles (GOE, GUE, and GSE) are completely determined by time-reversal symmetry, depending on a single parameter β = 1, 2, or 4, the Dyson index.
Since the early days of RMT, level spacing distributions, i.e. the distribution of the distance, s = ε i+1 − ε i > 0, between consecutive energy levels, ε i+1 , ε i , have proved to be a very useful and hence popular measure of spectral correlations in integrable and chaotic systems, i.e. a signature of quantum chaos. Indeed, for closed systems, spacings between uncorrelated levels display level clustering, while RMT statistics lead to level repulsion, with a characteristic power-law behavior of the spacing distribution, P (s) ∝ s β as s → 0, in the respective universality classes. Rather remarkably, the spacing distribution in the (universal) large-matrix-size limit is very well described by that obtained for 2 × 2 matrices, the Wigner surmise. Spacing distributions further allow the study of intermediate statistics, either with crossovers between Poisson and RMT statistics [48][49][50][51][52][53][54][55] or transitions between different RMT universality classes [50,[56][57][58][59]. Statistics of higher-order spacings (i.e. distance between kth-nearest neighbors) have also been considered over the years [60][61][62][63][64][65].
For nonhermitian systems, by a direct generalization [66] of the Berry-Tabor and Bohigas-Giannoni-Schmit conjectures to dissipative systems, we expect classically integrable systems and classically chaotic systems to follow Poisson and Ginibre level statistics, respectively. arXiv:1910.12784v1 [cond-mat.stat-mech] 28 Oct 2019 For random matrices from the Ginibre ensembles (i.e. matrices where all entries are iid Gaussian random variables) one finds cubic level repulsion, P (s) ∝ s 3 . Interestingly, all three Ginibre ensembles (GinOE, GinUE, and GinSE) have the same cubic level repulsion [66][67][68][69], independently of the Dyson index β. For those, a Wigner-like surmise, in terms of modified Bessel functions, has been recently proposed in Ref. [70], in which it was also shown that non-cubic level repulsion can exist in nonhermitian ensembles with different symmetries.
In order to compare theoretical predictions of RMT with actual measured or computed level sequences, one has to eliminate the dependence of the spacing distribution on the local mean spectral density, which is nonuniversal and system-dependent. This is achieved by a procedure known as unfolding [68,71], in which, in the case of a real spectrum, one changes from a sequence E j of levels to a new sequence e j = N (E j ), where N (x) is the level staircase measuring the mean number of levels below x. At the unfolded scale, the spacing distribution has a mean unit spacing and thus fluctuations can be uniformly compared across the spectrum. Unfolding is a nontrivial procedure since it requires an analytic expression (or accurate estimate) of the level density, which is not available in general. Furthermore, numerical unfolding sometimes proves ambiguous and numerically unreliable. In the case of a two-dimensional, i.e. complex, spectrum the situation is worse: there the unfolding is in principle ambiguous; even so, one can find a minimal prescription that guarantees uniform unfolded complex level density [69].
While spacing (and spacing ratio) distributions for real spectra are well understood [68,71,[92][93][94], and some results exist for spacings in complex spectra [66-70, 95, 96], two major shortcomings in the latter case remain to be addressed. On the one hand, to bypass the difficult and unreliable unfolding procedure, one is naturally lead to consider ratios of spacings in the complex plane; this remains, however, an open question. On the other hand, the existent studies on spacings in complex spectra focused solely on the distance, s > 0, between the complex eigenvalue and its nearest neighbor, neglecting the additional information contained in the angular (directional) correlations.
In this paper, we tackle both issues above by introduc-ing complex spacing ratios, as the ratio of the distance (taken as a complex number with magnitude and direction) from a given level to its nearest neighbor (NN) by the (complex) distance to the next-to-nearest neighbor (NNN); for a precise definition see Sec. II. Two comments are in order regarding these complex spacing ratios. First, when defining ratios for real spectra, level sequences are usually assumed to be ordered. There is, however, no global order in the complex plane, and hence all ratios which relied on the ordering have to be abandoned. The only remaining spacing ratio is, then indeed, the NN by NNN ratio, the modulus of which was introduced in Ref. [63] (and kth-nearest neighbor generalizations) for studies of real spectra. Second, this new spacing ratio (and not only its modulus!) can be as well defined for real spectra. It does not coincide with any of the aforementioned ratios; in particular, it adds a sign to the NN by NNN ratio of Ref. [63]. We emphasize that, while this sign might seem a minor difference in the case of real spectra, for complex spectra, the full angular dependence constitutes, arguably, the cleanest signature of dissipative quantum chaos. The paper is organized as follows. In Sec. II we define the complex spacing ratio, mention some of its qualitative features, point out the differences for integrable and chaotic spectra and state the key ideas behind our analytical results. In Sec. III we present exact analytical distributions and small-N surmises. In Sec. IV examples of application to different physical problems (driven spin-chains, nonhermitian many-body localization, and classical stochastic processes) are studied. We draw our conclusions in Sec. V. A detailed derivation of the analytical results is given in three appendices: the ratio distributions for uncorrelated random variables in d dimensions are computed in Appendix A; exact analytical distributions and small-size surmises are derived for hermitian random matrix ensembles in Appendix B, and for nonhermitian ensembles in Appendix C.

II. OVERVIEW AND MAIN RESULTS
Let the set {λ k } N k=1 be the spectrum of some hermitian or nonhermitian matrix. The levels λ k may, correspondingly, be real or complex. For each λ k , we find its NN (with respect to the distance in R or in C), λ NN k , and its NNN, λ NNN k , and define the (in general complex) ratio This definition is illustrated in Fig. 1-(a). We then seek the probability distribution function (N ) (z) of finding a spacing ratio with value z, which is defined either in the limit N → ∞, or, for a finite N , upon averaging over spectra of an ensemble of random matrices. If the spectrum is real, z ≡ r satisfies −1 ≤ r ≤ 1 and may not coincide with the ratio of consecutive spacings. If the spectrum is complex, z ≡ re iθ ≡ x + iy, with 0 ≤ r ≤ 1, and the distribution is not necessarily isotropic. We also consider the radial and angular marginal distributions, (r) = dθ r (r, θ) and (θ) = dr r (r, θ), respectively. We start by considering two paradigmatic cases: synthetic uncorrelated levels (corresponding to random diagonal matrices) and the Ginibre Ensembles. By natural extensions of the Berry-Tabor and Bohigas-Giannoni-Schmit conjectures, one expects integrable systems to have the same ratio statistics as uncorrelated levels and chaotic systems to follow Ginibre statistics. Due to the independence of levels in the synthetic spectrum, the presence of a reference level does not influence its two nearest neighbors and hence all ratios z have the same probability, which yields a flat distribution. In contrast, for random matrices, we expect the usual repulsion, with two immediate consequences. First, the ratio density should vanish at the origin; second, the repulsion should spread all the neighbors of the reference level evenly around it, leading to a suppression of the ratio density for small angles. Figure 1 shows the ratio density (z) in the complex plane for uncorrelated levels, (b), and GinUE matrices, (c), and confirms the expectations above. For uncorrelated levels the ratio is indeed flat inside the unit circle, i.e. Poi (z) = (1/π)Θ(1 − |z|), with Θ the Heaviside step-function. It immediately follows that the radial and angular marginal distributions are, respectively, Poi (θ) = 1/(2π) and Poi (r) = 2r, and thus cos θ = dθ cos θ Poi (θ) = 0. GinUE random matrices, on the contrary, have cubic level repulsion, GinUE (r) ∝ r 3 as r → 0 (note that one power of r comes from the area element on the plane), and the distribution shows some anisotropy, measured, for instance, by cos θ = dθ cos θ GinUE (θ) 0.24.
For a real (complex) spectrum, Fig. 2 (3) shows the distribution function of the level spacing ratio, z, both for uncorrelated levels and for GUE (GinUE) random matrices of different sizes as well as the radial (and angular) marginal distributions. Contrary to the case of consecutive spacings ratios, the distribution function for smallsize GUE or GinUE matrices, say with N = 3 or N = 4, does not qualitatively capture the large-N asymptotics, for large N (yellow histogram), there is a high probability of finding negative ratios, while for small N (red and blue dashed lines), the probability of positive r is higher. This small-N peak inversion can be understood as a boundary effect. For definiteness, consider matrices drawn from a hermitian ensemble. For N = 3, the sign of the ratios is completely fixed: the two levels at the edges must, by construction, have both neighbors to the same side and hence r > 0; the central level has one neighbor to each side and hence r < 0; it follows that the area below the peak for negative r is exactly 1/3 and for positive r is 2/3 (the analytical expressions below confirm this reasoning exactly). As N increases, the edge levels, which always have positive ratios, loose importance relative to the growing number of bulk levels, which tend to have negative ratios, and peak inversion follows. The argument for nonhermitian matrices is analogous: bulk levels favor large angles while boundary levels lead to small angles; at small N , boundary levels dominate, but they cannot compete in number with bulk levels at large N .
The strong N -dependence thus precludes any smallsize Wigner-like surmise using GinUE-drawn matrices. One of our main results is that these boundary effects can be overcome by using different ensembles with the same asymptotic large-N distribution. Figure 4 sketches the main idea of our approach. For a real spectrum, we obtain a surmise using the spacing ratios of the circular unitary ensemble (CUE) [68,97], whose spectrum lies along the unit circle, therefore avoiding boundary effects. Figure 2-(b) shows that the predictions of this method (solid red and blue lines) converge rapidly for  increasing N and give already a very good quantitative agreement for N = 3 and N = 4. The toric unitary ensemble (TUE), introduced in the next section, generalizes this idea for the case of a complex spectrum. Figs. 3-(g) and (h) show that the predictions obtained in this way for small N (solid and red lines) also qualitatively reproduce the large-N results.
A second main result of our work is to verify that these distributions do generalize the Berry-Tabor and Bohigas-Giannoni-Schmit conjectures to physical situations where the relevant operators have complex-valued spectra. By studying different physical models where nonhermitian matrices arise, we show that known integrable cases have a flat distribution of complex spacing ratios whereas generic cases, expected to be chaotic, conform to Ran-dom Matrix Theory predictions. Figures 6-(a)-(e) below illustrates our findings for a spin-1/2 chain, subject to boundary driving and/or bulk dissipation, modeled by Markovian Lindblad dynamics. The flat distribution of Fig. 6-(a), corresponds to a boundary driven XX chain with bulk dephasing, which is known to have an integrable Liouvilian. This contrasts with the nonintegrable cases, (b)-(e), where the distribution of complex spacing ratios is highly asymmetric and which are expected to reach the GinUE distribution in the thermodynamic limit. Similar results are reported in Sec. IV B for the case of non-unitary Hamiltonian dynamics, and in Sec. IV C for the spectrum of the Markov matrix describing the ASEP. This provides solid evidence that the complex level-spacing ratio distribution can be used to distinguish chaotic from integrable dynamics of operators with complex-valued spectra.

III. ANALYTICAL RESULTS: EXACT DISTRIBUTION FUNCTIONS AND SURMISES
In this section, we summarize our main analytical results regarding the complex spacing ratio distribution of independent levels and the small-N surmises obtained for the CUE and the TUE.
For independent levels the spacing ratios are isotropic. Therefore the only non-trivial distribution is that of r = |z|, which can be obtained analytically for d-dimensions (generalizing real, d = 1, and complex, d = 2, spectra). Furthermore, all joint-spacing distributions of more than one spacing factorize into single-spacing distribu-tionsP (s) and one can write the ratio distribution in terms ofP (s) only: At the unfolded scale, the d-dimensional single-spacing distributionP (s) is a Brody distribution [48], which recovers the standard exponential distribution for one-dimensional spectra. The ratio distribution in d dimensions, then follows. This shows that (after introducing a ddimensional volume element) the ratio distribution is, indeed, flat. For more detail on spacing ratios for uncorrelated random variables, and some generalizations, see Appendix A. We now address random matrix ensembles starting by the case of real spectra. The level-spacing ratio distribution function for N × N matrices drawn from arbitrary hermitian ensembles, (N ) (r), can be formally written as an (N −1)-fold integral over the joint eigenvalue distribution function [Eq. (B3)]. By specializing to the Gaussian ensembles, this quantity can be explicitly computed for small-size matrices, e.g. N = 3 [Eq. (B8)]. Other small sizes are still amenable to a brute-force evaluation of the integrals. However, we were not able to determine the complete asymptotic large-N distribution using this approach. Nonetheless, it can be employed to capture the scaling, (N →∞) GUE (r) ∝ r β , in the vicinity of r = 0. As shown in the last section, although larger values of N suppress the weight of boundary effects, the convergence towards the infinite-N limit is very slow. Convergence is much faster in the case of the circular ensembles (CE), where results for small-size matrices (N = 3, N = 4) from the CE already capture most of the features of the large-N asymptotics. Since for N → ∞, CE and GE have the same level-spacing ratio statistics, we can use CE small-size matrices as surmises for the GE large-N distribution. As for GE, for CE the level-spacing ratio distribution function for N × N matrices can be formally obtained in the form of an (N − 1)-fold integral. For N = 3, the ratio distribution reads which is evaluated in Eq. (B17), yielding a ratio of polynomials of r, whose explicit form is given in the Supplemental Material [98]. A similar expression was also obtained for N = 4 [Eq. (B19)]. For further details on hermitian ensembles, we refer the reader to Appendix B.
Finally, we turn to nonhermitian ensembles, considering, for simplicity, only the case β = 2. The general expression of the ratio distribution, for an arbitrary ensemble, is a 2(N − 1)-fold real integral over the ensemble's joint eigenvalue distribution [Eq. (C4)]. For the GinUE, the distribution for N = 3 can be computed explicitly [Eq. (C7)], but, again, it does not correctly describe the large-N asymptotics. The leading-order expansion in powers of r yields (N ) (r) ∝ r 3 , but is valid only around r = 0.
In order to eliminate boundary effects from a complex spectrum we consider the 2-dimensional analogue of the circular ensemble. This novel ensemble has eigenvalues equally distributed on the 2-dimensional (Clifford) torus, , which can be parameterized by two angles, ϑ ∈ (−π, π], ϕ ∈ (−π, π]. In analogy with the CUE, we dubbed it the Toric Unitary Ensemble (TUE). Therefore, P Setting N = 3, we compute a Wigner-like surmise for the complex spacing ratio distribution for nonhermitian random matrices, The integral of Eq. show that the convergence of the radial marginal distribution is similar to that of the real case: both N = 3 and N = 4 provide good approximations, the latter being already almost indistinguishable from large-N exact diagonalization data. The angular marginal distribution has a much slower convergence, especially near θ = ±π. Although the qualitative features are already captured for N = 3, quantitatively, one can still distinguish the discrepancies even for N = 5 in Fig. 3-(h), although the agreement does improve as N increases. For further details on nonhermitian ensembles, see Appendix C.

IV. PHYSICAL APPLICATIONS
We now determine the complex spacing ratio distribution of several different numerical examples of current interest. In Sec. IV A we consider the Lindbladian description of boundary-driven dissipative spin-chains, in Sec. IV B we address a nonhermitian Hamiltonian modeling many-body localization, and in Sec. IV C we study a classical stochastic process. A. Boundary-driven dissipative spin-chains A simple way of modeling quantum open systems is by employing a master equation approach to describe the dynamics of the system's reduced density matrix. When the environment is Markovian, this procedure substantially simplifies and the master equation acquires the Lindblad form where H is the Hamiltonian and W µ , with µ = 1, . . . , D, are called jump operators, modeling the systemenvironment interaction.
Here, we study the spectrum of a family of nonhermitian operators L for a well studied physical setup of a chain of spins-1/2. In the middle of the chain, the magnetization along z is conserved and the net role of the environment is to dephase the system, i.e. decrease off-diagonal amplitudes of the density matrix when written in the z-basis. At the two ends of the chain, the spin-magnetization can be injected or extracted at fixed rates. This model had been extensively used for studying non-equilibrium spin transport [5,8,99].

Model
We consider a chain of N spins 1/2 evolving in time by the action of a Lindblad-Liouvillian operator, given by Eq. (8), and schematically represented in Fig. 5-(a). H belongs to a family of next-to-nearest neighbor Heisen-berg XXZ Hamiltonians, with σ α the Pauli operators, α ∈ {x, y, z} and ∈ {1, 2 . . . , N }, and J (J ) the nearest (next-to-nearest) neighbor exchange coupling and z-axis anisotropy ∆ (∆ ). To model bulk dephasing and spin injection, we consider two types of incoherent jump processes (in total D = N + 4 of them): (i) bulk dephasing of all spins, (ii) amplitude damping (spin polarization) processes at the boundaries, Here, γ controls the dephasing rate and γ ± L/R controls the spin injection, (+), and extraction, (−), at the left, L, or right, R, ends of the chain. The model is, thus, characterized by the nine parameters J, J , ∆, ∆ , γ, γ ± L,R , which allows tuning its integrability or chaoticity. The Hilbert space is spanned by the states |s 1 , . . . , s N , with s = ±1. The space of density matrices-the Liouville space, K-in which L acts, is spanned by ||s 1 , . . . , s N ; s 1 , . . . , s N = |s 1 , . . . , s N ⊗ s 1 , . . . , s N | T . Using this notation, we formulate the spectral problem for the Liouvillian superoperator L in terms of a 4 N × 4 N matrix representation acting on a 4 N -dimensional density operators ρ ∈ K, The Note that, for M = N , each complex conjugate pair of eigenvalues of the Liouvillian is divided across two sectors of symmetric magnetization, i.e. if sector M contains the eigenvalue Λ, then sector |2N − M | contains the eigenvalue Λ * . The different sectors M must be analyzed separately because spectra corresponding to different conserved quantum numbers form independent level sequences that superimpose without interacting [71].

Numerical results
Numerical results were obtained by exactly diagonalizing the matrix representation of L, Eq. (12), for different chain length N and in specific sectors M . The largest system we diagonalized was N = 10 spins in the sector with magnetization M = 7, which corresponds to a 77520 × 77520 matrix. The following four cases of parameters were studied: • (Deph) Boundary driven XX chain with bulk dephasing. Numerical parameters chosen as J = 1, This model can be mapped onto the Fermi-Hubbard model with imaginary interaction U = iγ [5] and hence is Bethe-ansatz integrable.
• (A) Isotropic Heisenberg (XXX) chain with puresource/pure-sink driving and no dephasing. Numerical parameters chosen as The steady-state of this model is known to be integrable [8], but the bulk of the spectrum is likely not integrable.
• (B) XXX chain with arbitrary boundary-driving and no dephasing. Numerical parameters chosen as J = ∆ = 1, J = ∆ = 0, γ = 0, γ + L = 0.5, The bulk Hamiltonian of this model is integrable, but, by adding a generic boundary-driving, not even the steadystate is expected to be exactly-solvable.
• (C) XXZ chain with next-to-nearest neighbor interactions, arbitrary boundary-driving, and no dephasing. Numerical parameters chosen as J = J = 1, ∆ = 0.5, ∆ = 1.5, For this model, not even the bulk Hamiltonian is integrable.
We applied the procedure described at the beginning of Sec. II to compute the distribution of the complex spacing ratios for the five models depicted in Fig. 6. There is a striking difference between the integrable model (Deph), and the others, which are expected to be chaotic. The dephasing-XX model, Fig. 6-(a), displays a distribution similar to that of uncorrelated levels. Models B, C and RL, Figs. 6-(c), (d), (e), respectively, clearly conform to RMT statistics. Model A, Fig. 6-(b), on the other Table I. Single-number signatures of integrability/chaos for different Liouvillians: models Deph, A, B, C, and RL. They are compared with exact analytical results for uncorrelated random variables (labeled Poisson), numerical exact diagonalization of (10 4 × 10 4 ) random GinUE matrices, and TUE surmise estimates for N = 3, 4, 5 (subscripts denote matrix size) computed from Eq. (C9). The convergence of cos θ computed from the TUE surmises is much slower than that of r , as noted in the text.

Poisson
Deph

Single-number signatures
Next, we try to capture the main features of the distribution of complex spacing ratios through a reduced set of numbers, which we call single-number signatures. A popular single-number signature, used for the ratio of undirected spacings, is the degree of level repulsion α, i.e. the exponent describing the power-law behavior of the radial marginal distribution, (r) ∝ r α , as r → 0 or, equivalently, α = lim r→0 log (r)/ log r. For hermitian random matrices it is given by the Dyson index, α = β, for nonhermitian random matrices from the universality class of either GinOE, GinUE, or GinSE it is α = 3, while for real independent random variables it is α = 0, and for complex uncorrelated random variables α = 1. Although the degrees of repulsion α just stated can be easily checked against the numerical spectra and the above predictions confirmed, an actual computation of α for a given spectrum introduces a large relative error. An alternative measure of the radial distribution is given by its moments, for instance, the mean r . For independent random variables, we can compute exactly r = 2/3, while for GinUE matrices we numerically find r ≈ 0.74.
To measure the anisotropy of the angular marginal distribution, we consider cos θ , which is zero for a flat distribution and positive (negative) when small angles are enhanced (suppressed), in particular, cos θ ≈ −0.24 for large-N GinUE matrices.
We give the values of cos θ and r for the five Liouvillians in Table I (the spin-chain Liouvillian values are for N = 10, M = 7). From the radial measure r it is difficult to discern the integrability or chaoticity of the different models. Indeed, the values for all four models A, B, C, RL are within 3% of each other. On the contrary, as anticipated in Sec. II, the angular distribution offers a more sensitive signature. From the value of cos θ , the dephasing-XX model clearly supports Poisson statistics and models C and RL are very close to RMT statistics. Model B, which seemed also very close to RMT statistics from Fig. 6 and from the value of r here shows a more significant deviation. Finally, model A has a value of cos θ almost exactly halfway between uncorrelated levels and RMT statistics, attesting its intermediate behavior, at least for the sector dimensions considered.

Finite-size scaling
We now provide a finite-size analysis of the dephasing-XX model, which conforms to Poisson level statistics, model A, with intermediate statistics, and model C, with RMT statistics.
Dephasing-XX model.-We considered single-number signatures cos θ and r as a function of sector dimension in Fig. 7. Both clearly tend to the expected value for uncorrelated random variables (dashed line) as k NM increases. There is also a visible difference between sectors with even or odd M , with sectors of even M tending faster to the large-dimension universal limit. This aspect is also visible in Fig. 7-(d).
Model A.-We observed that (i) there is a smaller degree of level repulsion here than for fully chaotic systems (and which does not increase substantially when k NM grows by nearly two orders of magnitude) and (ii) some anisotropy is developing as k NM increases. In Fig. 8 we plot the two single-number signatures cos θ and r . While the anisotropy indeed grows (slowly) with k NM , the average of the radial marginal distribution is approximately flat. No difference between even M and odd M is visible in this case. Contrary to the dephasing-XX model above, for which the N = 10, M = 7 sector is already very close to the limiting Poisson statistics, the convergence of model A towards either Poisson or GinUE statistics is much slower. From these results, it is, therefore, inconclusive whether the model is tending very slowly to RMT statistics (as favored by Fig. 8-(a)) or if it follows some type of intermediate statistics. Considerably larger sector dimensions are, unfortunately, out of reach of current computational capabilities.
Model C.-Finally, we consider a chaotic Liouvillian, model C. Here, the universal limit of RMT statistics is quickly attained. Figure 8 depicts the two single-number signatures cos θ and r and confirms the fast convergence. For the largest sectors diagonalized, the results are already compatible, within their statistical errors, with the infinite-size limit.

B. Disordered open system and detection of many-body localized regime
After a quench, local observables of chaotic systems thermalize to values that can be predicted by a thermodynamic ensemble average [100]. However, in the presence of sufficiently strong disorder for a given system size [101], isolated quantum systems, even interacting ones, may fail to thermalize-a phenomenon dubbed many-body localization (MBL) [47,102]. Spectral properties in the many-body localized regime resemble that of integrable models. In fact, some proposals to model MBL rely on approximate locally conserved quantities [103]. Recently, numerical observation of MBL regime has also been reported for nonhermitian Hamiltonians [36]. Moreover, within the delocalized (ergodic) regime, we show that the complex spacing ratio distribution is able to distinguish between GinUE statistics and those of another symmetry class, AI † [25,70]. This firmly supports the claim of Ref. [36] that the model considered therein belongs to this symmetry (universality) class.

Model
We consider the model of Ref. [70] consisting of hardcore bosons on a one-dimensional lattice with N sites and periodic boundary conditions. Nonhermiticity arises due to an alternating on-site gain/loss terms. The (nonhermitian) Hamiltonian reads where b † j (b j ) is the creation (annihilation) operator of a hard-core boson at site j, n j = b † j b j is the particlenumber operator, J is the hopping strength, U gives short-range repulsion, γ measures the nonhermiticity, and the local disorder h j is uniformly distributed in [−h, h]. The Hamiltonian conserves particle number, hence we divide the Hilbert space into sectors of fixed particle number D. We decompose H = N D=0 H D , where each H D is a N D × N D matrix.
Applying the numerical procedure described at the start of Sec. II, we computed the distribution of the complex spacing ratios for the localized and delocalized While in the localized regime (h = 10) the finite-size scaling is consistent with a statistical signature of uncorrelated levels, the delocalized regime, even if clearly non-Poissonian, does not conform to GinUE statistics. Instead, it attains the values labeled by AI † , obtained by sampling random matrices from the AI † symmetry class [25,70].
These findings show that complex spacing ratios are not only effective in discriminating between localized and delocalized phases, but can also be used to distinguish between random matrix ensembles with different symmetries.

C. Classical stochastic process
Classical stochastic processes are widely used to model physical, chemical and biological systems. The solution for a classical stochastic process is obtained by specifying the continuous-time evolution of a probability vector of the system, P , governed by a Markov matrix M : ∂ t P (t) = M P (t), i.e. P (t) = exp{M t}P (t = 0). By conservation of probability, the columns of the Markov matrix must add up to zero. It then follows that the diagonal elements of M are fully-determined by the off-diagonal elements and we can write M jk = A jk − δ jk m A mk , with δ jk the Kronecker delta and A jj = 0. Among the most studied classical stochastic process are asymmetric simple exclusion processes (ASEP), used to study transport of interacting particles in one dimension. In the following, we analyze the complex-valued spectrum of the ma- trix M for integrable and non-integrable ASEP using the complex spacing ratio distribution. We show that while the first case follows Poisson statistics of uncorrelated levels, the second conforms to RMT predictions.

Model
Consider a set of hard-core classical particles on a Nsite ring with nearest neighbor hoppings. The hard-core condition reduces the dimension of configuration space to 2 N . Within each time interval dt, any particle can hop from site j to site j + 1 with probability p dt and from site j to site j − 1 with probability q dt. When p = q, this defines the ASEP [104][105][106][107]. To break integrability, we further consider a staggering of the hoping probabilities by requiring that the probability of hopping from odd to even sites (p 1 dt if hopping clockwise) is different from that of hopping from even to odd sites (p 2 dt), and similarly from anticlockwise jumps, with probabilities q 1 dt and q 2 dt, respectively. Finally, we admit the possibility of particles entering or leaving the system at site j = 1, with probabilities µ + dt and µ − dt, in each time interval. Assuming N to be even, the matrix A for this process is given by

Numerical results
We numerically diagonalized the Markov matrix M , described in the preceding section, for a ring with N = 16 sites (M is 65536 × 65536). For simplicity, we considered a totally asymmetric exclusion process (TASEP, q 1 = q 2 = 0), fixed p 2 = 1 and set µ + = µ − = 0.5. Since µ + , µ − = 0, particle conservation is broken and, hence, we do not restrict M to sectors of fixed particle number. We then considered two cases: non-staggered hopping, p 1 = p 2 = 1, for which the model is known to be integrable [104], and staggered hopping, p 1 = 0.2 = p 2 , which we expect to break integrability. This expectation is confirmed by the distribution of complex spacing ratios, see Fig. 11. The complex spacing ratio distribution for nonstaggered hopping, Fig. 11-(a) is approximately flat (the inhomogeneity can, as before, be related to finite-size effects). The distribution for staggered hopping, Fig. 11-(b) clearly presents level repulsion and suppression at small angles, with − cos θ = 0.2356 (24) and r = 0.7382 (7). Both effects are compatible with the Ginibre universality class (recall that, for 10 4 × 10 4 matrices from the GinUE, we found − cos θ = 0.24051(61) and r = 0.73810 (18)). These results show that the complex spacing ratio distribution is also capable of discriminating between integrable and nonintegrable classical stochastic processes.

V. CONCLUSIONS AND OUTLOOK
We introduced complex spacing ratios to analyze universal spectral features of nonhermitian systems (integrable and chaotic). We found that angular correlations between levels in dissipative systems provide a clean signature of quantum chaos: uncorrelated random variables, which describe integrable systems, have a flat, and hence isotropic, ratio distribution in the complex plane, while for RMT ensembles from the Ginibre universality class there is a suppression of small angles in the large-N limit. We also reencountered the familiar cubic level repulsion in the latter case.
Our results show that complex spacing ratios allow to clearly distinguish (known or conjectured) integrable systems from chaotic ones. Compelling numerical evidence for this claim has been given by a finite-size analysis of boundary-driven spin-chain Liouvillians and classical stochastic processes. Complex spacing ratios can also differentiate the many-body-localized regime from the delocalized regime in the nonhermitian disordered many-body systems. Furthermore, in the delocalized phase, single number signatures, − cos θ and r can also discriminate between Hamiltonians in different symmetry classes.
We provided surmises of the large-N complex spacing ratio distribution for GUE and GinUE ensembles. These were obtained for small matrices, with N = 3, 4, using the CUE and its two-dimensional generalizationthe Toric Unitary Ensemble-, which overcome the large finite-size effects observed for small-size GUE and GinUE matrices. Even so, the angular marginal distribution was found to have a somehow slow convergence towards the N → ∞ limit.
Due to their ability to unambiguously discriminate between regular and chaotic dynamics, without the need of unfolding, we expect complex spacing ratios to play an important role in future studies of dissipative quantum chaos and classical stochastic processes. Specifically, complex spacing ratio statistic can be used as a clean and simple empirical detector of integrability, as well as an order parameter characterizing ergodicity-breaking transitions in nonhermitian systems.
An interesting open question is whether complex spacing rations can be used to discriminate between symmetry classes other than the Ginibre and AI † , for instance, those introduced recently in Ref. [70].
Finally, the toric unitary ensemble introduced in Sec. II, modeling the Coulomb gas on the Clifford torus, also warrants further study. Besides analyzing the properties of random matrix realizations of this novel ensemble, it would be interesting to encounter physical systems for which the TUE arises naturally. Taking isotropy as a starting point, i.e. assuming that the distribution of complex spacing ratios only depends on the absolute value of the ratio, r, we now show that it is, indeed, flat for uncorrelated random variables (the Poisson spectrum). The independence of the levels simplifies the problem enough so that we are also able to exactly compute the more general ratio, r mk , of the distance to the m th -nearest neighbor (mNN) by the distance to the k th -nearest neighbor (kNN) (which reduces to the ratio discussed in the main text when m = 1, k = 2). We do all calculations for spectral points (levels) represented by vectors in d-dimensional Euclidean space. Real, complex, and quaternionic spectra correspond to d = 1, 2, 4, respectively, but our results also apply to other cases, say, uncorrelated random vectors in three-dimensional space.

Joint spacing distributions
By translational invariance, we can consider the level for which the ratio is being computed (the reference level) at the origin. To compute the probabilityP (s) ds of finding its NN at a distance s, we introduce the conditional probability g(s) ds of finding a level in [s, s + ds] given our reference level at the origin, and the probability H(s) = ∞ s ds P (s ) = 1 − s 0 ds P (s ) of having no level in [0, s] (the hole probability). By independence of the levels, the probability g(s) ds is actually independent of the presence of the reference level. For the NN to be at s we must verify that (i) there is a level at s and (ii) there are no levels in [0, s], whence we conclude that P (s) = g(s) H(s) . (A1) Noting that the hole probability is equal to 1 − F (s), where F (s) is the cumulative distribution ofP (s), we can equally well expressP (s) solely in terms of H(s),P (s) = − dH ds . Alternatively, we can also write g(s) as a function of H(s) only, g(s) = (A2) Analogously, the joint distribution of the first-kNN spacings iŝ It is worthwhile to note that we can express the whole hierarchy of joint probabilities solely in terms of the singlevariable functionsP (s), g(s) or H(s), whichever is easier to compute in a given situation. This is, of course, a particularity of independent random variables, and does not carry over to random matrix ensembles.
The distribution of the (absolute value) of the ratio r = s 1 /s 2 is given in terms of the joint distributionP (s 1 , s 2 ), and, hence, it is also completely determined by the single spacing distributionP (s) [Eq. (2)]: In the last line, we have expressed the ratio distribution solely in terms of the single spacing probability. It only remains to compute one ofP (s), which we do in d-dimensions in the next section. Likewise, the mNN by kNN ratio, r mk ≡ s m /s k , is defined in terms of the joint spacing distributionP (s 1 , . . . , s k ) and is fully determined by the single spacing distribution: mk (r mk ) = ds 1 · · · ds m · · · ds kP (s 1 , . . . , s m , . . . , s k ) δ r − s m s k = ds 1 · · · ds m−1 ds m+1 · · · ds kP (s 1 , . . . , s m−1 , rs k , s m+1 , . . . , s k ) (A4)

Uncorrelated random variables in d-dimensional space
We consider the spectrum to be composed by N iid random variables, supported in a d-dimensional ball of radius R. At a later point, we shall take the limits N, R → ∞ with constant mean density N R −d = 1. The probabilities g(s) ds and H(s) are then given by ratios of d-dimensional volumes V d (L) = π d/2 /Γ(d/2 + 1)L d , where L is a length.
To determine g(s), we note that any one of N − 1 levels can be the NN if it falls inside the interval [s, s + ds], whence it follows that Taking the limits N, R → ∞, we obtain immediately that g(s) ∝ s d−1 .
Regarding H(s), since all other (N − 2) levels are independent and must lie beyond a distance s, we have To be able to properly take the limits, we need to unfold the spectrum to a unit mean, i.e. we change variables to s = s/ s . Note that, in the computation of g(s), the unfolding would only give an overall constant, so we did not need it to proceed. Using Eqs. (A1), (A5), (A6), we havê which is correctly normalized, as it should be. We then have or, taking N → ∞ and using the asymptotic behavior of the Γ function, lim N →∞ N α Γ(N )/Γ(N + α) = 1, for any α ∈ C, s = Γ(1 + 1/d)N −1/d R.
In terms of the unfolded variable s, the hole probability reads Taking limits, H(s) = exp −Γ(1 + 1/d) d s d and the (unfolded) spacing distribution is given bŷ Note that, for d = 1, we recover the standard exponential distribution,P (s) = e −s . In d-dimensions, the spacing follows, instead, a Brody distribution [48]. The NN-and NNN-joint spacing distributionP (s 1 , s 2 ) can be written solely in terms of the single spacing distri-butionP (s) by inserting Eq. (A10) into Eq. (A2), Finally, the joint distribution of the kNN spacings, Eq. (A3), reads, in d dimensions:

Ratio distribution
We now turn to the ratio distributions. Henceforth, we always assume that we are at the unfolded scale and denote the spacings by s instead of s. Inserting Eq. (A10) into the last equality of Eq. (2), we obtain Eq. (4), The constraint enforced by the Θ-function implies that the distribution is supported in the d-dimensional unit ball, which we parameterize by the radial distance r and the (d − 1)-dimensional solid angle Ω d−1 . By recalling that of Ω d−1 , by using Eq. (4), and by noting that dΩ d−1 = S d−1 gives the area of the unit sphere in d dimensions and i.e. distribution is indeed flat since it is given by the inverse of the volume of its support. We next consider the distribution of the mNN by kNN ratio, of which the above result is a special case (m = 1, k = 2). Inserting Eq. (A10) into the last equality of Eq. (A4), we obtain Finally, note that, a flat distribution in the d-dimensional unit ball is only possible if (r mk ) ∝ (r mk ) d−1 for all r mk , as in Eq. (4). This implies m = 1 and k − m = 1 ⇒ k = 2. We thus see that a flat distribution is a peculiarity of the NN by NNN ratio and is not achieved by any other combination of m, k.

Gaussian ensembles
Equation (B3) is valid for an arbitrary hermitian ensemble. We now specialize for the case of the Gaussian ensembles, GO/U/SE, labeled by the Dyson index β. The joint eigenvalue distribution reads In terms of the variables of Eq. (B3), we have The integration in u is Gaussian and can be readily performed, yielding du exp We finally obtain the distribution of the ratio as an (N − 2)-fold integral: For small N (N = 3, 4) the integrals in Eq. (B7) can be computed exactly. Unfortunately, contrary to the consecutive spacings ratio, for NN by NNN ratios, the small-size expressions do not accurately describe the large-N asymptotics.
For N = 3, no s j -integrals exist in Eq. (B7). Furthermore, the r-dependence can be factored out of the v-integral and no integrals have to be performed at all: where the β-dependant normalization is N = 9/4 for β = 1, N = 27 √ 3/(2π) for β = 2 and N = 243 √ 3/(2π) for β = 4. The distribution of Eq. (B8) for β = 2 is plotted in Fig. 2-(a), in comparison with exact diagonalization results. For N = 4, we must perform an additional integral in s (here for β = 2), If we denote then Eq. (B9) reads which is evaluated (with the correct normalization) to (4) GUE (r) = 1 4π with the polynomials B (4) k given in the Supplemental Material [98]. We compare the analytical predictions for the ratio distribution for small-size matrices with numerical results from exact diagonalization of GUE-drawn random matrices in Fig. 13-(a) and (b). The agreement is perfect, which was to be expected since the computation is exact. These results are not, however, particularly useful in practice, since universality is only displayed for large-N . We thus turn to the case N → ∞.
When N → ∞, we can rewrite Eq. (B7), discarding all exponentials suppressed by 1/N : Although we cannot compute this integral exactly, note that its multiplying prefactor gives the exact distribution of r for r → 0. It also qualitatively describes the distribution for all r, albeit missing the exact heights of the peaks of positive and negative r. We thus obtain the loosely approximating distribution, At any rate, the absence of a term −r inside the denominator (which was killed by the limit N → ∞) completely distinguishes this result from the case N = 3, see Fig. 2 (the black line is the approximation of Eq. (B14)).

Circular ensembles
We discussed in Sec. II how the difference between the N = 3 and N → ∞ statistics is due to boundary effects. To eliminate them we should consider periodic boundary conditions, i.e. identify the ends of the spectrum. Hence, we consider the circular ensembles, whose spectrum is supported on the unit circle and whose joint eigenvalue distribution is Note that, although the eigenvalues are complex (e iφj ), they are fully described by real angles φ j ∈ (−π, π]. The spacing ratio is defined in terms of the real variables, r = (φ 2 − φ 1 )/(φ 3 − φ 1 ), i.e. we are measuring the spacings on the circle, not in the embedding space, C. By rotational invariance of the circle, we may set φ 1 = 0. We can rewrite the Vandermonde interaction as e iφj − e iφ k = sin β (|φ j − φ k | /2). The general result of Eq. (B3), applied to the circular ensembles, reads We now evaluate the preceding integral for N = 3 and N = 4, restricting ourselves to the complex case, β = 2. We have which yields, after normalization,   For N = 4, we have an additional integral in s to perform, Following the procedure leading to Eq. (B12), the integral is evaluated (with the correct normalization) to with the polynomials Q (4) k given in the Supplemental Material [98]. The distribution of Eq. (B17) is plotted in black in Fig. 13-(d), in comparison with numerical diagonalization of N = 3 CUE matrices, and in blue in Fig. 2-(b), in comparison with exact diagonalization of large-N GUE matrices. Although in the large-N limit we can give only approximate expressions for of the complex spacing ratio distribution, the small-size surmises computed for the circular ensembles describe very well the universal large-N asymptotics. Indeed, the distribution for the N = 4 CUE is already indistinguishable (to the naked eye) from the numerical N → ∞ results.
In terms of these new variables, the δ-function constraints (fixing the real and imaginary parts of z) are and the Θ-function constraints (requiring all λ n with n > 3 to be further away from λ 1 than λ 3 ) are The distribution function of z is again obtained by integrating the joint eigenvalue distribution multiplied by the constraints of Eqs. (C2) and (C3). Integrating over p and q using the δ-functions, we arrive at the distribution for z for an arbitrary nonhermitian ensemble: × P (N ) (u, u + sx − ty, u + s, u + a 1 , . . . , u + a N −3 ; v, v + tx + sy, v + t, v + b 1 , . . . , v + b N −3 ) . (C4)

Ginibre Unitary Ensemble
We now restrict ourselves to the GinUE (complex Gaussian iid entries, β = 2), whose joint eigenvalue distribution reads: Replacing the x j , y j by the variables of Eq. (C4) and performing the Gaussian integration over the two variables u, v, we arrive at the ratio distribution for the Ginibre ensemble, GinUE (x, y) ∝ Θ 1 − (x 2 + y 2 ) (x 2 + y 2 )(1 + x 2 + y 2 − 2x) × ds dt (s 2 + t 2 ) 4 exp −(s 2 + t 2 ) (1 + As before, the distribution for N = 3 follows from Eq. (C6) without the need to perform any integrals explicitly. Indeed, in polar coordinates x = r cos θ, y = r sin θ, we get, after normalization, As for the hermitian case, the leading order behavior (i.e. the first term in a power expansion in r) of the distribution for N → ∞ can be obtained without carrying out any integral. Although it does not give a good quantitative match, it captures the high (low) density at large (small) angles. By discarding all exponentials suppressed by 1/N from Eq. (C6), factoring out terms containing z from the s and t integrals, we obtain the prefactor, (N →∞) GinUE (r, θ) ≈ 12 π r 2 (1 + r 2 − 2r cos θ) (1 + r 2 ) 5 Θ(1 − r) . (C8)

Toric Unitary Ensemble
We now want to eliminate boundary effects from a complex spectrum by considering a 2-dimensional analog of the circular ensembles. Recall that the circular ensembles had eigenvalues on the unit circle S 1 . A possible generalization would be to consider eigenvalues on the sphere S 2 ⊂ R 3 , which would be provided by the Spherical Unitary Ensemble (SUE) [93,[108][109][110][111], of matrices A −1 B with both A, B GinUE matrices. However, while belonging to the same universality class as the GinUE, also for the SUE the convergence to the large-N limit is quite slow. Instead, we consider eigenvalues on the 2-dimensional (Clifford) torus T 2 = S 1 ×S 1 ⊂ S 3 ⊂ R 4 , which show a very fast convergence. We parameterize the torus by two angles ϑ ∈ (−π, π], ϕ ∈ (−π, π], with a generic point P ∈ T 2 given by P = (1/ √ 2) (cos ϑ, sin ϑ, cos ϕ, sin ϕ). In analogy with the CUE, we want to construct a flat joint eigenvalue distribution on T 2 , which we call the distribution of the toric unitary ensemble (TUE). Since the Clifford torus has no curvature, the distribution is simply given by the Vandermonde interaction, P (N ) TUE (ϑ 1 , . . . , ϑ N ; ϕ, . . . , ϕ N ) ∝ |∆ T 2 | 2 . ∆ T 2 is given by the distance between points in the embedding space parameterized by the eigenvalues. That is, if P j ∈ T 2 is parameterized by the angles (ϑ j , ϕ j ) then the Vandermonde interaction is |∆ T 2 | = j<k P j − P k R 4 . One can check that, with this reasoning, the usual Vandermonde terms for the Gaussian, Ginibre, circular and spherical ensembles coincide with, respectively, |∆ R | = j<k P j − P k R , |∆ R 2 | = j<k P j − P k R 2 , |∆ S 1 | = j<k P j − P k R 2 , |∆ S 2 | = j<k P j − P k R 3 (with P j in the respective embedding spaces). Using our parameterizantion of the torus and considering only β = 2, the Vandermonde interaction reads |∆ T 2 | 2 = j<k [2 − cos(ϑ j − ϑ k ) − cos(ϕ j − ϕ k )].
We can then write down the joint eigenvalue distribution for the toric unitary ensemble (TUE), given by Eq. (6) Having introduced the relevant joint eigenvalue distribution, the remaining procedure is straightforward. By rotational invariance in both factors S 1 , we can set ϑ 1 = 0 and ϕ 1 = 0. The complex spacing ratio is, accordingly, z = (ϑ 2 + iϕ 2 )/(ϑ 3 + iϕ 3 ). If we then insert Eq. (6) into the general ratio distribution, Eq. (C4), we obtain For N = 3, the double integral to be performed was given in Eq. (7), The integral of Eq. (7), and its generalizations for N = 4, 5, . . . , can be numerically integrated (the analytic expression is far too involved to be useful) and describe very well the large-N asymptotics of the GinUE universality class, see Figs. 3-(e)-(h).