Statistics of the Spectral Form Factor in the Self-Dual Kicked Ising Model

We compute the full probability distribution of the spectral form factor in the self-dual kicked Ising model by providing an exact lower bound for each moment and verifying numerically that the latter is saturated. We show that at large enough times the probability distribution agrees exactly with the prediction of Random Matrix Theory if one identifies the appropriate ensemble of random matrices. We find that this ensemble is not the circular orthogonal one - composed of symmetric random unitary matrices and associated with time-reversal-invariant evolution operators - but is an ensemble of random matrices on a more restricted symmetric space (depending on the parity of the number of sites this space is either ${Sp(N)/U(N)}$ or ${O(2N)/{O(N)\!\times\!O(N)}}$). Even if the latter ensembles yield the same averaged spectral form factor as the circular orthogonal ensemble they show substantially enhanced fluctuations. This behaviour is due to a recently identified additional anti-unitary symmetry of the self-dual kicked Ising model.


I. INTRODUCTION
The quantum chaos conjecture [1][2][3] states that a quantum system is chaotic if the correlations of its energy levels have the same structure as those of random hermitian matrices [4,5]. This "conjecture" originates from studies on single-particle quantum systems, where the aforementioned property can be connected to the conventional chaoticity of the system (i.e. sensitivity of system's trajectories to initial conditions) in the classical limit [6][7][8][9][10][11].
For quantum many-body systems with no well defined classical limit the quantum chaos conjecture can be taken as a definition of quantum chaos. Indeed, an extensive number of numerical studies (see, e.g., Refs. [12][13][14][15]) established that systems with random-matrix spectral correlations display many features that are intuitively connected to chaos. In particular, spectral correlations are a widespread diagnostic tool to test numerically whether a many-body system is expected to be ergodic. Until recently, however, the theoretical explanations of this phenomenon where extremely scarce: no analytical method was known to deduce the spectral correlations from the Hamiltonian of the system or from the time evolution operator.
The situation has changed drastically over the last few years, when a number of settings and methods have been proposed to derive analytically the spectral form factor (SFF) (i.e. the Fourier transform of the twopoint correlation function of energy levels). Specifically, Refs. [16,17] established random matrix spectral fluctuations in long-ranged (but non-mean-field) periodically driven spin chains. Further on, Refs. [18][19][20][21] demonstrated the emergence of random-matrix spectral correlations in periodically driven local random circuits, where the interactions are determined by random two-site gates acting on neighbouring sites and chosen (once and for all) at the beginning of the evolution. In particular, analytical results were provided in the limit of large local Hilbert space dimension. Finally, Ref. [22] provided an exact result for the spectral form factor in the self-dual kicked Ising model: a system of spin-1/2 variables which are interacting locally with an Ising Hamiltonian and are periodically "kicked" by a longitudinal magnetic field. The term "self-dual" indicates that the longitudinal field and the Ising coupling are set to specific values. The key property to obtain the exact result is that, at aforementioned specific values of the couplings, the problem can be formulated in terms of a transfer matrix "in space" (i.e. propagating in the spatial direction, rather than in the temporal one) which is unitary.
The spectral form factor alone, however, is not a sufficient evidence for claiming the chaoticity of a system. Indeed, to invoke the quantum chaos conjecture one needs all the spectral correlation functions, not just the two point one. The goal of this paper is to provide such a result in the case of the self-dual kicked Ising model. We will generalise the space-transfer-matrix method of Ref. [22] to find expressions for higher moments of the spectral form factor and use them to obtain rigorous lower bounds. Then, we will demonstrate numerically that the bounds are saturated.
The rest of the paper is laid out as follows. In Sec. II we introduce the model and the quantities of interest (i.e. the spectral form factor and its higher moments). In Sec. III we identify the ensembles of random matrices which is relevant for the self-dual kicked Ising model and provide a prediction for the higher cumulants of the spectral form factor. In Sec. IV we provide the aforementioned lower bounds on the higher cumulants and in Sec. V we show numerically that the bounds are saturated. Finally, Sec. VI contains our conclusions. Appendix A reports some details on the spectrum of the space transfer matrix for short (finite) times.

II. THE MODEL
We consider the self-dual kicked Ising model [22,23], described by the following time-dependent Hamiltonian where δ p (t) = ∞ m=−∞ δ(t − mτ ) is the periodic delta function and Here τ is time interval between two kicks, L denotes the volume of the system, 1 x is the identity operator in (C 2 ) ⊗x , {σ α j } α=x,y,z are Pauli matrices at position j, and we impose The parameter h = (h 1 , . . . , h L ) describes a position dependent longitudinal field measured in units of τ −1 . From now on τ is set to 1 to simplify the notation. The Floquet operator generated by (1) reads as In Floquet systems it is customary to introduce quasienergies {ϕ n } defined as the phases of the eigenvalues of the Floquet operator. The quasienergies take values in the interval [0, 2π] and their number N = 2 L is the dimension of the Hilbert space where (1) acts, namely To characterise the distribution of quasienergies (and especially the correlations among them) it is convenient to consider the SFF This quantity represents an efficient diagnostic tool able to tell apart chaotic (non-integrable) systems from integrable ones even in the thermodynamic limit (L → ∞). Indeed, the former are believed to show uncorrelated (Poisson distributed) quasienergies [24] while the latter to display quasienergies distributed as in random unitary matrices [12][13][14][15][16][17][18][19][20][21][22]. In the first case the SFF (7) is independent of time, while it shows a linear ramp in the second. Importantly, the probability distribution of the SFF does not become a delta function in the thermodynamic limit [25,26] (this property is commonly referred to as "non-self-averaging" property [26]). This means that, to have a meaningful comparison with the prediction of RMT, one has to study the probability distribution of the SFF over an ensemble of systems. The ensemble can be formed by considering similar systems with different numerical values of the parameters or the same system at different times. Here we follow Ref. [22] and consider the distribution of (7) in an ensemble formed by self-dual kicked Ising models (1) with random longitudinal fields. Specifically we assume that the longitudinal magnetic fields at different spatial points h j are independently distributed Gaussian variables with mean valueh and variance σ 2 > 0. Differently from Ref. [22], however, here we are interested in the thermodynamic limit of all moments of the distribution of |tr[U t KI [h]]| 2 not just in the average. Namely we consider where the symbol E h [·] denotes the average over the longitudinal fields In this language the thermodynamic limit of the SFF corresponds to K 1 (t).

III. PREDICTION OF RANDOM MATRIX THEORY
Before computing (8) in the self-dual kicked Ising model we compute the moments for an ensemble of random unitary matrices subject to the same constraints -or symmetries -as the Floquet operator U KI [h] (cf. Eq. (5)). Indeed, due to some special symmetries of U KI [h], such ensemble is not the "standard" circular orthogonal ensemble (COE) -composed of symmetric unitary matrices. To see that let us start by reviewing the symmetries of U KI [h].

A. Symmetries of the Time-Evolution Operator
To analyse the symmetries of (5) it is convenient to make the following basis transformation This transformation leaves (8) invariant and brings the operator in a manifestly symmetric form SinceŪ KI [h] is unitary and symmetric we immediately have where (·) * denotes complex conjugation in the computational basis (the standard Pauli basis where both matrices σ x and σ z are real) and C is the anti-unitary operator implementing it in the Hilbert space. This is the most obvious anti-unitary symmetry of the time-evolution operator and corresponds to the standard time-reversal symmetry T (with T 2 = 1).
As observed in Ref. [27], however, T is not the only anti-unitary symmetry ofŪ KI [h]. Indeed, defining and noting one readily finds Reshaping (11) and (17) we will now see that they correspond to the constraints on random matrix ensembles associated to two compact symmetric spaces [28,29] (two different symmetric spaces will correspond to even and odd L). To see this, we note that, permuting the computational basis, F y can be brought to one of the two following block-diagonal forms, depending on the parity of L where P 1 P T 1 = P 2 P T 2 = 1 N , N = N /2, and {s j } N j=1 is a specific string of +1s and −1s.

L odd
The matrix (18) is a non-singular real skew-symmetric matrix (i.e. Ω N = −Ω T N ) multiplied by i L , this means that defininĝ we have that (11) and (17) becomê Namely the unitary matrixÛ KI [h] is constrained to be symmetric and symplectic. The compact symmetric space characterised by this constraint is and corresponds to CI in Cartan's classification [28,29].
Note that here we denoted by As shown in [28,29] matrices belonging to this symmetric space can be parametrized by The matrix on the r.h.s. of (19), instead, can be written as the square of Moreover, it can be brought to the following diagonal form by means of an orthogonal transformation P 3 So that defininĝ we havê This means that the unitary matrixÛ KI [h] is constrained to be real orthogonal and fulfil (29). The compact symmetric space characterised by these constraints is and corresponds to BDI in Cartan's classification [28,29]. As shown in [28,29] the matrices in this symmetric space can be parametrized by The random matrix ensembles corresponding to the symmetric spaces S − (N ) and S + (N ) have been introduced in Refs. [28,29]: for both ensembles one finds that the quasienergies come in pairs of opposite values {ϕ j , −ϕ j } N j=1 . Moreover, from the probability measure induced by the Riemannian metric of the symmetric spaces one finds the following (joint) probability distri- where the proportionality constant is chosen to ensure that their integral over [0, π] N is one. Changing variables from ϕ j to x j = cos ϕ j ∈ [−1, 1], both (32) and (33) are brought into the so-called Jacobi-ensemble form [30] with β = 1 and, respectively, a = b = 0 for S − (N ) and a = b = −1 for S + (N ).

C. Thermodynamic Limit of the Moments
Let us now turn to the main objective of this section: computing the moments of SFF (8) where the matrix U KI [h] is replaced by a random matrix in U ∈ S ± (N ) and the average E h [·] is replaced by where +/− are respectively chosen for L even/odd. To compute the moments (8) it is convenient to find the full probability distribution of the linear statistics where ϕ j are the quasienergies of U and in the last step we used that they can only appear in complex conjugated pairs. An immediate consequence of this relation is that, as opposed to what happens in the COE, the random variable T t,L is real.
Since we are interested in the thermodynamic limit (L → ∞) we do not need to find the statistics of T t,L exactly: it is sufficient to find its leading behaviour for large N = 2 L−1 . This can be efficiently done using "log-gas methods" [5,30,31], i.e. studying the statistical mechanics of quasienergies through the formal analogy with a gas of charged particles in two dimensions (confined in a one-dimensional domain). Specifically, here we will follow the treatment of Ref. [30].
Let us start by considering the average of (36). First we note that, introducing the n-point function of the density of quasienergies (37) the average can be expressed as Then we take the thermodynamic limit: using Proposition 3.6.3 of Ref. [30] (see also Exercises 14.2) with β = 1 and This result agrees with the infinite L limit of the exact one-point function in the Jacobi ensemble (cf. Proposition 6.3.3 of Ref. [30]) at the two special values of interest here and, in particular, implies where mod(n, m) = n mod m is the mod function. Next, we consider the variance where ρ c ±,2 (ϕ 1 , ϕ 2 ) is the connected two point function Using Eq. 14.56 of Ref. [30] with a(cos θ) = 2 cos(tθ) and At this point, introducing the probability distribution of we can use Eq. 14.68 of Ref. [30] and find The probability distribution (45) produces the following central moments in the thermodynamic limit (46) and therefore (8) read as The result (47) is very different from the one found for U ∈ COE. Indeed, in the latter case the expression (36) is complex and Ref. [25] found the following joint distribution for its real and imaginary part (respectively x and y) in the thermodynamic limit This distribution yields We see that, even though (49) and (47) agree for n = 1 and t odd, they are generically very different. In particular the moments (47) are much larger that (49) indicating that the fluctuations in the ensembles S ± (N ) are larger than those in the COE.

IV. LOWER BOUND FROM THE SPACE-TRANSFER-MATRIX APPROACH
Equipped with the random matrix theory prediction (47) we can now move on to our main goal: computing the moments K n (t) in the self-dual kicked Ising model. In this section we will determine a rigorous lower bound for K n (t).
1. An illustration of the n-th moment of the spectral form factor Kn(t). The lattice depicts a system of L spins that are propagated to time t. T2n acts as a transfer matrix on 2n copies of the lattice. The average (E h j ) is performed over the longitudinal magnetic fields hj. The loops on the edges of the lattice indicate that we need to compute the trace of T2n to get the n-th moment of the spectral form factor.

A. Transfer Matrix in Space
To derive the lower bound we will follow Ref. [22] and use the transfer matrix in space. The starting point is the following identity, which holds for the self dual kicked Ising model [22,23] Here ε is a vector with t entries equal to one andŨ KI [h] takes the form (5) with the only difference that the size L is replaced by t in (2) and (3). Note that the trace on the right hand side of Eq. (50) is over H t = (C 2 ) ⊗t . Equation (50) can be used to rewrite the n-th moment of the SFF as follows with T 2n ∈ End(H ⊗2n t ) defined as By looking at the graphical representation in Fig. IV where we introduced Note that hereŨ is the transfer matrix in space at the average magnetic field, and is the magnetisation (in the α direction) for a chain of length t.

B. Trace of U t KI [h]
Before embarking on the analysis of Eq. (51) it is useful to look at a simpler observable that can be studied with the same method, namely Indeed, the RMT prediction for this quantity is nontrivial (cf. Eq. (40)) and offers a convenient opportunity for testing the quantum chaos conjecture. Moreover, performing the calculation in this simple example will best illustrate some of the main ideas. Considering (59) and using (50) we have where in this case the space-transfer matrix reads as The limit (60) can be computed as follows. First we observe that the eigenvalues of the transfer matrix T are at most of unit magnitude and, additionally, geometric and algebraic multiplicity of any eigenvalue with magnitude one coincide. This can be seen by using the relation and reasoning as in the proof of Property 1 of Ref. [22]. This observation implies that B(t) is given by the number of eigenvectors |A corresponding to unimodular eigenvalues.
Next, we observe that -because of Eq. (62) -all unimodular eigenvalues of T lie in the eigenspace of O 1,0 corresponding to eigenvalue one. Given the form of the operator O 1,0 , this means that all relevant eigenvectors |A must be in the kernel of the operator M z , i.e.
This relation allows us to conclude the analysis of odd times. Indeed, since in that case there can be no vectors in the kernel of M z (a spin-1/2 chain of odd length cannot have zero magnetisation), we find immediately find that B(t) vanishes. To find the result for even t we continue by acting on |A with T, this yields This equation, together with (63), implies whereŨ is defined as in (14) but with L replaced by t.
The first of these equations can be verified by using the identitiesŨ while the second follows from Eq. (65) and (64).
Since the operatorŨ squares to 1 t we have A state that satisfies equations (65) and (66) is directly identified as with P i,j = 1 2 1 + 1 2 α σ α i σ α j being the transposition of the spins on sites i and j. In particular, it is easy to verify that (70) fulfils (65) and (66) with Assuming that (70) is the only eigenvector of T corresponding to unit magnitude eigenvalues we have which agrees with the RMT prediction (40). Note that, for even values of L, Eq. (72) gives a lower bound for B(t). Indeed, given the general structure (69) of the eigenvalues one can immediately see that the contribution of each eigenvalue to B(t) is always positive for L even.

Numerical Checks
The prediction (72) can be checked by finding numerically all unimodular eigenvalues of T for short times. The results for times up to t = 25 are shown in Tab. I. No eigenvectors are found for odd t while for even t the only eigenvalue is the one given in Eq. (71). The only exceptions are for t = 6 and t = 10. In these two cases we find an additional unit-magnitude eigenvalue and its corresponding eigenvectors have been identified in Ref. [22] (cf. Eqs. (171) and (175) of the Supplemental Material). As no other additional eigenvector can be found for t > 10 we conjecture that the presence of (73) is a short-time fluke.

C. Higher Moments of the Spectral Form Factor
Let us now move on to the main objective of this section and consider the moments (51). The steps to determine a lower bound for these quantities are similar to the ones taken in the previous subsection. In particular, a relation analogue to Eq. (62) still holds with T and O 1,0 replaced by T 2n and O n,n , namely As a consequence, the eigenvalues of T 2n have again magnitude bounded by one and those with unit magnitude have coinciding algebraic and geometric multiplicity. Another aspect that is unchanged is that the eigenvectors corresponding to the eigenvalues with unit magnitude belong to the eigenspace of O n,n with eigenvalue one. This immediately leads to the following two conditions on the relevant (i.e. corresponding to unit-magnitude eigenvalues) eigenvectors of T 2n Reasoning along the lines of the previous subsection, one can readily prove that (75)-(76) are equivalent to where we definedŨ Again, usingŨ 2 = 1 t , we have e iϕ = ±1.
To find a set of eigenvectors {|A } fulfilling (77)-(78) is useful to follow Ref. [22] and introduce the a state-tooperator map. This is implemented as follows. First we consider the coefficients A i1,...,i2n of |A in the basis where {|i } is the computational basis of H t . Namely Then, we define the operator A n in End(H ⊗n t ) by means of the following matrix elements i 1 · · · i n |A n |j 1 · · · j n = A i1,...,in,j1,...,jn . (82) In this way we can express the conditions (77) and (78) as The first observation is that, even tough both +1 and −1 are possible eigenvalues ofŨ , it is reasonable to restrict ourself to the case of positive eigenvalues. Indeed, as we will see in the following, negative eigenvalues are expected to be rare and appear only for small times. Moreover, considering only positive eigenvalues produces a lower bound for (51) if we only focus on even lengths. For this reason, we get rid of the contribution of negative eigenvalues by averaging the results for even and odd lengths, i.e. we definē A set of eigenvectors with eigenvalue one can be determined by finding the number of all linearly-independent operators that commute with the set of operators {U n , M α,n }. This set can be found by observing that, as shown in Ref. [22], the elements of the dihedral group G t commute with the set of operators {U, M α }. The group G t is a symmetry group of a polygon with t vertices and its elements be expressed as with Π denoting the periodic shift for one site and R reflection. These operators are represented in End(H ⊗n t ) as where P i,j is the transposition. The number of linearly independent elements of this representation of the dihedral group is [22] The above facts imply that any operator written as commutes with {U n , M α,n }. This means that the number of operators commuting with {U n , M α,n } is at least number of elements of the dihedral group to the power n. There is, however, an additional combinatorial prefactor that one should take into account to attain a tighter lower bound. The combinatorial prefactor arises from an arbitrariness in the definition (82) of the operator A. Indeed, it is easy to see that defining with τ, σ ∈ S n permutations of n elements, leads to operators fulfilling (96)-(97) for any τ and σ. These operators are not all linearly independent: since the set of all operators B (cf. (89)) is invariant under permutations of the copies in the tensor product, only A 1σ can be independent. This leads to a combinatorial prefactor n!. Such a combinatorial prefactor leads to a lower bound on the higher moments of the SFF that agrees with the standard COE prediction.
To see this we first note that considering the unitary mapping withF y,n ≡F y ⊗ · · · ⊗F y n (93) andF y,n defined as in (13) with L replaced with t, the conditions (77)-(77) become M α,n ⊗ 1 tn + 1 tn ⊗ M α,n |A = 0 , Mapping these into relations for operators by means of the definition (91) (with A replaced by A ) we then find we find that it fulfils (96)-(97) for all σ ∈ S 2n .
Taking again into account the invariance of the set {B} under permutations of the copies in the tensor product and noting that the set is also invariant under transposition in each single copy we obtain the following combinatorial prefactor Together with this additional factor a lower bound for K n (t) can then be expressed as We see that for odd times larger than 5 this bound agrees with the RMT prediction (47) and, therefore, we expect it to be tight. For even times we can find additional operators fulfilling (96)-(97) by considering |ψ ψ| with |ψ given in (70). In particular we find the following additional solutions with a combinatorial prefactor of Taking into account also these solutions we have that the bound agrees with the RMT prediction (47) for all times larger that 6.

Numerical Checks
The arguments of this section can again be tested (for short times) by identifying numerically all eigenvectors of the space-transfer matrix that have eigenvalues equal to ±1. Here we present an analysis of the simplest nontrivial case, i.e. n = 2. By repeatedly applying T 4 to a random state and then projecting to different fixedmomentum subspaces (power method) we enumerated all its unimodular eigenvectors up to t = 7: the results are gathered in Tab. II.
The first point to note is that negative eigenvalues are less common than positive ones. For odd times we did not find any eigenvalue −1. The next observation is that, as expected, the number of eigenvectors is much bigger than the standard COE prediction.
However, since we can only investigate the short-time behaviour, we observe some short-time effects that we believe will disappear for larger times. In particular, we observe two main phenomena. First, the number of linearly independent vectors in some subspaces is smaller than expected because vectors are "not long enough". In other words, for short times the operators identified in the previous section are not all linearly independent. Second, for short even times there are some additional eigenstates (similarly to what happens for t = 6 and t = 10 in Sec. IV B). Since these special states seem to appear only for even times we can avoid this complication by considering only odd times. The first phenomenon, however, remains also there. An example can be readily observed for t = 3. In this case we find only 59 eigenvectors with eigenvalue +1 even tough the lower bound from Eq. (100) predicts at least 75 of them. A similar effect can be seen for t = 7 where we found 587 eigenvectors, whereas the expected lower bound is higher by one. On the other hand at t = 5 the number of eigenvectors matches the predicted lower bound.
To obtain more detailed information we note that T 4 commutes with the four translation operators and count how many linearly independent eigenvectors with unit-magnitude eigenvalue exist in each subspace with fixed four-quasi-momentum {k 1 , k 2 , k 3 , k 4 } (see Appendix A for more details). By analysing the resultsreported in the Tables III-VIII -we identify the following general structure 1. The relevant eigenvectors appear in sectors where four momenta can be arranged into two pairs. Each pair (k 1 , k 2 ) contains two equal momenta k 1 = k 2 , or two momenta in the relation k 1 = t − k 2 ≡ −k 2 .
2. The number of linearly independent vectors in a sector is the same as the number of ways in which four momenta can be grouped into two pairs. This means that we can get the degeneracies one or three in a typical sector. For example The total number of vectors in a sector is therefore always given by the product of two numbers: the number of all possible pairs and that of all possible permutations of the momenta.
3. When a sector has momenta k/2 or 0, one gets independent contributions from even and odd reflection eigenspaces.
For short times, however, some of the reflection eigenspaces can vanish, or be smaller than expected. For example, at t = 7 in the reflection odd part of the sector with all four momenta equal to zero, we obtain only 2 independent vectors instead of the expected three. The same problem occurs for t = 3 in almost all sectors. The number of sectors where this happens decreased when t increases and this problem is expected to disappear for larger times.
It is interesting to check if by applying the above principles we can calculate the final result for the number of eigenvectors. For (large enough) odd times the result is exactly 12t 2 , while for even times we get 12t 2 + 12t + 1 (see Appendix A). Both results agree with the lower bound (100) and with the RMT prediction (47).

V. MONTE-CARLO SIMULATIONS
In this section we present numerical evidence substantiating the tightness of the bound (100). Our numerical results are obtained by means of simple Monte-Carlo simulations based on direct time propagation with U KI [h] followed by an average over different configurations of the longitudinal magnetic fields h j .
The trace of U t KI [h] is computed by restricting the sum to a set R containing m random states of C N . The states |r ∈ R are obtained by producing and normalising vectors with independent and identically distributed complex Gaussian random variables. The number of states m can be much smaller than 2 L and we expect fluctuations of the order O (1/ √ m). For example, for n = 2 the trace is approximated by and r 1 = r 2 = r 3 = r 4 . The results are obtained for finite-length chains and consequently the thermodynamic limit behaviour can only be observed for times t < L. Fig. 2 reports the results of the Monte-Carlo simulations for K 1 (t), K 2 (t) and K 3 (t). As we see these results indicate that the first, the second and the third moment of the SFF grow with time as predicted by Eq. (100). The corrections for even times can not be seen due to the fluctuations in the Monte-Carlo method. For comparison we also plotted the results for the time-reversal invariant dual-unitary circuits with random gates. The Floquet propagator has the form described in Ref. [32] (equations (23) and (24)) with J = 0 and We see that, unlike for the self-dual kicked Ising, the moments agree with the COE predictions.
Finally, in order to see whether all eigenvectors are identified, in Figure 3  well for all times except for t = 6 and L even. This might indicate that some additional eigenvectors with eigenvalues +1 and −1 are not identified. Other causes of disagreement might be finite-size corrections in the Monte-Carlo simulation or fluctuations due to the finite number of realisations.

VI. CONCLUSIONS
In this paper we computed the statistics of the spectral form factor in the self-dual kicked Ising model. Our strategy has been to establish a rigorous lower bound on the higher moments (generalising the space transfer matrix method of Ref. [22]) and to check its saturation numerically (via Monte-Carlo simulations). We found that, even though the spectral form factor takes the standard COE form, the fluctuations are consistently higher. We explained this result by noting that, since the self dual kicked Ising model has two anti-unitary symmetries [27], the relevant random matrix ensemble is not the COE but is defined on a more restricted symmetric space. We found that this space is either Sp(N )/U (N ) or O(2N )/O(N )×O(N ) depending on the parity of the number of sites.
Our work suggests several possible directions for future research. An obvious one is to prove rigorously the findings of this paper in the spirit of Ref. [22]. Namely, devise a mathematical proof of the bound's saturation. Our numerical analysis of the short time behaviour suggests that such a proof is concretely within reach, at least in the case of odd times.
Moreover, it is interesting to apply the method adopted here to the study of the spectral-form-factor statistics in other systems. Our numerical results, together with recent compelling analytical evidence [32][33][34][35][36][37][38][39][40], suggest that dual-unitary circuits [32] provide a very convenient framework where these questions can be investigated analytically. Indeed, preliminary results indicate that all circuits in this class are characterised by a vanishing Thouless time, meaning that there is no characteristic time scale other than the Heisenberg time given by the dimension of the Hilbert space. In fact, they seem to provide an arena where one can generate many-different random matrix ensembles by including increasingly more anti-unitary symmetries in the local gates.
Finally, it is interesting to ask whether the method of this work can be successfully applied to "generic systems" with non-unitary space transfer matrix. There a meaningful comparison with RMT can only be performed in a finite volume due to a Thouless time increasing monotonically with the volume [16,19].

VII. ACKNOWLEDGMENTS
All authors have been supported by the EU Horizon 2020 program through the ERC Advanced Grant OMNES No. 694544, and by the Slovenian Research Agency (ARRS) under the Programme P1-0402. TP acknowledges a fruitful discussion with Nick Hunter-Jones in the preliminary stage of this work.