Spectral statistics of chaotic many-body systems

We derive a trace formula that expresses the level density of chaotic many-body systems as a smooth term plus a sum over contributions associated to solutions of the nonlinear Schr\"odinger (or Gross-Pitaevski) equation. Our formula applies to bosonic systems with discretised positions, such as the Bose-Hubbard model, in the semiclassical limit as well as in the limit where the number of particles is taken to infinity. We use the trace formula to investigate the spectral statistics of these systems, by studying interference between solutions of the nonlinear Schr\"odinger equation. We show that in the limits taken the statistics of fully chaotic many-particle systems becomes universal and agrees with predictions from the Wigner-Dyson ensembles of random matrix theory. The conditions for Wigner-Dyson statistics involve a gap in the spectrum of the Frobenius-Perron operator, leaving the possibility of different statistics for systems with weaker chaotic properties.


Introduction
Feynman's path integral (see e.g. [1]) provides a convenient way to represent the propagator of a quantum mechanical system, and an excellent starting point for semiclassical and related approximations.
Prime examples are van Vleck's approximation of the propagator of a quantum system as a sum over contributions of classical trajectories [2], and Gutzwiller's seminal work [3] relating the energy spectrum of chaotic single-particle systems to periodic classical trajectories. These semiclassical methods provide one of the foundations of the field of quantum chaos [4,5,6]. For a many-particle system identifying the semiclassical limit is less obvious. A promising approach is to consider the path integral in second quantisation, running over different choices for the macroscopic wave function parameterised by position and time. One can show that in the semiclassical limit and in the limit where the number of particles N is taken to infinity this path integral is dominated by stationary points of the action, corresponding to solutions of the nonlinear Schrödinger equation or Gross-Pitaevski equation. These solutions take on a role analogous to the one played by classical trajectories in the single-particle theory. However previous work usually focused either on studying the full problem in second quantisation or, as e.g. in nuclear dynamics [7], on one dominating solution to the Gross-Pitaevski equation [8,9]. This does not exhaust the power of the approximation. In particular keeping the sum over different solutions of the Gross-Pitaevski equation allows to account for crucial interference effects between such solutions. An approach where the semiclassical propagator of bosonic many-particle systems was used to study these effects was pioneered in [10] for coherent backscattering, see also [11] for applications to fermionic systems.
In the present paper we will focus on a further fundamental problem for which the interference between stationary points of second quantised path integrals is of vital importance: the statistics of the energy levels of many-body systems. To do so we will first derive an approximation of the level density in terms of stationary points of the action, and then study the interference between these points. To our knowledge the consequences of interference effects for many-body spectral statistics have not yet been investigated explicitly. We will see that this statistics depends crucially on the dynamics generated by the Gross-Pitaevski equation. If the dynamics is fully chaotic (in the sense to be specified below) the statistics in the limits considered becomes universal, and agrees with predictions from random-matrix theory (RMT). These predictions entail, for instance, a repulsion between the energy levels.
This universal behaviour mirrors the behaviour of chaotic single-particle systems studied in the semiclassical limit. For such systems spectral statistics faithful to random matrix theory was conjectured in [12], and a semiclassical explanation was developed in [13,14,15,16,17,18]. This explanation is based on Gutzwiller's trace formula [3,4,5,19] d(E) ≈d(E) + 1 π Re which expresses the level density as a smooth termd(E) plus fluctuations associated to classical periodic orbits p (with energy E) of the system. Here S p is the reduced action of the orbit given by S p = p(t) ·q(t)dt where q(t) denotes the vector of generalised coordinates and p(t) denotes the associated momentum. The amplitude depends on the primitive period T prim p , the so-called Maslov index µ p and the stability matrix M p relating deviations in the end of the orbit p to deviations in the beginning; I is a unit matrix. (Throughout this paper I will denote unit matrices with a subscript indicating their size if it is not clear from context.) Our first aim will be to generalise the trace formula to bosonic many-particle systems in second quantisation, with solutions of the Gross-Pitaevski equation taking the role of classical trajectories. We will then use this result to investigate spectral statistics. An observable characterising spectral statistics is the two-point correlation function of the level density d(E) = j δ(E − E j ) (where E j are the energy levels of the system). This correlation function is defined by where the average . . . is taken over an interval of E for whichd can be taken constant as well as over a small range of . Inserting the trace formula one obtains a double sum over solutions of the Gross-Pitaevski equation, and by taking into account interference between solutions we indeed recover statistics in agreement with RMT. More precisely this agreement holds for the statistics inside appropriate subspectra defined by the symmetries of the problem; for many-body systems we have at least one symmetry, particle number conservation, requiring to consider subspectra associated to a fixed particle number. Further refinements (to be discussed below) arise in case of geometrical symmetries. The precise ensemble to be chosen depends on the behaviour of the system under time reversal. The most frequent case involves systems invariant under a time-reversal operator squaring to 1. In this case (assuming that there are no further symmetries) one has to use Wigner's and Dyson's Gaussian Orthogonal Ensemble (GOE), i.e. predictions for spectral statistics are obtained by modelling the Hamiltonian through a real symmetric matrix, and then averaging over all possible such matrices with a Gaussian weight. In the absence of time-reversal invariance one averages instead over the ensemble of Hermitian matrices with a Gaussian weight, the Gaussian Unitary Ensemble (GUE).
A paradigmatic example for the systems to be considered is the Bose-Hubbard model, a model with L discrete sites (labelled by k = 0 . . . L−1) accommodating bosonic particles. Denoting the creation and annihilation operators for particles at these sites byâ † k andâ k the second-quantised Hamiltonian can be written aŝ describing hopping between the sites as well as interaction between particles on the same site. One can consider periodic boundary conditions and setâ L =â 0 ,â † L =â † 0 . More generally we are interested in discrete bosonic many-body Hamiltonians that are of the formĤ = kl h klâ † kâ l + klmn U klmnâ † kâ † lâ mân (4) or have even higher-order interactions; here real coefficients imply time-reversal invariance.
A crucial requirement is that the underlying 'classical dynamics' is chaotic. This dynamics is obtained by replacing the creation and annihilation operators by mutually complex conjugate time dependent variables ψ * k , ψ k where ψ = (ψ 0 , ψ 1 , . . .) can be interpreted as a macroscopic wave function and ψ k as its value at site k. AsN = kâ † kâ k is the particle number operator, the macroscopic wave function is normalised to have k |ψ k | 2 = N where N is the particle number. The associated analogue of Hamilton's equations for ψ k can be shown to read i ψ k = ∂H ∂ψ * k and takes the role of a discrete nonlinear Schrödinger equation.
We note that the 'classical' Hamiltonian entering into this equation depends on the ordering of operators. For the Bose-Hubbard model the normal-ordered Hamiltonian in (3) yields an interaction term U 2 k |ψ k | 4 ; however if the Hamiltonian is first brought to Weyl-ordered form (with all possible orderings of operators in a product contributing symmetrically) before replacing the operators by macroscopic wave functions one obtains . For the Bose-Hubbard model the dynamics generated by the discrete nonlinear Schrödinger equation has been found to be mainly chaotic (with chaotic regions of phase space dominating compared to regular ones) in the case of several sites and comparable hopping and interaction terms [20]. In the same regime, numerical studies suggest spectral statistics in line with the GOE [21,22]; see also [23] for fermionic systems.
To explain the observed faithfulness to RMT, we follow a semiclassical approach inspired by single-body spectral statistics as well as [10]. Our first aim is to derive a trace formula for second quantised Bose-Hubbard-like systems. This will be done in Sect. 2, after a brief reminder of the corresponding derivation for one-body chaotic systems. Special emphasis is placed on the treatment of conserved quantities. In Sect. 3 it is demonstrated how the obtained trace formula can be used in order to predict the spectral statistics of the system, especially the 2−point correlation function. We explain in detail how to generalise the approach previously used for one-body chaotic systems in first quantisation, and we discuss the range of validity of this approach. In Sect. 4 the consequences of discrete symmetries of chaotic many-body system, in particular for Bose Hubbard model are investigated. In Sect. 5 we present numerical results supporting our claims. Some more technical details of the derivation are put in two appendices.

Trace formula for single-body systems
In order to prepare our derivation of a trace formula for bosonic many-particle systems, we want to briefly review the calculation leading to the trace formula for single-body systems. We refer to [4,5,6,24,25] for further details.
The level density can be accessed from the trace of the time evolution operator where an infinitesimal positive imaginary part has been added to the energy to ensure convergence. The trace of the time evolution operator can itself be expressed as a path integral over phase space trajectories (q(t), p(t)) where the action weighting each path is determined by the classical Hamiltonian H as Now we are interested in an approximation of this path integral in the semiclassical limit, i.e., the limit of large quantum numbers implying that typical classical actions are much larger than ; formally this limit is often denoted by → 0. In the semiclasssical limit the path integral is dominated by stationary points of the action, corresponding to periodic orbits that satisfy Hamilton's equations of motion, and a stationary-phase approximation leads to a discrete sum over such periodic orbits. For chaotic dynamics, free from continuous symmetries, the periodic orbits are isolated but one has to take into account that each of them formally gives rise to a one-parameter family of stationary points distinguished by which phase-space point along the orbit is taken as the initial one associated to t = 0. The Laplace transform in (5) can then be carried out in a further stationary-phase approximation and one obtains the trace formula already shown earlier.
Here the determinant and the phase factor involving the Maslov index arise from the Hessian matrix of the action entering the stationary-phase approximation, and the factor T prim p results from integration over different choices of initial points. The summandd(E) can be derived from a careful treatment of the lower limit of the time integral (5).

Path integral for many-particle systems
We now want to generalise the trace formula to many-particle systems in second quantisation, for the case of bosonic particles at discrete sites. In second quantisation it is natural to work in a basis of coherent states, i.e., normalised joint eigenstates of all annihilatorsâ k with eigenvalues ψ k . Further differences from the single-particle setting arise from operator ordering as well as the conservation of the particle number.
Again we access the level density from the trace of the time evolution operator e − i Ĥ t using (5). For many-particle systems the latter trace is given by the path integral [26] tr over all macroscopic wave functions ψ(t ) that return to their initial value after time t.
Here the action is Eq. (9) is related to but somewhat simpler than the path integral for matrix elements of coherent state time evolution operator [27]; our use of tr e − i Ĥ t is motivated by [24,25].
In case of normal ordering the path integral can be derived by splitting the time interval t into J steps of width τ = t J and then using the result for the short-time propagator after integration over the macroscopic wave functions at all time steps this leads to a discrete path integral with the action and then to (10) after taking the limit J → ∞. To fix notation, throughout the paper j will be a time index, k will label degrees of freedom e.g. associated to sites, and bold vectors assemble all choices for k. j is defined modulo J and k is taken modulo L. With these conventions the integration measure in the path integral is given by jk . The integral thus obtained agrees with the phase-space formulation of the path integral in first quantisation if the macroscopic wave function is related to canonical coordinates and momenta through ψ k = 1 √ 2 (q k + ip k ) or ψ k = I k e −iθ k . We can now perform a semiclassical approximation valid in the limit → 0. Here ψ k is taken of the order 1 √ such that the first term in the action (10) becomes independent of ; the coefficients in the Hamiltonian are scaled following J/ → J, U/ 2 → U in such a way that the Hamiltonian is independent of as well (We will later also give an interpretation in terms of the limit N → ∞.) In our limit the integral is dominated by periodic solutions of the nonlinear Schrödinger equation i ψ k = ∂H ∂ψ * k (and its complex conjugate) following from the stationarity of R.

Particle number conservation
Importantly, however, our Hamiltonian has a continuous (gauge) symmetry w.r.t. multiplying all components of the macroscopic wave function with the same phase factor. A consequence of this symmetry is that as mentioned the total particle numberN is a conserved quantity commuting with the Hamiltonian. It is thus preferable to consider the density of levels forming the subspectrum associated to a fixed particle number.
To implement this restriction we subject the I k , θ k defined above to a canonical transformation. (See [28,9] for alternative approaches.) This transformation is chosen to lead to I k , θ k where I 0 = k I k = k |ψ k | 2 = N , the remaining I k are linear combinations of the I j , and the θ k are defined such that the overall transformation becomes canonical. If we choose the transformation in such a way that the range of possible θ k is limited to 2π it is also convenient to let ψ k = I k e −iθ k . As I 0 is a conserved quantity the corresponding canonical coordinate θ 0 must be absent from the Hamiltonian. The remaining variables with k ≥ 1 parameterise a reduced phase space associated to the particle number N . We are interested in the part of the spectrum associated to a given value N of the particle number. The associated level density can be formally written as with the Kronecker delta As one might expect, d N (E) is determined by solutions periodic in the reduced phase space. To see this formally we split the exponential e iφN from (14) into factors e iφN /J to be inserted into each of the short-time propagators (11). This gives an additional termR to be added to the action, leading tõ in the limit J → ∞. The phase now becomes stationary if where the addition only affects the dynamics of θ 0 , replacing it bẏ Hence, as anticipated, the requirement of periodicity is nontrivial only for the variables ψ k with k ≥ 1. I 0 is conserved and hence trivially periodic, and θ 0 is required to be periodic only after modifying the dynamics through the φ-dependent term in the action, but all possible choices of φ are integrated over.

Determinant of the Hessian matrix
We now have to determine the weight associated to each periodic solution. We recall that if the stationary points p of a given action R[x] (with x ∈ R n for now) are isolated and we are taking the limit → 0, the integral over e iR/ can be approximated by the following sum over contributions associated to stationary points, here µ p is the number of the negative eigenvalues of the Hessian matrix ∂ 2 Rp ∂x 2 for the stationary point p. If the stationary points are not isolated this leads to vanishing eigenvalues of the Hessian. In this case the integral over the directions orthogonal to the stationary-point manifold can still be computed using (19), but it has to be accompanied by an integral over the manifold itself.
The stationary points of (12) are not isolated. In particular continuous time shifts of a solution of the nonlinear Schrödinger equation lead to a different solution, and the same applies to simultaneous shifts of all θ j0 . This is a consequence of the two conserved quantities, for the energy conjugate to time and the particle number conjugate to the θ j0 . However as a first step it is still helpful to compute the (rescaled) Hessian involving derivatives w.r.t. all components of the macroscopic wave functions at all time steps. Adopting a complex notation we consider Using the discretised action (12) the derivatives are given by where I = I L is a unit matrix of dimension L. The Hessian w.r.t. these variables then assumes periodic block tridiagonal form where n = 2J and the blocks are matrices of size L × L given by Here we neglected corrections of order τ 2 arising from the fact that the arguments of the Hamiltonian in (21) are taken at slightly different times. We can now use a general formula for determinants of block tridiagonal matrices as in (22) that was derived in [29] using a transfer matrix approach, Here we have In Eqs. (24), (25) we have slightly modified the numbering of indices from [29] and taken the dimension of the block matrices as L. Due to n = 2J the sign factor in (24) is just equal to 1. To evaluateM we group the factors in (25) into pairs which using (23) can be simplified tõ One can now show that multiplication withM j maps small deviations (δψ * j , δψ j ) from a given solution of the nonlinear Schrödinger equation at time jτ to the resulting deviation at time (j + 1)τ . This follows immediately by linearizing around the equation and its complex conjugate. Hence the product (understood in the limit J → ∞) maps deviations in at time 0 to those at time t. The matrix P is given by evaluating its determinant and taking the continuum limit then leads to The factor (det P ) −1/2 arising from this term in (detH) −1/2 is known as the Solari-Kochetov (SK) phase; it is in line with previous work about the propagator in normal ordering [30,31,27].

Treatment of the conserved quantities
We now modify our treatment to take into account particle number and energy conservation. For notational convenience it is helpful to complement our earlier canonical transformation singling out the particle number by a further transformation affecting only the variables with k ≥ 1. This transformation is defined only in the vicinity of a periodic solution and leads to I 1 indicating the energy on the orbit, and θ 1 the time along the orbit. The remaining variables I k , θ k with k ≥ 2 indicate transverse deviations from this orbit. If desired they can again be turned into complex variables ψ k as above , and the vector formed by all ψ k with k ≥ 2 will be denoted by ψ ⊥ . The present transformation is analogous to the transformation to parallel and perpendicular coordinates in the derivation of the trace formula in first quantisation [3,4,5]. We note that it cannot be expanded to a transformation in the full phase space as this would imply integrability.
When determining the weight associated to periodic solutions, it is now crucial to take into account that shifting all time coordinates θ j1 by the same amount leads to a valid solution, and the same applies to coordinated changes of the variables θ j0 conjugate to the particle number. Hence there are two linearly independent ways in which continuous changes from one stationary point of the action lead to a different This requires suitable rescaling to make both quantities dimensionless and restrict the range of θ k to 2π.
one. As second derivatives of the action in the associated directions must necessarily be zero the matrixH defined above must have a two-fold eigenvalue zero.
To compare this to the behaviour ofM − I 2L we first of all observe thatM maps deviations of the variables associated to k = 0, 1 only to deviations associated to the same k. Written in terms of θ k , I k the stability matrix multiplies deviations (δθ k , δI k ) withM Here the diagonal elements indicate that the conserved quantities I k stay fixed and changes of θ k are translated into equal changes in the end. In the right upper element ∆θ k indicates the increase of θ k along the orbit, e.g. the period for k = 1. The coefficient b k takes into account that a change of the energy typically changes the period of the orbit, and similarly a change of the particle number typically changes the difference of the initial and final θ 0 . As a consequence of (32), To deal with these zeroes one can consider a perturbation of the Hamiltonian [25] that replaces one of the zero eigenvalues of the Hessian by a small value . In the correspondingM (k) a factor J then enters in the lower left corner, turning the determinant into − J b k . A brief account of these perturbative results for the present case is given in Appendix A. The determinant of the Hessian matrix H omitting the directions associated to conserved quantities can now be evaluated by considering perturbations for both k = 0 and k = 1 and dividing out the two factors , leading to Here M is defined in analogy to (27) and (29) but only w.r.t. the variables ψ ⊥ omitting k = 0, 1. M maps initial deviations of ψ ⊥ * , ψ ⊥ to the corresponding final values. This meaning is precisely equivalent to the stability matrix appearing in the conventional trace formula.
Our result for det H allows to evaluate (in a stationary phase approximation) the integral over all momenta and coordinates apart from θ j0 , θ j1 , as well as the fluctuations of θ j0 , θ j1 as j is varied.
It remains to consider the constant (in j) Fourier modes of θ j0 , θ j1 . Importantly, if we perform a discrete Fourier transform of θ jk and want the associated Jacobian determinant to be 1, the integration variable parameterising the constant mode has to be chosen as j θ jk √ J , i.e.

√
J times the average of the θ jk . More natural parameterisations of the constant Fourier modes, such as by the average of the θ jk , entail a Jacobian √ J for each k = 0, 1. These jointly cancel the J −1 from (det H) −1/2 .
Integration over the constant modes now leads to multiplication with the integration ranges. For θ 0j the range is 2π, and for the time coordinates the range is normally the period. However if the periodic solution consists of several repetitions of a shorter 'primitive' periodic solutions, all distinctive stationary points of the action can be accessed by integration over the period T prim p of the primitive solution only. Equating T p = T prim p if our solution does not consist of repetitions of a shorter one we thus obtain a factor T prim p in all cases. We still have to evaluate the integral over φ where theR from (16) equals φI 0 . Due to dR dφ = I 0 the stationary phase condition gives I 0 = N as expected. Now due to (18) a change of φ is equivalent to changing the range ∆θ 0 by an opposite amount. This implies that d 2R The b 0 thus obtained cancels with the one from (33).
Altogether the trace of the time evolution operator, restricted to a given value of the particle number, is thus approximated by the following sum over solutions p of the nonlinear Schrödinger equation that are periodic with period T p = t in our reduced phase space Here various factors from the integration measure, Eq. (19), Eq. (14), and the θ 0 integral have cancelled mutually. The result is in line with [24,25]. The Maslov index µ p counts the negative eigenvalues of the part of the Hessian matrix H associated to k = 2, 3, . . . (when brought to symmetric form involving derivatives w.r.t. I jk , θ jk ). The index µ p takes a similar role for k = 1; it is equal to 1 if b 1 < 0 and 0 otherwise. An analogous phase associated to k = 0 has already been cancelled by the stationary-phase approximation of (34). For notational convenience we also use the Maslov index to absorb the Solari-Kochetov phase [30,31,27] We note that like the stability matrix also the action can be taken within our reduced phase space. The only difference between the two is the k = 0 contribution to the first term in (10), which is given by however for integer N this term has no influence on the phase factor e iR/ . Moreover, as the Hamiltonian is independent of θ j0 , the change in the dynamics of θ j0 mentioned above does not affect any of the quantities entering the trace formula.

Level density
Finally the Laplace transform of (35) can be performed in a further stationary-phase approximation analogous to [3,4,5] (and in fact similar to the φ integral above). The phase in e i(Et+Rp)/ becomes stationary if E = − dRp dt . However as shown e.g. in [4] − dRp dt is precisely the energy of the solution with period t. Hence the result of the stationary phase approximation will be a sum over all solutions that are periodic in the sense above and have fixed energy E rather than fixed period. For these solutions Et then cancels with the second term in the action (10), hence the associated phase will be determined only by the first term S p = − i t 0 dt ψ * (t )·ψ(t ) also referred to as the reduced action. Given that p π 2 arising from the stationary-phase approximation combines nicely with the one from (35). Altogether we thus obtain the anticipated trace formula The summandd N (E) gives the smooth part of the level density and arises from the lower limit 0 of the time integral. It is proportional to the volume of the energy shell in our reduced phase space and given by the pertinent variant of the Weyl or Thomas-Fermi formulad This result (which is more meaningful if one avoids the canonical transformation in the beginning of subsection 2.5) is readily obtained using that for small times no time slicing is necessary. The action (12) then boils down to −H(ψ * 0 , ψ 0 )t and the result follows directly after Laplace transformation. For the semiclassical approach to spectral statistics it will not be necessary to evaluate the -often involved -integral ford N (E). As in first quantisation [32] subleading corrections to this result are to be expected.

Discussion
Noting that |ψ| 2 = N the present approximation is valid not only for → 0 but also in the limit N → ∞. In the latter case it is helpful to rescale ψ → √ N ψ to avoid changing the normalisation of the variables in the limit taken. If we want the hopping and interaction terms in a Bose-Hubbard like system to remain comparable we then have to adjust the corresponding coefficients in a way similar to the case → 0. This now entails J → J, N U → U . In case of the Bose-Hubbard model the resulting Hamiltonian entering the trace formula then satisfies with U and J constant, the normalisation |ψ| 2 = 1 and the large parameter N and the whole action (12) being proportional to N .
If the Hamiltonian is written in Weyl instead of normal ordering one obtains an analogous result however without a Solari-Kochetov phase. This is in line with [24] as well as results for the propagator in [27,33]. A derivation along the lines followed here will be given in Appendix B; it uses the analogue of the discretised action stated in [33] and the corresponding Hessian is again evaluated with the help of [29].
An alternative derivation of the trace formula can be based on the semiclassical approximation [27,33,9] of matrix elements ψ (f ) |e −iĤt/ |ψ (i) in a basis of coherent states. Interestingly, in this case ψ(t ) and ψ * (t ) first have to be treated as independent functions subject to the conditions ψ(0) = ψ (i) , ψ * (t) = ψ (f ) * only. The two functions become complex conjugate after evaluating the trace in a stationary-phase approximation, and the resulting trace formula coincides with the one obtained above.

Spectral statistics
We are now equipped to study spectral statistics. Inserting the trace formula into the definition of the two-point correlation function one obtains, as for single-particle systems, the double sum The correlation function is thus expressed in terms of periodic solutions of the nonlinear Schrödinger equation in a way that allows to keep track of crucial interference effects. The double sum can be evaluated in the same way as for chaotic single-particle systems. We now want to discuss this evaluation in more detail. In doing so we will emphasise the ingredients entering the calculation and check the conditions under which the reasoning for single-particle systems carries through to many-particle systems in second quantisation. For the details of the calculation for single-particle systems based on these ingredients we refer to the original literature quoted below as well as [4,34].

Conditions
The phase space of predominantly chaotic many-particle systems typically still has small stability islands. Hence it is important to stress that our theory describes the behaviour of states supported by the chaotic part of phase space. The spectral statistics is dominated by this contribution if the regular parts of phase space are small in comparison, as in the case of the Bose-Hubbard model in the regimes considered.
In the chaotic part of the phase space our treatment requires a gap in the spectrum of the Frobenius-Perron operator. This condition implies various weaker requirements such as ergodicity and hyperbolicity. The Frobenius-Perron operator describes the time evolution of classical phase space densities [4], leading from a density ρ 0 (x) at time 0 to the density ρ t (x) = d n x δ(x − Φ t (x ))ρ 0 (x ) at time t. Here Φ t (x ) gives the image of the phase space point x under classical time evolution over time t. The leading eigenvalue of this operator is 1, with the associated eigenfunction corresponding to a uniform density on the energy shell. The remaining eigenvalues can be written as e −γmt where γ m with Re γ m ≥ 0 are the Ruelle-Pollicott resonances. The system is said to have a spectral gap if the remaining Re γ m are bounded away from zero so that all modes associated to non-uniform phase space densities decay in time with at least a minimal rate. For many-particle systems the dynamics required to satisfy this condition is the one induced by the discrete nonlinear Schrödinger equation in the reduced phase space parameterised by I k , θ k (k = 1, 2, . . . , L − 1).
As a further condition we assume for now that our system has no further symmetries beyond the particle number conservation already taken into account, however we will extend our treatment to deal with discrete geometric symmetries at a later point.

Diagonal approximation
When evaluating (42) it is important to look for pairs of solutions whose (reduced) actions S p , S p are similar as this means that their contributions can interference constructively. In contrast, terms in (42) arising from pairs of orbits with large action differences oscillate rapidly as the energy is varied and are washed out by the energy average.
The simplest pairs of solutions with similar actions involve two solutions that are identical (apart from the slight difference due to the offset in their energy arguments). If we neglect the difference of the corresponding factors A p , A p and Taylor expand the exponent using that dSp dE gives the period T p , the contribution from identical solutions can be written as the single sum Under the condition of a spectral gap such sums over periodic orbits can be evaluated using a sum rule derived in the quantum chaos context by Hannay and Ozorio de Almeida [14]. In the notation used here it can be written as where the dots represent an arbitrary property of the solutions that depends only on their period T . This rule is a general statistical property of periodic solutions in systems with a spectral gap; it is very helpful to extract information from these solutions even in situations where it is difficult to determine the solutions individually. Eq. (44) implies that even very long orbits give important collective contributions: while the factors |A p | 2 associated to these orbits decrease with increasing period their number increases, and both effects approximately compensate. Using (44) the sum in (43) can be evaluated to give This result, originally derived in [13,14], is known as the diagonal approximation. It gives the first nontrivial term in an expansion of the two-point correlation function in 1 , after the leading term 1 present in (42). For time-reversal invariant systems we also have to keep track of mutually time reversed solutions, leading to a doubling of this result.

Encounters
Contributions of higher order in 1/ arise from pairs of orbits differing noticeably only in so-called encounters [16,17]. Inside these encounters two or more parts of the same orbit come close up to time reversal, and a partner orbit can then obtained by connecting the ends of these orbit parts differently and a subsequent adjustment to obtain a valid periodic orbit satisfying the equations of motion. For time-reversal invariant systems one also has to include the case where parts of an orbit are almost mutually time reversed. Two simplified sketches of such pairs of orbits are shown in Fig. 1.
The existence of such pairs of orbits requires hyperbolicity which follows from the existence of a spectral gap. In a hyperbolic system the possible directions at each point in phase space (apart from the direction of the flow and the direction of increasing energy) are spanned by pairs of stable and unstable directions. Deviations between (parts of) trajectories along the stable directions decrease exponentially in time, whereas deviations along stable directions increase exponentially but decrease for large negative times. This allows to 'change connections' inside an encounter, by constructing a part of an orbit p that moves away from one part of p and towards a different part of p; its deviation from the first part has to be unstable whereas the deviation from the second part is stable.
When deriving the contribution of such pairs of orbits for single-particle systems, the deviations between encountering parts of an orbit were measured in a system of coordinates associated to the stable and unstable directions [35,36,37]. This carries over in the present scenario as well. For example, if two parts come close we consider the points where these two parts pierce through a Poincaré surface of section in our reduced phase space. By transformation from the difference between the I k , θ k (k = 2, . . . , L−1) we then obtain L − 2 pairs of coordinates s k , u k characterising the deviations between the two orbit parts in the stable and unstable directions. If the parts are almost timereversed instead of close in phase space we instead have to consider the deviation of one part from the time-reversed of the other.
We can now determine the difference between the (reduced) actions of the partner orbits. Apart from the difference associated to the energy offset already seen in the diagonal approximation this has a further contribution accounting for the change of action due to the changed connections in the encounters. For an encounter of two parts the latter contribution can be written as [35,36,37] where the sum runs over pairs of associated stable and unstable coordinates. Generalisations to encounters involving more orbit parts are given in [17]. As in the phase of (42) the action difference is divided by , systematic contributions of constructively interfering orbits will have action differences of this scale and hence very close encounters.
As a further ingredient one needs to determine the probability that encounters with given separations between the stretches arise in a periodic orbit/solution. For the corresponding formula we refer to [17]. We stress that the only requirement for its validity is the existence of spectral gap. This is used to derive an 'ergodic' probability for different parts of an orbit to come close, as well as an estimate for the duration of an encounter, both of which enter the formula.
With all ingredients for the single-particle treatment remaining valid, the contributions of all 'diagrams' of orbits differing in encounters (such as those displayed in Fig. 1) remain unchanged. These contributions involve powers of 1/ with the power increasing for more complex diagrams. The sum over all diagrams can be performed using the combinatorial techniques discussed in [17]. For systems without time-reversal invariance the contributions cancel leaving only the diagonal approximation. For timereversal invariant systems one obtains a series expansion to all orders in 1/ . In either case the results agree with the non-oscillatory terms with coefficients c n in the random matrix prediction 2i n for n ≥ 4. We note that the treatment summarised here assumes that our system (in the reduced phase space) has no further symmetries as these would lead to additional pairs of orbits with similar or identical actions, such as mutually reflected orbits in a system with reflection symmetry.
The oscillatory terms and systems with discrete geometric symmetries will be discussed later.
As in case of normal ordering the trace formula is modified to include the Solari-Kochetov (SK) phase, we have to check that this does not affect the results arising in the present approach. To do so, we note that the contributions to the double sum only involve differences between phases associated to closeby partner orbits. Due to ψ ∝ 1 √ the SK phase itself is of an order independent of . As the SK phase is time-reversal invariant the phases of two partner orbits can only differ due to the slight changes inside encounters. However for the encounters relevant for spectral statistics the encounters become very close in the semiclassical limit, meaning that the difference between the SK phases becomes negligible.

Oscillatory contributions
The random-matrix prediction (47) for R( ) also involves oscillatory contributions proportional to Re 1 (i ) n e 2i . To access these contributions a more careful semiclassical approximation is needed [38,18]. In this approximation the level density is accessed from spectral determinants via An approximation for the spectral determinant on the level of the trace formula is Here the sum is taken over sets of orbits Γ with n Γ elements and cumulative reduced action S Γ ; the amplitude F Γ depends on the stability and the Maslov indices of the contributing orbits andN (E) is the smooth approximation for the number of energy levels below E. However using the spectral determinant allows to incorporate further quantum mechanical information, in particular the fact that due toĤ being Hermitian det(E −Ĥ) has to be real for real arguments E. As shown in [38] this leads to an approximation for the spectral determinant where the contributions from sets of orbits with cumulative periods larger than half of the Heisenberg tine T H = 2π d are replaced by the complex conjugate of the contributions from sets with cumulative periods below this threshold. This 'Riemann-Siegel lookalike formula' readily generalises to the manybody systems under consideration as its key ingredient, the semiclassical approximation for the trace of the resolvent (E −Ĥ) −1 , can be accessed from the trace of the propagator as above; essentially one only has to omit the restriction to the imaginary part in (13). Again periodicity is required only in the reduced phase space.
Using the improved approximation of the level density as well as the same general ideas as outlined above it is possible to resolve oscillatory contributions as well. As each level density brings in two determinants, evaluating the two-point correlation function then requires to study interference between quadruplets of (possibly empty) sets of orbits. Hence one also needs to take into account contributions where, say, after changing connections inside encounters an orbit is broken into two orbits with similar cumulative action. A treatment of these more involved correlations shows that for chaotic quantum systems R( ) fully agrees with the predictions from RMT [18].

Systems without a spectral gap
Interestingly, there are chaotic systems that do not have a spectral gap but still satisfy some of the ingredients of our calculation, such as hyperbolicity which is required for the existence of orbit pairs differing in encounters. For these systems many of the techniques sketched here are applicable, but the final result of Wigner-Dyson statistics does not carry over as e.g. the sum rule (44) becomes invalid. An interesting question for future work is whether chaotic many-body systems without a spectral gap could be faithful to RMT ensembles that incorporate more information about the problem at hand and reduce to Wigner-Dyson statistics in important regimes, e.g. the embedded many-body ensembles [39,40] or ensembles sensitive to the spacial structure of the problem. An example for a situation where the spacial structure is important is the continuum limit with diverging number of sites; additional orbit correlations relevant in this case were identified in [41].

Symmetries
In a system with additional discrete symmetries one needs to consider the spectral statistics inside subspectra determined by the symmetry group. A prime example is the discrete (disorder-free) Bose-Hubbard model with periodic boundary conditions as introduced above. Its symmetry group is the dihedral group, consisting of the discrete translation ψ k → ψ k+1 and its iterates; the reflection ψ k → ψ L−1−k ; and combinations of translations and reflection that can be viewed here as reflection about a different centre. More formally the symmetry group for this model is the dihedral group D L . Here D n stands for the dihedral group of order 2n. The spectrum of the Bose-Hubbard model decomposes into subspectra labelled by the eigenvalues e 2πiκ/L of the discrete translation operator where κ = 0, . . . , L − 1 denotes the quasi-momentum. The spectra with κ = 0 and (if integer) κ = L/2 decompose further into components even and odd under reflection, whereas the remaining subspectra come in energy-degenerate pairs κ, L − κ related by reflection. Altogether we obtain L + 1 subspectra for odd L and L + 2 subspectra for even L. A representation-theoretic justification of this decomposition will be given below and we refer to [42] for ways to implement the decomposition numerically.
The level density associated to each subspectrum is obtained using a trace formula where the classical orbits, or in the present context the solutions of the Gross-Pitaevski equation, are required to be periodic (in the sense above) inside a fundamental domain of the system [43]. (See also [44].) For the case at hand this domain can be defined e.g. by converting ψ k to ψ k for a fixed choice of ψ 0 and then imposing certain conditions on ψ k . Demanding that Reψ k is smallest for k = 0 guarantees that applying the translation operator to a point inside the fundamental domain leads to a point outside. The same applies to most other elements of the symmetry group but in order to guarantee that reflection about the 0-th site leads outside the fundamental domain we need an additional requirement such as Reψ 1 < Reψ L−1 .
Each subspectrum α (not distinguishing between degenerate subspectra) can now be associated to an irreducible representation of the symmetry group, and the corresponding trace formula [43] has a form similar to (1), The only difference from (1) apart from the restriction to the fundamental domain is the additional factor χ α (g p ). Here g p is the group element relating the initial and final point of the orbit p if the orbit is considered in the full phase space as opposed to the fundamental domain. The character χ α (g p ) is the trace of the matrix representing g p in the representation α. As in [43] the derivation of (50) requires that a projection operator on the part of the Hilbert space associated to our subspectrum is inserted inside the trace taken over the Hilbert space; doing this in our many-particle calculations starting from (13) leads to exactly the same modifications as observed in [43] for single-particle systems.
To apply (50) to the Bose-Hubbard model we need the irreducible representations of the dihedral group [45]. The representation of the translation operator must have an L-fold power I; its two-dimensional representations are therefore be given by The reflection operator is then represented by and the representations of all other group elements can be found using that the matrix representation of a product of symmetry operations must be the product of the matrix representations. The general theory of symmetries in quantum mechanics (see e.g. [46]) implies that the energy eigenstates associated to this representation can be grouped into pairs (representable as vectors) with identical energy. One can show that for the two-dimensional representation at hand this leads precisely to the energy-degenerate subspectra associated to eigenvalues of the translation operator e 2πiκ/L (κ = 0, L 2 ) as mentioned above.
However if we have κ = 0 or, assuming even L, κ = L 2 the matrix (51) becomes diagonal and the two-dimensional representations defined above become reducible.
Instead of using these representations the associated subspectra therefore have to be described by one-dimensional representations. These represent the translation operator either by 1 or in case of even L alternatively by -1, whereas the reflection operator can be represented by either 1 or −1 regardless of L. Again these spectra precisely correspond to what we said above about the cases κ = 0, L 2 . In general the consideration of discrete symmetries can change the appropriate RMT ensemble compared to the one expected based on the time-reversal properties alone [47]. However using Eq. (50) one can show that this does not happen for subspectra associated to representations by real matrices [48,49]. As the aforementioned representations of the dihedral group are real the spectral statistics of all its subspectra is thus in line with the GOE as observed numerically in [22] as well as in the next section.

Numerical results
To support our results numerically and complement the previous numerical studies we have performed a numerical analysis of the chaotic properties of both the quantum and the classical Bose-Hubbard model. Our results support Wigner-Dyson statistics as well as mostly classical chaotic dynamics in the regime where the hopping and interaction terms are comparable.
For the quantum model we are interested in the spectral statistics as discussed above. We use a statistical observable slightly more convenient for computations than R( ), namely the normalised distribution P (r) of ratios between subsequent level spacings; if the ordered quantum levels are denoted by E n , these ratios are given by The ratio distribution is especially suited for our purposes as there is no requirement to explicitly evaluate the average level density. A random matrix prediction for P (r) was obtained in [50] by considering 3 × 3 random matrices; for the GOE it reads In particular for r → 0 one can see that P (r) ∝ r, which is due to level repulsion. Notice that for large r the distribution has a fat tail in contrast to the level spacing distribution: P (r) ∝ r −3 for r 1. As for the density of nearest-neighbour spacings the results for large matrices are expected to be very similar in value but analytically more complicated. For comparison, we mention that the ratio distribution for an integrable system (assuming that the energy levels are independent) is given by We determined the spectrum of the quantum Bose-Hubbard model via exact diagonalisation. For each of the subspectra described in section 4 we computed the histogram of r to get a numerical estimate of the ratio distribution P num (r). Examples of such numerical histograms are displayed in Fig. 2. Then we determined the L 1 norm (i.e. the integral over the absolute value) of the difference P num (r) − P RM T (r). Finally this difference was averaged over all different subspectra. Fig. 3 top shows the norm for the case L = 5 and N = 15 as well as N = 25. We obtain good agreement for the case that the interaction term (of order U N 2 ) is comparable to or slightly larger than the hopping term (of order JN ). The agreement improves as N increases, which is reassuring as our theoretical derivation of Wigner-Dyson statistics holds in the limit N → ∞ (or equivalently in the semiclassical limit). For the classical Bose-Hubbard model we introduce a numerically determined measure of ergodicity, which indicates whether the classical limit is mostly chaotic.
To determine this measure we considered the Hamiltonian where we identified ψ L = ψ 0 and ψ * L = ψ * 0 . The classical dynamics is given by If an initial point (ψ * 0 , ψ 0 ) is given in phase space, Eq. (56) allows to compute the unique trajectory (ψ * (t), ψ(t)) such that (ψ * (0), ψ(0)) = (ψ * 0 , ψ 0 ). It may be useful to add that as each point is part of a classical trajectory the point (ψ * 0 , ψ 0 ) is parameterised by 2L real numbers as the coefficients of ψ * 0 are the complex conjugates of the coefficients of ψ 0 .
As the energy E and the particle number N are conserved, the trajectory can only access a subspace in phase space, denoted by M N,E , which is the constant energy surface associated to E in the Fock space with N particles. Numerically a triangulation of M N,E is performed by picking random points on this surface. These points are generated using Halton sequences, which are commonly used to explore high dimensional spaces. First a random point is sampled on the hypersphere ||ψ|| 2 = N . If its energy coincides with E it is kept and stored as a point in M N,E . Otherwise we continue the sampling. In order to cover substantially M N,E for typical values of E, U N/J we chose to pick 2 × 10 6 random points. Then the mean distance between any two of these points is calculated, let us call it δ. This gives a typical distance scale for the constant energy surface at a given energy E. This typical distance is used to define balls of radius δ/10 around each points. After having determined a cover of M N,E the trajectory starting from (ψ * 0 , ψ 0 ) is computed for a sequence of increasing (final) times. We took large enough final values of the propagation time so that the number of visited balls does not vary significantly after this time.
Finally a measure for the ergodicity of a trajectory is defined by the ratio between the number of visited balls and the total number of balls. In order to have generic values, this ergodicity measure is averaged over different choices of initial conditions (ψ * 0 , ψ 0 ). Our results are displayed in Fig. 3 bottom. By definition our measure is a real number between 0 and 1. 1 indicates perfect ergodicity. The closer it is to 1 the more ergodic the trajectory is. While not a rigorous check of the conditions for random matrix statistics, our measure indicates whether the system is mostly ergodic or has a substantially mixed phase space with large stability islands. Fig. 3 bottom shows that, for the same range of U N/J, the classical dynamics is predominantly ergodic and the quantum system agrees closely with RMT prediction.
Our results are compatible with those obtained through a different approach and a modified version of the model in [51]. There the authors chose a different normalisation of ψ and claimed that their version of the classical Bose-Hubbard model becomes more and more classical when increasing U/J. This is line with the left part of Fig. 3 bottom.
we choose δH(θ k ) periodic in θ k with the period coinciding with the range of values covered by θ k .
We first investigate how the perturbation modifies the vanishing eigenvalue of the Hessian H. The corresponding eigenvector is associated to a simultaneous shift of θ jk for our k at all time steps j. It is thus convenient to rewrite the Hessian in terms of derivatives w.r.t. I jk , θ jk ; this also absorbs the divisor in (21) and it has the advantage that the Hessian becomes Hermitian. The perturbation of the Hessian δH has as its only non-vanishing matrix elements the derivatives − ∂ 2 δH(θ jk ) ∂θ 2 jk τ . In the same coordinate system the eigenvector e associated to the vanishing eigenvalue has identical entries for the J components associated to the θ jk associated to our k and the remaining components are zero. First-order perturbation theory then yields a perturbed eigenvalue e · δHe e · e = − 1 J To study the effect on the stability matrix we first consider the matrixM (k) j mapping a deviation of (θ jk , I jk ) to the resulting deviation after a time step of size τ . We obtain, e.g. by converting Eq. (27)  Here the lower left entry vanishes for the unperturbed Hamiltonian but is moved away from zero by the perturbation δH(θ k ). As a consequence of this perturbation the product given in (32)

Appendix B. Weyl ordering
Finally we want to show explicitly that in line with [27,24,25,9,33] the SK phase does not arise in case of Weyl ordering. Using the Hamiltonian in Weyl ordering the short-time propagator can be written as [33] ψ j |e − i Ĥ τ |ψ j−1 ≈ 2 L d[w j , w * j ] exp − 2|w j | 2 + 2ψ * j · w j + 2ψ j−1 · w * j −|ψ j | 2 /2 − |ψ j−1 | 2 /2 − ψ * j · ψ j−1 − i H(w * j , w j )τ (B.1) involving an integration over a pair of complex conjugate variables w j , w * j that were not required in normal ordering. The action in the discretised coherent-state path integral then turns into [33] and the path integral runs over w j , w * j as well as ψ j , ψ * j . The action becomes stationary under variation of ψ j , ψ * j if w j = ψ j + ψ j−1 2 (B.3) and its complex conjugate hold. This suggests that in the limit of fine discretisation the w j should be interpreted as macroscopic wave functions halfway between two discretisation steps. Furthermore stationarity w.r.t. variation of w j , w * j implies i w j − ψ j−1 τ /2 = ∂H(w * j , w j ) ∂w * j (B.4) and its complex conjugate, which is the appropriate discretisation of the nonlinear Schrödinger equation i ψ = ∂H ∂ψ * for the time interval τ /2. After eliminating w j with the help of (B.3) a short calculation shows that the combination of macroscopic wave functions in the action (B.2) agrees with the simpler form in (12).
We now have to evaluate the Hessian of R and its determinant. In analogy to subsection 2.4 we start with the Hessian containing all second derivatives. If we first write the derivatives w.r.t. all w * j , w j and then those w.r.t. all ψ * j , ψ j the Hessian assumes the form i D −2E −2E T E + E T where D is block-diagonal with the blocks (whereM j is defined as in (27) but with w j , w * j taking the role of ψ j , ψ * j ) and again ignoring terms of quadratic and higher order in τ in the entries we obtain As usual the stationary-phase approximation requires the inverse square root of this determinant, containing a factor 2 −J L . This factor is cancelled by the 2 J L in the integration measure for w j , w * j . As anticipated, we thus obtain the same result as with normal ordering apart from the SK phase. The remaining steps, in particular taking into account the conservation laws, carry over directly from the case of normal ordering discussed in the main text.