Level statistics of the one-dimensional ionic Hubbard model

In this work we analyze the spectral level statistics of the one-dimensional ionic Hubbard model, the Hubbard model with an alternating on-site potential. In particular, we focus on the statistics of the gap ratios between consecutive energy levels. This quantity is often used in order to signal whether a many-body system is integrable or chaotic. A chaotic system has typically the statistics of a Gaussian ensemble of random matrices while the spectral properties of the integrable system follow a Poisson statistics. We find that whereas the Hubbard model without alternating potential is known to be integrable and its spectral properties follow a Poissonian statistics, the presence of an alternating potential causes a drastic change in the spectral properties which resemble the one of a Gaussian ensemble of random matrices. However, to uncover this behavior one has to separately consider the blocks of all symmetries of the ionic Hubbard model.


I. INTRODUCTION
In the past years a significant effort has been devoted to the understanding of the dynamics of quantum systems. In contrast to the equilibrium physics at low temperature which is typically dominated by the low energy properties of a quantum system, the non-equilibrium physics relies often on properties of an arbitrary part of the spectrum. A central question is the determination of the conditions under which an isolated quantum many-body system thermalizes or fails to thermalize [1][2][3][4][5]. In this context it is generally believed that a generic chaotic system is best suited to exhibit thermalization towards a suitable statistical mechanics ensemble [6].
Recent studies have put forward several examples of nonequilibrium phenomena, which provide examples of systems failing to thermalize, such as many-body localization [4,5,[7][8][9], and Hilbert space fragmentation [10,11]. These systems feature a large number of (almost) conserved quantities and are thus not chaotic. The phenomenon of quantum many-body scarring [12][13][14][15] highlights that even in overall chaotic systems tiny subspaces can be found, which fail to thermalize. This scenario is especially relevant if the initial conditions lie in this subspace [16].
A powerful probe of the properties of a many-body quantum system in the context of quantum chaos are the universal properties of the level statistics. There exist two cornerstone limits: the Poisson statistics and the statistics stemming from random matrix theory [17]. The random matrix theory statistics is conjectured to hold for generic chaotic systems [18,19], where the energy levels repel each other. The exact class of random matrices is determined by the symmetries present in the system. In contrast, the Poisson statistics is found for quantum integrable systems.
The level statistics of different many-body quantum systems have been classified as for example of the prime examples of the Hubbard model [20] and variants of it [21], or its bosonic counterpart [22], a kicked-parameter model of spinless fermions [23], or more recently a family of Sachdev-Ye-Kitaev models [24].
The Hubbard model is one of the most studied models in condensed matter physics, since it is one of the simplest models containing the competition between the kinetic term and the interaction term of fermions. A large amount of work has been devoted to study the properties of the Hubbard model [25] and in particular of its one-dimensional version [26,27]. The one-dimensional model has the particularity that it is Bethe-Ansatz integrable which enabled a lot of detailed studies [26]. In agreement with the conjecture mentioned before, the level statistics of the one-dimensional Hubbard model was found to follow the Poisson statistics [20], while e.g. the twodimensional Hubbard model at low filling exhibits the statistics of Gaussian orthogonal ensemble (GOE) [28].
In this work, we will investigate the spectral properties of a generalization of the one-dimensional Hubbard model by adding an alternating local potential term. Such an ionic Hubbard model has been realized in cold atomic gases using a superlattice potential [29,30] and proposed to describe GeSe [31]. Further, the ionic Hubbard model has been devised in order to study the physics occurring at neutral-ionic phase transitions as they occur in solids for example, ionic to neutral transitions in organic charge-transfer solids [32,33] and at ferroelectric transitions in perovskites [34]. An alternating local potential is also at the heart of the staggered fermion formulation of massive fermions in quantum field theories on a lattice [35], and is of relevance to ongoing efforts to simulate quantum field theories.
Such an alternating potential is widely believed to break the integrability of the Hubbard model. However, recently, numerical results pointing towards a Poisson statistics of the levels [21] were reported. The tension with the general expectation motivated us to investigate the level statistics properties of this model in more detail. Indeed we find that after all the symmetries of the Hamiltonian are taken into account, the level statistics are clearly of GOE nature in a wide range of parameters considered.
In Sec. II we describe the one-dimensional ionic Hubbard model considered in this work. In Sec. III we present the level statistics approach we use to show the chaotic character of the ionic Hubbard model. In order to employ this procedure we first need to identify all symmetries of the model, which we present in Sec. IV. We determine the eigenenergies of the ionic Hubbard model by performing complete numerical exact diagonalization using the identified symmetries. Our results regarding the spectral properties for generic and half filling are presented in Sec. V.

II. MODEL
The Hubbard model is one of the standard models of condensed matter physics. It is the simplest model which describes the competition between the motion of two species of fermions, called spin up and spin down fermions, in a periodic lattice potential and their on-site interaction. The onedimensional Hubbard model is given by the Hamiltonian Here c j,σ are the fermionic operators describing a fermion with spin σ =↑, ↓ on site j. J is the tunneling amplitude of the fermions from one lattice site to the neighboring ones, and U the strength of the on-site interaction between the different fermionic species. L denotes the length of the chain, which is chosen to be an even number. The total particle number of the fermions is denoted by N . We note that we focus on even N , therefore the total spin is integer. We use in the following periodic boundary conditions, i.e. we identify c 1,σ . Depending on the parameters chosen, the one-dimensional Hubbard model has many interesting phases [26,27]. These reach from a metallic/Luttinger liquid phase (U = 0, and U > 0 away of half filling) over a band insulating phase for the totally filled case, a Mott insulating phase at half filling for U > 0 to a superfluid phase for the attractive interaction U < 0. Here we mentioned only the phases occuring for a spin-balanced situation.
One experimentally relevant extension of the Hubbard model is the so-called ionic Hubbard model which has an additional alternating local potential. Its Hamiltonian is given by where η gives the strength of the alternating potential. The ground state phase diagram of the ionic Hubbard model has attracted a lot of interest  and many interesting phases have been pointed out. In particular, the very interesting phase of a bond order arises. These results have been obtained using either approximative methods or numerical methods because of the common believe that the ionic Hubbard model for generic couplings (i.e. away from the known integrable or free parameter sets) is not Bethe-Ansatz solvable.

III. LEVEL STATISTICS
In this section we would like to summarize some properties of the spectra of quantum many-body models. It has been conjectured that the spectral statistics of quantum manybody systems displays either universal features described by random matrix theory (RMT) for chaotic quantum systems [18,19,58,59] (if all known symmetries are removed), or follows Poisson statistics for quantum integrable systems. There exist numerous examples in the literature that show the usefulness of the spectral analysis in order to obtain information on the chaoticity or integrability of many-body quantum models [9,20,[22][23][24][60][61][62][63].
In order to quantify the spectral properties of a model one important quantity is the distance between adjacent manybody eigenvalues where {E n } are the eigenvalues of the Hamiltonian sorted in an ascending order. In the case of integrable models, for which one can find an extensive number of conserved quantities, the distribution of the level spacings should follow a Poisson distribution where ∆ is the mean level spacing. In contrast, for a chaotic model the underlying symmetries, as the time-reversal symmetry, determine with which random matrix ensemble the model shares the same universal features. Similar to Hubbard model the ionic Hubbard Hamiltonian considered here (Eq. 2) is invariant under time reversal operator T and has rotational symmetry which leads to the Gaussian orthogonal ensemble (GOE) in random matrix theory [59]. For the GOE the level spacing distribution has the Wigner-Dyson form Numerically, often an alternative way is employed in order to characterize the spectral properties. One considers the behavior of the gap ratio for consecutive levels defined by [62] r n = min (δ n , δ n+1 ) max (δ n , δ n+1 ) .
This approach has the advantage of bypassing the need to compute the density of states. The computation of the density of states would imply a procedure for the unfolding of the spectrum, which can lead to inaccurate numerical results [65].
The probability distribution of the consecutive gap ratios for the Poisson statistics is given by and for the GOE statistics by [66] P GOE (r) = 27 4 r + r 2 (1 + r + r 2 ) 5/2 .
One can quantify the proximity of the computed distribution of gap ratios to Poisson or GOE by using the mean value r or the value P (r = 0). The expected values for the Poisson and GOE distributions are given in Table I The comparison of the spectral properties of a quantum chaotic Hamiltonian with the ones of the RMT is typically performed within each symmetry sector of the symmetries present in the Hamiltonian. These symmetries comprise for example spatial translations and reflections, but also more special symmetries as particle hole transformations. In the case in which a few symmetries are known for the Hamiltonian, one typically block diagonalizes the Hamiltonian in order to perform the analysis of the level statistics. However, in certain situations the symmetries might not be known or the block diagonalization procedure may be impractical or not possible. Several work investigated different situations where additional symmetries are present [67][68][69][70].
In particular, Giraud and coworkers [67] present analytical results for the distributions of random matrices comprised of several independent symmetry blocks, corresponding to the presence of a discrete symmetry. In particular, we make use of the results of Ref. [67] for the GOE distribution for m = 2 or 3 blocks of equal size, to show the chaotic nature of the halffilled ionic Hubbard model. The mean value r and the value P (r = 0) for the GOE distributions with m = 2 or 3 are given in Table I and the distributions are plotted in Fig. 7. The analytical expression of the distributions can be found in the Supplementary Material of Ref. [67]. One can observe that by having multiple symmetry sectors the value P (r = 0) becomes finite and increases for larger m (see Table I). The mean value r is already closer to the value expected for the Poisson distribution, Eq. (7), than the one for the GOE distribution, Eq. (8), even for just m = 2 or 3 blocks. Whereas the distribution of consecutive level spacings for m = 2 shows very pronounced leveling off at low values as r → 0, the distribution of m = 3 already shows a behavior very similar to the Poisson distribution. By taking a very large number of symmetry sectors, m → ∞, one will recover the Poisson distribution [67]. Thus, extreme care has to be taken in the numerical analysis of the spectral properties of many body quantum model, since missing just a few symmetries the numerical results for a chaotic system could resemble an integrable one. We demonstrate this explicitly for the ionic Hubbard model in the Appendix A.
We note that in the following we mainly make use of the results for the GOE distributions for which the m symmetry blocks are of equal size. However, one can perform the same analysis also for the case of unequal blocks [67], as we use it in Appendix A. In the case of m = 2, as the ratio of the sizes of the two blocks goes from 0 to 1 the values r and P (r = 0) interpolate between the values for m = 1 and m = 2.

IV. SYMMETRIES
In order to study the level statistics of the ionic Hubbard model we need to identify the discrete symmetries existing in the model. In the following we list the symmetries (beside the time-reversal symmetry discussed already in the previous section) occurring in the ionic Hubbard model with periodic boundary conditions and comment on how we take them into account.
a. Gauge symmetry: The global U (1) gauge symmetry where φ is a real number, leads to the conservation of the particle number N σ with spin σ and consequently also of the total particle number N = N ↑ + N ↓ . This conservation is directly implemented in the numerical exact diagonalization method. b. Spin rotation symmetry: The SU (2) spin rotation symmetry, well known from the Hubbard model, also exists in the presence of an alternating potential. The spin operators are defined by where σ α are the Pauli matrices. These spin operators represent a SU(2) Lie algebra and generate the rotations in spin space. As all spin operators commute with the Hamiltonian of the ionic Hubbard model in Eq.
(2), [H, S α ] = 0. Thus, the Hamiltonian is rotationally invariant in spin space which corresponds to the conservation of the total spin S. The symmetry sectors corresponding to this symmetry are labeled by |s, m z , where S 2 |s, m z = s(s + 1) |s, m z and S z |s, m z = m z |s, m z (whereh = 1). The conservation of the S z component is directly implemented in the numerical exact diagonalization method. However, the symmetry sectors corresponding to the absolute spin S are reconstructed from the numerical results. For more explanation on how this is performed see Appendix B. c. Translational symmetry: The lattice in the ionic Hubbard model has a two site unit cell. Therefore, a translational symmetry with respect to a translation by an even number of lattice sites exists, i.e.
The corresponding conserved (dimensionless) momentum k determine the eigenvalues exp(ik) of the generator T x .
The possible momentum values are k = 2π L/2 j with j = 0, 1, . . . , L 2 − 1, for a chain with an even number L of sites. Thus, each symmetry sector is characterized by its dimensionless momentum k.
d. Parity (reflection) symmetry: In the ionic Hubbard model there exists a reflection symmetry given by the transformation This transformation also leaves the Hamiltonian invariant.
Here, the reflection center is chosen as the first site of the chain. Due to the translational symmetry, the reflection center can be chosen on any site of the chain (but not on a bond center). The eigenvalues corresponding to the reflection symmetry are p = ±1. Note that all symmetry generators besides (T x , P ) commute. As the reflection symmetry changes the direction of the momentum, it does not commute with the translational symmetry in general, but only if the momentum is k = 0 or k = π. The momentum k = π only exists if L 4 is an integer number. Therefore, we cannot diagonalize the matrix with respect to both the translation and reflection symmetry beside in the special case of k = 0 or k = π. We choose to use the translation symmetry first and the reflection symmetry only for the special case of k = 0, π.
Combining all symmetries, the individual symmetry blocks are thus characterized by (N, s, m z , k, p). The last quantum number is only present for k = 0, π. Since the level statistics are independent of the respective symmetry sectors, if not stated otherwise, the distributions are calculated within all computable sectors individually and are only afterwards combined to decrease the statistical error of the distributions.
e. Particle-hole like symmetry: In the ionic Hubbard model an additional symmetry is present at half-filling. This additional symmetry resembles a particle hole symmetry, here between two sublattices, i.e. with the transformation The ionic Hubbard Hamiltonian at half filling is invariant under this transformation. We discuss in the result section how this symmetry influences the properties of the level statistics. f. Symmetries used in exact diagonalization In the exact diagonalization code we build subsectors for fixed N , m z , k and reflection quantum number, if applicable. It is tedious to combine the particle hole-like symmetry with the other symmetries in the existing code, and we did therefore not implement it. However, in the subsequent analysis we are able to detect the presence of this symmetry in the resulting level statistics distributions and to account for it in the analysis, and additionally we found a way to explicitly break the particle hole-like symmetry by alternating Hubbard on-site interactions, while retaining all other symmetries.
In the following we diagonalize individual blocks of sizes up to ∼ 1.27 × 10 4 . For the averages over all the symmetry blocks up to ∼ 3.84 × 10 5 eigenvalues were taken into account.

V. RESULTS
In this section we show our results on the spectral properties of the ionic Hubbard model and an extension using alternating interactions. We start discussing the results at quarter filling (N/L = 1/2) and at filling N/L = 2/3 as typical fillings (Sec. V A), only the standard symmetries (Sec. IV (a)-(d)) are present for this case. We recover the GOE like behavior which points towards the chaotic character of the ionic Hubbard model. Additionally, we discuss in Sec. V B the case of half filling, which is special due to the extra particle hole like symmetry (Sec. IV (e)). We show that also in this case, the results are in agreement with a GOE like behavior occurring within the different symmetry sectors.

A. Away from half filling
We start considering the ionic Hubbard model, Eq. (2), at quarter filling, i.e. N/L = 1/2. We compute the energy levels of the Hamiltonian in each block of the symmetries discussed in Sec. IV and consider the level distribution separately within these blocks. The data obtained from different symmetry sectors at quarter filling is then assembled in a single histogram. In the representation of the histograms the chosen size of the bins plays a role. We choose the bin size in order to minimize the statistical fluctuations and still obtain a good representation of the distribution. Since the distributions vary depending on the Hamiltonian parameters, we adapted the bin size correspondingly.
In Fig. 1(a) we present the distribution of the ratios of consecutive level-spacings averaged over all symmetry blocks for the standard Hubbard model, i.e. η/J = 0, as a benchmark of our procedure. For the standard Hubbard model, the ratios of consecutive level spacings are known to follow Poissonian level statistics [20]. For the chosen intermediate value of the interaction, U/J = 2, we find a distribution which follows nicely the Poisson prediction, as can be seen from the good overall agreement between the numerical histogram and the Poisson prediction for P (r), including the correct limiting behavior as r → 0 + . Furthermore the mean value r = 0.386 of the numerical data, which does not depend on the number of bins, agrees with the Poisson distribution value in the first three digits.
Once we switch on the alternating potential of the ionic Hubbard model, the behavior is changed drastically. In Fig. 1(b) the distribution P (r) from the numerical data is shown for a value of the alternating potential of η/J = 1 following closely the distribution expected for the GOE random matrix theory ensemble. The behavior for small r, where P (r) clearly starts at 0 and rises linearly for small values of r, is in stark contrast with the Poisson behavior and is a consequence of the level repulsion of chaotic systems. The mean value r = 0.515 is also rather close to the GOE prediction 0.536.
In order to investigate in more detail the behavior as a function of the amplitude of the alternating potential η/J and the dependence on U/J, we use the average value of the ratio of the consecutive level-spacings r . For the Poissonian distribution the mean value of this ratio is given by r P = 0.386 and for the GOE the mean value is r GOE = 0.536 (see Table I). The dependence on η/J for fixed U/J is shown in Fig. 2 (a-c) for different interaction strength U/J = 2, 4. As discussed above, at η = 0 the value of r is very close to r P of the Poisson distribution. We observe that for the smallest non-zero values of η we consider and for all considered interactions, the average value of r increases very rapidly to almost the value expected for a GOE ensemble. A very good agreement with the expected value of the GOE is found in particular for 0.5 < ∼ η/J < ∼ 2. Since, in principle one expects an infinitesimal value of η in the thermodynamic limit to follow the GOE like distributions, we investigate the transition from the Poisson distribution at η = 0 to the GOE distribution for η > 0 more carefully and study the numerical distributions at very small values of η. Fig. 3 shows the distribution of level spacing for four different values of η for a finite system of size L = 12, at quarter filling for U/J = 2. The main difference between the Poisson distribution and the GOE like distribution is the finite value arising in the Poisson distribution at r = 0 and its suppression due to the level repulsion in the GOE like ensemble. This shows    The GOE distribution is expected for any finite value of the alternating potential (η/J = 0) for sufficiently large system sizes. We see in Fig. 4 that the suppression of the distribution at low values of r is very sensible to system size. Whereas for L = 12 still the distribution for low r lies considerably above the GOE distribution, this deviation is already reduced for L = 16 and the distribution lies much closer to the GOE distribution. This hints that already at small values of the alternating potential the level distribution follows the predictions from the GOE if sufficiently large system sizes are considered.
Also for large values of approximately η/J > ∼ 2 the numerically found average value of r drops below the value of the GOE ensemble (cf. Fig. 2). We attribute the deviations from the expected GOE ensemble to finite size effects as discussed for small values of η/J. This is supported by the results obtained for larger system sizes in Fig. 2 moving towards the GOE predictions with increasing system size both at small and large η/J. The deviations between the different system lengths are in particular large for low and large values of the alternating potential, where the GOE expectation is not yet met.
The dependence of the average value of r on the interaction strength is shown in Fig. 5 for two different value of η. A drastic rise is clearly seen from the non-interacting integrable case at U = 0 which lies close to the predictions for the Poisson distribution to a finite interaction strength. Already at interaction strength of about U/J ≈ 1 an almost steady value is reached for η/J = 1, which agrees well with the predictions for a GOE ensemble even for a smaller system size considered [see Fig. 5(a)]. At larger value of the alternating potential, η/J = 2 similar behavior is observed but a larger system size of L = 16 is needed to reach the expected GOE value due to finite size effects as discussed above [see Fig. 5 At large values of the interaction we would expect again that deviations from the GOE ensemble arise due to the large separations occurring in the energy scales.
In order to show that the discussed behavior is generic for different fillings, we show in Fig. 6 the dependence of the average value of r on both the alternating potential η and the interaction strength U for N/L = 2/3. We find that the shown behavior resembles very much the behavior of the quarter filling. The only difference is that the deviations from the GOE value at a larger values of the alternating potential are smaller than for quarter filling.
We therefore conclude this section, that the energy spectra of the ionic Hubbard model at a generic filling (here N/L = 1/2 and N/L = 2/3) nicely follow the expected behavior of a GOE ensemble and therefore point towards a chaotic nature of the model.  For the case of half filling N/L = 1, the obtained distributions change drastically in comparison to what we described   Table I).
for quarter filling, as seen in Fig. 7. This motivated us to investigate this filling in more detail. From our symmetry considerations, see Sec. IV, we know that in the case of half filling an additional symmetry, Eq. (13), is present in the ionic Hubbard model. As it was pointed out in Ref. [67] such an additional symmetry changes the behavior seen in the spectra drastically, if one does not consider the different symmetry blocks independently. As seen in Sec. III depending on how many symmetries blocks are considered together, the distribution of the consecutive level spacing changes from the GOE like distribution to a more Poisson like distribution [67]. In Fig. 7 the distribution of the consecutive level spacings is shown for the case of half filling, where we have separated the spectrum into the blocks of the symmetries, Sec. IV (a)-(d), but not into the additional particle hole symmetry, Sec. IV (e). Note that resulting distribution neither follows a GOE like distribution nor a Poisson like distribution, but lies somewhere in between. These distributions are genuinely in between GOE and Poissonian, and do not converge to either of them with larger system sizes, unlike the crossover behavior discussed in the previous section. Using the predictions from Ref. [67] we confront the obtained distribution to the ones which correspond to the presence of two (m = 2) or three (m = 3) equally sized symmetry blocks (see Sec. III). The numerical results lie very close to the results for a GOE RMT with two symmetry blocks present, as expected from our symmetry analysis. This is further supported considering the mean value of r and its dependence on the value of the alternating potential shown in Fig. 8 (a). As for the case of quarter filling, at L = 12 intermediate value of 0.5 < ∼ η/J < ∼ 2 lead to an excellent agreement with the predicted value of the one predicted for the distribution GOE m = 2. Larger deviations occur at larger values of η. We attribute these again to finite size effects, as we can show that the longer system sizes considered strongly approach the expected mean value compared to the smaller system sizes. The average value of r for the system size L = 12 already agrees nicely up to η/J ≈ 2 to the expected value fo the distribution of the GOE with m = 2.
In Fig. 8 (b) the average value and its dependence on different values of U/J is shown. For U/J > 0.5, the obtained main value fluctuates very closely around the one predicted for the distribution GOE m = 2 and is distinct from the one predicted for the distribution with m = 3 equal symmetry blocks.
To conclude our investigation at half-filling, let us note, that the form of the distribution drastically changes and cannot be distinguished from the Poisson distribution, if the level statistics is not separated with respect to the symmetries, Sec. IV (a)-(d), as shown in Appendix A.
Thus, our entire results indicate that the ionic Hubbard model shows similar properties as the GOE like ensemble, if all known symmetries are considered.
C. Half filling -breaking the particle-hole symmetry In order to further corroborate our finding of the importance of the particle-hole like symmetry at half filling, we introduce an additional term in the Hamiltonian which explicitly breaks this symmetry. This is achieved by substituting the homogeneous on-site interaction with an interaction using an additional alternation in the on-site interaction term: n 2j,↑ n 2j,↓ + U 1 L/2 j=1 n 2j+1,↑ n 2j+1,↓ (14) This term breaks the particle-hole like symmetry, but leaves all the other discussed symmetries (a)-(d), such as the reflection, translation and spin rotation symmetry invariant. As shown in Fig. 9 the presence of such a symmetry breaking term leads to a distribution which follows very clearly the GOE distribution for a single (m = 1) symmetry block. In order to investigate this behavior in more detail we show in Fig. 10 the average value of r with the symmetry breaking term, i.e. with the difference in the interaction strength U 1 − U 0 . A very steep rise of the average value of r is seen at small amplitudes of the symmetry breaking term, but already at (U 1 −U 0 )/J ≈ 0.2, the average value stabilizes close to the expected value of the GOE m = 1 ensemble. We expect that this rise becomes more and more steep with increasing system size signaling a direct change of the behavior from the m = 2 situation for the presence of the particle-hole like symmetry to the standard GOE (m = 1) behavior when this additional symmetry is explicitly broken.   Table I). System size N = 10, L = 10 considered symmetry sector mz = 0 combining all s and k, p symmetry sectors.

VI. CONCLUSION
In this work, we have investigated the energy level statistics of the ionic Hubbard model and a generalization which breaks the particle-hole like symmetry at half filling. We find that in general the ionic Hubbard model and its generalization show, when all known symmetries are considered, the spectral features of a GOE ensemble. This implies the chaotic behavior of the level statistics of the ionic Hubbard model. Following the conjecture that the GOE like statistics indicates nonintegrability of a many-body Hamiltonian [20], to which to our knowledge no generic counterexamples are known, our results indicate the non-integrability of the ionic Hubbard model and its generalization. Our findings are in agreement with the general belief that the ionic Hubbard model is, for generic parameters, a non-integrable model. In this appendix, we show how neglecting the separation into the different symmetry sectors the in principle GOE like distribution of the consecutive level spacings approaches the Poisson-like distribution. In particular, we focus on half fill-ing, for parameters close to the ones considered in Hosseinzadeh et al. [21] where the Poisson statistics was identified.
To emphasize the effects of the separation into the different symmetry sectors, we present in Fig. A1 the distribution of the consecutive level spacings at half filling, but only taking a fixed S z sector, and not a fixed sector of the total spin S. In this case the expected distribution for a chaotic system would be a GOE with m = 10 symmetry blocks [67]. The number of blocks results from the five symmetry sectors due to the total spin, which are of unequal size, which should furthermore be split in half due to the particle-hole symmetry present at half filling (see Sec. V B). The sizes of the symmetry sectors are: s = 0, 3880 states; s = 1, 5940 states; s = 2, 2475 states; s = 3, 385 states; s = 4, 20 states. We can notice in Fig. A1 that this GOE distribution with m = 10 is very close to a Poisson distribution due to the large number of symmetry blocks. Our numerical data is in agreement with this distribution and thus, also lies very close to the Poisson distribution. This again underlines the importance of the separation of the spectra into the different symmetry sectors.

Appendix B: The implementation of the spin rotation symmetries
As explained in the main text the quantum numbers associated to the spin rotation are the total spin S with quantum number s and the spin in z-direction S z with quantum number m z (which varies from −s to s) with the corresponding eigenvector in each symmetry sector labeled as |s, m z . In our computation only the quantum number of S z is implemented but one can take the advantage of the degeneracy of the eigenvalues for a fixed total spin s and different m z 's. The eigenvectors |s, m z and |s, m z are connected to each other with the ladder operators S + and S − . As the Hamiltonian commutes with the ladder operators, the eigenvectors |s, m z and |s, m z are degenerate. Thus, knowing all eigenvalues for different m z 's one can pick the ones for the desired quantum numbers s. We start with the symmetry block with largest possible value of m z (for N spin-half particles this value is N/2) which corresponds only to one total spin s = m z and the solution is clear. In the next symmetry sector with smaller quantum number m z = s − 1 there exist the eigenvalues for two quantum numbers, s and s − 1, the ones with total spin s are the degenerate ones with the previous sector we looked at and the remaining eigenvalues correspond to the symmetry sector with total spin s − 1. One can decrease the quantum number m z to zero and and with this procedure find all eigenvalues labeled with both quantum numbers s and m z . Note, that the eigenvalues corresponding to the negative spin in zdirection (m z < 0 ) for a fixed total spin s are degenerate with the ones with its positive value (|m z |).