Chaos in the Bose-Hubbard model and random two-body Hamiltonians

We investigate the chaotic phase of the Bose-Hubbard model [L. Pausch et al, Phys. Rev. Lett. 126, 150601 (2021)] in relation to the bosonic embedded random matrix ensemble, which mirrors the dominant few-body nature of many-particle interactions, and hence the Fock space sparsity of quantum many-body systems. The energy dependence of the chaotic regime is well described by the bosonic embedded ensemble, which also reproduces the Bose-Hubbard chaotic eigenvector features, quantified by the expectation value and eigenstate-to-eigenstate fluctuations of fractal dimensions. Despite this agreement, in terms of the fractal dimension distribution, these two models depart from each other and from the Gaussian orthogonal ensemble as Hilbert space grows. These results provide further evidence of a way to discriminate among different many-body Hamiltonians in the chaotic regime.


I. INTRODUCTION
Chaotic behavior emerges frequently in both classical and quantum physics [1][2][3][4][5][6]. Classically, the defining property of chaos is the exponential sensitivity of dynamics on the initial conditions [3], giving rise to bundles of wandering trajectories in phase space. Chaos is hence intimately related to ergodicity, which is a feature of just few physical systems and typically difficult to verify [1][2][3]. For quantum systems with a strictly chaotic classical limit, statistical properties of the energy spectrum were shown to be in agreement with the (intrinsically ergodic [7]) Gaussian ensembles of random matrix theory (RMT) [4,[8][9][10][11]. RMT has therefore set the baseline to define chaos in quantum systems: Spectral statistics are commonly used to distinguish chaotic from integrable and localized regimes [12][13][14][15][16][17][18][19][20], as are RMT predictions for eigenvector features [21,22] and other observables [23].
Given that most dynamical systems have a mixed rather than purely chaotic phase space [24][25][26][27][28][29], deviations from RMT would be expected [30]. Additionally, the standard Gaussian ensembles, described by dense matrices, fail to capture a characteristic feature of many-body Hamiltonians, namely the structural Fock-space sparsity induced by the typical few-body nature of interactions. To account for the latter, further random matrix ensembles, so-called embedded ensembles, have been introduced [31][32][33][34][35][36] with the aim of obtaining a more accurate description of the universal properties of the chaotic regime of many-body quantum systems.
In a recent publication [63], we have investigated the chaotic and non-chaotic phases of one-dimensional interacting bosons described by the Bose-Hubbard Hamiltonian (BHH). Our study revealed a remarkable agreement with the predictions of the Gaussian orthogonal ensemble (GOE) of RMT, at the level of spectral statistics and the typical value and fluctuations of fractal dimensions for chaotic eigenvectors. Most interestingly, we found that, despite this agreement, both models become ever more distinguishable in terms of the fractal dimension distribution as the dimensionality of Hilbert space grows. This establishes an appealing connection with the current problem of certification of distinctive features of complex quantum systems [64][65][66][67][68], and suggests the enticing prospect of having an accessible figure of merit that could unveil the specifics of the underlying many-body Hamiltonian in the chaotic regime, which is otherwise mainly determined by universal features. We nonetheless also contemplated the possibility that such deviation could be accounted for by the more refined embedded ensembles.
The purpose of the present work is to provide a detailed analysis of the chaotic properties of BHH and of the bosonic embedded Gaussian orthogonal random matrix ensemble (EGOE), and to perform a comparison between both models to come to a definite answer to the question posed above.
The remainder of the manuscript is organized as follows. In sections II and III, we present the physical models and the theoretical tools used to characterize their chaotic behaviour. In section IV, we discuss fingerprints of chaos in the BHH spectrum and eigenstates as a function of appropriately scaled energy and interaction strength. The properties of EGOE and the comparison between both models are first presented in section V. Building upon these results, we investigate, in section VI, the scaling of the chaos markers and of the fractal dimension distributions with Hilbert space dimension, correlating the behaviours of GOE, EGOE and BHH. We summarize our findings and conclude in section VII.

A. Bose-Hubbard Hamiltonian
The one-dimensional Bose-Hubbard Hamiltonian [69][70][71][72] of bosons on lattice sites is the sum of a one-particle tunneling observable and a local two-particle interaction, = tun + int , where Here, † , are standard bosonic creation and annihilation operators associated with Wannier orbitals localized at each lattice site. tun describes particle tunneling between neighbouring sites with tunneling strength , and int provides repulsive on-site interactions with interaction strength > 0. The sum in the tunneling Hamiltonian runs from = 1 to = − 1 for hard-wall boundary conditions (HWBCs) and from = 1 to = for periodic boundary conditions (PBCs), where we identify ( †) +1 := ( †) 1 . For both boundary conditions exhibits reflection symmetry about the center of the chain, ↦ → + 1 − , and for PBCs it is also invariant under translations, ↦ → ( + ) mod . The underlying Hilbert space can therefore be decomposed into a symmetric (even parity, = +1) and an antisymmetric (odd parity, = −1) subspace with respect to reflection symmetry for HWBCs, while for PBCs the Hilbert space decomposes into subspaces characterized by the total quasimomentum = 0, 1, . . . , − 1 and, in the case of = 0 or = /2, by further reflection symmetry.
Both int and tun , which constitute the respective limits of for vanishing tunneling ( = 0) or interaction ( = 0), are individually integrable, with eigenstates | = | 1 , . . . , and |˜ = |˜1, . . . ,˜ , respectively, which are uniquely identified by the eigenvalues of the number operators = † and = † . Natural bases of the BHH are hence symmetrized versions (according to reflection and/or translation symmetry) of the interaction Fock basis {| } and of the tunneling Fock basis {|˜ }, that we refer to as the interaction basis and the tunneling basis, respectively. Note that these two bases are conjugated in the following sense: The eigenstates of in the limit → 0 ( → 0) are localized in the interaction (tunneling) basis but highly delocalized in the conjugated basis.
The non-integrability of resulting from the combination of interaction and tunneling leads to the emergence of a chaotic phase that leaves an imprint in the spectral and eigenvector properties [14,18,52,53,63].
Due to particle number conservation, [ , ] = 0 with = = ˜, we restrict our analysis to subspaces of fixed . In particular, we consider filling factor one, = , if not stated differently. As → ∞, the BHH reaches its classical limit [18,28,73], whose dynamics is governed only by the scaled energy / 2 and the scaled tunneling strength = / .
As we demonstrated in Ref. [63], is the parameter that controls the appearance of the chaotic phase in the quantum system at unit density, even for a moderate number of particles. We obtain full spectra of the BHH numerically using exact diagonalization, and also calculate eigenstates and eigenenergies around chosen target energies [74][75][76]. These calculations are performed for several irreducible symmetry-induced subspaces, for different boundary conditions and in both natural bases. We analyze spectral and eigenvector features for different as functions of the scaled energy and . Since tun ∼ = 2 and int ∼ 2 for 1 and arbitrary , note that max − min ∼ 2 for fixed and large . Hence, scales the energy in the same way as / 2 , thus effectively providing the relevant energy scale in the classical limit.

B. Bosonic embedded ensemble
The Gaussian Orthogonal Ensemble (GOE), the typical benchmark of quantum chaos in time-reversal invariant systems such as the BHH, consists of fully dense matrices, in which each basis state is connected to every other. In a many-body system, however, many-particle basis states are typically coupled via few-body processes, giving rise to a characteristic sparsity in the Hamiltonian's matrix representation.
To describe such structural constraint of many-body systems, we consider the bosonic two-body GOE embedded ensemble (EGOE), which is defined by the Hamiltonian where † , are bosonic creation and annihilation operators of arbitrary orthogonal single-particle states. (1) and (2) are independent GOE random matrices in the space of one and two particles, respectively: Their matrix elements are independent Gaussian random variables with zero mean and variances var( (1) ) = 1 + , var( (2) , ) = 1 + [6,35]. Note that the addition of the single-and two-particle terms leads to correlations in the matrix elements of EGOE , e.g., its + −1 [77] diagonal entries in the corresponding Fock basis of indistinguishable bosons, , † † /(1 + ), are defined by just independent numbers (1) and ( + 1)/2 independent numbers (2) , . The parameter tunes the strength of the interaction part and can be used to induce the emergence of chaotic behavior [35].
In order to access the same Hilbert space dimensions as for the symmetry-induced subspaces of BHH we demand that EGOE obeys reflection symmetry as BHH. This imposes additional constraints on (1) and (2) , (2) that ensure a block diagonal structure with respect to parity .
By definition, EGOE captures the maximal Fock space connectivity allowed by at most two-particle operators. The sparsity of this model (i.e., its zero matrix element structure) is therefore contained in the matrix of the BHH in any basis. Figure 1 illustrates the matrix structure of EGOE and of BHH for HWBCs in both natural bases: In the interaction basis, different Fock states may be coupled by one-particle nearestneighbour tunneling only [Eq. (1)], and hence the BHH matrix is even sparser than EGOE. In the tunneling basis, however, the two-particle interaction [Eq. (5)] becomes strongly off-diagonal, and the BHH matrix structure closely resembles that of EGOE.
Like for the BHH, the spectral and eigenvector features of EGOE matrices will be numerically investigated using exact diagonalization, implementing an ensemble (disorder) average over 100 realizations of EGOE , unless otherwise stated.

A. Level spacing ratios
A convenient tool to characterize the short-range properties of a quantum system's energy spectrum is given by the level spacing ratios , defined as [15,16,78] = min where = +1 − is the th level spacing. In contradistinction to the distribution of the bare level spacings, the calculation of does not require any-potentially intricate [79,80]-unfolding procedure. For Gaussian random matrix ensembles, an analytic approximation to the probability density of is known [78], which, for the particular case of GOE, yields the mean level spacing ratio GOE = 4−2 √ 3 ≈ 0.536, in good agreement with GOE ≈ 0.5307 obtained from large-scale numerics [78].
For a state | = | in a Hilbert space of dimension N spanned by the orthonormal basis {| }, the -moments, for ∈ R + , are defined as = | | 2 . As N → ∞, the -moments of a state that is equally distributed over the full basis scale as = N 1− , while they become independent of N for a state that lives on a fixed finite set of basis elements. For generic states, one may expect the N → ∞ scaling to be of the form ∼ N ( We stress that the GFDs are basis-dependent quantities by definition, and will typically differ among different bases.
In order to study localization in finite-dimensional Hilbert spaces, we define finite-size generalized fractal dimensions which are rescaled versions of the Rényi entropies = − ln /( − 1) [98][99][100][101][102]. The GFDs are then found as the limits = lim N→∞˜. We focus on the special values As N grows, these fractal dimensions govern the scaling of exp( ), where is the information entropy, of the inverse participation ratio IPR ≡ 2 = | | 4 [103], and of the maximum intensity of | , respectively. We note that˜is a monotonically decreasing function of [104], i.e.,˜1 ≥˜2 ≥ ∞ for any state.

IV. CHAOS IN THE BOSE-HUBBARD MODEL
We now unveil the chaotic phase of the BHH by investigating the features of its spectrum and eigenstates, as functions of scaled interaction strength [Eq. (7)] and scaled energy [Eq. (8)].

A. Localization properties of individual eigenstates
In Fig. 2, we plot the fractal dimensions˜( = 1, 2, ∞) of all BHH eigenstates, in both natural bases, as functions of their scaled eigenenergies for four different values of the parameter = / . The subspace considered is that of parity = −1 and total quasimomentum = 0, for = 12 particles on = 12 sites with PBCs (N = 55898). Furthermore, we depict the density of states (DOS) and the distributions of˜over all eigenstates.
The DOS is dominated by discrete peaks at small = 0.0028 (i.e., prevailing interaction), changes to a seemingly Gaussian shape for intermediate = 0.0252 and = 0.2759, and reaches a structure featuring several peaks on top of the Gaussian shape at large = 2.5166 (i.e., dominating tunneling). Upon increasing , the DOS shifts from the low energy part of the spectrum (with its maximum at ≈ 0.2) to an almost symmetric shape around = 0.5.
The behavior of the DOS for extreme values of is inherited from the limiting Hamiltonians. For dominant interaction ( 1), is close to int , whose eigenenergies form degenerate manifolds at integer multiples of . The multiplicity of the manifolds, and hence the height of the associated DOS peaks is given by the number of permutations of the on-site occupation numbers { 1 , . . . } that give rise to the same eigenenergies. This number is largest, ∼ !, when most are different. Each particle of such states interacts only with − 1 particles, and the corresponding energy lies in the lower half of the spectrum. In contrast, high-energy levels correspond to states with − particles on a single site, , and their multiplicity is of order ∼ −1 !. In the opposite limit ( 1), is close to tun , whose spectrum is the sum of identical single-particle spectra. According to the central limit theorem, the DOS should then converge to a Gaussian. Deviations from this shape are expected due to the discrete nature of the single-particle spectrum for finite .
The fractal dimensions also show a clear development with that is qualitatively the same for all˜considered. For small , the values of˜for states emerging from the same degenerate manifold (for = 0) are widely spread irrespective of the basis under consideration. In this limit, according to Sec. II A, the eigenstates will tend to be localized in the interaction basis, and hence exhibit lower fractal dimensions than in the tunneling basis, where the eigenstates are more delocalized. Note, however, that full localization (i.e.,˜→ 0) in the interaction basis for small finite cannot occur, since the actual eigenstates as → 0 follow from the diagonalization of the perturbation tun on the degenerate manifolds, mixing many Fock states (an ever increasing number with ) of the same degenerate manifold.
For intermediate values, the˜distribution for states that are close in energy becomes very narrow in both bases. This narrowing is enhanced at the bulk of the spectrum-hence anchored to the position of the DOS maximum-, which broadens and shifts to higher energies upon increasing (compare = 0.0252 and = 0.2759). In fact, in the interaction basis the distribution of fractal dimensions becomes significantly flat, a behaviour which does not distinctively appear in the tunneling basis. Nevertheless, the emergence of a very narrow distribution in the bulk of the spectrum for intermediate seems to be a basis-independent feature.
For the largest shown, = 2.5166, the distribution ofb roadens again. In the interaction basis the fractal dimensions denote a delocalization tendency of the eigenstates (albeit reaching slightly lower˜values than for intermediate = 0.2759), whereas in the tuneling basis they tend towards small values without reaching full localization, in agreement with the conjugated nature of both bases discussed earlier and with mixing many degenerate tunneling Fock states in the actual eigenstates [105].

B. Energy-resolved spectral statistics and fractal dimensions
We now proceed to investigate systematically the relation between spectral and eigenstate structure changes as functions of and . To this end, we characterize statistically the distribution of energy levels and fractal dimensions over small energy windows. Figure 3 shows the evolution with and of mean level spacing ratio , and mean ˜ , variance var ˜ , and absolute value of skewness skew ˜ = ˜− ˜ 3 /var ˜ 3/2 of fractal dimensions for = 1, 2, ∞ and = = 12 with PBCs [subspace of = 0 and = +1 (N = 56822)]. To obtain these statistical quantifiers, the axis is divided into 100 bins of equal width, and the and˜distributions are built from the eigenvalues and eigenstates that fall into each bin at fixed . Additionally, Fig. 3(b) shows averaged over the innermost % of the energy levels as a function of , with ranging from 100 to 40, such that the influence of the edges of the spectrum is progressively reduced.
The energy-resolved level spacing ratio [ Fig. 3(a)] signals spectral chaos within a slightly bent oval region for 10 −2 1 and 0.1 0.9, where approaches the prediction of random-matrix theory, GOE = 0.5307. The emergence of chaotic spectral statistics depends strongly on the energy. In particular, the lower threshold shifts to larger (i.e., stronger tunneling) upon increasing the energy, and the width of the chaotic region (in terms of ) is reduced towards the edges of the spectrum. The shift of the center of the chaotic phase across the -plane can be correlated with the drift of the DOS maximum towards larger upon increasing observed in Fig. 2 and indicated by the green dash-dotted line in Fig. 3 (a).
The spectrum-averaged level spacing ratio in Fig. 3(b) also shows agreement with GOE for 10 −2 1, regardless of whether the average is taken over the full spectrum or just the bulk energy levels. Hence, spectral chaos occurs not only for a wide range of energies, but also for a large majority of the levels. The transition into the ≈ GOE plateau is rather sharp and becomes even sharper when averaging over a smaller percentage of the inner spectrum. This feature is most notably visible at the plateau's lower edge, and reflects the dependence of the onset of chaos on energy discussed above. Far outside the chaotic region and deep within its center, however, averaging over different percentages of the inner spectrum has only a weak influence on , which shows that, for these regimes, most levels belong to energy regions of similar spectral statistics.
As observed in Figs. 3(c)-(e), the boundary of the spectrally chaotic region correlates with a pronounced increase of the mean fractal dimensions ˜ , i.e., with a distinct delocalization tendency of the states in Fock space. In the interaction basis, ˜ grows rapidly from, e.g., ˜1 ≈ 0.5 to ˜1 0.9 in a small, -dependent range which follows the left boundary of the chaotic region. In the tunneling basis, however, given the conjugated character of the bases, the GFD surge unveils most clearly the upper boundary of the chaotic phase.
Most interestingly, the variation of the finite-size GFDs for close-in-energy eigenstates, captured by var ˜ , registers a pronounced depression in a parameter domain that correlates remarkably with the spectrally chaotic region. Although this variance shows some slight quantitative differences between the interaction and the tunneling basis in Fig. 3(c)-(e), in both cases it exhibits the same qualitative behaviour: It decreases dramatically by several orders of magnitude when entering the chaotic phase. This figure of merit thus appears as a rather sensitive probe which can serve as a basis-independent quantifier of chaos, according to our results.
The behaviour of mean and variance of˜reveals a clear connection between spectral chaos and eigenstate structure: When the spectrum indicates chaos in the sense of RMT, the involved eigenstates share a similar structure in Fock space, characterized by a marked delocalization tendency in both natural bases.
Like the variance, the skewness of the˜distribution also exhibits a qualitatively basis independent behaviour: In the -parameter space, skew ˜ undergoes a sharp increase (by at least one order of magnitude) at the boundaries of the chaotic region. This reflects that upon tuning not all eigenstates with similar energies flow into the chaotic region (i.e., toward higher˜values) at the same speed, which gives rise to an asymmetric tail in the˜distribution, and hence to an enhanced skewness. Within the chaotic region, on the other hand, the fractal dimension distribution is markedly symmetric for = 1, 2: At the center of the chaotic phase skew ˜1 For = ∞, however, the reduction of the skewness in the chaotic region is far less pronounced, and the corresponding distribution there remains distinctively asymmetric as we show in more detail in Sec. VI.
As can be seen in Figs. 3(c)-(e), the three quantifiers of eigenstate structure considered, ˜ , var ˜ , and skew ˜ , exhibit the same qualitative behaviour for increasing , albeit the contrast of the chaotic phase is reduced, e.g., the minimum of var ˜ becomes less deep for larger (note that the same color scale applies to all values). Since˜∞ is determined by the largest intensity of each state [Eq. (18)], already the statistics over one specific intensity in Fock space reveals the emergence of chaos.

V. EGOE VERSUS BOSE-HUBBARD
In the following, we investigate how BHH relates to EGOE [Eqs. (9)-(11)] within the chaotic phase. This will be done first as a function of energy, for a fixed value of within the BHH's chaotic phase, and subsequently for fixed target energies changing .

A. Comparison of energy dependent features of chaos
In order to compare both models, from now on we fix = 1 in Eq. (9) for EGOE, since for this value the chaotic behaviour is fully developed, meaning that the studied figures of merit do not change upon further increase in , for any value of the energy , see appendix A. For BHH the scaled tunneling strength is set to = 0.1905, which is indicated by vertical dashed lines in Fig. 3 and is deep within the chaotic phase. Figure 4 provides an overview of the spectral and eigenstate structural features, quantified in terms of the finite-size GFDs, of the two Hamiltonians as functions of , for the = −1 Fock subspace of the = = 10 system (N = 46126). Mean level spacing ratios and moments of˜distributions are obtained after discretizing the energy axis, as explained above, and for EGOE they follow from an additional average over 100 ensemble realizations. It must be stressed that the ensemble-averaged EGOE spectral and eigenstate properties are symmetric around = 0.5 by definition [Eqs. (9)- (11)], while such symmetry does in general not hold for the BHH (it emerges only for sufficiently large ). For the value considered, the BHH's DOS starts approaching such symmetric behaviour, which enables the model comparison for the same . Nonetheless, the BHH features may naturally exhibit an asymmetry at the spectrum's tails which cannot be quantitatively captured by EGOE.
In contrast to GOE, whose DOS is dictated by the semi-circle law [6], the EGOE and the BHH energy level distributions [ Figs. 4(a) and (b), respectively] are found to exhibit a comparable Gaussian-like shape, centered around = 0.5 ( = 0.45) for EGOE (BHH). Close inspection reveals clear deviations from Gaussianity in the DOS that are in both cases accurately described by an Edgeworth expansion up to second order [cf. blue versus black lines in panels (a) and (b)] [35], in which a Gaussian is multiplied by a polynomial determined by the third and fourth central moments. Such expansion has also been found to apply to fermionic embedded ensembles [35].
The chaotic phase in both models is revealed by the agreement of the mean level spacing ratio with the GOE value in the range 0.2 0.8, flanked by strong fluctuations of at the spectral edges [ Fig. 4(c) and (d)]. Note that the overall smoother behavior of for EGOE is a consequence of the additional ensemble average. The analysis of shows that EGOE captures the absence of chaos at the spectral edges in qualitative agreement with the BHH. This is in great contrast to GOE, whose spectral (as well as eigenvector) properties do not carry any dependence on energy [6,10,107], and suggests that sparse Fock space connectivity must be responsible for the features observed at the spectrum's tails.
Within the spectrally chaotic region, for the three cases analyzed (EGOE and BHH in both natural bases), ˜ rise to their maximum values, and simultaneously var ˜ achieve their minima after dropping by several orders of magnitude, as shown in Fig. 4 This behaviour is accompanied by a levelling of var ˜ around very low values that agree best between EGOE and BHH in the interaction basis.
As discussed in the previous section, the skewness of the finite-size GFDs was particularly sensitive to the boundaries of the chaotic phase. This sensitivity is manifestly shown as a function of the energy in panels (k)-(m). In both models, the absolute skewness for = 1 and = 2 is suppressed within the chaotic region, with comparable values in all cases, and maximal at its borders, which for BHH also correlate with a clear increase in the uncertainties. In contrast to = 1, 2, the skewness for = ∞ is not suppressed within the chaotic region but rather acquires a common value skew ˜∞ ≈ 0.7 in all three cases, indicating the persistent asymmetry of the˜∞ distribution.
Although EGOE captures the major qualitative features of the finite-size GFDs' behaviour in the BHH, it must be noticed that the agreement is manifestly better with the analysis in the interaction basis, despite the fact that the EGOE Fock space connectivity is far closer to the one of the BHH in the tunneling basis [ Fig. 1]. In BHH, the correlations among the nonzero connectivity amplitudes in Hilbert space are much stronger than for EGOE, since BHH is defined by the parameters and only, while EGOE depends on all independent entries of the matrices (1) and (2) . Our findings thus show that these correlations also play a significant role. Figure 5 shows the evolution of ˜1 , var ˜1 and skew ˜1 for the BHH as functions of around specific target energies ranging from 0.2 to 0.8 and system sizes = ∈ [7,13]. These statistical quantifiers are obtained from the 100 eigenstates closest to the considered energy.

B. Chaos around specific target energies
In the interaction basis, for all system sizes and values, ˜1 registers a sharp surge upon increasing that tends to develop into a plateau as grows. For each system size, the onset and width of the plateau depend on energy, in accordance with the chaotic phase qualitatively identified in Fig. 3. The emergence of the plateau is mirrored by a synchronous drop in the variance of˜1 by several orders of magnitude, towards a flat minimum which becomes ever lower with increasing . This behaviour is also accompanied by a distinctive suppression of the skewness of the˜1 distribution. As discussed above, the emergence of chaos can be most clearly observed in the behaviour of var ˜1 . The chaotic phase is more prominent [lower var ˜1 and wider plateaux] for ≤ 0.6 for which there is always a range of values where the corresponding falls within the bulk of the spectrum, as can be observed in Fig. 3(a) from the evolution of the inner 40% of energy levels as a function of . Energy = 0.8 can be considered to belong to the tail of the spectrum for all , and nonetheless chaos leaves a distinct imprint in the variance of˜1.
For increasing system size, the left edges of the plateaux show a tendency to converge to an independent value (although depending on ), which confirms that the semiclassically inspired [18,28,73]  transition into the chaotic phase. Remarkably, upon increasing the plateaux' right termination point shifts towards larger [108], indicating that the chaotic region might extend to arbitrarily large (arbitrarily small interactions) in the thermodynamic limit (with a discontinuity at the integrable point = ∞). This latter observation agrees with results for fermions [109], spin systems and hard-core bosons [110,111], where the critical perturbation strength to drive an integrable system into chaos was found to decay linearly with system size.
In order to compare the BHH chaotic behaviour at fixed against EGOE, we must take into account the -dependent shift of the spectrum bulk for BHH. In other words, the identification of energy bulk and tails for both models is in general not a one to one correspondence in : The EGOE maximum DOS is always located at = 0.5, while this is not the case for BHH, even in its chaotic phase. We will take the position of the DOS maximum as the bulk centre. Therefore, once we have located the chaotic regime in at fixed for BHH [mainly identified by the behaviour of var ˜1 ], we analyze the distance of from the bulk centre for such range, and find the corresponding energies in EGOE (which we refer to as EGOE in the caption of Fig. 5) which lie at the same distance from its bulk centre (see appendix A). In this way we compare on an equal footing for both models the influence of energy deviation from the bulk centre on the considered markers of chaos.
As seen in Fig. 5, the GFD figures of merit within the BHH chaotic phase for all considered agree remarkably well with the EGOE values for the corresponding system size. This shows that the change of the BHH chaotic phase across the energy spectrum is correctly captured by EGOE, once the above subtlety in the energy comparison is taken into account.
In contrast to the interaction basis, for the system sizes studied, the formation of ˜1 plateaux in the tunneling basis is not so clear, and arguably only for = 0.2 can they be seen in agreement with the EGOE value. Despite this, and most interestingly, the development of the chaotic phase for increasing is unmistakably exposed for all by the behaviour of the variance and the skewness of˜1, although agreement with EGOE is not observed for > 0.2. This reflects the natural basis dependence of chaotic eigenvector features, which might also translate into different scaling behaviours with system size.

VI. SCALING WITH HILBERT SPACE DIMENSION
We now analyse the scaling behaviour of the GFDs for BHH and EGOE with Hilbert space dimension N , in comparison to analytical results for GOE, and determine whether or not the models approach each other at the level of the full GFD distribution within the chaotic regime. Figure 6 shows the dependence on N of ˜ and var ˜ of both models for three pairs of ( , ). We restrict to ≤ 0.6 for which, as discussed above, the BHH chaotic phase is more prominent. The corresponding is fixed to ensure that var ˜1 for BHH reaches minimal values either in the tunneling [(0. the following analytical results for GOE [63,108,112]:

A. Scaling of mean and variance
var ˜1 GOE = var D GOE where ℎ = =1 −1 = ∫ 1 0 (1 − )/(1 − ) d is the harmonic number, (1) is the first derivative of the digamma function [113], and the corresponding quantities for = ∞ are obtained using for integer . The asymptotic behaviour that ensues is var which also applies toD GOE 2 , and var ˜∞ GOE ∼ ln −4 N .
As observed in Fig. 6, ˜ and var ˜ are very close for both random ensembles and almost coincide for the largest accessible N . We can then conclude that the EGOE values essentially follow the GOE behavior, in particular that the EGOE dominant finite-size corrections are given by Eqs. (24)- (27).
In the tunneling (for = 0.2) and interaction basis (for = 0.4, 0.6), ˜ and var ˜ approach rather quickly, for all and all symmetry-induced subspaces, the corresponding (E)GOE values when increasing N . Such convergence can also be inferred in the interaction basis at = 0.2, albeit much slower. For these cases, the numerical evidence indicates that ˜ [var ˜ ] takes the ergodic value 1 [0] in the thermodynamic limit, and does so according to Eqs. (24)- (27).
In contrast, in the tunneling basis at = 0.4 and 0.6, both statistical quantifiers show pronounced quantitative deviations from the (E)GOE trajectories (as also observed in Fig. 5), to the point that a linear extrapolation in 1/ln N of the available ˜ data actually yields lim N→∞ ˜ < 1, indicating that the dominant finite size correction of BHH in this case must have a different dependence on N than (E)GOE. This observation clearly manifests that the development of the chaotic phase in the eigenvector features and the scaling trajectory followed towards ergodicity depend strongly on the basis considered.

B. Comparison of the full distributions
In Ref. [63] we have shown that, while ˜ and var ˜ of BHH eigenstates converge towards GOE, both models depart from each other at the level of the full GFD distributions with increasing N . The corresponding analysis for BHH and EGOE will determine whether the two-body Hamiltonian structure is what sets BHH apart from GOE.
As the left panels in Fig. 7 reveal, the GFD distributions of BHH and EGOE drift apart as N grows. In the right panels, we quantify the distance between the distributions by the relative difference of the means, and by the square root of the Kullback-Leibler divergence (information entropy) [114,115] (P, Q) = Both distance measures increase with Hilbert space dimension for BHH versus either GOE or EGOE for all considered, although the growth is less pronounced for = ∞. Note that the EGOE distribution is closer to BHH than GOE, and hence the distance (EGOE-BHH) is always bounded from above by (GOE-BHH). For N ≥ 10 6 the mean˜1 ,2 of BHH differs from GOE already by 10 , and, by extrapolation, the distance from EGOE, although smaller, will be of comparable magnitude. We therefore conclude that the statistical departure between BHH and GOE in terms of the GFD distributions, cannot be solely explained by the two-body nature of BHH and the consequent Fock space sparsity, since these features are also present in EGOE. The observed deviation must then be ascribed to the very specific nature of the BHH on-site interactions and nearest-neighbour tunneling, and the existing correlations among the Hamiltonian matrix elements.
Interestingly, an increasing trend of and for = 1, 2 is also observed for EGOE versus GOE, whereas for = ∞ the GFD distributions of both models appear to approach each other in the accessible N range. Hence, despite the quantitative agreement exhibited by ˜ and var ˜ , random two-body Hamiltonians become more distinguishable from full GOE in terms of P (˜1 ,2 ) when approaching the thermodynamic limit.
The comparison of the numerical results for BHH (with respect to the interaction basis, at the center of the spectrum) and EGOE with the analytical GOE predictions will help us understand the origin of the departure of P (˜1 ,2 ) in terms of the finite-size corrections to ˜ . As discussed in section VI A, the numerical results suggest that the functional dependence of the dominant correction for all models is given by Eqs. (24)- (27). Consequently, since √︃ var ˜1 ,2 BHH ∼ 1/ √ N ln N and 1 − ˜1 ,2 ∼ 1/ln N for all models, we can single out three main possibilities leading to three different scaling behaviours of : (i) the coefficient of the dominant correction differs between GOE and the other models, which leads to 1,2 ∼ √ N [indicated by solid lines in the right panels of Fig. 7], (ii) the dominant correction is exactly the same in all models, and BHH and EGOE carry subleading terms in ˜1 ,2 decaying slower than 1/ √ N ln N , which entails an increasing behaviour of 1,2 slower than √ N , and (iii) the dominant correction is exactly the same in all models, and the second term in BHH and EGOE decays faster than 1/ √ N ln N , which would imply a decreasing behaviour of 1,2 , and hence an approach of the distributions.
The tendencies shown in the right panels of Fig. 7 suggest that case (ii) holds. This implies that asymptotically 1 − ˜1 ,2 = 1 /ln N , with the same 1 for the three models. Therefore, the increasing distinguishability of the distributions cannot be attributed to the so-called 'weak ergodicity' [93]. In fact, as described above, a common 1 is still compatible with either a departure or a convergence of the distributions, a behaviour that is determined by the functional dependence of the second finite N correction.
A similar argumentation applies to = ∞: Were the coefficient of the dominant correction in ˜∞ to differ between GOE and the other models, one should observe ∞ ∼ ln(ln N ) ln N . In this case however, given the slow logarithmic increase, it is more challenging to draw a definite conclusion about BHH. We also note that the second and third GOE corrections decay slower than √︃ var ˜∞ ∼ ln −2 N [63], which implies that an increasing behaviour of ∞ is still compatible with BHH having the first three finite-size corrections in ˜∞ dictated by GOE. For EGOE on the other hand, the decreasing behaviour of ∞ indicates that both random ensembles do share at least all correction terms that decay slower than ln −2 N .

VII. CONCLUSIONS
We have scrutinized the properties of the chaotic phase of the Bose-Hubbard Hamiltonian (BHH) via spectral and eigenvector features in relation to the bosonic embedded Gaussian orthogonal random-matrix ensemble (EGOE). The emergence of chaos at the level of individual eigenstates as well as statistically over energy intervals for both natural bases of BHH has been analyzed in detail. We have confirmed that the fluctuations of the fractal dimensions for close-in-energy eigenstates are an efficient, and qualitatively basis independent, probe of chaos. Furthermore, the boundary of the chaotic region is signalled by a strongly skewed distribution of the fractal dimensions for close-in-energy eigenstates.
The analysis of the chaotic phase for EGOE in terms of spectral statistics and eigenstate fractal dimensions reveals the same dependence on energy as for BHH, capturing in particular the disappearance of chaos at the spectral edges, in contrast to GOE. We therefore attribute this feature to the sparse Fock-space connectivity of the Hamiltonian. EGOE describes remarkably well the typical value and the fluctuation of fractal dimensions of BHH eigenvectors in the interaction basis at different energies from the centre of the spectrum's bulk. However, the agreement between EGOE and BHH in the tunneling basis is more restricted, despite the close matrix structure of both Hamiltonians.
Scaling for increasing Hilbert space dimension N confirms the convergence to ergodicity of the BHH eigenvector characteristics in the chaotic region, and also exposes the strong basis dependence of such process for certain parameters. The scaling analyses suggest that the coefficient of the dominant asymptotic ln −1 N correction that governs the convergence of the mean fractal dimensions ˜1 ,2 to 1, for BHH in the interaction basis and for EGOE, coincides with the one for GOE. Notwithstanding this agreement, the three models become more distinguishable from each other in terms of the fractal dimension distributions as Hilbert space dimension grows, in a similar way as we found previously for BHH and GOE [63]. This behaviour shows that the distinctive statistics of BHH in the chaotic regime cannot be traced back solely to the two-body nature of the Hamiltonian and the associated Fock-space sparsity, but must be connected to specific correlations among the Hamiltonian's matrix elements. These results provide further evidence of a way to discriminate among different many-body Hamiltonians in the chaotic regime, where coarse-grained dynamics are governed by common universal features.
The emergence of the chaotic phase in Hamiltonian (9) as a function of is depicted in Fig. 8 in terms of ˜1 and var ˜1 , obtained from 100 eigenstates for different and averaged over 100 ensemble realizations. While for 0.1, ˜1 clearly increases and var ˜1 decreases with , both quantities essentially converge to a constant value for 0.5 and all considered. We have checked that this holds true also for˜2 ,∞ . Hence our choice of = 1 in the main text. Furthermore, the behaviour observed for and 1 − is very similar, in agreement with the symmetry of the ensemble averaged properties of EGOE around = 0.5 according to its definition, Eqs (9)- (11) [see also Fig. 4].
The chaotic properties of BHH and EGOE cannot be compared around the same value, due to the shift of the bulk of the , and alternatively derived from the condition to bound the same percentage of the total spectrum from above (squares). Dashed lines and thick shades on the axes highlight the intervals associated with the BHH chaotic regions in Fig. 5 and the corresponding EGOE ranges. spectrum for BHH with . In order to enable the comparison, we determine a correspondence between BHH and EGOE in two different ways: First, one may demand that BHH and EGOE be located at the same distance from the corresponding bulk centre, which we identify with the position of the DOS maximum, * BHH ( ) and * EGOE = 0.5, respectively. Then, for chosen BHH and , we have EGOE ( BHH , ) = BHH − * BHH ( ) + 0.5.
Alternatively, the energy equivalence could be established by requiring that EGOE and BHH have the same position in terms of percentage of the spectrum, i.e., that they bound from above the same fraction of the associated spectra. Figure 9 shows the energy EGOE as a function of for four different BHH values according to the two methods described, highlighting the ranges of EGOE corresponding to the intervals associated with the chaotic region of BHH in Fig. 5. Even though the two definitions of EGOE differ for some , they both yield a very similar EGOE interval for the chaotic region of BHH at the considered BHH . Note that due to the discussed symmetry of EGOE, the EGOE intervals for BHH = 0.2 and 0.4 can be truncated at 0.5.