Many-body spectral statistics of relativistic quantum billiard systems

In the ﬁeld of quantum chaos, spectral statistics is one of the most extensively investigated characteristics. Despite a large body of existing literature, the effects of many-body interactions on the spectral statistics of relativistic quantum systems remain poorly understood. Treating electron-electron interactions through the one-orbital mean-ﬁeld Hubbard model, we address this fundamental issue using graphene billiards with the geometric shape of a circular sector as prototypical systems. Our approach is to consider the two characteristically different cases where the statistics are Poisson and Gaussian orthogonal ensemble (GOE) so the corresponding classical dynamics are typically integrable and chaotic, respectively, and to systematically investigate how the statistics change as the Hubbard interaction strength increases from zero. We ﬁnd that, for energies near the Dirac point, the Hubbard interactions have a signiﬁcant effect on the spectral statistics. Regardless of the type of spectral statistics to begin with, increasing the Hubbard interaction strength up to a critical value causes the statistics to approach GOE, rendering more applicable the random matrix theory. As the interaction strength increases beyond the critical value, the statistics evolve toward Poisson, due to the emergence of an energy gap rendering the quasiparticles massive. We also ﬁnd that the energy levels above and below the Dirac point can exhibit different statistics, and the many-body interactions have little effect on the statistics for levels far from the Dirac point. These results reveal the intriguing interplay between many-body interactions and the spectral statistics, which we develop a physical picture to understand.


I. INTRODUCTION
Quantum chaos is the interdisciplinary field investigating the quantum manifestations or signatures of chaos that arise in the classical limit [1][2][3].In real-world quantum systems, many-body interactions are ubiquitous, and their effects on measurable quantities need to be identified and understood.In quantum chaos, the vast majority of studies in the existing literature was carried out in the single-particle framework without considering many-body interactions.Nonetheless, there have been works on the effects of many-body interactions on the energy level spacing statistics in nonrelativistic quantum systems described by the Schrödinger equation [4][5][6][7][8][9][10][11][12][13][14].Many-body interactions have also been studied in quantum thermalization [15,16] using the spin-chain models [17], and there were sporadic works on the many-body effects through a mean-field potential in quantum billiard systems [10,18] where the classical orbits play a nontrivial role in the spectral statistics.Despite these works, the effects of many-body interactions on basic quantities in the field of relativistic quantum chaos [19][20][21][22] have remained poorly understood, which represents a knowledge gap in basic physics.This problem also has practical implications.Starting from the separation of graphene two decades ago [23][24][25][26][27], the important role of relativistic quantum mechanics in solid-state systems has been recognized, as the low-energy excitations in graphene and a plethora of other two-dimensional (2D) Dirac materials are governed by the 2D Dirac equation.Understanding the interplay among many-body interactions, relativistic quantum mechanics, and classical dynamics thus has both basic and applied values.
In the single-particle framework in which the many-body interactions are completely ignored, the universal behaviors of the spectral statistics of nonrelativistic quantum systems with respect to the generic classical dynamics have been well documented and understood.For fully chaotic systems with time-reversal symmetry in the classical limit, the level spacing statistics are equivalent to those of random matrices from the Gaussian orthogonal ensemble (GOE) [28][29][30].For classically integrable dynamics, the spectral statistics are Poisson [31], due to the independence of quantum numbers associated with different degrees of freedom and the consequent loss of any correlation between adjacent energy levels.Previous works also revealed that the correspondences between the nature of the classical dynamics and the GOE or Poisson statistics hold for relativistic quantum systems [19,[32][33][34][35][36][37].However, exceptions to the correspondences were recently uncovered where, for circular sector billiards filled with graphene, the spectra of energy levels close to the Dirac point follow GOE [38], despite the classically dynamics being integrable.This is surprising because the conventional wisdom from the studies of nonrelativistic quantum billiards stipulates that the level spacing statistics should be Poisson in this case.The physical reason for the emergence of GOE statistics in circular sector graphene billiard systems can be attributed to the complicated boundary conditions of the spinor wave functions that break the integrability of the corresponding classical dynamics [38].
In this work, we investigate the effects of electron-electron interactions on the spectral properties of relativistic quantum systems, using circular sector graphene billiards as a prototypical model.Our concrete technical approach is to consider two distinct types of circular sectors whose spectral statistics are Poisson and GOE, respectively, in the absence of manybody interactions [38], and to systematically investigate if and how characteristic changes in the statistics occur as the electron-electron interaction strength increases from zero.To accomplish this task, an essential requirement is to incorporate the many-body interactions into the graphene billiard system.We use the Hubbard model [39], which was proposed for understanding the transition between conducting and insulating phases [40].In the model, the electron-electron interactions are approximated by a screened Coulomb potential, which results in a conventional tight-binding model obeying the Pauli exclusion principle [41,42].Despite its introduction more than half a century ago, the Hubbard model and its extensions remain an active research area in condensed matter physics due to its high relevance to frontier areas such as heavy-fermion materials [43,44], high-temperature superconductors [45][46][47], and edge magnetism in graphene [48][49][50][51][52].The Hubbard model represents a theoretical and computation paradigm to gain physical insights into many-body interactions in various physical systems.
Our work and the main contributions can be described as follows.We calculate the full set of eigenenergies and eigenstates for circular sector graphene billiards by diagonalizing the mean-field Hubbard Hamiltonian through an iterative approach and obtain the spectral statistics in different energy ranges.In general, the spectral properties depend on the Hubbard interactions.To identify the nature of the spectral statistics, we take the representative eigenwave functions in the energy range of interest and examine their regularity, where regular and irregular eigenstates typically correspond to Poisson and GOE statistics, respectively.A common situation is that the spectral properties lie somewhere in between Poisson and GOE statistics.In this case, we employ an interpolation fit by assuming that the Hamiltonian matrix is a composition, with appropriate weights, of the matrices that generate Poisson and GOE statistics [53,54].For nonrelativistic quantum systems, previous works [4,10] indicated that strengthening the electron-electron interactions makes the Hamiltonian more complicated, and their spectra align more closely with the random-matrix ensemble statistics.However, we find that, for graphene billiards in most energy ranges, the spectral statistics are insensitive to the Hubbard interaction but with an exception: for energies close to the Dirac point, the Hubbard interactions tend to drive the spectral statistics toward the GOE, regardless of the nature of the initial statistics in the absence of such interactions.This continual evolution toward GOE holds until the Hubbard interaction strength reaches a critical point of transition from ferromagnetic to antiferromagnetic order [25,55].A surprising phenomenon is that, as the interaction strength increases further from the critical value, due to the emergence of a gap and a metamorphosis in the behavior of the quasiparticles, the spectral statistics tend to approach Poisson.Taken together, as schematically illustrated in Fig. 1, when the interaction strength increases from zero, three distinct regimes can arise: (1) a weak interaction regime in which distinct spectral statistics exist, (2) an intermediate interaction regime (about the critical point) in which the statistics are universally GOE, and (3) a strong interaction regime in which the statistics are universally Poisson.We articulate a physical picture based on the energy band structure to heuristically understand these behaviors.The transition toward one universal class (GOE) and then to another (Poisson) as the many-body interaction strength increases, regardless of the nature of the spectral statistics to begin with, is particularly striking.Our results elucidate the interplay between many-body interactions and the spectral statistics of relativistic quantum systems and represent a contribution not only to the field of relativistic quantum chaos but, more broadly, to fundamental physics.

A. Mean-field Hubbard Hamiltonian
To incorporate electron-electron interactions into the quantum dynamics of a sector-shaped graphene billiard, we use the one-orbital mean-field Hubbard model [18,56,57].The model Hamiltonian has two parts: a π -orbital tight-binding Hamiltonian H TB for a single particle and the Hubbard potential H U that describes the repulsive electron-electron Coulomb interactions: where the summation is over all pairs of nearest neighboring sites i, j , the index σ denotes spin, and c † i,σ (c j,σ ) is the creation (annihilation) operator at site i ( j) for spin σ .The nearest neighbor hopping energy is t = 2.8 eV [25].The electronic band structure of graphene can be described by the tight-binding model.The electrons interact with each other via a screened Coulomb potential, which can be approximated by the Hubbard potential that is introduced through the short-range repulsive onsite Coulomb interaction strength U FIG. 1. Main result of the paper: effects of many-body Hubbard interactions on the spectral statistics of circular sector graphene billiards for energies near the Dirac point.Schematically shown is the characteristic evolution of the energy level spacing distribution as the Hubbard interaction strength U increases from zero.The starting point is U = 0, at which the level spacing statistics can be either Poisson or Gaussian orthogonal ensemble (GOE), where the inset shows the corresponding spatial distribution of the eigenstates, which is typically regular for Poisson and irregular for GOE.As U increases toward a critical point U c , the level spacing distributions approach GOE, regardless of the nature of the initial distribution at U = 0.As U increases from U c , the statistics evolve toward universally Poisson.as where σ is the opposite spin to σ stipulated by the Pauli exclusion principle, and n i,σ = c † i,σ c i,σ is the number operator.Due to the exponential growth of the Hilbert space dimension with the number N of electrons, the diagonalization and analysis become computationally difficult even for moderate system sizes with only tens of atoms.A solution is to take advantage of the mean-field approximation [58][59][60][61][62][63]: where the mean-field Hamiltonian describes the situation in which a spin-σ electron at site i interacts with the average spin-σ electron population n i, σ at the same site and vice versa.This mean-field Hubbard model is like those deduced from the first-principle or quantum Monte Carlo calculations, suggesting that the mean-field approximation is suitable for graphene [48,59,60,64,65], especially in the weak-coupling regime [48,59].Another issue of possible concern is the range of the physically meaningful value of U in Eq. (3).At present there is no consensus on the actual values of U for graphene [56], due to the lack of experiments with magnetic graphene systems which would allow U to be estimated more accurately [57].Previous discussions indicated that the value of U can range from 2.8 to 8.4 eV for graphene quantum dots [55,61,66].It has been known that, when the onsite Coulomb repulsion potential approaches the critical value U c /t ∼ 2.23, the honeycomb lattice Hamiltonian will experience a symmetry-breaking process with a phase transition from ferromagnetic order (below the critical value) to antiferromagnetic order (above the critical value) [25,55,67].
For our theoretical study herein, we employ the honeycomb lattice model and treat U as a control parameter in the range of U/t ∈ [0, 3] to uncover and quantify the effects of the electron-electron interaction strength on the spectral statistics, especially before and after the transition point U c .As will be shown in Sec.III, about the phase transition point, the spectral statistics are close to GOE.
To solve the eigenvalues and eigenvectors from the Hamiltonian in Eq. ( 4), we employ an iterative procedure [18] and assume that the graphene system is at half-filling of π electrons, corresponding to electronically neutral graphene, at zero temperature [56].This half-filling is the commonly exploited situation in the literature [48,56,[68][69][70], as it exhibits interesting phenomena such as Mott insulating behavior and antiferromagnetic order [71,72].Assigning an initial value of n i, σ to site i, we write the Hamiltonian in Eq. ( 4) for spin-σ electrons in a matrix form, where the eigenenergies and eigenstates can be solved using direct diagonalization.From the eigenstates ψ α,σ associated with spin σ , the mean occupation number at site i can be obtained as With n i,σ , the eigenenergies and eigenstates for spin σ can be obtained directly from the corresponding Hamiltonian in Eq. (4), from which n i, σ can be derived similarly, completing one iteration.The process repeats until n i,σ and n i, σ reach a steady state, yielding spin-resolved single-particle spectra [56,59,73].

B. Circular sector graphene billiards
A circular sector graphene billiard can be cut from a graphene sheet, where it was found previously that, for most concentric angles, the spectral statistics are close to GOE for energy levels about the Dirac point in the absence of electron-electron interactions [38].(It is worth noting that, for the corresponding Schrödinger billiards, the statistics are Poisson [38].)There are a few exceptions.For example, for the 60 • graphene sector with perfect zigzag or armchair radial boundaries, the spectral statistics are Poisson.To be concrete, we consider one example from each category: (1) a 15 For both systems, the iterative procedure can obtain the corresponding spin-resolved singleparticle spectrum of N energy levels [56,59,73]; thus, we can have a sufficiently large number of energy levels, e.g., 500 or 1000, for different energy ranges to generate accurate spectral statistics.

C. Spectral statistics
We divide the spectra about the Dirac point into several components and unfold each separately.For eigenenergies be the smoothed function of the counting function N c (E ) for the number of eigenenergies below E .The unfolded spectra are given by E u j ≡ N c (E j ) , and the level spacing is S j = E u j+1 − E u j .The following results have been well established in the field of traditional (nonrelativistic) quantum chaos [2].If the classical system is integrable, the level spacing follows the Poisson distribution: When the classical dynamics are fully chaotic and do not possess any geometric symmetry, GOE statistics arise: An intermediate type of statistics between Poisson and GOE, the so-called semi-Poisson statistics [75], can also arise: With the level spacing distribution P(S), the accumulated distribution can be obtained as Two additional relevant quantities are the number variance 2 (L) and the spectral rigidity 3 (L), where L is the number of mean level spacing [2].
To have a reference for the spectral statistics, we adopt a parameter-dependent random matrix model that interpolates between random matrices exhibiting Poisson and GOE statistics [53,54]: where H 0 belongs to a diagonal matrix of random Poisson numbers, H 1 is a random matrix with GOE statistics, and the variances of the matrix elements of H 0 and H 1 are chosen such that their eigenvalues have the same mean spacing.As the number of levels in each group in this paper is within 600, for matrices H 0 and H 1 , the size is also set as 600.As the control parameter λ changes from zero to infinity, the spectral fluctuations of H (λ) exhibit a transition from Poisson to GOE.Since the spectral properties are already close to GOE for λ 1, we set the maximum value of λ to be six.To estimate the parameter λ for a given spectrum, we generate a series of spectra from the random matrix model Eq. ( 5) with λ values ranging from 0 to 6 at the incremental step of 0.01 and calculate the relevant functions 2 (L) and 3 (L).
We then calculate the variance between these functions from the given spectral data and the standard references from the random matrix model.Those with the minimum variance are deemed proper.In our calculation, 2 (L) and 3 (L) typically account for long-range correlations, but 2 (L) would be too sensitive when L is large and may yield spurious results.To overcome this difficulty, in the calculation of λ, we constrain L in the range (0, 2) for 2 (L) and L ∈ (0, 10) for 3 (L).

III. RESULTS
Examining the spectral statistics in different energy ranges, we find that the Hubbard electron-electron interactions will affect the spectral statistics but only for energies near the Dirac point.Particularly, since the Dirac point (or the gap) is around N/2, in the following, we shall focus on the energy levels in the range [N/2 − 1000, N/2 + 1000] with a total of ∼2000 levels and investigate how the spectral statistics change as the Hubbard interaction increases.Our calculations reveal a simplification: the spectral statistics for spin-up and spin-down energy spectra are indistinguishable.It thus suffices to focus on the spectral statistics of the eigenenergies for one spin direction, e.g., the spin-down spectra.

A. Energy range and spin for spectral statistics
To determine a proper energy range for analysis, we compare the level spacing statistics between two cases: one without electron-electron interaction (U = 0) and the other with a sufficiently large Hubbard strength (U = 2t) for the 60 • AA graphene billiard, as shown in Fig. 3, where the insets display the density of states with the blue shaded region indicating the energy range used to calculate the level spacing statistics.The left and right columns of Fig. 3 are for U = 0 and 2t, respectively.Figures 3(a) and 3(e) are for the energy levels about the Dirac point, where the level-spacing distribution is Poisson for U = 0 and GOE for U = 2t.Figures 3(b) and 3(f) show the statistics in an energy range slightly further from the Dirac point for U = 0 and 2t, respectively, and the corresponding results in an energy range far below the band edge are shown in Figs.3(c) and 3(g).In these four cases, the statistics are GOE.Figures 3(d) and 3(h) show the statistics at the band edge so that the quasiparticles are described by the Schrödinger equation with a modified effective mass, where the level spacing statistics are Poisson for both U = 0 and 2t.From these results, we conclude that the Hubbard interactions affect the level spacing statistics but only for energies around the Dirac point, where a change from Poisson [Fig.3(a)] to GOE [Fig.3(e)] occurs as the interaction strength changes from U = 0 to 2t.Away from the Dirac point, the statistics are either GOE or Poisson, regardless of the presence of the Hubbard interactions.Note that the slight deviation of the result in Fig. 3(a) from the Poisson statistics is due to the finite size of the system [38].In the following, we shall focus on the eigenenergies near the Dirac point.
To distinguish the effect of spin orientation, we choose the spin-resolved energy spectra in the ranges (E (N/2), E (N/2 + 500)) and (E (N/2 − 500), E (N/2)) for spin-up and spindown spectra, respectively, for the two graphene sector billiards 15 • ZM and 60 • AA with U = 2t and a case of U = 1t for 60 • AA, where E (N/2) = 0 is the energy of the Dirac point.We disregard the edge states in the vicinity of the Dirac point.Figure 4 shows the spectral statistics, i.e., the distribution of the unfolded level spacing P(S), the cumulative distribution I (S) = S 0 P(S )dS , the number variance 2 (L), and the spectral rigidity 3 (L) for the three cases and the two spin orientations.It can be seen that P(S), I (S), and 3 (L) are practically indistinguishable for the spin-up and spin-down spectra, with only a small discrepancy appearing in the 2 (L) statistics.This can be understood from the mean-field Hubbard Hamiltonian in Eq. ( 4), where the physical properties associated with the two spin orientations are symmetric to each other.Figure 5 shows the best-fit parameter λ as a function of the relative electron-electron interaction strength U/t.In all cases, for U = 0 so that the energies are all close to the Dirac point, λ is small and ∼0.5, indicating Poisson-like spectral statistics.Note that the λ values for regions 3 and 4 are slightly larger than those for regions 1 and 2 because they are further away from the Dirac point than regions 1 and 2. This behavior is consistent with those in Figs.3(a (1,2) and regions (3,4) is that, the former reaches a peak about U/t = 2.1 with λ values >1, indicating a transition to GOE-like spectral statistics.As U increases further, λ reduces to ∼0.5, signifying a return to Poisson-like spectral statistics.For regions 3 and 4, λ increases approximately monotonically with U , indicating a monotonous transition from Poisson-like to GOE-like statistics.The behaviors of λ extracted from 2 (L) and 3 (L) are mostly consistent with each other, except a few cases due to the different L ranges about U/t = 3, as shown in Figs.5(c) and 5(d).The different responses of λ between regions 1 (2) and 3 (4) give a rough indicator of the energy range that can be deemed as the neighborhood of the Dirac point.
The transition in the spectral statistics is further elaborated in Fig. 6 which shows I (S) − I GOE (S) and 2 (L) for different U values.In the plot of I (S) − I GOE (S), three theoretical curves are also presented as a reference: the corresponding quantity for Poisson, the intermediate semi-Poisson, and GOE (the horizontal line at zero).As U/t increases from 0 to 2.1, starting from close to Poisson, the curves for the data shrink continually, passing the semi-Poisson curve, and approach the zero horizontal line for GOE.However, as U increases further, the data curves in the I (S) − I GOE (S) plot expand and deviate away from the zero horizontal line, i.e., returning to Poisson.Note that the semi-Poisson curve is included just as a reference, as it represents only an intermediate case in the change of the parameter U .The plots of 2 (L) show similar features, i.e., the data curves decrease from close to Poisson, and arrive at a curve that is close to GOE at U/t = 2.1 but then increase and pass the semi-Poisson line.
The distinct spectral statistics for different values of the Hubbard interaction strength U are also corroborated by the eigenwave functions, as shown by four different cases in Fig. 7, where six eigenstates around the Dirac point are displayed for each case.For U = 0, the spectral statistics are close to Poisson, and most of the eigenwave functions exhibit regular patterns, as exemplified in panels 1-4 and 6 in Fig. 7(a), but some eigenstates exhibit irregular patterns, as illustrated in panel 5, which are not uncommon when the spectral statistics follow those of random matrix ensembles.These irregular states are the main cause of the deviation of the spectral statistics from Poisson.For U = 2t, as shown in Fig. 7(b), most eigenstates are irregular, which is in accordance with the observed GOE-like statistics demonstrated in Fig. 6 for the same U value.This suggests that, as the Hubbard interaction increases from zero, most of the eigenstates change from regular to irregular, leading to the observed change from Poisson-like to GOE-like spectral statistics.
For U = 3t, most states become regular again, as shown in Fig. 7(c), due to the onset of antiferromagnetic order leading to the effective staggered potential.In addition, some states are like the eigenstates of the nonrelativistic quantum billiards described by the Schrödinger equation with the Dirichlet boundary condition, e.g., ψ S,mn ∼ sin mπφ φ 0 J mπ/φ 0 (k mn ρ) with m, n = 1, 2, • • • and k mn being the eigenwave numbers.For example, the quantum number sets for panels 1 and 2 are (m, n) = (6, 1) and (2, 7), respectively.These are different from the eigenstates close to the Dirac point in a graphene billiard without Hubbard interactions, as shown in Fig. 7(a).This indicates that, while the spectral statistics return to Poisson-like in the large U regime, the eigenstates of the system do not return to those in the U = 0 case but undergo a metamorphosis in that they become nonrelativistic quantum states as in a nonrelativistic quantum billiard corresponding to the case about the band edge, e.g., E /t ∼ ±3 in the U = 0 graphene billiard.Indeed, about U/t 2.2, a gap begins to open and becomes larger as U increases further [76], leading to nonrelativistic quantum behaviors, especially at the bottom of the band when the effective wave number is small.For regions 3 and 4, the general tendency of λ is to increase with U/t, from <1 to >1, as shown in Figs.5(c) and 5(d), signifying a transition of the spectral properties from Poissonlike to GOE-like statistics.For U = 3t, in contrast to the behaviors in regions 1 and 2, the spectral statistics do not return to Poisson but are even closer to GOE.This can also be corroborated from the wave function patterns, as exemplified in Fig. 7(d), which are apparently irregular.

C. Spectral statistics of 15 • ZM graphene billiard
We have seen from Sec. III B that, for the 60 • AA graphene billiard, when Hubbard interactions are not present, the spectral statistics are close to Poisson, making the changes in the spectral statistics significant as U increases.However, for most graphene sector billiards, the spectral statistics are already GOE-like for eigenenergies about the Dirac point even for U = 0, so the expectation is that, as U increases, the changes in the spectral statistics, if any, would be insignificant.To test if this expectation holds, we consider the 15 • ZM graphene billiard and examine the responses of the spectral statistics to tuning on the Hubbard interaction.As will be shown below, the interplay between many-body interactions and the spectral statistics can be interestingly nontrivial with unexpected features.
In regions 1 and 2 close to the Dirac point, as shown respectively in Figs.8(a) and 8(b), for U = 0, the spectral statistics are approximately GOE, so λ ≈ 2. A previous work demonstrated that, for larger systems, the statistics are perfectly GOE [38].As U increases, λ reaches the maximum value around U/t = 2.2 despite fluctuations.As U increases further to U/t = 3, λ decreases.For region 1 (below the Dirac point), the value of λ becomes <1, leading to semi-Poisson statistics.However, for region 2 (above the Dirac point), the value of λ decreases to ∼0, giving rise to Poisson statistics.The changes occur continuously with U .These results are surprising because the expectation is that the spectral behaviors for energies above and below the Dirac point should be the same.However, under strong Hubbard interactions, the system can be in the In regions 3 and 4, as shown in Figs.8(c) and 8(d), the values of λ are all >1 for all values of U , indicating that the spectral properties are all close to GOE.This is consistent with the results in Figs. 3 and 5 that higher energies or stronger electron-electron interactions typically lead to GOE-like spectral statistics.Nevertheless, as U increases, with oscillations, λ reaches the maximum around U = 2t according to both the 2 (L) and 3 (L) statistics.
Figure 9 displays the integrated level spacing distribution I (S) − I GOE (S) and the number variance 2 (L) for different U values for regions 1 and 2. For most U values, the I (S) − I GOE (S) curves are ∼0, indicating GOE-like spectral statistics.For U approaching 3t, the curves in region 1 rise slightly and approach the semi-Poisson curve, but those in region 2 rise significantly to approach Poisson.Similar behaviors occur in the 2 (L) curves.
To understand the distinct behaviors in regions 1 and 2 for U = 3t, we display the representative eigenwave functions from the two regions in Fig. 10, where Figs.10(a 1 ) and 10(a 2 ) show the lowest nonlocalized eigenstates with respect to the Dirac point from both sides.The eigenstates are similar and resemble the first eigenstate in the corresponding Schrödinger billiard.Note that, for such a large value U , the Dirac point has already been destroyed and a gap appears-here, the Dirac point is only used as a reference point for regions 1 and 2 with FIG. 9. Deviation of the accumulated distribution from Gaussian orthogonal ensemble (GOE) and the number variance for the 15 • ZM graphene billiard.Shown are I (S) − I GOE (S) and 2 (L) for regions 1 (blue) and 2 (red) for U values ranging from 0 to 3t.The circles are for region 1 (red curves), and the squares are for region 2 (blue curves).The coordinate of U/t is not linear for a clear comparison with the corresponding Poisson (solid), semi-Poisson (dash-dotted), and GOE (dashed) curves for U/t = 3. the highly localized states neglected.Figures 10(a 1 )-10(l 1 ) and 10(a 2 )-10(l 2 ) show 12 representative eigenstates for each region in a span of 440 eigenstates that are sampled at the interval of 40 between two neighboring panels.Apparently, the eigenstates in the two regions are quite distinct.In region 1, only the first four eigenstates are regular, and most others are irregular, as the spectral statistics are semi-Poisson.In contrast, in region 2, the first eight eigenstates are regular, and there are only a few irregular states, which are consistent with the Poisson spectral statistics in this region.

IV. HEURISTIC UNDERSTANDING OF MANY-BODY SPECTRAL STATISTICS
The many-body interaction mediated transition in the spectral properties observed from the 60 • AA and 15 • ZM graphene billiard systems can be understood through the mean-field Hamiltonian in Eq. ( 4).Recall that the eigenenergy spectrum is obtained by iteratively solving the eigenenergies of Eq. ( 4) when a steady state is reached.The mean-field Hamiltonian matrix can be written as where H TB is the hopping term (H TB,i j = −t if sites i and j are nearest neighbors, and H TB,i j = 0 otherwise), and H U σ is a diagonal matrix associated with the Hubbard interaction term for spin-σ states that can effectively be regarded as an onsite potential: Without loss of generality, we examine the average density distribution for spin-up electrons, i.e., n i,↑ .Below the critical point U c , we have n i,↑ = 0.5 for most atoms in the billiard region.Significant deviations from the average occur mostly on the atoms near the billiard boundary, which are then spin polarized (up or down).An example is presented in Fig. 11(a) for U ∼ 0.5t, where the values of n i,↑ for both A and B atoms are shown.The corresponding density distribution in the physical space is shown in Fig. 11(d).For U above the critical point, most atoms are polarized.If atom i belongs to the sublattice A (B), we have n i,↑ > 0.5 ( n i,↑ < 0.5).From the relation H U ↓,ii = U n i,↑ , we have that the onsite potential for spin-down electrons is positive for one sublattice and negative for the other sublattice with a shift U/2, which is effectively a staggered potential and generates a mass term mσ z for the effective Dirac Hamiltonian.The effective mass m depends on U and becomes larger as U increases.When U exceeds a critical value, the original Dirac point is destroyed, opening an energy gap.Since the energy range of interest is close to the Dirac point, we have h2 k 2 m, so the collective modes are now described by the Schrödinger equation, Specifically, we assume different but uniform onsite potentials U A and U B for A and B atoms (U A > U B ).About the original Dirac point (with a gap), the dispersion relation between E ↓ for the spin-down states and the wave vector is [77] where H 12 is given by the plane-wave assumption and the overlap integral t between the wave functions centered at atom A and one of its nearest neighbors B under the lattice potential can be determined through where a = 2.46 Å is the lattice constant.For |H 12 | (U A − U B ), Eq. ( 6) can be simplified to Close to the Dirac point, we have , which approximates the dispersion relation of a Schrödinger particle.
On the zigzag boundary or a zigzag segment, some A (B) atoms have a different value of n i,↑ than the values of most of the other A (B) atoms inside the billiard, as shown in Figs.11(d)-11(f).This results in the emergence of correlation among the energy levels, leading to a deviation from the Poisson toward GOE-like statistics.As shown in Fig. 11(b), if the zigzag boundaries are long, n i,↑ is generally not symmetric for A and B atoms, e.g., one type of sublattice atom can have a broader distribution of the density values than the other sublattice atoms.The atoms with extraordinary density values are mostly located at the boundaries, as illustrated in Figs.11(d)-11(f).For example, in our convention, the most outward atoms at the long bottom zigzag boundary of the 15 • ZM billiard are of type B. For most of B atoms in the inner part of the billiard, for U = 3t, we have n i,↑ ≈ 0.27.However, for the A atoms, we have n i,↑ ≈ 0.72.At the long bottom zigzag boundary [Fig.11(e)], the values of n i,↑ on the B atoms can be as high as 0.82, even higher than those of typical A atoms.On the contrary, the nearby A atoms have a markedly smaller value, i.e., ∼0.33 instead of 0.72, but still larger than the typical values for B atoms.These lead to abnormal onsite potentials on the zigzag boundaries and further result in longer-range arclike structures around them.In the inner domain of the billiard, a dual-value density distribution arises: 0.72 for A and 0.27 for B. Since the onsite potential can be regarded as the effective mass, the deviations from the values of the inner atoms effectively form a mass disorder.The amount and amplitude of the disorder are much larger for U B than those for U A , making E ↓,− more random than E ↓,+ .This is consistent with the results that, in energy region 1 (energy levels below the Dirac point, corresponding to E ↓,− and B atoms), the spectral statistic deviates further from Poisson than that in region 2 (energy levels above the Dirac point, corresponding to E ↓,+ and A atoms).

V. DISCUSSION
Uncovering and understanding the effects of many-body interactions on the spectral statistics has long been recognized as a fundamental problem in quantum chaos, but previous works focused on nonrelativistic, Schrödinger quantum systems [4][5][6][7][8][9][10][11][12][13][14].How many-body interactions affect the energy level statistics in relativistic quantum systems remains largely unexplored.This paper represents a step forward in this direction.Utilizing circular sector graphene billiards as a prototypical system and treating the electron-electron interactions through the Hubbard formalism, we systematically investigate the evolution of the spectral properties as the many-body Hubbard interaction strength increases from zero.The graphene billiard system is experimentally feasible and can be tuned to exhibit relativistic quantum or nonrelativistic quantum characteristics.Especially for energy levels about the Dirac point, the system is in the relativistic quantum regime, while for energies close to the band edge, the system becomes of the Schrödinger type.
The contribution of this paper is that, as the many-body Hubbard interaction strength U is tuned, the level spacing distribution in the relativistic quantum regime can exhibit a surprising transition scenario.Specifically, we have taken two distinct systems that exhibit characteristically different spectral statistics in the absence of many-body interactions: the 60 • AA and 15 • ZM circular sector graphene billiards, where the level spacing statistics are of the Poisson type for the former and GOE type for the latter for U = 0.For the 60 • AA billiard, the following transition scenario occurs.As U increases from zero, a transition from Poisson to GOE occurs below the critical value U c , but as U increases further, a transition back to Poisson occurs.The first transition (i.e., the one from Poisson to GOE for U < U c ) occurs gradually, with the spectral statistics most close to GOE slightly below the critical point.For the 15 • ZM graphene billiard, as U increases from zero, the spectral statistics remain GOE-like and become maximally GOE (as characterized by the quantity λ) near the critical point.As U increases further, a transition to Poisson-like statistics occurs.The transitions can be visually confirmed by examining the spatial patterns of the eigenstates, as the Poisson and GOE statistics typically correspond to regular and irregular patterns, respectively.As summarized in Fig. 1, the key finding of this paper is the following striking phenomenon: regardless of the nature of the spectral statistics to begin with, as the Hubbard interaction strength increases from zero toward a critical point, GOE statistics dominate but, passing through the critical point as the interaction strength increases further, a transition to Poisson statistics occurs for both of the circular sector billiards.
In general, Hubbard interactions will make the system more complex so that the spectral statistics would be GOE.Previous works revealed a sharp transition to GOE with even weak Hubbard interactions at low filling [8].For the 60 • AA billiard system studied in this paper, at half-filling, the transition in the spectral statistics from Poisson to GOE occurs gradually.Furthermore, for strong Hubbard interactions (e.g., U = 3t) and for systems with long zigzag boundaries, the effective potential of the atoms on sublattices A and B are not symmetric.Along the long zigzag boundaries, the outermost B atoms mostly have low effective potentials, but some B atoms can have a drastically large potential accompanied by a sharp reduction in the potential value for nearby A atoms, as described in Sec.IV.This uneven effective potential distribution for atoms on the two sublattices of graphene is the reason for the observed distinct transition behaviors in the statistics for the spectrum below and above the gap at large U .
The main finding of this paper, i.e., the transitions to GOE and then to Poisson level spacing statistics, was obtained using the mean-field Hubbard model.A question is whether the transitions are not an artifact of the mean-field approximation.Indeed, this approximation is valid in the weak interaction regime with relatively small values of the electron-electron interaction strength U , as previously demonstrated through a comparison with the results from density function theory [59] or from the quantum Monte Carlo method [67].In the strong interaction regime, the mean-field model may break down.However, there is a physical justification for the validity of our finding.The uncovered transitions to GOE and then to Poisson level spacing statistics hinge on two factors: (1) the occurrence of the phase transition in the honeycomb lattice from ferromagnetic to antiferromagnetic order and (2) the opening of an energy gap.For example, a transition in the spectral statistics takes place when U is about U c ∼ 2.1t where the statistics are effectively GOE, as demonstrated in Figs.5(a for the 15 • ZM billiard.This value of U c is close to the phase transition point 2.23t where the honeycomb lattice transitions from ferromagnetic to antiferromagnetic order [67].Beyond the critical value U c , a gap opens, and the effective onsite potentials for the A and B atoms in the graphene lattice begin to deviate from each other.The difference becomes larger as U increases further, contributing effectively a mass term to the Dirac Hamiltonian and rendering the relativistic quantum quasiparticles on the lattice massive.As a result, in the energy range close to the Dirac point where the relative wave number k is small, the massive Dirac equation can be approximated by the Schrödinger equation.When this happens, a circular sector quantum billiard, regardless of the angle, becomes integrable with the Poisson spectral statistics.Summarizing it, the transition in the spectral statistics uncovered in this work will be valid insofar as there is a ferromagnetic-to-antiferromagnetic transition which has been revealed in previous works [25,67,68,78,79].Our finding thus is not an artifact of the mean-field model but a robust physical phenomenon. Another issue concerns the possible role of many-body localization (MBL) in the spectral statistics transitions.MBL, a many-body phenomenon in quantum systems, hinders the system from reaching a thermal equilibrium [80].Previous works on spin-chain models revealed that MBL can lead to Poisson statistics for the energy spectrum [81,82], which is consistent with the results on graphene billiards where the localization of quantum states can generate Poisson-like spectral statistics [33,36,83] due to the breaking of the correlation between the states [31].(Spectral correlation due to, e.g., chaotic dynamics, generally results in GOE-like spectral statistics [28][29][30].)However, we found that the trend toward Poisson for large U in this work is not caused by any localized states but due to the transition from massless to massive Dirac quasiparticles rendered by the opening of an energy gap.We have calculated the participation ratio that characterizes the degree of localization of the eigenstates [84] and found no systematic change in this ratio associated with the transition across U c , even in the large U regime toward the Poisson statistics.
Finally, it was conjectured that, for systems under a metalinsulator phase transition, the spectral statistics should follow semi-Poisson close to the transition point [75,85].In our study, semi-Poisson statistics can also arise when the Hubbard interaction strength falls in an appropriate interval (i.e., U ∼ t to 1.5t), but the statistics are not related to the critical states.

APPENDIX: SPECTRAL PROPERTIES OF ADDITIONAL CIRCULAR SECTOR GRAPHENE BILLIARDS
To further demonstrate that the energy levels above and below the Dirac point can follow distinct spectral statistics as observed from the 15 • ZM billiard at U = 3t, we examine several additional sector billiards with different angles and boundaries.As shown in Fig. 12, the 2 (L) statistics of the 45 • ZM billiard is most close to that of the 15 • ZM billiard treated in the main text, while the 2 (L) statistics of the 45 • AM billiard is close to that of the 60 • AA billiard.These results are consistent with the analysis in the main text.

FIG. 2 .
FIG. 2. Circular sector graphene billiards.The concentric angle is (a) 15 • and (b) 60 • .The two nonequivalent atoms in a unit cell of the graphene lattice, namely, A and B atoms, are marked as blue and orange circles, respectively.The painting of the lattice uses pybinding [74].
• sector with one radial edge being zigzag and the other radial edge being a mixture of zigzag and armchair [denoted as 15 • ZM, Fig. 2(a)], and (2) a 60 • sector with both radial edges being armchair [denoted as 60 • AA, Fig. 2(b)], where the respective spectral statistics are GOE and Poisson for energy levels about the Dirac point in the absence of Hubbard interactions.The parameter values of the two sector billiard systems are as follows.The 15 • ZM billiard has N = 40 328 atoms, with size ∼90 × 25 nm.The 60 • AA billiard has N = 42 540 atoms, with size ∼46 × 40 nm.

FIG. 4 .
FIG. 4. Comparison of spectral statistics for spin-up and spindown spectra for the 15 • ZM and 60 • AA graphene sector billiards with eigenenergies close to the Dirac point.(a) The nearest neighbor level spacing distribution P(S), (b) the cumulative distribution I (S) subtracting that of the Poisson distribution I P (S), (c) the number variance 2 (L), and (d) the spectral rigidity 3 (L).For spin-up and spin-down spectra, the energy ranges are (E (N/2), E (N/2 + 500)) and (E (N/2 − 500), E (N/2)), respectively.The solid, dashdotted, and dashed curves are the random matrix results for Poisson, semi-Poisson, and Gaussian orthogonal ensemble (GOE) statistics, respectively.

FIG. 5 .
FIG. 5. Determination of the random matrix parameter λ as an indicator to interpolate between Poisson and Gaussian orthogonal ensemble (GOE).Shown is the best-fitted parameter λ obtained through a comparison between the calculated spectra with the standard 2 (L) (red up triangles) and 3 (L) (blue circles) statistics from random matrices constructed from Eq. (5) vs U/t for the 60 • AA graphene billiard.(a)-(d) Results for regions 1-4, respectively, which are marked in the insets of (a) for U = 0 and (d) for U/t = 2.3.
FIG.6.Deviation of the accumulated distribution from Gaussian orthogonal ensemble (GOE) and the number variance for the 60 • AA graphene billiard.Shown are I (S) − I GOE (S) and 2 (L) for region 1 under different electron-electron interaction strength ranging from 0 to 3t.Here, the coordinates U/t are not linear, as they are normalized to the corresponding Poisson (solid), semi-Poisson (dash-dotted), and GOE (dashed) curves at U/t = 3.

FIG. 8 .
FIG. 8.The random matrix parameter λ interpolating between Poisson and Gaussian orthogonal ensemble (GOE) for the 15 • ZM graphene billiard.Legends are the same as those in Fig. 5.

FIG. 10 .
FIG. 10.Distributions of wave function intensity |ψ j | 2 for the 15 • ZM sector graphene billiard for U = 3t.(a 1 )-(l 1 ) |ψ j | 2 of the 20 138th to the 19 698th eigenstates with a step of 40 eigenstates between two adjacent panels, whose eigenenergies belong to the energy range of region 1. (a 2 )-(l 2 ) |ψ j | 2 for the 20 199th to the 20 639th eigenstates in region 2. Brighter colors correspond to larger values of |ψ j | 2 .The scales in different panels are chosen differently for better visualization.

FIG. 11 .
FIG. 11.Average electron density n i,↑ for spin-up states associated with A and B lattices in an increasing order.(a) 15 • ZM with U = 0.5t, (b) 15 • ZM with U = 3t, and (c) 60 • AA with U = 3t.(d)-(f) Spatial distribution of n i,↑ corresponding to (a)-(c), respectively, where a magnification of a small fragment of edges marked by the small red rectangle in each panel is shown.Only those atoms whose n i,↑ values deviate from the average electron density in the billiard are plotted.
FIG. 12. Spectral properties of additional circular sector graphene billiards.Shown are the 2 (L) statistics for different circular sector graphene billiards for U = 3t: (a) 15 • AM, (b) 60 • ZZ, (c) 45 • AM, and (d) 45 • ZM, where the notations A, M, and Z stand for armchair, mixed, and zigzag boundaries, respectively.A zoom-in close to the tip is shown as the insets.The orange circles and blue squares represent the data of regions 1 and 2, respectively.

ACKNOWLEDGMENTS
This paper was supported by the National Natural Science Foundation of China under Grants No. 12175090, No. 11775101, No. 12247101, and No. 12047501, and by the 111 Project under Grant No. B20063.X.C. acknowledges financial support from the China Scholarship Council No. 202006180040 and thanks B. Dietz for discussions and the suggestion to use Eq.(5) to interpolate between Poisson and GOE statistics.The efforts at Arizona State University were supported by the Air Force Office of Scientific Research through Grant No. FA9550-21-1-0186.