Quantum chaos and ensemble inequivalence of quantum long-range Ising chains

We use large-scale exact diagonalization to study the quantum Ising chain in a transverse field with long-range power-law interactions decaying with exponent $\alpha$. We numerically study various probes for quantum chaos and eigenstate thermalization {on} the level of eigenvalues and eigenstates. The level-spacing statistics yields a clear sign towards a Wigner-Dyson distribution and therefore towards quantum chaos across all values of $\alpha>0$. Yet, for $\alpha<1$ we find that the microcanonical entropy is nonconvex. This is due to the fact that the spectrum is organized in energetically separated multiplets for $\alpha<1$. While quantum chaotic behaviour develops within the individual multiplets, many multiplets don't overlap and don't mix with each other, as we analytically and numerically argue. Our findings suggest that a small fraction of the multiplets could persist at low energies for $\alpha\ll 1$ even for large $N$, giving rise to ensemble inequivalence.


I. INTRODUCTION
Thermalization in classical Hamiltonian systems is well understood in terms of chaotic dynamics and the related essentially ergodic exploration of the phase space [1][2][3]. From the quantum point of view the physical mechanism is quite different, with the eigenstates of the Hamiltonian behaving similar to the eigenstates of a random matrix with the additional property that they appear thermal from the point of view of local measurements. This is the paradigm of eigenstate thermalization (ETH) introduced in Refs. [4][5][6][7]. In general there is correspondence between classical and quantum thermalization [5,[8][9][10][11][12][13], but due to the different physical mechanism there can be cases where quantization breaks ergodicity, as for manybody localization (see [15] for a review) and many-body dynamical localization [7,[16][17][18][19]21].
In quantum short-range thermalizing systems there are three strictly related properties. First of all eigenstate thermalization, that's to say that almost all the excited eigenstates locally behave equal to the microcanonical or thermal density matrix [22]. So, expectation values of local observables equal the corresponding microcanonical ones, up to fluctuations vanishing in the thermodynamic limit. This property is strictly related to a second one: quantum chaos [22]. Quantum chaos means that the spectrum of the Hamiltonian behaves essentially as the one of a random matrix [88] and this occurs typically for many-body nonintegrable models [23] and for Hamiltonians obtained quantizing classical chaotic systems [12]. Hamiltonians show in general eigenstate thermalization together with quantum chaos and behave as random matrices [22] (with some caveats [14]). This fact gives rise to random eigenstates which look locally thermal as appropriate for ETH. A third property relevant in thermalized short-range interacting systems is additivity and ensemble equivalence which are strictly related to a convex microcanonical entropy [37].
An interesting question is if what is the relation between quantum chaos, ETH and ensemble equivalence in quantum systems with long-range interactions. In the classical case, for instance, the thermalization behavior is very different in the case of short-and long-range interactions. For classical systems with short-range interactions, any nonlinear Hamiltonian with more than two degrees of freedom and no conservation law beyond energy gives rise to chaos, essentially ergodic dynamics [1] and ensemble equivalence [2]. In the long-range case the situation is very different. A central aspect of long-range classical systems is the inequivalence of canonical and microcanonical ensemble due to the lack of additivity of the Hamiltonian [24,36,37]. This implies that the dynamics does not lead to a simple thermalization behavior, even in presence of chaos. One can see an effectively regular behavior dominated by one or few degrees of freedom [24][25][26][27][28][29] which has been exploited to obtain a classical Hamiltonian time crystal [30].
Although ensemble inequivalence for the exactlysolvable infinite-range anisotropic quantum Heisenberg model has been studied in [31][32][33], the relation between quantum chaos and ensemble equivalence in generic interacting quantum long-range systems has not yet been explored. We fill here this gap focusing on a long-range ferromagnetic Ising spin-1/2 chain model. Similar models have been already studied. One very well studied case is the Ising model with infinite-range interactions (the so called Lipkin-Meshkov-Glick model) which is known to be integrable [34,35,39]. It is also known that the isotropic Heisenberg chain with power-law interactions with exponent α = 2 is integrable [40,41] as well as some anisotropic spin-chain models with α = 2 [42][43][44]. Spin chains with disorder and power-law interactions are known to undergo a transition between a many-bodylocalized-like and an ergodic phase [45][46][47][48][49][50][51].
In our work we focus on this same model and widely extend the ETH and quantum chaos analysis by using exact diagonalization and exploring a wide range of α and energies. The main question is the relation between eigenstate thermalization, quantum chaos and convex microcanonical entropy. For α < 1 we find a very interesting behavior.
On the one hand the level spacing statistics gives a clear answer pointing towards a random-matrix Wigner-Dyson form. This is valid for any value of 0 < α < ∞, but for the region around α ≈ 2 for weak transverse fields around, hinting to possible the vicinity of some integrable point.
On the other hand ETH indicators (eigenstate expectations and eigenstate half-system entanglement entropies) yield a much less clear perspective for finite system sizes, in particular, for α < 1. We find that the permutation symmetry, which is only exact at α = 0, leaves behind a strong fingerprint in many ETH indicators: The α = 0 symmetry-protected multiplets in the energy spectrum represent a relatively rigid structure for 0 < α < 1. They affect the eigenstate quantities and forbid them a smooth ETH dependence on energy, in contrast with short-range interacting quantum-chaotic systems [22].
These multiplets have another important consequence: The microcanonical entropy becomes a nonconvex function of energy, which in the thermodynamic limit excludes ensemble equivalence in a thermodynamic sense. We provide an analytical argument for the rigidity of the multiplets for large but finite N when α < 1. For α 1 we observe that some of the multiplets at low energies persist also for large N . As a consequence, we argue that the system doesn't obey ensemble equivalence.
These observations on the multiplet structure seem to contradict our findings for the level spacing statistics. These results are reconciled by what we call a partial spectral quantum chaos. The states in individual multiplets, which are separated in energy with respect to each other, mix in a quantum chaotic fashion, whereas the multiplets don't yet mix among each other for the accessible system sizes. Each multiplet in the bulk of the spectrum behaves as a separate random matrix leading to a overall Wigner-Dyson level statistics. This is a significant result: each multiplet behaves a random matrix from a spectral point of view, so its spectrum tends to a smooth continuum for N → ∞. This is in contrast to integrable long-range system whose spectrum has been claimed to be pure point also in the thermodynamic limit [62].
We emphasize again that we expect the multiplet structure to be most rigid at low energy densities, which might have important consequences for the absence of thermalization observed in low-energy quenches [53,61].
The paper is organized as follows. In Sec. II we define the model Hamiltonian. In Sec. III we study the quantum chaos properties at the level of the spectrum. We show a generalized tendency towards a Wigner-Dyson levelspacing statistics for increasing system size. In Sec. IV we discuss an analytical argument based on the randommatrix behavior of each multiplet. We show that the spectral multiplet width increases linearly in α, in agreement with numerics, and that part of the multiplets persist in the large-N limit, for low energies and α 1. In Sec. V we better discuss the multiplet spectral structure for small α and finite N and study the corresponding nonconvex behavior of the microcanonical entropy related to ensemble inequivalence. In Sec. VI we study the broken symmetry edge (the energy density below which there is Z 2 symmetry breaking) and find a different behavior in the canonical and microcanonical ensemble, although there are too strong finite-size effects to allow to make statements on ensemble inequivalence. We study also the eigenstate properties by considering the eigenstate expectation values of a local operator, the longitudinal nearest-neighbour correlation (Sec. VII), and of the halfsystem entanglement entropies of the eigenstates (Appendix A).
In Appendix B we discuss the Hilbert-Schmidt distance of the α > 0 Hamiltonian from the α = 0 Hamiltonian, showing its linearity in the limit α → 0. This fact, together with the random-matrix assumption, allows us to explain the linearity in α of the multiplet spectral width in Sec. IV.

II. MODEL HAMILTONIAN
In this work we study the ferromagnetic long-range interacting quantum Ising chain in a transverse field: Here, σ α i with α = x, y, z denotes the Pauli matrices at lattice site i = 1, . . . , N with N the system size. We use periodic boundary conditions implemented through the definition [57] in order to preserve extensivity of the Hamiltonian.
We use exact diagonalization. We largely exploit the translation, inversion and Z 2 (σ z i → −σ z i ) symmetries of the model in order to restrict to an invariant subspace of the Hamiltonian. In most of the text we restrict to the subspace fully symmetric under all the symmetries of the Hamiltonian. We call this Hamiltonian eigenspace H S and we define it as the zero-momentum sector subspace even with respect to inversion and Z 2 symmetry. For future convenience we define N S ≡ dim H S . In Sec. VI we are interested in the spectral pairing properties of the model, which requires to consider both Z 2 symmetry sectors: We consider here the zero-momentum sector subspace even only with respect to inversion. We denote the eigenstates of the Hamiltonian |ϕ µ and the corresponding eigenenergies E µ (taken in increasing order), while always specifying which subspace we are considering.
In the limit α → ∞ the model in Eq. (1) reduces to the nearest-neighbour quantum Ising chain. This model is integrable and undergoes a quantum phase transition: Its ground state breaks the Z 2 symmetry for h < 1 [78,79]. For any finite system size, the ground state is doubly degenerate made up by the two states symmetric and antisymmetric under the global Z 2 symmetry, with a splitting exponentially small in the system size. The states in the doublet show long-range order and the doublet becomes degenerate in the thermodynamic limit, giving rise to symmetry breaking.
In the limit α = 0, on the opposite, Eq. (1) reduces to the Lipkin-Meshkov-Glick model. This model is also integrable, thanks to the full permutation symmetry, and it shows a symmetry-broken phase for h < 1. In contrast to the α → ∞ case, all the spectrum up to an extensive energy N e * is organized in doublets with exponentially small splitting and the corresponding eigenstates have long range order [34,35,39]. Due to the full permutation symmetry, the Hilbert space is factorized in a number of invariant subspaces, differently transforming under the permutation symmetries [35]. The number of these subspaces is exponential in N , and many of them have the same level structure. This gives rise to massively degenerate multiplets, whose levels belong to different symmetry sectors, a property which will be quite relevant in the following.
For α = 0, the number of distinct multiplets is set by the possible distinct simultaneous eigenstates of the square total spinŜ 2 = 1 4 ( jσ ) 2 and the total spin z componentŜ z = 1 2 j σ z j . This is a consequence of the total-spin conservation and the permutation symmetry of the Hamiltonian [34]. The square total spin has eigenvalues S(S + 1) with S going from S = 0 to S = N/2 and for each value of S the total z component can acquire 2S + 1 values. Assuming N from now on even -so that S assumes only integer values -the number of multiplets is Q = N/2 S=0 (2S + 1) = (N/2 + 1) 2 . For α = 0 each multiplet is degenerate with degeneracy g(S) given only by S and N through the formula [35] g (S) In the remainder of the paper we consider the case of intermediate α.

III. QUANTUM CHAOS AND LEVEL SPACING STATISTICS
First, we study the quantum chaos properties focusing on the level spacing statistics. The model in Eq. (1) is integrable for the limits α = 0 (infinite-range case) and α → ∞ (nearest-neighbour case). We now aim at exploring the behavior at intermediate α.
For concreteness, we don't scan extensively across the transverse fields, but rather focus on two representative values h = 0.1 and h = 0.5. In Fig. 1 we investigate the spectral properties of the model as a function of α upon varying the system size N . Specifically, we plot the average level spacing ratio, r (introduced in [81]), which is a central probe for quantum chaos and is defined as With the time-reversal symmetry properties of our Hamiltonian, a value r = r WD 0.5295 would be associated with a fully quantum-chaotic random-matrixlike behavior given by the Gaussian Orthogonal Ensemble (GOE) and a Wigner-Dyson distribution for the level spacings [88]. On the opposite, a value r = r P 0.386 is known to be related to a Poisson distribution of the level spacings, which implies integrable behavior [93].
Before considering the behavior for large α (Sec. III A) and α 1 (Sec. III B), and the associated tendency towards quantum chaos for increasing N , let us say something about the strong minimum at α = 2 appearing in Fig. 1(a). It suggests a behavior closer to integrability (and the corresponding Poisson value) which persists at least up to N = 22. It is important to remind that there are spin models with power-law interactions decaying with α = 2 that are integrable, such as the longrange isotropic Heisenberg chain [40] or other anisotropic long-range models [42][43][44]. It could be an interesting question for future research to investigate whether this phenomenon is related to the proximity to an integrable point.

A. Large α
For large α we see in Fig. 1 that there is a crossover towards the Poisson value r P . At some larger value of α there is another crossover towards a value even smaller than Poisson. This behavior of r is a finite-size effect due to the proximity of the integrable α → ∞ point. The spectrum becomes quantum chaotic in the thermodynamic limit: As we are going to show, the crossover towards Poisson shifts to large α for increasing N .
We can argue this shift towards integrability as follows. In a free-fermion model (corresponding to our α → ∞ case), any arbitrarily small integrability-breaking next-nearest-neighbour interaction restores thermalization in the thermodynamic limit [75,76]. Similarly, in our case, for α 1, the next-nearest neighbour terms are the stronger ones breaking the integrability of the nearest-neighbour α → ∞ model. For increasing N , the next-nearest-neighbour terms become at some point large enough compared to the level spacings, and the model becomes quantum chaotic Let us now roughly estimate the crossover scale at which the system becomes quantum chaotic for α 1, by comparing the next-nearest neighbour interaction term with the relevant gap ∆ of the integrable nearest-neighbour model. The next-nearest neighbour term is of order V ∼ J/(N (α)2 α ). We can understand the relevant gap of the nearest-neighbour model, moving to its fermionic representation via the Jordan-Wigner transformation [94]. In this representation, the nearest-neighbour model is integrable and its excitations are fermionic quasiparticles [79,95] We have k ∈ [0, π] and, for finite system size N , k can take only N discrete equally spaced values. In the fermionic representation the next-nearest-neighbour term becomes a fourfermion term which induces inelastic scattering between the fermionic quasiparticles. If momenta k 1 and k 2 We can roughly estimate ∆ by taking twice the bandwidth of k and dividing it by N , the number of allowed equally-spaced k values. We find Imposing that V ∆, one finds that quantum chaotic behavior is obeyed for α α * . We evaluate α * numerically, and find that α * asymptotically increases as log 2 N (see Fig. 2). So, for N → ∞ there is quantum chaos for all values of h.

B. The role of multiplets for α 1
For α 1 r is close to the Wigner-Dyson value ( Fig. 1). Therefore, our numerics suggests that the integrable behavior at α = 0 [34] is unstable to a small perturbation in α which breaks the full permutation symmetry at α = 0.
As we have already discussed in Sec. II, the multiplets at α = 0 do not correspond to a given permutation symmetry class, but contain states belonging to different invariant subspaces, differently transforming under permutation. There are many subspaces with the same energy levels inside [35]. When perturbation symmetry is broken by α 1, the degenerate states inside each multiplet can mix and so all the subspaces are mixed by the Hamiltonian. This leads to quantum chaos, as we are going to argue.
Since there is no gap to protect the subspaces from mixing, this change happens abruptly as soon as α > 0 and the multiplet degeneracy is lifted. We can see an example of that in Fig. 3. We plot E µ versus µ/N S for h = 0.1 and two values of α, α = 0 and α = 0.15. For α = 0 there are many degenerate multiplets at all energies, as we can see in the magnifying insets. For α = 0.15 the multiplets merge into a smooth continuum at large energy (right inset) but can be still well identified at low energy (left inset). The organization of the spectrum in multiplets for small α is also evident in the eigenstate expectation of local observables (Sec. VII) and the halfsystem entanglement entropy of these eigenstates (Appendix A).
This multiplet structure is apparently in contrast with the average level spacing ratio being close to the Wigner-Dyson value. In order to explain this apparent contradiction, we notice that the number of gaps among multiplets is much smaller than the total number of states. The number of discontinuity points scales as the number of distinct multiplets at α = 0, which scales as N (N + 1)/2 (see Sec. V), while the number of states equals N S which is exponential in N . So, if each of the multiplets behaves separately as a random matrix, the overall average level spacing ratio is Wigner Dyson in the large N limit. This is exactly what happens, as we show in detail in the next section.

IV. RANDOM-MATRIX BEHAVIOR AND MULTIPLET SPECTRAL WIDTH FOR α < 1
Goal of this section is to argue that each multiplet broadens by an amount proportional to α. This numerically verified statement relies on the Hamiltonian projected to a multiplet subspace behaving like a GOE random matrix, as we argue in Sec. IV A. The main implication is that the total multiplet width is linear in N and much smaller than the total spectral width for α 1. This result has important consequences for the rigidity of part of the multiplet structure in the large-N limit, as we clarify in Sec. IV B.
A. Width of a single multiplet Let us focus on ∆Ĥ(α, N ) =Ĥ (α) −Ĥ (0) , the difference of the two Hamiltonians at α and at α = 0. We choose the basis |i of eigenstates ofĤ (0) such that the matrix elements H S denoting the energy of the multiplet with spin S at α = 0. Then we consider the square root of the quadratic average of the matrix elements of ∆Ĥ(α, N ), defined in the following way N in the denominator is the number of nonvanishing matrix elements of ∆Ĥ(α, N ). In order to quantify it we recall that ∆Ĥ(α, N ) is a sum of terms of the form σ z j σ z l . Under a global rotation, σ z j σ z l transforms like the sum of a scalar and a tensor, i.e. an object with spin 2. Thus, by Wigner-Eckart theorem [112], and by the rules of spin addition, we have that, if |S, i is a state with spin S, then σ z j σ z l |S, i is a superposition of states whose spin is in the set {S − 2, S − 1, S, S + 1, S + 2}. Considering that in each spin-S sector there are 2S + 1 multiplets, and that ∆Ĥ(α, N ) commutes with the total spin along z, we can therefore evaluate N as The numerator in Eq. (5) is the Hilbert-Schmidt norm of ∆Ĥ(α, N ), whose symbol is ∆H(α, N ) HS . As we show in Appendix B, the scaling behavior of this norm is where K > 0 is a numerical factor. We emphasize that K is order 1 for the values of α < 1 we are considering (see Appendix B). dim H= 2 N is the dimension of the full Hilbert space. (Restricting to the fully even subspace will only modify dim H and g(S) by a factor 1/N , leaving Eq. (7) and our conclusions unchanged.) We assume now that: (i) the gaps separating each multiplet from the neighbouring ones are much larger than the matrix elements coupling it to them; (ii) when we restrict to a multiplet, the spectrum resembles that of a random matrix from the GOE ensemble. We might expect the second assumption to hold on the one hand due to our results on quantum chaos and on the other hand since the projection onto a multiplet is an highly nonlocal operation that will destroy any locality -or sparsitystructure from H (α) . When these assumptions hold, the eigenvalue spectrum in each multiplet resembles Wigner's semicircle law [88,89], and the multiplet spectral width is given by with the multiplet-degeneracy g(S) given in Eq. (2), and N in Eq. (6). We emphasize that averaging the square matrix elements over all the Hilbert space does not contradict the fact that each multiplet separately behaves as a random matrix, as long as assumption (i) is valid and there is no mixing between multiplets. Eq. (7) tells us that our assumption of random-matrix behavior inside a multiplet gives rise to the prediction of a w(N, S) linear in α. We can numerically verify that this is exactly what happens for multiplets in the bulk of the spectrum (see Fig. 4). So, each multiplet separately behaves as a random matrix and all together give rise to the Wigner-Dyson statistics. Near the edges of the spectrum the behavior is probably different, but states near the spectral edges are a small fraction, vanishing in the limit of large N .

B. Total multiplet width and spectral rigidity
In order to better understand the rigidity of the multiplets upon increasing system size N , we now consider the total multiplet width [90] We evaluate this quantity using Eqs. (7) [97] and (2) and numerically compute the factorials using the Lanczos formula [98]. We see that W (N ) increases linearly in N [see inset of Fig. 5(a)] with a slope obtained from a linear fit β W = 0.9.
In order to understand if the majority of the multiplets overlaps for large N , or if there is a significant fraction of them which survives, we need to compare W (N ) with the total spectral width ∆E(N ) ≡ max µ (E µ ) − min µ (E µ ), So, when α < β∆ Kβ W , the total multiplet width W (N ) is asymptotically smaller than the total spectral width ∆E(N ). In particular, when α 1 [more precisely, α min(1, β∆ Kβ W )], we expect that the spectral structure seen in Fig. 3 persists for larger system size, with a multiplet structure visible at low energy densities. When α 1 we have W (N ) ∆E(N ) for large N and we expect that some multiplets persist.
Looking at Fig. 3 (see also Figs. 6 and 7), we see that the persisting multiplets lie at low energy densities. The rigidity of these multiplets, and the related ensemble inequivalence, are likely behind the effective nonergodic behavior and the persistent longitudinal magnetization appearing in low-energy quenches [53,61] for α < 2.

V. NONCONVEX MICROCANONICAL ENTROPY AND ENSEMBLE INEQUIVALENCE
The spectrum being organized in multiplets gives rise to a nonconvex microcanonical entropy, with many maxima, one per each multiplet. As we have seen above, for α 1, part of the multiplets persists for very large N . A nonconvex microcanonical entropy in this limit gives rise to ensemble inequivalence, as it happens in classical long-range systems [37].
In order to visually show how the presence of multiplets gives rise to a nonconvex microcanonical entropy, let us numerically evaluate the microcanonical entropy S th (E) in a case of finite N . To define the entropy, we start from the density of states We average it over an energy shell (we divide the energy spectrum in N Shell equal energy shells and mark the energy-shell average as · · · Shell ) and we define S th (E) = ln ρ Shell (E) (for each shell, E is the middle-point energy and we take k B = 1). We show our results in Figs. 6 and 7.
In Fig. 6 (a) we plot S th (E) versus the energy density E/N for α = 0.05, h = 0.1 and two system sizes. At low and intermediate energy densities, we clearly see the peaks corresponding each to a multiplet and we do not see a strong tendency for them to disappear for increasing system size. We can see something similar for α = 0.25, h = 0.1 [ Fig. 6 (b)] where the low and intermediate energy density multiplet structure becomes more evident for increasing system size. So, multiplets strongly affect the dynamics for finite system sizes giving rise to a nonconvex microcanonical entropy. For α < 1 we clearly see the same nonconvex structure for both h = 0.1 and h = 0.5 (Fig. 7). We remark that each peak corresponds to a multiplet, an object with many levels giving rise to a smooth random-matrix continuum for N → ∞. So each peak is something physical, very different from the spikes appearing at finite size in the density of states of the short-range Ising model, when a energy shell smaller than the finite-size gaps between the eigenenergies is considered.
In the plots in Fig. 6 we notice that at the lowest energy densities we have only few levels in the multiplets and there are significant gaps separating the multiplets. The first two or three multiplets survive even at larger α, as we can see in the density-of-states plots of

VI. SPECTRAL PAIRING AND BROKEN SYMMETRY EDGE
It is well known that the long-range quantum Ising chain exhibits a symmetry-breaking transition at nonzero temperature as soon as α < 2. [106] The corresponding microcanonical or even single-eigenstate properties have, however, not been explored extensively, except the notable Ref. [74]. Here we study the long-range order of the eigenstates which gives rise to Z 2 symmetry breaking in the thermodynamic limit. In particular, we want to quantify whether for α = 0 there are states with longrange order at finite excitation energy density and to estimate the critical energy density e * below which the eigenstates break the symmetry in the thermodynamic limit (e * is called broken symmetry edge [35]). The existence of the broken-symmetry edge is well known for the case α = 0 [35], h < 1, but it is not explored in detail for α = 0. We are going to compare this quantity with the corresponding canonical one and show that the two differ from each other for the accessible α ≤ 1.5 values.
For the microcanonical analysis, we need both the two Z 2 symmetry sectors. Therefore, we restrict to the sub- space corresponding to the zero-momentum sector and even only with respect to inversion. We target the single eigenstates and study the energy gaps between nearby states: If there is symmetry breaking in the thermodynamic limit, the eigenstates must appear in quasidegenerate doublets, which become degenerate in the thermodynamic limit (the splitting is exponentially small in the system size). We make use of this property to determine the broken-symmetry edge. We consider the splitting inside pairs of nearby eigenenergies ∆ (1) n = E 2n −E 2n−1 , (n is an integer number labeling the eigenvalues in increasing order) and the gap between nearby pairs, evaluated as the difference of next-nearest neighbor eigenenergies ∆ (2) n = E 2n+1 − E 2n−1 . If we are in presence a quasidegenerate doublet (E 2n−1 and E 2n belong to the same doublet), ∆ (1) n should be much smaller than ∆ (2) n and the ratio ∆ (1) n /∆ (2) n should scale to 0 with the system size. It is convenient to average ∆ (1) n and ∆ (2) n on energy shells, in order to reduce fluctuations. We define the N Shell energy shells as in Sec. V and we consider the ratio n Shell (E) ∆ (2) n Shell (E) (10) of the averages over the energy shells ∆ (1) n Shell (E) and ∆ (2) n Shell (E). We term D(E) as the relative splitting and plot it versus E/N for different system sizes in Fig. 8. We consider h = 0.1 and two values of α, α = 0.05 [ Fig. 8. (a)] and α = 0.5 [ Fig. 8. (b)]. For the first value of α the spectrum is organized in multiplets for the system sizes we have access to, while for the second it does not. For α = 0.5 we can see that the curves for different N clearly cross: There is a value of E/N below which D(E) decreases with the system size and above which increases. This is exactly what one would expect for a broken-symmetry edge, and we take this crossing point as an estimate for the broken symmetry edge, with an errorbar given by the mesh in E.
In contrast to the α = 0.5 case, for α = 0.05 we do not see any crossing as smooth as this one [ Fig. 8 (a)]. For this value of α and these system sizes, the dynamics is strongly affected by the above-discussed multiplets. A noisy behavior appears in Fig. 8 (a) and does not allow us to clearly give an estimate for e * . We will estimate the broken symmetry edge only for those values of α and N where we do not see a noisy multiplet structure in the crossing region.
We plot the resulting microcanonical e * versus α in Fig. 9 for h = 0.1 and h = 0.5 with the label "Micro". We obtain it considering the crossing of the relative-splitting curves for N = 20 and N = 22 and for α = 0 we take the theoretical value e * = −h found in [35]. We can reliably estimate e * with our method up to α = 1.5. Above that value larger system sizes are needed.
We compare it with the canonical broken-symmetry edge labeled as "Canonical" in Fig. 9. The latter is evaluated considering the Binder cumulant, a measure of Z 2 symmetry breaking particularly effective in the canonical ensemble [38]. DefiningŜ q z ≡ jσ z j q , the Binder cumulant is given by , where · · · th is the thermal canonical expectation. Varying the temperature, both B and the corresponding energy density e = Ĥ th /N vary. We plot B versus e for a set of parameters and two different values of N in Fig. 10. The canonical symmetry breaking threshold is estimated as the crossing between these two curves, in a way similar to what done in [61]. Here the thermal canonical expectations · · · th are obtained by evolving in imaginary time a purified infinite-temperature state [108,109]. The imaginary-time evolution is performed through the TDVP algorithm [110,111].
The canonical e * versus α (Fig. 9) shows a strong dependence on N and so that the canonical e * increases if we take the crossing of curves for larger N : The difference with the microcanonical value increases. This fact suggests ensemble inequivalence, but finite-size effects are too strong for making a precise statement. Moreover, considering that the ground-state is at e GS −1), Fig. 8 gives us the nontrivial conclusion that for α ≤ 1.5 the system shows Z 2 symmetry breaking at finite excitation energy densities. So, there is a finite fraction of the energy-spectrum width where the eigenstates show long-range order, similarly to the α = 0 and the disordered case. This is in agreement with the findings of [53,61], where the long-time dynamics supports longtime magnetization in the range α ≤ 1.5 and beyond.

VII. ETH PROPERTIES
After having studied in detail spectral properties, we now take a step further and aim to study eigenstatethermalization properties. For concreteness, we consider the longitudinal nearest-neighbour correlation operator as a representative for local observables. We focus on the properties of the eigenstate expectation values G µ ≡ ϕ µ |Ĝ|ϕ µ . We expect that the same behaviour occurs for any local observable. As we show in Appendix A also the entanglement entropy (involving half of the system size) shows a similar behaviour. We consider the scatter plots of G µ versus E µ in Fig. 11. Most importantly, these expectation values as a function of energy don't always exhibit a smooth dependence with small fluctuations, as expected in a system obeying ETH [6] even though the level spacing ratio Eq. (3) is close to Wigner-Dyson. The finite-size effects are too strong, mainly related to the spectrum being organized in multiplets for α < 1, and no quantitative extrapolation to larger size is possible. Nevertheless we see a lack of correspondence between quantum chaos and ETH, in contrast with short-range interacting systems.
The most noteworthy case is α = 0.05 [ Fig. 11 (a) and (b)] where we see many almost vertical lines, as many as the multiplets. Each of these lines is a continuous curve, as if ETH was to hold just within a multiplet but not across them. As we have argued in Sec. IV, when N is increased, part of the multiplets should survive, and then this behaviour should persist. What we see in Fig. 11 (a) and (b) is nevertheless strongly affected by finite size effects.
Another interesting case is provided by α = 0.5 [ Fig. 11 (c) and (d)]. For h = 0.1 [panel (c)] we can see a qualitatively different behavior at large and small energy. In the center of the spectrum we observe a quite smooth curve with some small fluctuations, which appears as a prototypical example of a system obeying ETH. Overall, for these small system sizes, this doesn't seem to follow the predictions by ETH.
For larger α [α = 1.5 in Fig. 11 (e), (f) and α = 2 in Fig. 11 (g), (h)] we see a fully developed ETH behavior for h = 0.5: very smooth curves with noise at the edges of the spectrum [panels (f) and (h)]. On the opposite, for h = 0.1 [panels (e) and (g)], the situation is not at all ETH, in close correspondence with the average level spacing ratio being different from Wigner-Dyson [ Fig. 1 (a)]. In particular, the case α = 2 is very regular-like with some scattered points between the horizontal lines suggesting a stronger mixing at larger system sizes.

VIII. CONCLUSION
In conclusion we have considered the long-range Ising model with power-law interactions and used exact diagonalization to study the relation between quantum chaos, eigenstate thermalization and convexity of the microcanonical entropy. For small α we have remarkably found that the level spacing distribution is Wigner Dyson but this does not reflect a full-random-matrix behaving Hamiltonian.
The reason comes from the strong effect of the α = 0 integrable point, where the Hilbert space decomposes into many identical subspaces with the same energy levels, due to the full permutation symmetry. Even an in- finitesimal α > 0 mixes the degenerate levels belonging to different subspaces; the resulting spectrum is organized in multiplets and we argue that multiplets in the bulk of the spectrum separately behave as a random matrices, with a negligible role of the spectral edges. Due to the strong effect of multiplets, this Wigner-Dyson spectral statistics appears in association with anomalous thermalization properties.
The randommatrix behavior of the multiplets suggests that part of the multiplets persists at large N and α < 1. This holds in particular in the α 1 limit. So, also at large N there are multiplets, and they give rise to a nonconvex microcanonical entropy as a function of energy, implying ensemble inequivalence [37]. From the numerics, we expect that the multiplets persisting at large N lie at low energy densities; they are probably involved in the persistent magnetization, which has been observed in the low-energy dynamics of this model [53,61].
We further analyse the eigenstate thermalization properties and we see that at small α the local observable eigenstate expectation values and the corresponding half-system entanglement entropies do not organize into smooth curves as a function of the energy, as one should naively expect from quantum chaotic behavior in the Wigner-Dyson level spacing statistics. In contrast to short-range interacting systems the spectrum is organized in multiplets and there is no simple ETH behavior. Quantitative probes (see Appendix A) suggest that the curves become smoother for increasing system sizes and we cannot tell if this is due to the ETH being obeyed better and better inside the multiples or to the fact that the multiplets at large energy densities tend to merge.
We remark that our exact diagonalization results show a persisting nonergodic behavior for h = 0.1 and α around the value α ≈ 2. This is a suggestive result because there are other long-range models with α = 2 which are integrable, but the system sizes we have access to do not allow to state if this effect persists in the thermodynamic limit. Nevertheless, a nonchaotic behavior for N = 22 is already remarkable and might suggest at least the proximity of an integrable point. In all the other cases we see an ergodic behavior.
Perspectives of future work will focus on the connection between the dynamical phase transition in α undergone by this model [53,56,61] and the corresponding lowenergy confinement-deconfinement transition [57]. Another direction of research will be to study the relation between quantum chaos in sectors of the Hilbert space and ensemble inequivalence in models with Hilbert space fragmentation [113]. ETH properties of eigenstates can be explored also by means of the entanglement entropy. This is not a local object because it involves correlations extending up to a distance N/2, but eigenstate thermalization has been proved valid for subsystems up to this size [80]. Considering an eigenstate |ϕ µ , and decomposing the system in two parts A and B in physical real space, we define (A1) Specifically, we focus on the half-system entanglement entropy S In Fig. 12 we show the scatter plots of the entanglement entropy S (µ) N/2 [defined in Eq. (A1)] versus the corresponding eigenstate energy E µ . ETH is strictly related to these curves looking "smooth", as appropriate for microcanonical entropy [80]. Let us first discuss this point qualitatively. We consider a small value of α, α = 0.05 [panels (a), (c)]. The S (µ) N/2 versus E µ look like smooth curves, as in the ETH case, only if we restrict inside the multiplets. This result fits with the average level spacing ratio being Wigner-Dyson for these small values of α (Sec. III) and each multiplet behaving separately as a random matrix (Sec. IV). The nonconvex entanglement entropy of these curves corresponds to a nonconvex microcanonical entropy and to ensemble inequivalence (see Sec. V).
Increasing α the multiplet structure disappears, first at higher, then at lower energy densities, as one can see in Fig. 12 (a) and (c) already for α = 0.5 and α = 0.75. The scatter plot for α = 2 and h = 0.1 [ Fig. 12 (b)] is remarkable. Here the scatter plot looks fuzzy and loses the smoothness typical of ETH. For this value of h, α = 2 corresponds to a minimum in the level spacing ratio (see Fig. 1 (a) ).
Let us move to quantify the smoothness of the entanglement-entropy curves. Considering S (µ) N/2 , we wish to characterize its eigenstate to eigenstate fluctuations. In ETH these fluctuations should be smaller compared to other contexts, because S (µ) N/2 should resemble the microcanonical curve, smooth in E µ . In order to quantify the fluctuations we consider Here, |ϕ µ and |ϕ µ+1 are "nearby eigenstates" [81] with the E µ and E µ+1 in increasing order. (unique for α > 0 and inside H S , where there are no degeneracies.) A quantity similar to M was introduced in [81] in the disordered Heisenberg chain taking instead of S (µ) N/2 the local magnetizations. In case of a system obeying ETH, M is expected to exhibit a rapid decay upon increasing system size N .
We plot M versus N in Fig. 13. We compare with the case of the α → ∞ (nearest-neighbour) Ising model in transverse field in Fig. 13 (b) and (c). The nearestneighbour model is integrable [95], and, consistently with that, the value M stays more or less constant with the size N . On the opposite, in the long-range model Eq. (1), M clearly decreases with N for most of the considered values of α. We emphasize that this occurs for the small values of α, but we cannot tell if this is due to the entanglement-entropy curves getting smoother inside the multiplets or to the fact that the multiplets tend to merge with each other for increasing N . We see that there is a close correspondence between the decay of M with N and the Wigner-Dyson value of the level spacing ratio (see Fig. 1). Indeed, the only conditions where we see something different from a decrease of M with N in Fig. 13 correspond to values of α where the average level spacing ratio has not yet attained the Wigner-Dyson value. This is true for α = 8 [ Fig. 13 (b), (c)] and, as we have argued in Sec. III, this is most probably a finite-size effect. This is also true for h = 0.1 and α = 2, 2.25 [ Fig. 13. (b)]. The effect is very strong for α = 2, again suggesting a connection with the integrability of other α = 2 long-range spin chain models.
Another quantitative analysis relevant for the study of ETH is the comparison with the Page value. ETH eigenstates with the largest entanglement are expected to approach the so-called Page value [100] upon increasing the system size N (the Page value corresponds to the entanglement entropy of a fully-random state [99]). We want to quantitatively probe this fact and consider the following two quantities introduced in [20]. The first one is defined as The rationale is that the logarithm overweights the smallest values of the argument and the high-entropy statescorresponding to the smallest values of the difference in the argument -give the strongest contribution to the average. If the highest-entropy states tend to the Page value, Λ S (N ) takes more and more negative values. In order to define the second quantity, we need to first define the integer number 1 ≤ µ * ≤ dim H S as the value of µ such that the quantity |S N/2 | is minimum over µ. Restricting the average of the entanglement entropy to states around the energy E µ * , we focus on the highest entropy states, the ones nearest to the Page value. More formally, if we term the width of the energy spectrum as ∆E(N ) = max µ (E µ ) − min µ (E µ ), we restrict the sum to the states with eigenenergy E µ ∈ [E µ * − f 2 ∆E(N ), E µ * + f 2 ∆E(N )] (call their number N f ). In this way we can define (A4) We choose f = 0.2, so that the sum is restricted around the state with entropy nearest to the Page value, that's to say to the infinite-temperature value. If Λ S (N ) and (S (Page) N/2 − S N/2 f )/N get smaller, the system becomes more ETH.
We report the results for Λ S (N ) versus α for different values of N in Fig. 14 (a), (c), and those for (S (Page) N/2 − S N/2 f )/N in Fig. 14 (b), (d). The steady decrease with N for h = 0.5 suggest a tendency to ETH for increasing system size. The largest-α crossing point between curves with nearby values of N tends to shift right for increasing N . The increase in N for large α is therefore a finite-size effect. Results for h = 0.1, on the opposite, are not that conclusive. Although the behavior at small and large α is similar to the h = 0.5 case, we find an interval of α (α ∈ [1, 1.5]) where both the considered quantities seem to saturate with N . Quite remarkably, in this interval of α the average level spacing ratio is significantly different from the Wigner-Dyson value [see Fig. 1 (a)] and probably finite-size effects are too strong.
Appendix B: Hilbert-Schmidt distance from the infinite-range model The Hilbert-Schmidt distance is an operator distance used in quantum information [91,92] and is defined by the norm Ô HS = Tr Ô †Ô . We are going to show that the Hilbert-Schmidt distance of the Hamiltonian at α > 0 from the infinite-range Hamiltonian at α = 0 increases linearly with α when α is small.
We consider the Hamiltonian Eq. (1), and we want to quantify the Hilbert-Schmidt distance ofĤ (α) from its infinite-range α = 0 counterpartĤ (0) . We define the distance as Taking the trace, all term but the first one vanish, so that d(α, N ) = 2 N/2 N i,j, i =j We numerically compute this quantity for various values of N and report it versus α in Fig. 15. We clearly see that it increases linearly in α for small α.
We strongly remark that, for α < 1, d(α, N )/2 N/2 fast saturates to a constant when N is increased. This point is crucial: The fact that d(α, N )/2 N/2 is asymptotically constant with N is at the root of our argument in Sec. IV. This result can be seen in Fig. 15 and can also be analytically checked in the large-N limit, by using translational invariance and writing approximately