Long-range Ising chains: eigenstate thermalization and symmetry breaking of excited states
LLong-range Ising chains: eigenstate thermalization and symmetry breaking of excitedstates
Angelo Russomanno, Michele Fava, and Markus Heyl Max-Planck-Institut f¨ur Physik Komplexer Systeme,N¨othnitzer Straße 38, D-01187, Dresden, Germany Rudolf Peierls Centre for Theoretical Physics, Clarendon Laboratory,University of Oxford, Oxford OX1 3PU, United Kingdom
We use large-scale exact diagonalization to study the quantum Ising chain in a transverse fieldwith long-range power-law interactions decaying with exponent α . Analyzing various eigenstate andeigenvalue properties, we find numerical evidence for ergodic behavior in the thermodynamic limitfor α > i. e. for the slightest breaking of the permutation symmetry at α = 0. Considering anexcited-state fidelity susceptibility, an energy-resolved average level-spacing ratio and the eigenstateexpectations of observables, we observe that a behavior consistent with eigenstate thermalizationfirst emerges at high energy densities for finite system sizes, as soon as α >
0. We argue thatergodicity moves towards lower energy densities for increasing system sizes. While we argue thesystem to be ergodic for any α >
0, we also find a peculiar behaviour near α = 2 suggesting theproximity to a yet unknown integrable point. We further study the symmetry-breaking properties ofthe eigenstates. We argue that for weak transverse fields the eigenstates break the Z symmetry, andshow long-range order, at finite excitation energy densities for all the values of α we can technicallyaddress ( α ≤ . I. INTRODUCTION
Thermalization in classical Hamiltonian systems is wellunderstood in terms of chaotic dynamics and the relatedessentially ergodic exploration of the phase space [1–3].From the quantum point of view the physical mech-anism is quite different, involving the eigenstates ofthe Hamiltonian being fully random strongly entangledstates which appear thermal from the point of view oflocal measurements. This is the paradigm of eigenstatethermalization (ETH) introduced in Refs. [4–7].In many cases one can identify a correspondence be-tween classical and quantum thermalization. This canbe seen from a variety of different arguments, includingBerry’s random-wave conjecture for the energy eigenfunc-tions [5, 8–11], the analogy between quantized chaoticsystems and random matrix theory [12], and the semi-classical periodic orbit expansion, assuming a certain ran-domness for the periodic orbits [13]. In summary, whena classically chaotic Hamiltonian is quantized, it givesrise to a Hamiltonian matrix which is essentially a ran-dom matrix (with some caveats [14]) and its eigenstateshave exactly the required properties for local thermaliza-tion. This correspondence is nevertheless highly non triv-ial, because quantum effects can give rise to ergodicity-breaking phenomena with no analog in the classical do-main.This has been very evident in the research of the lastyears in the context of nearest-neighbour many-bodysystems. From a classical point of view, any nearest-neighbour nonlinear Hamiltonian system with more thantwo degrees of freedom and no conservation law beyondenergy gives rise to chaos and essentially ergodic dynam-ics [1]. From the quantum point of view the situation is different: There are plenty of examples of ergodicitybreaking in many-body quantum systems like many bodylocalization (see [15] for a review), many body dynamicallocalization [16–19], and other forms of regular behaviourin quantum many-body nonlinear systems [7, 20–23].Comparatively less attention has been provided to thequantum properties of models with long-range interac-tions. In the classical case the behaviour of these sys-tems is quite remarkable. At high energy densities, thechaotic fraction of the phase space can become vanishingin the thermodynamic limit [1], giving rise to an effec-tively regular behaviour dominated by one or few degreesof freedom [24–29] which has been exploited to obtain aclassical Hamiltonian time crystal [30].The situation in the quantum regime has been stud-ied in the case of long-range spin-1 / α , the exponent ofthe power-law decaying interaction. One very well stud-ied case is the Ising model with infinite-range interac-tions ( α = 0 – the so called Lipkin-Meshkov-Glick model)which is known to be integrable [31–33]. It is also knownthat the isotropic Heisenberg chain with power-law inter-actions with exponent α = 2 is integrable [34, 35] as wellas some anisotropic spin-chain models with α = 2 [36–38]. Spin chains with disorder and power-law interactionsare known to undergo a transition between a regularmany-body-localized-like and an ergodic phase [39–45].Quite remarkably, these models show a regular many-body-localized behaviour when the power-law interac-tions decay fast, with an exponent α above a certainthreshold [44, 45].The situation is less clear when there is no disorder andone considers a clean long-range interacting spin mod-els with power law exponent. Although these models a r X i v : . [ c ond - m a t . s t a t - m ec h ] D ec have been extensively studied in the context of quantumquenches [46–56], and their dynamics has attracted a lotof experimental interest [40, 57–63], an analysis of thethermalization properties of the eigenstates is generallylacking. A significant exception is [64] which showed er-godicity at low energies for α = 1 . / α and energies.The dynamics of this model has been intensively studied,mostly in connection with the persistence of long-rangeorder in the asymptotic state of the dynamics [33, 46–49, 51–55], for different values of α and small transversefield, but it is not known if this asymptotic state is ther-mal.Based on an analysis of the eigenstate and spectralproperties, we find numerical evidence that for 0 < α (cid:28) α = 0 the model isperfectly integrable. Any infinitesimally small α > α = 0. The α = 0integrability is indeed a quite peculiar one: There aremany invariant subspaces transforming differently underthe permutation symmetries and many of these subspaceshave identical level structure [32]. So, because of degen-eracy, even the smallest α > α →
0, as we show in Appendix Aby discussing the behaviour of an operator distance inthis limit.The integrable point at α = 0 affects spectral structureat 0 < α (cid:28)
1: The spectrum is organized in multipletsseparated by gaps. At α = 0 these multiplets are per-fectly degenerate due to the full permutation symmetry,when α (cid:28) α (cid:29) α → ∞ , in agreement withthe general picture of Ref. [65, 66]. For intermediate val-ues of α the model behaves ergodically, with few excep-tions. A quite remarkable one is an ergodicity breakingoccurring for a weak transverse field h = 0 . α around2 and system sizes up to N = 22. This appears rathersuggestive in view of the fact that the Heisenberg modelat α = 2 is exactly integrable [34], but our numerics doesnot allow us to settle whether this behaviour persists inthe thermodynamic limit. Nevertheless we cannot ex-clude the proximity to an integrable point.We also systematically explore ergodicity in the modelacross the many-body spectrum. For finite-size systemsergodicity first emerges for high energy densities. We ar-gue that the ergodic region extends towards lower energy densities upon increasing system size (this can be seen inthe results for the excited-states fidelity of Sec. III D).Nevertheless, our numerics can reach too small systemsizes for providing a definite proof.Then, we systematically explore the Z -symmetry-breaking properties at the level of single eigenstates. Thepresence of long-range order at low energies for this modelwas already demonstrated at α = 1 . Z -symmetry breaking in the thermodynamic limit)with spectral pairing of eigenvalues. When α ≤ .
5, fortransverse fields h = 0 . h = 0 .
5, we find a thresholdin the energy density below which the excited states showlong-range order (and symmetry breaking in the thermo-dynamic limit). This threshold is the broken symmetryedge e ∗ ; We are able to estimate its value and study itsdependence on α . We find that the broken symmetryedge depends continuously on α and we find that it isabove the ground-state value. This implies that there islong-range order at finite excitation energy density: fora finite fraction of the energy spectral width the eigen-states show long-range order, like the α = 0 and the dis-ordered case. This implies that the low-energy quench ofan initial state polarized along z gives rise to a magne-tization persisting for long times in the thermodynamiclimit, which confirms the findings of [47, 55].The paper is organized as follows. In Sec. II we de-fine the model Hamiltonian. In Sec. III we study thethermalization properties. First, we consider the spectralproperties, focusing on the level spacing ratio averagedover all the energy spectrum (Sec. III A), and then on itsversion averaged over the energy bins (Sec. III A 1), inorder to explore ergodicity across different energy densi-ties. In Sec. III A 2 we discuss the fact that the spectrumis organized in multiplets for small α and small energydensity and study the behaviour of the density of states.In Sec. III B we study ergodicity by means of the en-tanglement entropy properties of the eigenstates and inSec. III C we do the same by considering the eigenstateexpectation values of a local operator, the longitudinalnearest-neighbour correlation. In Sec. III D we consideran excited-state fidelity susceptibility. They suggest aregular behaviour at small energy density which, at leastfor h = 0 .
1, is confirmed by the level spacing ratio av-eraged over the energy bins. In Sec. IV we discuss thespectral-pairing properties and find the broken-symmetryedge as a function of α . In Appendix A we discuss theHilbert-Schmidt distance of the α > α = 0 Hamiltonian, showing its continuity in thelimit α →
0. In Sec. V we draw our conclusions.
II. MODEL HAMILTONIAN AND RELEVANTQUANTITIES
In this work we study the long-range interacting quan-tum Ising chain in a transverse field:ˆ H = − JN ( α ) N (cid:88) i,j, i (cid:54) = j ˆ σ zi ˆ σ zj D αi, j + h N (cid:88) i =1 ˆ σ xi . (1)Here, σ αi with α = x, y, z denotes the Pauli matrices atlattice site i = 1 , . . . , N with N the system size. We useperiodic boundary conditions implemented through thedefinition [51] D i, j ≡ min[ | i − j | , N − | i − j | ]; we definethe Kac factor [67] N ( α ) ≡ N − (cid:80) i,j, i (cid:54) = j D αi, j in order topreserve extensivity of the Hamiltonian throughout thefull spectrum.Throughout the paper we will use exact diagonaliza-tion. We will largely exploit the translation, inversionand Z (ˆ σ zi → − ˆ σ zi ) symmetries of the model in orderto restrict to an invariant subspace of the Hamiltonian.In Section III we will restrict to the subspace fully sym-metric under all the symmetries of the Hamiltonian. Wecall this Hamiltonian eigenspace H S and we define it asthe subspace corresponding to the zero-momentum sectorand also symmetric with respect to inversion and Z sym-metry. For future convenience we define N S ≡ dim H S .In Section IV we will be interested in the spectral pair-ing properties of the model, which requires to considerboth Z symmetry sectors. Therefore we will considerthe subspace corresponding to the zero-momentum sec-tor and even only with respect to inversion. Throughoutthe paper, we will denote the eigenstates of the Hamilto-nian | ϕ µ (cid:105) and the corresponding eigenenergies E µ (takenin increasing order), while always specifying which sub-space we are considering.In the limit α → ∞ the model in Eq. (1) reduces to thenearest-neighbour quantum Ising chain. This model is in-tegrable and undergoes a quantum phase transition: Itsground state breaks the Z symmetry for h < Z symmetry, with a split-ting exponentially small in the system size. The statesin the doublet show long-range order and the doublet be-comes degenerate in the thermodynamic limit, giving riseto symmetry breaking.In the limit α = 0, on the opposite, Eq. (1) reducesto the Lipkin-Meshkov-Glick model. This model is alsointegrable, thanks to the full permutation symmetry, andit shows a symmetry-broken phase for h <
1. In contrastto the α → ∞ case, all the spectrum up to an extensiveenergy N e ∗ is organized in doublets with exponentiallysmall splitting and the corresponding eigenstates havelong range order and break the Z symmetry in the ther-modynamic limit [31–33]. Due to the full permutationsymmetry, the Hilbert space is factorized in a number ofinvariant subspaces, differently transforming under the permutation symmetries [32]. The number of these sub-spaces is exponential in N , and many of them have thesame level structure. This gives rise to massively degen-erate multiplets, whose levels belong to different symme-try sectors, a property which will be quite relevant in thefollowing.In the remainder of the paper we will consider the caseof intermediate α , trying to characterize for which modelparameters the dynamics is ergodic or regular and, fur-thermore, which fraction of the spectrum is organized insymmetry-breaking doublets.For the study of the eigenstate properties we will con-sider the longitudinal nearest-neighbour correlation op-erator: ˆ G = 1 N N (cid:88) j =1 ˆ σ zj ˆ σ zj +1 , (2)as a representative for local observables. We will focuson the properties of the eigenstate expectation values G µ ≡ (cid:104) ϕ µ | ˆ G| ϕ µ (cid:105) . (3)Another quantity we consider is the half-chain entangle-ment entropy for the eigenstates. This is not a localobject because it involves correlations extending up to adistance N/
2, but quite remarkably eigenstate thermal-ization has been proved valid for reduced density matricesrelated to subsystems up to this size [70]. Considering aneigenstate | ϕ µ (cid:105) , and decomposing the system in two parts A and B in physical real space, we define S ( µ ) A = − Tr[ˆ ρ A log ˆ ρ A ] with ˆ ρ A = Tr B [ | ϕ µ (cid:105) (cid:104) ϕ µ | ] . (4)Specifically, we will focus on the half-system entangle-ment entropy S ( µ ) N/ taking each bipartition made up of N/ G µ and S ( µ ) A are equal to their microcanonicalvalue at energy E µ , up to relative fluctuations decreas-ing with the system size. (The microcanonical equivalentof S ( µ ) N/ is the thermodynamical entropy for half of thesystem.)In order to characterize the eigenstate to eigenstatefluctuations of our observables, we will consider M ( O ) ≡ N S − N S − (cid:88) µ =1 |O µ +1 − O µ | , (5)with O µ ( O µ = S ( µ ) N/ or O µ = G µ ) restricted to thesubspace H S . Here, | ϕ µ (cid:105) and | ϕ µ +1 (cid:105) are “nearby eigen-states” [71] with the E µ and E µ +1 in increasing order . M ( O ) was introduced in [71] taking as observables thelocal magnetizations in the disordered Heisenberg chain.In case of a system obeying ETH, M ( O ) is expectedto exhibit a rapid decay upon increasing system size N .Practically, this becomes manifest as a reduction of fluc-tuations as N is increased when plotting O µ as a func-tion of E µ showing convergence to the smooth functionas expected from ETH. In contrast with that, in case ofmany-body localized systems, M ( O ) is expected to tendto a non-vanishing constant for increasing N , as verifiedin [71].It is, however, not fully clear how this behaves forclean nonergodic systems, which appear in some cases inthe present work and correspond to a non Wigner-Dysonlevel spacing ratio (see Fig. 1). Actually, we find two dif-ferent behaviours for M ( S N/ ) and M ( G ). On one side M ( S N/ ) shows a clear correspondence with the averagelevel spacing ratio behaviour: Wigner-Dyson correspondsto a decay in N , otherwise there is a saturation or evenan increase in N (see Sec. III B). So, for M ( S N/ ), thesituation is quite similar to the case of the crossover fromETH and many-body localization considered in [71].On the other side M ( G ) always shows a clear decaywith N (see Sec. III C). So, in order to distinguish ergodicfrom ergodicity breaking, we need a finer analysis. Wehave to find the form of the decay behaviour in the ETHcase and compare with it the actual behaviour of M ( G ).We now derive the ETH behaviour and show the finalresult in Eq. (11). We start considering that in the ETHcase one knows that [72–76] for a local operator (in thiscase we take as local operator G ) |G µ +1 − G µ | ∼ e −S ( E µ +1 ) / (6)where S ( E µ − ) is the thermodynamical entropy at en-ergy E µ − . So we see that in the ETH case M ( G ) ∼ N S − (cid:80) N S − µ =1 e −S ( E µ +1 ) / . We can estimate this quan-tity by proceeding as in [20, 77]. By using the density ofstates ρ ( E ) [Eq. (15)] we can write M ETH ( G ) ∼ (cid:82) ρ ( E )e −S ( E ) / d E (cid:82) d Eρ ( E ) . (7)Considering that ρ ( E ) = e S ( E ) and introducing the re-duced variables s ≡ S /N and (cid:15) ≡ E/N we can write M ETH ( G ) ∼ (cid:82) e Ns ( (cid:15) ) / d E (cid:82) e Ns ( (cid:15) ) d E . (8)Expanding s ( (cid:15) ) around the maximum (corresponding tothe T = ∞ value) as s ( ε ) = s ( ε ∞ ) − ( ε − ε ∞ ) W (9)and applying the saddle-point approximation in the limit N (cid:29) M ETH ( G ) ∼ √ − Ns ( (cid:15) ∞ ) / . (10)Because e Ns ( (cid:15) ∞ ) ∼ N S [73], we find M ETH ( G ) ∼ N / S . (11) This is the ETH value with whom we will compare the M ( G ) results in order to assess quantitatively whetherthe system obeys ETH or not.We remark that it is crucial to compute M ( O ) in sub-spaces accounting for the symmetries of the consideredmodel which avoids degeneracies in the spectrum andtherefore eigenstates can be put in increasing order ina unique way. Therefore, in the case α = 0, where thespectrum is organized in massively degenerate multiplets, M ( O ) has no definite meaning in the subspace H S , be-cause there is not a unique way of ordering the eigen-states in each multiplet. On the opposite, for α > α → ∞ , the spectrum inside H S is non-degenerate andwe can evaluate M ( O ). In the next Section we are go-ing to use this and other quantities to probe eigenstatethermalization in our model. III. ERGODICITY PROPERTIES OFEIGENSTATES AND EIGENENERGIES
In the following we aim at exploring the ergodicityproperties of the long-range interacting quantum Isingchain. For that purpose we will study the system’s spec-tral and eigenstate properties using large-scale exact di-agonalization.
A. Spectral properties
The model in Eq. (1) is integrable for the limits α = 0(infinite-range case) and α → ∞ (nearest-neighbourcase). We now aim at exploring the behaviour at inter-mediate α , which has not yet been extensively studied.As we will show, we identify various numerical evidencefor ergodicity for all the intermediate values of α . Forconcreteness, we don’t scan extensively across the trans-verse fields, but rather focus on the two representativevalues h = 0 . h = 0 .
5. In Fig. 1 we investigate thespectral properties of the model as a function of α uponvarying the system size N . Specifically, we plot the aver-age level spacing ratio, r (introduced in [71]), which hasturned into a central probe for quantum ergodicity andis defined as r = 1 N S − N S − (cid:88) µ =1 min( E µ +2 − E µ +1 , E µ +1 − E µ )max( E µ +2 − E µ +1 , E µ +1 − E µ ) . (12)Let us start by considering the limit α (cid:28)
1. Wesee in Fig. 1 that in this range of α r is close to theWigner-Dyson value r WD = 0 . α = 0is unstable to a small perturbation in α . Thus, breakingthe full permutation symmetry at α = 0 turns the systeminto an ergodic one.As we have already discussed in Sec. II, the multipletsat α = 0 do not correspond to a given permutation sym-metry class, but instead contain states which belong todifferent invariant subspaces, differently transforming un-der permutation. Indeed, there are many identical sub-spaces, with the same energy levels [32]. The momentperturbation symmetry is weakly broken by α (cid:28)
1, thestates inside each multiplet can mix and so all the sub-spaces are mixed by the Hamiltonian. This is a dramaticchange of the Hilbert space leading to ergodicity and,since there is no gap to protect the subspaces from mix-ing, this change is likely to happen abruptly as soon as α >
0. We remark that this fact is peculiar for infinite-range models with permutation symmetry leading to alarge number of identical Hilbert subspaces with a simi-lar form of the Hamiltonian. These subspaces are degen-erate with each other and even the slightest breaking ofthis symmetry mixes all of them. This is very differentfrom other cases, where the breaking of the degeneracy isnot associated with the breaking of a symmetry as strongas the full permutation one [20].We remark that there is no discontinuity in the Hamil-tonian operator when α becomes different from 0. Inorder to appreciate this fact, we consider the Hilbert-Schmidt distance, an operator distance defined by thenorm (cid:13)(cid:13)(cid:13) ˆ O (cid:13)(cid:13)(cid:13) HS = (cid:114) Tr (cid:16) ˆ O † ˆ O (cid:17) . This distance has foundapplications in quantum information [79], namely as en-tanglement witness [80]. As we show in Appendix A,the Hilbert-Schmidt distance of the Hamiltonian at α > α = 0 increaseslinearly with α when α is small. The situation is ex-actly the same as when a perturbation δ · ˆ V is applied:the Hilbert-Schmidt distance of the perturbed Hamilto-nian from the nonperturbed one linearly increases with δ . Therefore, it should be possible to understand thesmall- α spectral properties using perturbation theory in α , and that will be the subject of further studies.For h = 0 . α = 2which does not appear for h = 0 .
5. For h = 0 . α around 2 there is also a peculiar behaviour of all theother quantities probing eigenstate thermalization, as weare going to see in the next subsections. It is impor-tant to remind that there are spin models with power-law interactions decaying with α = 2 that are integrable,such as the long-range isotropic Heisenberg chain [34]or other long-range models anisotropic like ours [36–38].Assessing whether this ergodicity breaking persists in thethermodynamic limit at this particular parameter regimewould, however, require an analysis of even larger sys-tems which goes beyond what we can currently achieve.Nevertheless the persistence of ergodicity breaking up to N = 22 is quite remarkable and a valuable future ex-ploration would investigate whether this phenomenon isrelated to the proximity to an integrable point.For large α we see in Fig. 1 that the curves tend to-wards the Poisson value r P (cid:39) . α → ∞ integrable nearest-neighbour model. We can argue that this behaviour ofthe average level spacing ratio is an effect of the proxim-ity of the integrable α → ∞ point and that the systembecomes ergodic in the thermodynamic limit. First, wesee a clear increase of r with N for α ≤
6. Furthermore,Refs. [65, 66] consider a free-fermion model — as is thecase for the α → ∞ point — and show that an arbitrarilysmall integrability-breaking next-nearest-neighbour in-teraction restores thermalization in the thermodynamiclimit. Similarly, in our case, for α (cid:29)
1, the next-nearestneighbour terms are the stronger ones breaking the inte-grability of the nearest-neighbour α → ∞ model.In order to roughly estimate at which point the sys-tem becomes ergodic for α (cid:29)
1, we focus on the be-haviour of these terms. We now argue that it is suffi-cient to compare the next-nearest neighbour interactionterm with the relevant gap ∆ of the integrable nearest-neighbour model. The next-nearest neighbour term is oforder V ∼ J/ ( N ( α )2 α ). In order to understand the rel-evant gap of the nearest-neighbour model, we have totransform it in the fermionic representation using theJordan-Wigner transformation [82]. In this representa-tion, the nearest-neighbour model is integrable and itsexcitations are fermionic quasiparticles [69, 83] with en-ergy (cid:15) k = N ( α ) (cid:113) ( J cos k − h ) + J sin k . We have k ∈ [0 , π ] and, for finite system size N , k we can takeonly N discrete equally spaced values. In the fermionicrepresentation the next-nearest-neighbour term becomesa four-fermion term which induces inelastic scattering be-tween the fermionic quasiparticles. If momenta k and k go into momenta k , k + k − k , the relevant gap is∆ = (cid:15) k + (cid:15) k + k − k − (cid:15) k − (cid:15) k . We can provide a roughestimate of ∆ by taking twice the bandwidth of (cid:15) µ anddividing it by N , the number of allowed equally-spaced k values. We find ∆ ∼ hN · N ( α ) (13)Imposing that V (cid:38) ∆, one finds the condition for theergodic behaviour as α (cid:46) α ∗ ≡ log N + log (cid:18) J h (cid:19) (14)Applying this formula, one finds α ∗ (cid:39) . N = 22 and h = 0 . α ∗ (cid:39) . N = 22 and h = 0 .
5. In both cases,confirm Fig. 1, for α = α ∗ , r starts deviating from theergodic Wigner-Dyson value. The estimate seems indeedreasonable, because α ∗ → ∞ when N → ∞ as one ex-pects from the discussion above. Moreover, α ∗ decreaseswith increasing h , in agreement with the qualitative ob-servations one can do from Fig. 1. (a)(b)FIG. 1. Average level spacing ratio versus α . We consider h = 0 . h = 0 .
1. Energy-bin averaged level spacing
On a general level, ergodicity (or its breaking) can bespecific to the considered energy densities [84]. In orderto spectrally explore ergodicity, we will therefore con-sider averaging the level spacing ratio not just over thefull spectrum, but rather over bins of states at differentenergy densities [84]. We order the levels in the sense ofincreasing energy and divide them in N bin bins with thesame number of levels in each. Each bin has a differentextension in energy, as the density of states is not con-stant as a function of energy (see Sec. III A 2 for details).While the bins at the edges of the spectrum thereforerepresent a larger energy range due to smaller density ofstates, this procedure ensures sufficient statistics in ev-ery bin. To each bin we associate a energy E bin choosingthe middle point of the energy extension. In each bin weevaluate the average level spacing ratio, restricting thesum in Eq. (12) to the states inside the bin. We termthis quantity (cid:104) r (cid:105) bin .We plot (cid:104) r (cid:105) bin versus α and E bin /N for two values of h in Fig. 2. For the case h = 0 . N = 22 [panel (b)],when α (cid:46)
5, we see that there are fluctuations in E bin but (a)(b)FIG. 2. (cid:104) r (cid:105) bin versus α and E bin /N . [Panel (a)] h = 0 . h = 0 . N = 22 and N bin = 40. (cid:104) r (cid:105) bin (cid:39) r WD at all energy densities. So, the ergodicityis quite uniform in energy. Also the finite-size regularbehaviour at large α is quite uniform in energy. Thecrossover between ergodicity and regularity, instead, doesnot happen uniformly throughout the spectrum, but wecan identify energy intervals which retain ergodicity atmuch larger α than the rest of the spectrum. For h = 0 . (cid:104) r (cid:105) bin nearer to the Poisson value at small energy densityfor α ∈ [0 . , α = 2, there are alsoother intervals of energy density where the level spacingratio averaged over bins is Poisson-like. This fact shouldbe related to the minimum in r seen around α = 2 inFig. 1 (a).
2. Multiplets and density of states
In this section we argue that, for small α values, wesee a multiplet structure in the eigenstates at small andintermediate energy. Indeed, at α = 0 the spectrum isorganized in massively degenerate multiplets, due to thefull permutation symmetry. For α > α (cid:28)
1, themultiplet structure persists at low energy (the degeneracyinside each multiplet is of course broken). In Fig. 3 weshow an example of this fact. We plot E µ versus µ/ N S for h = 0 . α , α = 0 and α = 0 . α = 0 there are many degenerate multiplets at all FIG. 3. Plot of E µ versus µ/ N S for h = 0 . N = 22 and twodifferent values of α . energies, as we can see in the magnifying insets. For α = 0 .
15 the multiplets still exist at low energy (leftinset) and merge into a smooth continuum at large energy(right inset).It is important to estimate how the energy densitythreshold between these two different behaviours changeswith N . One possibility is that it moves towards theground-state value for increasing N and that for N → ∞ the multiplet structure disappears, as we know occursfor the Bose-Hubbard model [20]. To further study theseaspects of the multiplet structure we now consider thedensity of states defined as ρ ( E ) = 1 N S (cid:88) µ δ ( E − E µ ) . (15)In Fig. 4 (a) we plot the density of states versus theenergy density for α = 0 . h = 0 . α = 0 . h = 0 . α , aswe can see in the density-of-states plots of Fig. 5, bothfor h = 0 . h = 0 . α >
0. So, we see that at the lowestenergy densities the spectrum is made by two or threewell-separated multiplets with few states (they are notcaptured by (cid:104) r (cid:105) bin in Fig. 2 because it has not enough (a)(b)FIG. 4. Density of states ρ ( E ) versus E/N for different valuesof N . [Panel (a)] α = 0 . h = 0 .
1. [Panel (b)] α = 0 . h = 0 . resolution at small energy densities). These multipletspossibly survive even for larger system sizes giving rise toan effective nonergodic behaviour in quenches involvingsmall excitation energy densities. This is probably thesituation occurring in [47, 55] for α < B. Half-system entanglement entropies
After having explored the spectral properties of theconsidered long-range Ising chains we now aim to targettheir entanglement properties. In Fig. 6 we show the en-tanglement entropy S ( µ ) N/ [defined in Eq. (4)] versus thecorresponding eigenstate energy E µ . Let us start consid-ering first a small value of α , α = 0 .
05 [panels (a), (c)].Quite remarkably, inside each multiplet, the S ( µ ) N/ versus E µ look like smooth curves, as appropriate for a thermal-izing system where the eigenstate entropies correspond tothe microcanonical values [70]. This result nicely fits withthe average level spacing ratio being Wigner-Dyson forthese small values of α . We are going to quantitativelyconfirm this qualitative observation with the M ( S N/ )analysis below. Nevertheless, in Fig. 6 (a) and (c) we do (a)(b)FIG. 5. Density of states ρ ( E ) versus E/N for different valuesof α and h . System size N = 22. h = 0 . h = 0 . not see just a single smooth curve, as in the usual ETHsystems, but many disconnected smooth curves, becauseof the strong effect of the multiplets. For the systemsizes we have access to ( N ≤
22) the multiplet structureis very evident and affects a lot our numerics (as we aregoing to better specify in Sec. IV).Increasing α the multiplet structure disappears, firstat higher, then at lower energy densities, as one can seein Fig. 6 (a) and (c) already for α = 0 . α = 0 . α = 2 and h = 0 . M ( S N/ ) analysis below. For thisvalue of h , α = 2 corresponds to a minimum in the levelspacing ratio (see Fig. 1 (a) ). Still, we do not know ifthis is a finite size effect or it is a sign of an integrable-likebehaviour persisting also at larger system sizes.At this point we get more quantitative in orderto probe ETH with these quantities and we consider M ( S N/ ) we have defined in Eq. (5). We plot M ( S N/ )versus N in Fig. 7. We compare with the case of the α →∞ Ising model in transverse field labeled with “Shortrange” in Fig. 7 (b) and (c). The nearest-neighbour model is integrable [83], then breaks ergodicity and, con-sistently with that, the value M ( S N/ ) stays more or lessconstant with the size N . On the opposite, in the long-range model Eq. (1), M ( S N/ ) clearly decreases with N for most of the considered values of α . We see that thereis a close correspondence between the decay of M ( S N/ )with N and the average level spacing ratio having at-tained the Wigner-Dyson value (see Fig. 1). Indeed, theonly conditions where we see something different froma decrease of M ( S N/ ) with N in Fig. 7 correspond tovalues of α where the average level spacing ratio has notyet attained the Wigner-Dyson value. This is true for α = 8 [Fig. 7 (b), (c)] and, as we have argued above, thisis most probably a finite-size effect. This is also true for h = 0 . α = 2 , .
25 [Fig. 7. (b)]. The effect is verystrong for α = 2, again suggesting a connection with theintegrability of other α = 2 long-range spin chain models.
1. High-entropy states
In the following we aim to study the entanglementproperties of eigenstates with high entropy. Ergodiceigenstates with the largest entanglement are expectedto approach the so-called Page value [86] upon increas-ing the system size N (the Page value corresponds tothe entanglement entropy of a fully-random state [85]).As probes we consider the following two quantities intro-duced in [20]. The first one is defined asΛ S ( N ) = 1 N S (cid:88) µ log (cid:16) | S (Page) N/ − S ( µ ) N/ | (cid:17) . (16)The rationale is that the logarithm overweights the small-est values of the argument. Because the high-entropystates correspond to the smallest values of the differencein the argument, they give the strongest contribution tothis average. If the highest-entropy states tend to thePage value, Λ S ( N ) takes more and more negative values.In order to define the second quantity, we need to firstdefine the integer number 1 ≤ µ ∗ ≤ dim N S as the valueof µ such that the quantity | S (Page) N/ − S ( µ ∗ ) N/ | is minimumover µ . Restricting the average of the entanglement en-tropy to states around the energy E µ ∗ , we focus on thehighest entropy states, the ones nearest to the Page value.More formally, if we term the width of the energy spec-trum as D = E max − E min , we restrict the sum to thestates with eigenenergy E µ ∈ [ E µ ∗ − f D/ , E µ ∗ + f D/ N ). In this way we can define (cid:10) S N/ (cid:11) f = 1 N (cid:88) µ s . t . E µ ∈ [ E µ ∗ − fD/ ,E µ ∗ + fD/ S ( µ ) N/ . (17)We choose f = 0 .
2, so that the sum is restricted aroundthe state with entropy nearest to the Page value, that’sto say to the infinite-temperature value. If Λ S ( N ) and( S (Page) N/ − (cid:10) S N/ (cid:11) f ) /N get smaller, the system becomesindeed more ergodic. (a) (b)(c) (d)FIG. 6. Half-system entanglement entropy of the eigenstates, scatter plot of S µN/ versus E µ for different values of the parameters.We consider N = 20 and h = 0 . h = 0 . N = 20, the value corresponding to a fully random state [85]. We report the results for Λ S ( N ) versus J for dif-ferent values of N in Fig. 8 (a), (c), and those for( S (Page) N/ − (cid:10) S N/ (cid:11) f ) /N in Fig. 8 (b), (d). Results for h = 0 . N for large α but this is most probably a finite-size effect.Indeed the largest- α crossing point between curves withnearby values of N tends to shift right for increasing N ,suggesting a slow tendency towards ergodicity. Resultsfor h = 0 .
1, on the opposite, are not that conclusive. Al-though the behaviour at small and large α is similar tothe h = 0 . α ( α ∈ [1 , . N . Quite remarkably, in this interval of α the aver-age level spacing ratio is significantly different from theWigner-Dyson value [see Fig. 1 (a)]. We do not know ifthis is a finite-size effect and further research is neededto clarify this point. C. Longitudinal nearest-neighbour correlations
Now we aim to study the behaviour of the expectationvalues G µ of the longitudinal nearest-neighbour correla- tion in eigenstates | ϕ µ (cid:105) [see Eq. (3)]. We first considerthe scatter plots of G µ versus E µ in Fig. 9. Most impor-tantly, these expectation values as a function of energydon’t always exhibit a smooth dependence with smallfluctuations, as expected in a system obeying ETH [6]even though the level spacing ratio Eq. (12) is close toWigner-Dyson and therefore ergodic.The most noteworthy case is α = 0 .
05 [Fig. 9 (a)and (b)] where we see many almost vertical lines, as manyas the multiplets. Each of these lines is a continuouscurve, so that the ETH paradigm seems respected justwithin a multiplet but not across them. Nevertheless,the average level spacing ratio is close to the Wigner-Dyson value. In order to understand this fact, we noticethat the number of discontinuity points (giving rise tothe large gaps among multiplets) is much smaller thanthe total number of states with vanishing fraction in thethermodynamic limit. Indeed, the number of discontinu-ity points scales as the number of distinct multiplets at α = 0, which is quadratic in N [32], while the number ofstates equals dim H S which is exponential in N .Another interesting case is provided by α = 0 . h = 0 . M ( S N / ) Nh = 0.1 α = 0.05 α = 0.25 α = 0.5 α = 0.75 α = 1.0 α = 1.5 (a) M ( S N / ) Nh = 0.1 α = 2.0 α = 2.25 α = 3.0 α = 8.0 α ➝ ∞ (b) -0.2-0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 10 12 14 16 18 20 M ( S N / ) N h = 0.5 α = 0.05 α = 0.5 α = 1.0 α = 2.0 α = 4.0 α = 8.0 α ➝ ∞ (c)FIG. 7. Plot of M ( S N/ ) versus N for different values of α and h and, for comparison, the Ising integrable α → ∞ nearest-neighbour model. energy. In the center of the spectrum we observe a quitesmooth curve with some small fluctuations, which ap-pears as a prototypical example of a system obeyingETH. The situation is very different at low energies.There we see very well separated multiplets. There isalso some noise in between suggesting a stronger mixingfor larger system sizes. Overall, this doesn’t seem to fol-low the predictions by ETH.. This different behaviour atlow and high energies for α = 0 . h = 0 . (cid:104) r (cid:105) bin (Sec. III A 1). For these sameparameters, Fig. 2 (a) shows a change of behaviour fromPoisson to Wigner Dyson at energy density (cid:39) − . G µ occurs roughly atenergy density ∼ − . /
20 = − .
125 [Fig. 9 (c)]. On theopposite, the case α = 0 . h = 0 . (cid:104) r (cid:105) bin [Fig. 2 (b)].For larger α [ α = 1 . α = 2 inFig. 9 (g), (h)] we see a fully developed ETH behaviourfor h = 0 .
5: very smooth curves with noise at the edgesof the spectrum [panels (f) and (h)]. On the opposite,for h = 0 . α = 2 is very regular-like. There issome noise suggesting a stronger mixing at larger systemsizes. Our numerics does not allow us to tell if this mix-ing develops into a full ergodicity in the thermodynamiclimit.Similarly to what we have done for the half-system en-tanglement entropy in Fig. 7, we consider the dependenceof M ( G ) on N S [Eq. (5)] in the double-logarithmic plotsof Fig. 10. For h = 0 . N S (and then with N ) corresponding toregularity in the large-size limit only for α = 2 , .
25. Wesee a decrease for all the other considered values of α . For α < ∼ / N γS , but γ is clearly different from the ETH value 1 / α → ∞ case and α ≥ .
25 we see that the curves tendto deviate from the ETH behaviour as the system sizeincreases. The slope decreases in absolute value for in-creasing N . For even larger N we expect a behaviourclearly different from 1 / N / S for the nearest-neighbourcase and (at some point) a recover of ergodicity for thefinite values of α .For h = 0 . / N / S ETH case only for α ≤
1, consistently with the ergodicity of most of theeigenstates [Fig. 2 (b)]. For α = 4 , α → ∞ thedecay is clearly slower. This marks integrability for thenearest-neighbour case and a regular behaviour persistingfor quite large system sizes for the finite values of α . D. Excited-state fidelity susceptibilities
While we have been concentrating up to now on explor-ing the properties of the long-range Ising chain by meansof the spectra and observables, we now aim to quan-tify ergodicity by studying properties of eigenstates. For1 (a) (b) -1.2-1-0.8-0.6-0.4-0.2 0 0.2 0.4 0.6 0.8 0 1 2 3 4 5 6 7 8 Λ S ( N ) α h = 0.5N = 16N = 18N = 20 (c) ( S P age N / - < S N / > f ) / N α h = 0.5N = 16N = 18N = 20 (d)FIG. 8. Plot of the quantities Λ S ( N ) [Eq. (16) – panels (a), (c)] and ( S (Page) N/ − (cid:10) S N/ (cid:11) f ) /N [Eq. (16) – panels (b), (d) – f = 0 . J for different values of N . h = 0 . h = 0 . N = 22 we consider 14000randomly chosen eigenstates, in the other cases all the spectrum. that purpose we consider here the overlaps between corre-sponding eigenstates at nearby h . Such overlaps are verywell known in the context of quantum information [87],it has been used to characterize ground-state quantumphase transitions [88] and it has been used in the contextof random matrices to define the so-called excited-statefidelity susceptibility [89], which we will now also con-sider in this section. In order to define it, we fix N andproceed as follows. We fix α and h and compute alleigenstates. Then we keep the same α , consider h + δh and evaluate again all eigenstates. Finally, we evaluatean overlap between two corresponding eigenstates in thefollowing way: we fix the eigenstate | ϕ µ ( α, h ) (cid:105) at α and h and choose the eigenstate | ϕ µ ( α, h + δh ) (cid:105) at α , h + δh providing the maximum squared overlap with | ϕ µ ( α, h ) (cid:105) .(This procedure was already used in the context of peri-odically driven dynamics with Floquet states in order tofind their adiabatic continuation [90, 91].) The excited-state fidelity susceptibility is then defined as the quantity χ µ ( α, h ) ≡ − h log | (cid:104) ϕ µ ( α, h ) | ϕ µ ( α, h + δh ) (cid:105) | . (18)We easily see that the larger χ µ ( α, h ), the more the two states are different. The closer χ µ ( α, h ) to 0, the moresimilar are the two states.Considering the behaviour of the excited-state fidelitysusceptibility, first of all we notice a clear discontinuityin its behaviour when α moves from 0 to a value largerthan 0. We can see this fact in Fig. 11 where we see amarked and generalized decrease moving from α = 0 to α = 0 .
05. The point is that for α > δh . Then the eigenstates for h and those for h + δh become more similar to each other than in the α = 0 case and the excited-state fidelity susceptibilitydecreases.For α larger than some value, the multiplets clearlyseen in Fig. 11 disappear. We consider an example ofthat in Fig. 12 (b) where we consider α = 0 . (cid:104)· · · (cid:105) Shell ) in order toreduce the noise. Differently from the energy bins dis-cussed above, the energy shells are defined dividing theenergy-spectrum width in N Shell equal intervals. Each in-terval is a shell and we associate to it E Shell , the energyat its middle point. We rescale the fidelity susceptibilityby N , (cid:104) χ ( α, h ) (cid:105) Shell /N . In Fig. 12 (b) we see a quite2 (a)(b) (c)(d)(e)(f) (g)(h)FIG. 9. Scatter plots of G µ versus E µ for different values of the parameters. We consider N = 20. different behaviour at small and large energy densities E Shell /N . In the first regime, (cid:104) χ ( α, h ) (cid:105) Shell /N is verynear to 0 and essentially stays there as the system sizeincreases. In the second regime, it increases with thesystem size. The increase seems linear, as appropriatefor a fully chaotic system where the eigenstates behaveas random states and totally change moving from onevalue of h to a very nearby one. The boundary betweenthe two regimes seems to slightly shift left as the sys-tem size is increased, so that the ergodicity slowly ex-tends towards smaller energy densities, although largersystem sizes would be required in order to quantitativelyextrapolate to the thermodynamic limit. In the lower in-set we plot δ Shell χ ( α, h ) /N , the standard deviation of theexcited-state fidelity susceptibility over the energy shell,versus the corresponding energy density. We see that thestandard deviation is large in the ergodic regime. Notethat while for local for local observables fluctuations areexpected to be weak according to ETH, it is exactly oppo- site for the non-local excited-state fidelity susceptibility.We plot for comparison also the case α = 0 in Fig. 12 (a).Here, the situation is different due to the multiplets withno clear tendency as compared to the α = 0 . (cid:104) χ ( α, h ) (cid:105) Shell versus α and energydensity E Shell /N . We consider only α > α = 0 point. The presence of anergodic energy band is clearly seen up to α = 2. Com-paring the plot in Fig. 13 (a) with the density-of-statesplots in Fig. 5, we see that there is no relation betweenthe behaviour of the excited-state fidelity susceptibilityand the persistence of the multiplet structure.It is important to compare these results withthe energy-bin averaged level spacing ratio (cid:104) r (cid:105) bin (Sec. III A 1). Especially for h = 0 . (cid:104) χ ( α, h ) (cid:105) Shell /N [Fig. 2 (b)] is quite different from (cid:104) r (cid:105) bin [Fig. 13 (b)]: the level spacing ratio is near to Wigner-3 M ( G ) N S h = 0.1 α = 0.05 α = 0.25 α = 0.5 α = 0.75 α = 1.0 α = 1.51/ N S1/2 (a) M ( G ) N S h = 0.1 α = 2.0 α = 2.25 α = 3.0 α = 3.0 α ➝ ∞ N S1/2 (b) M ( G ) N S h = 0.5 α = 0.05 α = 0.5 α = 1.0 α = 2.0 α = 4.0 α = 8.0 α ➝ ∞ N S1/2 (c)FIG. 10. Double-logarithmic plot of M ( G ) versus N S fordifferent values of α and h . (Here the corresponding valuesof N go from N = 14 to the extreme left to N = 20 tothe extreme right.) For comparison, we plot also the Isingintegrable α → ∞ nearest-neighbour model case and the ETHtrend ∼ / N / S of Eq. (11). Dyson for all the considered energy densities. Considerthat the bins at low energy densities are much largerthan the shells due to the very small density of states, sothe lowest-energy bin corresponds to many energy shells.Moreover, in the lowest-energy bin the predominant roleis played by the states at larger energy density whichare larger in number. So the level spacing ratio averagedover bins has a smaller resolution at small energy densi-ties. Notwithstanding this fact, (cid:104) r (cid:105) bin shows ergodicity FIG. 11. Some excited-state fidelity susceptibilities for N =20, h = 0 . δh = 0 . up to energy densities smaller than what the excited-statefidelity susceptibility suggests. Moreover, the extensionin α of the large-fidelity-susceptibility region for h = 0 . α where r takes the Wigner-Dyson value Fig. 1 (b). We have shownthat the latter interval extends everywhere in the ther-modynamic limit and it is quite possible that also theformer does. For h = 0 .
1, in contrast, the situation ismore complex: In the interval α ∈ [0 . ,
2] the behaviourof (cid:104) r (cid:105) bin [Fig. 13 (a)] follows quite closely the behaviourof the excited-state fidelity susceptibility [Fig. 2 (a)]. So,in this case we can quite confidently state that there isa low-energy density ergodicity breaking at finite systemsize, strictly related to the small value of r around α = 2in Fig. 1 (a). IV. SPECTRAL PAIRING
It is well known that the long-range quantum Isingchain can exhibit a symmetry-breaking transition atnonzero temperature as soon as α <
2. [92] The corre-sponding microcanonical or even single-eigenstate prop-erties have, however, not been explored extensively, ex-cept the notable Ref. [64] which focuses on low energiesand α = 1 .
5. In the following we aim to fill this gapby systematically studying the long-range order of theeigenstates which gives rise to Z symmetry breakingin the thermodynamic limit. In particular, we want toquantify whether for α (cid:54) = 0 there are states with long-range order at finite excitation energy density and toestimate the critical energy density e ∗ below which theeigenstates break the symmetry in the thermodynamiclimit ( e ∗ is called broken symmetry edge [32]). The exis-tence of the broken-symmetry edge is well known for thecase α = 0 [32], h <
1, but it is not explored in detail for α (cid:54) = 0.In order to perform this analysis, we need to ac-cess both of the two Z symmetry sectors. Therefore,4 < χ ( α , h ) > S he ll / N E Shell /N α = 0, h = 0.1, δ h = 0.01 N = 20N = 18N = 16 (a) < χ ( α , h ) > S he ll / N E Shell /N α = 0.5, h = 0.1, δ h = 0.01 N = 20 N = 18 N = 16 δ S he ll χ ( α , h ) (b)FIG. 12. Some excited-state fidelity susceptibilities averagedover the energy shell versus the corresponding energy densityfor different system sizes and h = 0 . δh = 0 . α = 0[panel (a)] and α = 0 . α = 0 . N Shell = 50. we restrict to the subspace corresponding to the zero-momentum sector and even only with respect to inver-sion. We target the single eigenstates and study theenergy gaps between nearby states: If there is symme-try breaking in the thermodynamic limit, the eigenstatesmust appear in quasidegenerate doublets, which becomedegenerate in the thermodynamic limit (the splitting isexponentially small in the system size). We make useof this property to determine the broken-symmetry edge.For that purpose we consider the splitting inside pairs ofnearby eigenenergies∆ (1) n = E n − E n − , (19)( n is an integer number labeling the eigenvalues in in-creasing order) and the gap between nearby pairs, eval-uated as the difference of next-nearest neighbor eigenen-ergies ∆ (2) n = E n +1 − E n − . (20)The rationale is therefore to understand up to which en-ergy the spectrum is organized in quasidegenerate dou-blets (with energies inside the doublet E n − and E n ). If < χ ( α , h = 0.1)> Shell
0 0.5 1 1.5 2 2.5 3 α -1-0.8-0.6-0.4-0.2 0 0.2 0.4 0.6 0.8 E S he ll / N (a) < χ ( α , h = 0.5)> Shell
0 0.5 1 1.5 2 2.5 3 α -1.2-1-0.8-0.6-0.4-0.2 0 0.2 0.4 0.6 0.8 E S he ll / N (b)FIG. 13. (cid:104) χ ( α, h ) (cid:105) Shell /N versus α and E Shell /N for h = 0 . h = 0 . N Shell = 50, N = 20. we are in presence of such a doublet, ∆ (1) n should be muchsmaller than ∆ (2) n . In particular, the ratio ∆ (1) n / ∆ (2) n should scale to 0 with the system size. In order to com-pare the two ∆, it is convenient to average ∆ (1) n and ∆ (2) n on energy shells, to reduce the nearby-eigenstate fluctu-ations. We define the N Shell energy shells, and associateto each of them an energy value E Shell , as we have donein Sec. III D. We consider the averages over the energyshells (cid:104) ∆ (1) n (cid:105) Shell ( E Shell ) and (cid:104) ∆ (2) n (cid:105) Shell ( E Shell ) and theirratio D ( E Shell ) = (cid:104) ∆ (1) n (cid:105) Shell ( E Shell ) (cid:104) ∆ (2) n (cid:105) Shell ( E Shell ) (21)We term D ( E Shell ) as the relative splitting and plot it ver-sus E Shell /N for different system sizes in Fig. 14. We con-sider h = 0 . α , α = 0 .
05 [Fig. 14. (a)]and α = 0 . α the spectrum is organized in multiplets for the systemsizes we have access to, while for the second it does not.For α = 0 . N clearly cross: There is a value of E Shell /N below which D ( E Shell ) decreases with the system size and above whichincreases. This is exactly what one would expect for abroken-symmetry edge, and we take this crossing pointas an estimate for the broken symmetry edge, with anerrorbar given by the mesh in E Shell .In contrast to the α = 0 . α = 0 .
05 we donot see any crossing as smooth as this one [Fig. 14 (a)].The point is that for this value of α and these system5sizes, the dynamics is strongly affected by the multipletswe have discussed above. The multiplets give rise to thenoisy behaviour which appears in Fig. 14 (a) and doesnot allow us to clearly give an estimate for e ∗ . We willprovide an estimate for the broken symmetry edge onlyfor those values of α and N where we do not see a noisymultiplet structure in the crossing region.We plot the resulting e ∗ versus α in Fig. 15 for h = 0 . h = 0 .
5, obtained considering the crossing of thecurves for N = 20 and N = 22. For α = 0 we take thetheoretical value e ∗ = − h found in [32]. For h = 0 . e ∗ stays constant inside the errorbars up to α = 0 .
6. We can reliably estimate e ∗ with our methodup to α = 1 .
5. Above that value results are unclear andlarger system sizes are needed even to clearly state if abroken symmetry edge exists. Nevertheless, consideringthat the ground-state is at e GS (cid:39) − α ≤ . Z symmetry breaking at finite excitation energy densi-ties. So, there is a finite fraction of the energy-spectrumwidth where the eigenstates show long-range order, sim-ilarly to the α = 0 and the disordered case. This is inagreement with the findings of [47, 55], where the long-time dynamics supports long-time magnetization in therange of α we can study and beyond. -12 -10 -8 -6 D ( E S he ll ) E Shell /N α =0.05 N = 18N = 20N = 22 (a) -12 -10 -8 -6 D ( E S he ll ) E Shell /N α =0.5 N = 18N = 20N = 22 (b)FIG. 14. D ( E Shell ) versus E Shell /N for α = 0 .
05 [panel (a)]and α = 0 . h = 0 . N Shell = 50. -0.9-0.8-0.7-0.6-0.5-0.4-0.3-0.2-0.1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 e * α h = 0.1h = 0.5 FIG. 15. e ∗ versus α for two values of h <
1. For α (cid:54) = 0 e ∗ has been estimated from the crossing between the curves for N = 20 and N = 22 of D ( E Shell ) versus E Shell /N . N Shell =50. For α = 0 we use the theoretical value e ∗ = − h providedby [32]. V. CONCLUSION
In conclusion we have considered the long-range Isingmodel with power-law interactions and used exact diago-nalization to study the ergodicity and the Z symmetry-breaking properties of the Hamiltonian eigenstates at fi-nite energy densities. Through an analysis of the averagelevel spacing ratio, the properties of the half-chain entan-glement entropy, and the longitudinal nearest-neighbourcorrelation operator in the Hamiltonian eigenstates, wehave found evidence that the system shows ergodicityfor all α >
0, as soon as the full permutation symmetrypresent at α = 0 is broken. For α = 0 the Hilbert spacedecomposes into many identical subspaces with the sameenergy levels. We argue that even an infinitesimal α > α = 0 integrablepoint nevertheless plays an important role in the regime α (cid:28) α = 0 the spectrum is organized in degenerate multipletsdue to the full permutation symmetry, which survive atlow and intermediate energies for a finite system size at0 < α (cid:28) (cid:104) r (cid:105) bin suggest, however, that even-tually the full spectrum becomes ergodic, at least for h = 0 .
5. For h = 0 . (cid:104) r (cid:105) bin is significantly smaller than the ergodicWigner-Dyson value, in close correspondence with thesmall values of the susceptibility and the non-ETH be-haviour of the longitudinal nearest-neighbour correlation.It is, however, also clear that with such small transversefields we naturally operate in close proximity to an inte-grable point implying larger finite-size effects, it is overalllikely that in the thermodynamic limit ergodicity is re-stored also in this limit, but our numerics do not providea definitive proof.Our findings may be relevant for those quench-dynamics investigations which focus on probing lower-energy densities and show a regular-like behaviour forobservables [47, 55]. The results of the density-of-statesanalysis show low-energy few-states multiplets well sep-arated by gaps and persisting for increasing system size.Our exact diagonalizations show a persisting noner-godic behaviour for h = 0 . α around the value of α ≈
2. This is a very suggestive result because there areother long-range models with α = 2 which are integrable,but the system sizes we have access to do not allow tostate if this effect persists in the thermodynamic limit.Nevertheless, a nonintegrable behaviour for N = 22 is al-ready remarkable and might suggest the proximity of anintegrable point. In all the other cases we see an ergodicbehaviour.Then we have moved to inquire the Z symmetry-breaking properties of the excited Hamiltonian eigen-states, in connection with the low-energy quench dynam-ics. We have found a energy-density threshold e ∗ (brokensymmetry edge) below which the eigenstates show long-range order. We can clearly find it up to α = 1 . − α = 1 .
5. Thisimplies that a finite fraction of the energy spectrum pro-vides eigenstates which show long-range order. Thesestates are at low energy, and then a low-energy quenchof an initial state polarized along z should give rise toa macroscopic magnetization persisting for long timesin the thermodynamic limit, in line with the findingsof [47, 55]. In our analysis, we exploit that eigenstatesbreaking the Z symmetry in the thermodynamic limitcorrespond to eigenvalues organized in quasidegeneratedoublets with splitting exponentially small in the systemsize. Below e ∗ the doublet splitting decreases relativelyto the next-nearest-neighbor eigenenergies gap, above e ∗ the relative splitting increases.Perspectives of future work include the study of this model at thermal equilibrium in the canonical ensemble,in order to understand the threshold value of α belowwhich the broken symmetry edge and its canonical coun-terpart start differing from each other. For α = 0 isquite easy to see that they are different using the re-sults of [32] and [93]; this is a particular case of thegeneral fact that in long-range system microcanonicaland canonical ensemble give rise to different results [29].Another possibility would be to inquire ergodicity underthe breaking of the full permutation symmetry puttingthe long-range α > α = 0 infinite-range sys-tems, like the infinite-range Bose-Hubbard [33, 94], thegeneralized Jaynes-Cummings model [33], and the quan-tum clock models with infinite-range interactions [95].We will also use perturbation theory in α to understandthe spectral properties at small α , as the results for theHilbert-Schmidt distance suggest to do and explore howour results are related to the prethermalization dynamicsobserved in [57]. ACKNOWLEDGMENTS
We thank M. Dalmonte, R, Khasseh, S. Pappalardi,F. M. Surace and especially R. Fazio for fruitful discus-sions. A. R. warmly thanks D. Rossini and A. Tomadinfor the access to the GOLDRAKE cluster where partof the numerical work for this project was performed.This project has received funding from the EuropeanResearch Council (ERC) under the European Union’sHorizon 2020 research and innovation programme (grantagreement No. 853443), and M. H. further acknowledgessupport by the Deutsche Forschungsgemeinschaft (DFG)via the Gottfried Wilhelm Leibniz Prize program.
Appendix A: Hilbert-Schmidt distance frominfinite-range model
Considering the Hamiltonian Eq. (1), we explicitlywrite its dependence on α and h , ˆ H = ˆ H ( α, h ). We wantto quantify the Hilbert-Schmidt distance of ˆ H ( α, h ) fromits infinite-range α = 0 counterpart ˆ H (0 , h ). We definethe distance as d ( α, N ) = (cid:107) ∆ H ( α, N ) (cid:107) HS = (cid:115) Tr (cid:20)(cid:16) ∆ ˆ H ( α, N ) (cid:17) (cid:21) , (A1)with ∆ ˆ H ( α, N ) ≡ ˆ H ( α, h ) − ˆ H (0 , h ) independent of h .Note that for an Hermitian operator ˆ O with eigenvalues λ j , (cid:107) O (cid:107) HS = (cid:113)(cid:80) j λ j .To compute d , we write∆ ˆ H ( α, L ) = N (cid:88) i,j, i (cid:54) = j J (cid:48) i,j ( α ) σ zi σ zj , (A2)7 FIG. 16. d ( α, N ) versus α for different values of N . Noticethe linear increase with α . where J (cid:48) i,j = 1 N ( α ) D αi, j − N (0) . (A3) Then (cid:104) ∆ ˆ H ( α, N ) (cid:105) = N (cid:88) i,j, i (cid:54) = j (cid:2) J (cid:48) i,j ( α ) (cid:3) + N (cid:88) distinct i,j,k ( · · · )ˆ σ zi ˆ σ zj + N (cid:88) distinct i,j,k,l ( · · · )ˆ σ zi ˆ σ zj ˆ σ zµ ˆ σ zl . (A4)Taking the trace, all term but the first one vanish, sothat d ( α, N ) = (cid:118)(cid:117)(cid:117)(cid:116) N (cid:88) i,j, i (cid:54) = j (cid:2) J (cid:48) i,j ( α ) (cid:3) . (A5)We numerically compute this quantity for various valuesof N and report it versus α in Fig. 16. We clearly seethat it increases linearly in α for small α . [1] A. Lichtenberg and M. Lieberman, Regular and ChaoticMotion (Springer, 1992).[2] A. Vulpiani, M. Falcioni, and P. Castiglione,
Chaos andCoarse Graining in Statistical Mechanics (CambridgeUniversity Press, 2008).[3] M. V. Berry, in
Topics in Nonlinear Mechanics , Vol. 46,edited by S. Jorna (Am.Inst.Ph., 1978) pp. 16–120.[4] J. M. Deutsch, Phys. Rev. A , 2046 (1991).[5] M. Srednicki, Phys. Rev. E , 888 (1994).[6] M. Rigol, V. Dunjko, and M. Olshanii, Nature , 854(2008).[7] T. c. v. Prosen, Phys. Rev. E , 3949 (1999).[8] M. V. Berry, Journal of Physics A: Mathematical andGeneral , 2083 (1977).[9] P. Pechukas, Phys. Rev. Lett. , 943 (1983).[10] M. Feingold and A. Peres, Phys. Rev. A , 591 (1986).[11] T. Prosen, Annals of Physics , 115 (1994).[12] O. Bohigas, M. J. Giannoni, and C. Schmit, Phys. Rev.Lett. , 1 (1984).[13] B. Eckhardt and J. Main, Phys. Rev. Lett. , 2300(1995).[14] M. Haque, P. A. McClarty, and I. M. Khaymovich,“Entanglement of mid-spectrum eigenstates of chaoticmany-body systems – deviation from random ensem-bles,” (2020), arXiv:2008.12782 [cond-mat.stat-mech].[15] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev.Mod. Phys. , 021001 (2019).[16] E. B. Rozenbaum and V. Galitski, Physical Review B ,064303 (2017).[17] S. Notarnicola, F. Iemini, D. Rossini, R. Fazio, A. Silva,and A. Russomanno, Phys. Rev. E , 022202 (2018).[18] C. Rylands, E. Rozenbaum, V. Galitski, and R. Konik,Phys. Rev. Lett. , 155302 (2020). [19] M. Fava, R. Fazio, and A. Russomanno, Phys. Rev. B , 064302 (2020).[20] A. Russomanno, M. Fava, and R. Fazio, Physical ReviewB (2020), 10.1103/physrevb.102.144302.[21] T. c. v. Prosen, Phys. Rev. Lett. , 1808 (1998).[22] G. Biroli, C. Kollath, and A. M. L¨auchli, Phys. Rev.Lett. , 250401 (2010).[23] C. Kollath, A. M. L¨auchli, and E. Altman, Physicalreview letters , 180601 (2007).[24] M. Antoni and S. Ruffo, Phys. Rev. E , 2361 (1995).[25] M.-C. Firpo, Phys. Rev. E , 6599 (1998).[26] V. Latora, A. Rapisarda, and S. Ruffo, Phys. Rev. Lett. , 692 (1998).[27] C. Anteneodo and C. Tsallis, Phys. Rev. Lett. , 5313(1998).[28] T. M. R. Filho, A. E. Santana, M. A. Amato, andA. Figueiredo, Phys. Rev. E , 032133 (2014).[29] A. Campa, T. Dauxois, D. Fanelli, and S. Ruffo, Physicsof Long-Range Interacting Systems (Oxford, 2014).[30] R. Khasseh, R. Fazio, S. Ruffo, and A. Russomanno,Phys. Rev. Lett. , 184301 (2019).[31] H. Lipkin, N. Meshkov, and A. Glick, Nuclear Physics , 188 (1965).[32] G. Mazza and M. Fabrizio, Phys. Rev. B , 184303(2012).[33] B. Sciolla and G. Biroli, J. Stat. Mech.: Theor. and Ex-per. , P11003 (2011).[34] Z. N. C. Ha and F. D. M. Haldane, Phys. Rev. B ,12459 (1993).[35] F. D. M. Haldane, in Correlation Effects in Low-Dimensional Electron Systems , edited by A. Okiji andN. Kawakami (Springer Berlin Heidelberg, Berlin, Hei-delberg, 1994) pp. 3–20. [36] D. Uglov, “The trigonometric counterpart of the haldaneshastry model,” (1995), arXiv:hep-th/9508145 [hep-th].[37] I. Sechin and A. Zotov, Physics Letters B , 1 (2018).[38] J. Lamers, Phys. Rev. B , 214416 (2018).[39] P. Hauke and M. Heyl, Phys. Rev. B , 134204 (2015),arXiv:1410.1491 [cond-mat.dis-nn].[40] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess,P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, NaturePhysics , 907–911 (2016).[41] A. Burin, Annalen der Physik (2017),10.1002/andp.201600292.[42] S. Roy and D. E. Logan, SciPost Physics (2019),10.21468/scipostphys.7.4.042.[43] N. Y. Yao, C. R. Laumann, S. Gopalakrishnan, M. Knap,M. M¨uller, E. A. Demler, and M. D. Lukin, Phys. Rev.Lett. , 243002 (2014).[44] K. S. Tikhonov and A. D. Mirlin, , 214205 (2018).[45] A. L. Burin, Phys. Rev. B , 094202 (2015).[46] G. Piccitto and A. Silva, Journal of Statistical Mechanics:Theory and Experiment , 094017 (2019).[47] B. ˇZunkoviˇc, M. Heyl, M. Knap, and A. Silva, Phys.Rev. Lett. , 130601 (2018).[48] A. Lerose, B. ˇZunkoviˇc, J. Marino, A. Gambassi, andA. Silva, Phys. Rev. B , 045128 (2019).[49] S. Pappalardi, A. Russomanno, B. ˇZunkoviˇc, F. Iemini,A. Silva, and R. Fazio, Phys. Rev. B , 134303 (2018).[50] J. C. Halimeh, V. Zauner-Stauber, I. P. McCulloch, I. deVega, U. Schollw¨ock, and M. Kastner, Phys. Rev. B ,024302 (2017), arXiv:1610.01468 [cond-mat.quant-gas].[51] F. Liu, R. Lundgren, P. Titum, G. Pagano, J. Zhang,C. Monroe, and A. V. Gorshkov, Phys. Rev. Lett. ,150601 (2019).[52] A. Y. Guo, M. C. Tran, A. M. Childs, A. V.Gorshkov, and Z.-X. Gong, “Signaling and scram-bling with strongly long-range interactions,” (2019),arXiv:1906.02662 [quant-ph].[53] A. Lerose, B. ˇZunkoviˇc, A. Silva, and A. Gambassi, Phys.Rev. B , 121112 (2019).[54] R. Verdel, F. Liu, S. Whitsitt, A. V. Gorshkov, andM. Heyl, “Real-time dynamics of string breaking inquantum spin chains,” (2019), arXiv:1911.11382 [cond-mat.stat-mech].[55] R. Khasseh, A. Russomanno, M. Schmitt, M. Heyl, andR. Fazio, Phys. Rev. B , 014303 (2020).[56] D. J. Luitz and Y. Bar Lev, Phys. Rev. A , 010105(2019).[57] B. Neyenhuis, J. Smith, A. C. Lee, J. Zhang, P. Richerme,P. W. Hess, Z. X. Gong, A. V. Gorshkov, and C. Monroe,“Observation of prethermalization in long-range interact-ing spin chains,” (2016), arXiv:1608.00681 [quant-ph].[58] W. L. Tan, P. Becker, F. Liu, G. Pagano, K. S.Collins, A. De, L. Feng, H. B. Kaplan, A. Kyprianidis,R. Lundgren, W. Morong, S. Whitsitt, A. V. Gorshkov,and C. Monroe, “Observation of domain wall confine-ment and dynamics in a quantum simulator,” (2019),arXiv:1912.11117 [quant-ph].[59] C. Monroe, W. C. Campbell, L. M. Duan, Z. X. Gong,A. V. Gorshkov, P. Hess, R. Islam, K. Kim, N. Linke,G. Pagano, P. Richerme, C. Senko, and N. Y. Yao, “Pro-grammable quantum simulations of spin systems withtrapped ions,” (2020), arXiv:1912.07845 [quant-ph].[60] J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis,P. Becker, H. Kaplan, A. V. Gorshkov, Z.-X. Gong, andC. Monroe, Nature , 601–604 (2017). [61] P. Jurcevic, H. Shen, P. Hauke, C. Maier, T. Brydges,C. Hempel, B. Lanyon, M. Heyl, R. Blatt, and C. Roos,Physical Review Letters (2017), 10.1103/phys-revlett.119.080501.[62] T. Brydges, A. Elben, P. Jurcevic, B. Vermersch,C. Maier, B. P. Lanyon, P. Zoller, R. Blatt, and C. F.Roos, Science , 260–263 (2019).[63] B. P. Lanyon, C. Maier, M. Holz¨apfel, T. Baumgratz,C. Hempel, P. Jurcevic, I. Dhand, A. S. Buyskikh, A. J.Daley, M. Cramer, and et al., Nature Physics ,1158–1162 (2017).[64] K. R. Fratus and M. Srednicki, “Eigenstate thermaliza-tion and spontaneous symmetry breaking in the one-dimensional transverse-field ising model with power-lawinteractions,” (2017), arXiv:1611.03992 [cond-mat.stat-mech].[65] B. Bertini, F. H. L. Essler, S. Groha, and N. J. Robinson,Phys. Rev. Lett. , 180601 (2015).[66] B. Bertini, F. H. L. Essler, S. Groha, and N. J. Robinson,Phys. Rev. B , 245117 (2016).[67] M. Kac, J. Math. Phys. (NY) , 216 (1963).[68] S. Sachdev, Quantum Phase Transitions (CambridgeUniversity Press, 2011).[69] G. B. Mbeng, A. Russomanno, and G. E. Santoro,“The quantum ising chain for beginners,” (2020),arXiv:2009.09208 [quant-ph].[70] D. J. Luitz, Phys. Rev. B , 134201 (2016).[71] A. Pal and D. A. Huse, Phys. Rev. B , 174411 (2010).[72] M. Srednicki, Journal of Physics A: Mathematical andGeneral , 1163 (1999).[73] V. Khemani, C. R. Laumann, and A. Chandran, Phys.Rev. B , 161101 (2019).[74] I. Mondragon-Shem, M. G. Vavilov, and I. Martin, “Thefate of quantum many-body scars in the presence of dis-order,” (2020), arXiv:2010.10535 [cond-mat.quant-gas].[75] M. Srednicki, Journal of Physics A: Mathematical andGeneral , L75 (1996).[76] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn,and Z. Papi´c, Phys. Rev. B , 155134 (2018).[77] Y. Huang, Nuclear Physics B , 594 (2019).[78] F. Haake, Quantum Signatures of Chaos (Springer-Verlag New York, Inc., Secaucus, NJ, USA, 2006)Chap. 7, pp. 263–274.[79] P. J. Coles, M. Cerezo, and L. Cincio, Phys. Rev. A ,022103 (2019).[80] P. Pandya, O. Sakarya, and M. Wie´sniak, Phys. Rev. A , 012409 (2020).[81] M. V. Berry and M. Tabor, Proc. Roy. Soc. A , 375(1977).[82] E. Lieb, T. Schultz, and D. Mattis, Annals of Physics , 407 (1961).[83] P. Pfeuty, Annals of Physics , 79 (1970).[84] D. J. Luitz, N. Laflorencie, and F. Alet, Physical ReviewB , 081103 (2015).[85] We numerically evaluate the Page value by constructing afully random state | ψ (cid:105) in H S , evaluating the correspond-ing S N/ , and then averaging over N rand randomness real-izations. Calling | S (cid:105) the symmetrized spin configurations(which are a basis of H S ), we impose (cid:104) S | ψ (cid:105) = √ N S e − iφ S ,where φ S is a random variable uniformly distributed in[0 , π ]. The procedure is very similar to the one usedin [19, 20]. In the plots we have used N rand = 2000.[86] D. N. Page, Phys. Rev. Lett. , 1291 (1993). [87] M. Nielsen and I. L. Chuang, Quantum Computationand Quantum Information (Cambridge University Press,2000).[88] P. Zanardi and N. Paunkovi´c, Phys. Rev. E , 031123(2006).[89] P. Sierant, A. Maksymov, M. Ku´s, and J. Za-krzewski, Physical Review E (2019), 10.1103/phys-reve.99.050102.[90] A. Russomanno, B.-e. Friedman, and E. G. Dalla Torre,Phys. Rev. B , 045422 (2017).[91] A. Russomanno and E. G. D. Torre, EPL (EurophysicsLetters) , 30006 (2016). [92] A. Dutta and J. K. Bhattacharjee, Phys. Rev. B ,184106 (2001).[93] J. Wilms, J. Vidal, F. Verstraete, and S. Dusuel, Journalof Statistical Mechanics: Theory and Experiment ,P01023 (2012).[94] B. Sciolla and G. Biroli, Phys. Rev. Lett. , 220401(2010).[95] A. Offei-Danso, F. M. Surace, F. Iemini, A. Russomanno,and R. Fazio, Journal of Statistical Mechanics: Theoryand Experiment2020