Off-diagonal matrix elements of local operators in many-body quantum systems
OOff-diagonal matrix elements of local operators in many-body quantum systems
Wouter Beugeling, Roderich Moessner, and Masudul Haque
Max-Planck-Institut f¨ur Physik komplexer Systeme, N¨othnitzer Straße 38, 01187 Dresden, Germany (Dated: September 13, 2018)In the time evolution of isolated quantum systems out of equilibrium, local observables generally relax to along-time asymptotic value, governed by the expectation values (diagonal matrix elements) of the correspondingoperator in the eigenstates of the system. The temporal fluctuations around this value, response to furtherperturbations, and the relaxation toward this asymptotic value, are all determined by the off-diagonal matrixelements. Motivated by this non-equilibrium role, we present generic statistical properties of off-diagonal matrixelements of local observables in two families of interacting many-body systems with local interactions. Sinceintegrability (or lack thereof) is an important ingredient in the relaxation process, we analyze models that canbe continuously tuned to integrability. We show that, for generic non-integrable systems, the distribution of off-diagonal matrix elements is a gaussian centered at zero. As one approaches integrability, the peak around zerobecomes sharper, so that the distribution is approximately a combination of two gaussians. We characterize theproximity to integrability through the deviation of this distribution from a gaussian shape. We also determinethe scaling dependence on system size of the average magnitude of off-diagonal matrix elements.
PACS numbers: 05.30.-d,05.70.Ln,75.10.Pq
I. INTRODUCTION
The topic of non-equilibrium dynamics of thermally iso-lated quantum systems has enjoyed a resurgence of interest,partly because of experimental progress with cold atoms. Anisolated system has no relaxation mechanism toward the low-lying parts of the many-body spectrum. As a result, the prop-erties of eigenstates far from the edges of the spectrum maybe more important for a non-equilibrium experiment than thelow-energy parts of the spectrum, which is the traditional fo-cus of interest of many-body quantum theory.A key question in the non-equilibrium dynamics of isolatedquantum systems is the thermalization or relaxation of a sys-tem prepared far out of equilibrium and subject to a time-independent Hamiltonian. The value (if any) to which localobservables relax is determined by the diagonal matrix ele-ments A αα = (cid:104) ψ α | ˆ A | ψ α (cid:105) of the corresponding operator ˆ A inthe eigenstates | ψ α (cid:105) . The eigenstate thermalization hypothesis(ETH) [1–5] proposes that the mechanism for the thermaliza-tion of non-integrable (“chaotic”) systems is the smoothnessof A αα as a function of eigenenergies E α . Accordingly, diag-onal matrix elements of local operators have been the subjectof several studies [3, 6–16].Off-diagonal matrix elements, A αβ = (cid:104) ψ α | ˆ A | ψ β (cid:105) , pro-vide further information about the time evolution (cid:104) A (cid:105) ( t ) ofobservables. In any finite system initially prepared in a com-bination of many eigenstates, there will be residual tempo-ral fluctuations around the long-time average. These tempo-ral fluctuations have been the subject of several recent studies,both numerical [3, 17–22] and analytical [23–25]. The magni-tude of these fluctuations is determined by | A αβ | , weighted,of course, by the weights of the eigenstates in the non-equilibrium initial state. Autocorrelation functions (unequal-time correlators), interesting on their own and appearing inthe formulation of fluctuation-dissipation relations in the “re-laxed” state a long time after a quench [26], also are givenin terms of | A αβ | . Finally, the details of the temporal ap-proach to the final relaxed value are also determined by the off-diagonal matrix elements of the corresponding operator[17, 27]. The approach toward the final value has been cal-culated in some model systems [28–32], although the connec-tion to off-diagonal matrix elements has not been explored indetail.The (statistical) properties of off-diagonal matrix elementsof local operators, A αβ , are thus related to a range of temporalproperties of contemporary interest. In this work, we providea statistical study of these objects. We use Hamiltonians thatcan be tuned between integrable limits, and provide scalinganalyses as a function of system size. We thus study whathappens to the distributions of A αβ as a function of distancefrom integrability, as well as how the thermodynamic limit isapproached.Some statistical aspects of off-diagonal matrix elements A αβ have appeared in Ref. [26] in the context of a non-equilibrium fluctuation-dissipation relation, and in Ref. [12].The aim of the present paper is to focus directly on the A αβ ina manner independent of quench protocol and provide a thor-ough study of their statistical properties.In the time evolution (cid:104) A (cid:105) ( t ) , each matrix element A αβ con-tributes with a frequency equal to the eigenvalue difference E β − E α [26, 33, 34]. In many quenches of physical interest,the initial occupancies are confined to a small energy win-dow (e.g., [3, 35, 36]), yet involve many eigenstates [37]. Wetherefore pay particular attention to the behavior of the typicalvalues of A αβ for small E β − E α . At large frequencies, the av-erage | A αβ | falls off fast, exponentially or super-exponentiallywith E β − E α .We pay special attention to the proximity to integrability,since it is well-appreciated that the relaxation behavior ofchaotic or generic systems is quite different from systems sub-ject to integrability [6, 20, 31, 38–49] or to (many-body) lo-calization [22, 34, 41, 50–53]. We identify signatures of the A αβ typical to the integrable, close-to-integrable, and nonin-tegrable cases. Close to integrability, we show that the matrix | A αβ | has a block-like or banded structure as a function of theenergy difference (frequency) E β − E α , which is visible as a r X i v : . [ c ond - m a t . s t a t - m ec h ] J a n oscillatory behavior in the frequency-dependence of average | A αβ | values.We show that the distribution of the matrix elements inany small frequency window is peaked around zero, havinga near-gaussian form for generic non-integrable systems (cf.Ref. [12]). At or near integrability, there is a stronger peakaround zero, i.e., the probability distribution is a mixture oftwo gaussian-like curves with unequal widths. This differ-ence appears to be a basic distinction between generic (non-integrable) and integrable systems. We demonstrate how theproximity to integrability can be quantitatively characterizedthrough the shape of the distribution of A αβ values, e.g.,through the size dependence of the kurtosis of this distribu-tion.We find that the scaling behavior of the average value of | A αβ | is D − in terms of the Hilbert-space dimension D .The values of | A αβ | at low frequencies tend to be larger thanfor the generic matrix elements, but the scaling follows D − as well. The scaling analysis is analogous to studies of the di-agonal matrix elements A αα and related quantities as a func-tion of system size, performed, e.g., in Refs. [10, 12, 14, 54–56]. As for the diagonal fluctuations [13, 14], we can con-struct plausibility arguments based on an assumption of quasi-randomness of the vector coefficients of the energy eigen-states. As such assumptions are difficult to prove rigorously,we emphasize, as in Ref. [14], that such arguments are in-herently heuristic and that extensive, multi-system, numericalanalysis is required to establish scaling laws; this paper pro-vides such data.The size dependence of A αβ ’s is related to the size depen-dence of the magnitude of the temporal fluctuations aroundthe long-time average [19–21]. The D − scaling is consistentwith the exponential dependence of the long-time fluctuationson the system size [21].This paper is structured as follows. In Sec. II, we introduceour models: the XXZ ladder and the Bose-Hubbard chain.In Sec. III, we introduce the frequency-resolved average ofthe off-diagonal matrix elements. In Sec. IV, we analyze thedistribution of values of A αβ , characterizing how a mixeddistribution (with two components having different widths)emerges close to integrability. Sec. V provides a scaling anal-ysis of the size-dependence of the average values of | A αβ | ,focusing on the low-frequency matrix elements. Sections III,IV, and V show results for the XXZ ladder. We support thegenerality of these results by presenting corresponding datafor the Bose-Hubbard chain in Sec. VI. In the appendices,we provide details of the relationship between time evolution (cid:104) A (cid:105) ( t ) and the matrix elements A αβ , and about our quantifi-cation of the non-gaussian distributions. II. MODELS AND OBSERVABLES
We use two families of Hamiltonians, each of which can betuned to integrable points. Both have been used in our previ-ous work on diagonal matrix elements [14]. Because we areinterested in generic properties of matrix elements, we takecare to avoid spurious symmetries in our model systems. In the spirit of many thermalization studies using spin mod-els [12, 31, 34, 35, 46, 55, 57, 58], our first tunable model willbe the spin- Heisenberg XXZ ladder with the geometry in-troduced in Ref. [14]. One ladder leg has an extra site com-pared to the other. There are thus L = 2 p + 1 sites, with p rungs between the legs. This geometry avoids reflection sym-metries. We have nearest-neighbor Heisenberg couplings h i,j = (cid:0) S + i S − j + S − i S + j (cid:1) + ∆ S zi S zj , (1)with S ± i = S xi ± i S yi , where S µi ( µ = x, y, z ) are the spinoperators, and i, j denote the nearest-neighbor site pairs. Theanisotropy parameter ∆ is kept away from special values like and ± , in order to avoid SU (2) symmetry or special solv-able points; we use ∆ = 0 . . The Hamiltonian of the systemis H = H + λH , where H = p − (cid:88) i =1 h i,i +1 + p (cid:88) i = p h i,i +1 and H = p (cid:88) i =1 h i,i + p (2)are the intrachain (leg) and the interchain (rung) coupling, re-spectively. The rung coupling is multiplied by λ , which actsas a tuning parameter. The xy coupling along the ladder legssets the units of energy and frequency. For λ = 0 , the chainsare uncoupled and the model is integrable. For finite values of λ , the system is non-integrable. In the limit of large λ , wherethe rung couplings dominate, there is another integrable limit.The effect of varying λ on the fluctuations of diagonal matrixelements has been studied in detail in Ref. [14].The number N ↑ of up spins is a conserved quantity. Theanalysis can therefore be constrained to a fixed- N ↑ sector. Thedimension of the Hilbert space of the ( L, N ↑ ) sector is equal tothe binomial coefficient D = (cid:0) LN ↑ (cid:1) . In order to study scaling,we use a sequence of system sizes with almost constant fillingfraction. We present data for a sequence of systems with near-zero magnetization (near half filling), by choosing L = 2 p +1 and N ↑ = p for integer p .Discussion of thermalization generally concerns local ob-servables. We present data for S z and S z S zp +2 , which serveas representative examples of single-site and two-site opera-tors.The second tunable Hamiltonian is the Bose-Hubbardmodel, widely used in studies of thermalization [8, 59–66].We use the Bose-Hubbard Hamiltonian on an L -site chain,with an extra term at an edge site killing reflection symmetry,as in Ref. [14]: H BH = L − (cid:88) i =1 (cid:16) b † i b i +1 + b † i +1 b i (cid:17) + λ (cid:32)(cid:88) i b † i b † i b i b i + H ∆ (cid:33) (3)where b i is the creation operator at site i and H ∆ =∆ b † b † b b with ∆ = 0 . is a small perturbation to the interac-tion term at the first site. The system is integrable in the λ → and λ → ∞ limits, and nonintegrable for intermediate values.We show results for the sector of unit filling fraction, i.e., thenumber of bosons N b = L . This choice provides the samesequence of Hilbert-space sizes [for bosons, D = (cid:0) L + N b − N b (cid:1) ] ˆ A = S z S zp +2 (a) λ = 0 . − − E β − − E α (b) λ = 0 . − − − − E α (c) λ = 5 − − − −
10 0 10 E α . . . . | A αβ | FIG. 1. (Color online) Matrix structure of | A αβ | as function of E α and E β for the observable ˆ A = S z S zp +2 . The diagonal matrix elementsare ignored. The white bands near the edges are regions without eigenvalues. The dashed square indicates the central half of the energy range,i.e., (cid:2) E min + E max , E min + E max (cid:3) in each direction: this is the “bulk” of the spectrum on which we focus our analysis. The systemsize is ( L, N ↑ ) = (13 , ; the Hilbert space dimension is D = 1716 . The unit of energy is set by the xy coupling along the ladder legs. as the one given by L = 2 N ↑ + 1 for the XXZ ladder. Typicallocal observables in the study of this model include n i = b † i b i , b † i b i +1 + b † i +1 b i and n i n i +1 . III. FREQUENCY-RESOLVED AVERAGE MATRIXELEMENTS
In Fig. 1, we visualize through a density plot the structureof the matrix | A αβ | as a function of energies E α and E β , us-ing the rung correlator ˆ A = S z S zp +2 of the XXZ ladder asobservable. The diagonal matrix elements are not considered.The structure of darker bands parallel to the main diagonalsuggests that the magnitude of the | A αβ | depends roughly onthe difference E α − E β . Thus the energies ( E α , E β ) ratherthan the indices ( α, β ) are natural coordinates for this plot (cf.Refs. [17, 29]).To consider the | A αβ | from finite-size data as acontinuous function of frequency, we “smooth out” | A αβ | δ ( ω − ( E β − E α )) as a function of ω , by aver-aging the values of | A αβ | with E α − E β in the frequencywindow [ ω − ∆ ω, ω + ∆ ω ] , S A ( ω, ∆ ω ) ≡ N ω, ∆ ω (cid:88) α,βα (cid:54) = βE α − E β ∈ [ ω − ∆ ω,ω +∆ ω ] | A αβ | , (4)where ˜ N ω, ∆ ω is the number of state pairs satisfying E α − E β ∈ [ ω − ∆ ω, ω + ∆ ω ] . The frequency-window width ω is chosen such that the interval contains sufficientlymany pairs of states. We restrict ourselves to positive ω , since A αβ = A βα for hermitian observables. The quantity S A ( ω ) is closely related to fluctuations around the asymptotic valueto which (cid:104) A ( t ) (cid:105) relaxes a long time after a quantum quench(Appendix A). The quantity S A ( ω, ∆ ω ) is the standard devi-ation of the distribution formed by the A αβ in the frequencywindow. In the large-system limit, the number of states ˜ N ω, ∆ ω in thewindow can be approximated as ˜ N ω, ∆ ω ≈ ω ˜ ρ ( E ) , where ˜ ρ ( ω ) is the density of pairs , i.e., the density of values E α − E β .The density of pairs is defined as the autocorrelation integral ˜ ρ ( ω ) = (cid:90) ρ ( E ) ρ ( E − ω ) dE (5)of the density of states ρ ( E ) with itself. We note that thedensity of pairs does not show signatures of the level-spacingstatistics, because the density of states is considered on acoarser resolution than that of individual eigenvalues. The be-havior of ˜ ρ ( ω ) is shown in the bottom row of Figure 2 for thespin ladder system, for different values of the λ parameter.To distinguish frequency regimes, we define a “typical” fre-quency scale ω , as the root-mean-square of all possible fre-quencies: ω = var( E α − E β ) = 1 D (cid:88) α,β ( E α − E β ) = 2 var( E α ) . (6)In Fig. 2, the values ω are indicated by markers on the hori-zontal (frequency) axes.The frequency dependence of S A ( ω ) is shown in Fig. 2 forthe observables ˆ A = S z and ˆ A = S z S zp +2 in the XXZ laddermodel. In addition, through the shading in the top two panels,we indicate the distribution of the values of | A αβ | in each fre-quency window. The value of ∆ ω = 0 . used in this figureis a compromise between being sufficiently small to resolvethe details, and having sufficiently many state pairs for goodstatistics. In the cases of Fig. 2, the number of state pairs inthe window [ ω − ∆ ω, ω + ∆ ω ] is ∼ for ω (cid:46) ω .At high frequencies ( ω (cid:29) ω ), S A ( ω ) decreases as a func-tion of ω . The decrease is rapid; we have found this to begenerally exponential or super-exponential ( ∼ gaussian) with ω ; the details vary with the observable and appear to be non-universal. ˆ A = S z ˆ A = S z S zp +2 (a) λ = 0 . . . . . . S A , A α β . . . . S A , A α β ˜ ρ ( a r b . un it s ) ω S z S z S zp +2 (b) λ = 0 . ω S z S z S zp +2 (c) λ = 5 ω FIG. 2. (Color online) Frequency-resolved analysis of off-diagonalelements for the ( L, N ↑ ) = (13 , ladder. The shadings showthe distribution of the values of | A αβ | for each frequency window;darker shading indicates more occurrences of respective | A αβ | val-ues. The curves show S A as a function of the frequency ω . In thebottom panel, we indicate the density of pairs ˜ ρ for the central partof the spectrum (solid) and for the full spectrum (dashed), in arbi-trary units. The marker on the horizontal axis points to the typicalfrequency ω . The shading and solid curves are all results using thecentral part of the spectrum. The units of frequency are the same asthe units used for energy in Fig. 1. At medium frequencies, S A ( ω ) typically shows severalpeaks. The oscillatory behavior is more pronounced near in-tegrability, i.e., for small and large λ . We observe typicalsmall- λ behavior in Fig. 2(a): The quantity S A shows short-scale oscillations, while the density of pairs ˜ ρ ( ω ) is smooth.We conjecture that the oscillatory behavior in near-integrablesystems is due to the Hilbert space being decomposable intomany subspaces weakly coupled by the Hamiltonian. When-ever α and β are in different subspaces, A αβ ≈ .At large λ , the system splits into weakly coupled subspaceswhich are in addition separated in energy, as evidenced bythe block-like structure in Fig. 1(c). Thus, the peaks of S A are accompanied by those in the density of pairs ˜ ρ ( ω ) . Theblocks are separated by energy ∼ λ , which can be understoodfrom treating the system as uncoupled dimers in the λ → ∞ limit. These are also the approximate frequencies at whichpeaks can be seen in Fig. 2(c). IV. DISTRIBUTION OF OFF-DIAGONAL MATRIXELEMENTS
Having described the variance S A ( ω ) of the distribution ofthe values of A αβ , we now look at the full distribution. λ = 0 . λ = 0 . λ = 5 − . − .
05 0 .
00 0 .
05 0 . PDF E< . (a) − . . d e n s it y ( a r b . un it s ) − . − .
05 0 .
00 0 .
05 0 . PDF E< . (b) − . . − . − . − . − . . . . . . PDF E< . (c) − . . − . − .
05 0 .
00 0 .
05 0 . PDF E =0 . (d) − . . d e n s it y ( a r b . un it s ) A αβ − . − .
05 0 .
00 0 .
05 0 . PDF E =0 . (e) − . . A αβ − . − .
05 0 .
00 0 .
05 0 . PDF E =1 . (f) − . . A αβ FIG. 3. (Color online) Histograms (shaded area) of the off-diagonalelements A αβ with (a–c) E β − E α ∈ (0 , . ω ) and (d–f) E β − E α ∈ [0 . ω , . ω + δω ) , where δω = 0 . . The solid curve is agaussian fit and the dashed curve is a fit of a mixture of two gaussians,Eq. (8). The observable is ˆ A = S z S zp +2 and the system size is ( L, N ↑ ) = (13 , . The number of state pairs in these histogramcomputations ranges from to . A. Shapes of the distributions of A αβ In Fig. 2, we have shown using shading densities thefrequency-resolved distributions of values of | A αβ | . A fea-ture visible already in the density plots is that the distribu-tions are more strongly weighted near zero (near the horizon-tal axis) near integrability. This feature will be explored anddescribed in more detail below.In Fig. 3, we show the distributions of A αβ values, in twodifferent frequency windows. The top panels show the low-frequency regime (cutoff frequency ω max = 0 . ω ). Thebottom panels focus on a frequency window around . ω .Only the states in the central part of the spectrum (within thedashed square region in Fig. 1) are considered.The distributions are seen to be very nearly symmetricaround zero. Of course, the signs of individual A αβ val-ues are not meaningful since every eigenstate carries an ar-bitrary phase. However, from N (cid:29) eigenstates, one obtains N (cid:29) N matrix elements; so the overall shape of the distri-bution (roughly equal number of positive and negative values)cannot be altered by the choice of phases for the eigenstates.The solid curves are gaussian fits determined by the vari-ance of the A αβ , centered at . Far from integrability, this isseen to be a very good description. However, near integrabil-ity the distribution has a sharper peak than a gaussian, and ap-pears to be a mixture of two near-gaussian distributions withdifferent widths. This appears to be a fundamental distinctionbetween (near-)integrable and generic systems.We do not currently have a complete explanation for theextra peak in near-integrable systems, but we conjecture thefollowing mechanism which provides some intuition. In theintegrable case, there are many conserved quantities. The en-ergy eigenstates can be grouped into subspaces or symmetrysectors by the eigenvalues (“quantum numbers”) of the opera-tors corresponding to these conserved quantities. An approxi-mate version of this statement is true close to, but not at, inte-grability. An operator ˆ A corresponding to a local observable,when acting on an eigenstate | ψ α (cid:105) , changes the eigenstateonly locally, i.e., slightly. The resulting wavevector ˆ A | ψ α (cid:105) will thus be likely to have larger overlap with eigenstates hav-ing the same quantum numbers as | ψ α (cid:105) , and much smalleroverlaps with eigenstates having different quantum numbersfrom those of | ψ α (cid:105) . In other words, A αβ is close to zero when-ever α and β belong to different subspaces. Of course, | ψ α (cid:105) and | ψ β (cid:105) are orthogonal even if they belong to the same sub-space, so that the off-diagonal matrix element of a local op-erator is small anyway for large system sizes. The argumentis that, when they belong to different sectors, the states differadditionally by having different quantum numbers, not onlyby being orthogonal, and this should make the inter-subspacematrix elements statistically much smaller than intra-subspacematrix elements.This line of reasoning intuitively connects to the idea thatintegrablility makes a system “non-ergodic”. However, theargument is difficult to make rigorous. It is easy to con-struct special operators that connect different subspaces, e.g.,if (cid:80) j S zj is a conserved quantity in a spin Hamiltonian, thelocal operator S + j will connect different subspaces. How-ever, a generic operator is expected not to have such spe-cial relationships with many of the conserved quantities, sincemost conserved quantities in integrable lattice models haverather complicated form when expressed in terms of spa-tially local operators. Although the explanation providedby this “inter-subspace versus intra-subspace” perspective re-mains only heuristic at this stage, our data for multiple systemdemonstrates that near-integrable systems indeed have a sub-stantial number of extremely small matrix elements.In summary, numerical observations on the families of sys-tems (XXZ ladder, Bose-Hubbard chain) investigated in thiswork indicate that, in quasi-integrable cases, the studied localobservables tend to respect the symmetries that are dynami-cally conserved at exact integrability. We expect this behav-ior to be generic at (near-) integrability for local observablesin these types of models. Moreover, we conjecture that thisbehavior may also be generic in the class of integrable many-body systems at large, and for a wide class of local observ-ables.The gaussian shape of the distributions for generic non-integrable points can be explained heuristically by invokingthe central limit theorem. Writing c ( α ) γ ≡ (cid:104) φ γ | ψ α (cid:105) in terms ofthe eigenstates | ψ α (cid:105) of the Hamiltonian and | φ γ (cid:105) of A (witheigenvalues a γ ), we can write the matrix elements as A αβ = (cid:88) γ c ( α ) ∗ γ c ( β ) γ a γ . (7)For non-integrable systems, the summands c ( α ) ∗ γ c ( β ) γ a γ maybe expected to behave like quasi-independent random vari-ables. The central limit theorem then implies the gaussiandistribution of A αβ . As in Ref. [14], we stress that the ran-domness and independence of coefficients is a hypothesis anddifficult to prove rigorously. This is in the same spirit as argu-ments for scaling behaviors of diagonal matrix elements or ofinverse participation ratios based on similar randomness as- sumptions [13, 14]. The physical intuition for such random-ness assumptions is that an eigenstate in the middle of thespectrum of a generic system is so complex that the coeffi-cients behave as random and independent variables for manypurposes. B. Quantifying the distribution shapes
In order to characterize the nature of the distributions atsmall and large λ , we fit the numerically obtained histogramsto the sum of two gaussian distributions, defined as g ( A ) = an σ ( A ) + (1 − a ) n σ ( A ) , (8)where n σ i ( A ) is the gaussian distribution with variance σ i and zero mean, and ≤ a ≤ . There are three fit parameters, a , σ , and σ (with σ < σ ). Two parameters are determinedby equating the variance σ = aσ + (1 − a ) σ and excesskurtosis k = κ − a (1 − a )( σ − σ ) /σ of g ( A ) to thatof the data. We then perform a least-squares fit of the cumu-lative density function of the data to solve for the remainingdegree of freedom a . (See Appendix B for details.)The resulting distributions g ( A ) are plotted in Fig. 3 asdashed curves. The two-gaussian form works very well forsmall λ , and reasonably well for large λ . The discrepancy inFig. 3(f) may be simply due to the lack of sufficient data pointsto provide good statistics for these particular parameters.In Fig. 4, we show data related to this two-component de-scription ( σ , , σ , κ ), for the observables S z and S z S zp +2 in the ladder system. The two standard deviations gener-ally become equal at intermediate λ (the ratio σ /σ dropsto near unity), indicating that a single-gaussian descriptionworks well away from integrability. In (d), we show the kur-tosis κ of the distribution, used as an input for the fit. Thekurtosis is close to (the kurtosis value of the gaussian distri-bution) in the intermediate regime, again showing that a sin-gle gaussian is a good description for the distribution of A αβ values in generic systems. The kurtosis is significantly largerthan 3 as one approaches the integrable points, signifying astronger central peak than that of a single gaussian.In Figs. 4(e) and (f), we provide a scaling analysis by plot-ting σ /σ and κ as a function of the Hilbert-space dimension D for the observable ˆ A = S z S zp +2 . In the non-integrableregime, the values remain near σ /σ ≈ and κ ≈ asthe sizes are increased. For λ = 0 , the kurtosis κ increasesaway from 3 with larger D , indicating that the central peakgets stronger relative to the larger gaussian as the systemsize increases. This is consistent with our explanation of thetwo-component structure in terms of symmetry sectors: thenumber of eigenstate pairs belonging to different subspace in-creases faster with D compared to the number of eigenstatepairs within the same symmetry subspace.Also noteworthy is the behavior at the near-integrable point λ = 0 . : the data shows convergence with increasing D to-ward the non-integrable values σ /σ = 1 and κ = 3 . Inparticular κ shows non-monotonic behavior: it first increaseslike in the integrable case, and only beyond a certain size startsdecreasing back toward the single-gaussian value κ = 3 . This − − λ . . . . . . . . S z (a) . σ , σ , σ .
01 0 . − − λ . . . . . . . . (b) S z S zp +2 . σ , σ , σ .
01 0 . − − λ σ / σ (c) S z S z S zp +2 σ / σ .
01 0 . λ − − λ (d) κ .
01 0 . λ D S z S zp +2 (e) σ / σ D D (f) S z S zp +2 κ Dλ : .
05 0 . FIG. 4. (Color online) Characteristics of the distribution of A αβ .(a,b) Standard deviations σ , σ , and σ of the “inner” gaussian ofEq. (8), the full distribution, and the “outer” gaussian, in increasingorder. The thicker curve is σ . System size is ( L, N ↑ ) = (15 , .(c,d) Ratio σ /σ and kurtosis κ . The dashed horizontal lines arethe values for the gaussian distribution ( σ /σ = 1 , κ = 3 ). Weshow results for ˆ A = S z in red (squares) and for ˆ A = S z S zp +2 inblue (circles). (e,f) σ /σ and κ as a function of the Hilbert-spacedimension D for several values of λ . In (e), data for the smallestsystem size for λ = 5 is absent — the procedure does not yield asolution for σ i due to the low density of states. is a manifestation of the phenomenon that, near but not ex-actly at integrability, the system size needs to be large to showgeneric non-integrable behavior [14]. V. SCALING ANALYSIS
In this section we analyze the system-size dependence ofthe average magnitudes of | A αβ | , which corresponds to thewidths of the distributions studied in the previous section.The average value of | A αβ | close to the diagonal in thecentral part of the spectrum (omitting the lowest and highest of the energy range, as indicated by the dashed squaresin Fig. 1) is given by γ = 1˜ N ∼ (cid:88) α,βα (cid:54) = β | E β − E α |≤ ω max | A αβ | = S A (0 , ω max ) . (9)Here, ∼ (cid:80) denotes summation over the relevant state pairs: Itincludes all α and β within the bulk of the spectrum with ˆ A = S z − − − (a) γ low − − − γ − − − (b) γ all l o w a ll (c) γ l o w / γ a ll ˆ A = S z S zp +2 − − − (d) γ low − − − γ D − − − (e) γ all D l o w a ll (f) γ l o w / γ a ll Dλ : .
05 0 . FIG. 5. (Color online) System-size scaling analysis of average | A αβ | through the quantities (a) γ low , (b) γ all , and (c) γ low /γ all for the observable ˆ A = S z , for several values of λ . The respectiveresults for ˆ A = S z S zp +2 are shown in panels (d–f). The dotted linesin (a), (b), (d) and (e) are γ = 1 /D . α (cid:54) = β and with | E α − E β | ≤ ω max , where ω max acts asthe frequency cutoff.The quantity γ depends on the cutoff frequency ω max . Weconsider two values of ω max . First, we define a low-frequencymeasure, γ low = γ ( ω max = 0 . ω ) , where ω is the “typicalfrequency” [Eq. (6)]. Second, we define γ all = γ ( ω max →∞ ) including all state pairs within the bulk of the spectrum(dashed square in Figure 1).In Fig. 5, we plot the quantities γ low , γ all , and the ratio γ low /γ all as a function of Hilbert-space size D for severalvalues of λ , in the top row for the observable ˆ A = S z andin the lower row for ˆ A = S z S zp +2 . Both γ low and γ all showa power-law behavior, ∝ D − . the scaling is almost exact for γ all . This scaling behavior is consistent with the scaling of thetemporal fluctuations being exponential in L , as observed inRef. [21].The D − scaling for non-integrable systems can be ex-plained by using the central limit theorem invoked in the pre-vious section to explain the gaussian form of the distributionof A αβ values. From Eq. (7), we interpret A αβ as the aver-age of the random variables X γ ≡ Dc ( α ) ∗ γ c ( β ) γ a γ . Assuming c ( α ) γ and c ( β ) γ to be independent random variables, each withvariance /D due to normalization of the eigenfunctions, therandom variables X γ can be argued to be independent and tohave D -independent variance, var( X γ ) ∼ . The central limittheorem then states that the variance of A αβ (i.e., the averageof | A αβ | ) scales as var( X γ ) /D ∼ /D . As in Ref. [14], theargument relies or difficult-to-prove randomness assumptions.The D − scaling can be more directly understood by esti-mating the average value of all | A αβ | , including the edgesof the spectrum and the diagonal elements, which is equal to Tr( A ) /D . For local observables, Tr( A ) ∝ D . (In fact, forthe two observables in Fig. 5, Tr( A ) = D exactly.) The scal-ing of the average as ∝ D − immediately follows. Of course,both γ all and γ low are slightly different from Tr( A ) /D . For γ all , the states outside the central part and the diagonal ele-ments are not included, as opposed to Tr( A ) /D where theyare included. Nevertheless, in Fig. 5(b,e), γ all (data points)follows Tr( A ) /D = 1 /D (dotted line) very closely, for all λ . This shows that the contribution from the diagonal ele-ments and from the edge states are negligible.In Fig. 5(a,d), γ low shows approximate ∝ D − scaling. Themagnitudes are generally larger than /D for larger D , reflect-ing the fact that the low-frequency | A αβ | are on average largerthan other off-diagonal matrix elements (as seen previously inFigs. 1 and 2). This is also reflected, Fig. 5(c,f), in the ratio γ low /γ all . The ratio > for larger sizes. The effect is mostprominent for large λ , which reflects the very large concen-tration near the diagonal seen in Figs. 1(c) and 2(c). The ra-tios γ low /γ all increase with system size. It is conceivable thatthese ratios will converge to a constant at larger D , so that γ low also converges to a ∝ D − dependence. The availabledata hints at such behavior, but the available systems sizes areinsufficient to make a definitive statement.The scaling of γ low deviates from the ∝ D − especiallyfor smaller systems and close to integrability. This behavior isreminiscent of the fluctuations of A αα close to (but not exactlyat) an integrable point, where the scaling deviates from D − / for intermediate sizes but converges to D − / as the systemsize is increased [14]. VI. BOSE-HUBBARD CHAIN
To evaluate the generality of the results presented in previ-ous sections with the XXZ ladder system, we present in thissection a summary of analogous data for the Bose-Hubbardchain, Eq. (3). We show data for the observable ˆ A = b † b + b † b .The frequency-resolved analysis of the matrix elements A αβ is performed for the values λ = 1 , typical for the noninte-grable regime, and λ = 0 . and close to the two integrablelimits. The results, in Figs. 6(a)–(c), are qualitatively similarto the ones for the XXZ model in Fig. 2.In Figs. 6(d) and (e), we analyze the distribution of the val-ues A αβ for low frequencies, by fitting to the two-componentdistribution as described in Sect. IV. The ratio σ /σ andthe kurtosis κ are high ( (cid:29) and (cid:29) , respectively) at ornear integrability. In the nonintegrable regime (representedby λ = 1 ), both quantities are close to the values appropriatefor a gaussian distribution ( and , respectively).In Figs. 6(d)–(h) we show data for λ = 10 − as a sub-stitute for the exact integrable point λ = 0 , because thestrong oscillations at λ = 0 make our procedure for extract-ing σ , (Appendix B) unreliable. At accessible sizes, the λ = 10 − data indeed shows size-dependence characteristicof integrable points: increase of κ > with increasing sys-tem size. At some very large system size, κ ( D ) is expectedto decrease again. Such nonmonotonic behavior is a signatureof proximity to integrability. The non-monotonic behavior is ˆ A = b † b + b † b λ = 0 . (a) . . S A , A α β ω E α E β − − λ = 1 (b) ω E α E β λ = 10 (c) ω E α E β
50 705070 D (d) σ / σ D D (e) κ Dλ : .
001 0 . − − − (f) γ low − − − γ D − − − (g) γ all D (h) γ l o w / γ a ll D FIG. 6. (Color online) Bose-Hubbard model; observable ˆ A = b † b + b † b . (a,b,c) Frequency-dependence shown through S A , asin Fig. 2. Insets show fragments of the density plot of | A αβ | as func-tion of E α and E β , as in Fig. 1. The system size is ( L, N b ) = (7 , .(d,e) Analysis of the shape of the distribution through σ /σ and kur-tosis κ , as in Fig. 4. (f,g,h) Average matrix element γ low and γ all forlow and all frequencies, and their ratio, as in Fig. 5. In (f,g), thedotted lines are γ = 4 /D . visible at available system sizes for the λ = 0 . data.The root-mean-squared γ all [Fig. 6(g)] of the matrix ele-ments | A αβ | without frequency cutoff shows a D − scaling,the values being close to Tr( A ) /D = 4 /D [67]. The scal-ing of γ low is not equally clear at these sizes. The erratic be-havior close to integrability is possibly due to the presence ofmany very sharp peaks in S A , especially at low frequencies[see Fig. 6(a)]. VII. DISCUSSION
Motivated by the importance of off-diagonal matrix ele-ments ( A αβ ) of local operators in the physics of time evo-lution after a quantum quench, we have provided a detailedstudy of the statistical properties of such matrix elements, forsystems with short-range interactions. Data on off-diagonalmatrix elements have appeared in the non-equilibrium litera-ture (e.g. [3, 12, 26]); the present work extends such work toprovide a systematic account of these objects. We have cho-sen multiple observables and families of Hamiltonians, andhave thus been able to extract general features. We have alsoelucidated the role of proximity to integrability as well as theapproach to the thermodynamic limit.The distribution of values of A αβ is gaussian for genericsystems, but deviates in a particular way (stronger peak atzero, or mixture of two gaussian-like distributions) as one ap-proaches integrability. We have used this to formulate a quan-titative characterization of proximity to integrability, throughthe kurtosis κ of the distribution. We find κ ∼ for non-integrable (generic or chaotic) systems, and a larger κ that increases with system size for integrable systems. This dis-tinction makes it possible to graphically represent our idea,formulated in Ref. [14], that distance from integrability can becharacterized by a length scale —for near-integrable systems,the size-dependence κ ( D ) shows an initial increase followedby a decrease beyond a certain size. This size D where κ ( D ) is maximal characterizes the proximity to integrability, and in-creases as one approaches integrability. These properties mayprove to be useful signatures for proximity to integrability, inthe sense that they can be determined straightforwardly fromFourier transforms of the time evolution after a quench, whichmay be feasible to observe in experiments.The average magnitude of the matrix elements, S ( ω ) or γ ,determines the magnitude of temporal fluctuations of (cid:104) A (cid:105) ( t ) after a quantum quench. The scaling of this quantity for non-integrable systems, ∼ /D , is consistent with the scaling oftemporal fluctuations known from the literature [21]. We alsofind that the low-frequency average is higher than the aver-age over all frequencies, γ low > γ all (Figs. 5,6), reflecting theoverall decrease of S ( ω ) with increasing frequency (Fig. 2).This suggests that, for quenches in nonintegrable systems,low-frequency contributions are likely to dominate in the timeevolution, regardless of whether or not the initial conditionsare very local in energy.The ∼ /D scaling can be argued from the central limit the-orem assuming wavefunction coefficients of non-integrableHamiltonians to be effectively random. This is a recurringassumption in this field (e.g., [13, 14, 21]), usually withoutrigorous proof, but with a similar status as the ETH, namely,as a plausible hypothesis, the validity of which has to be es-tablished by numerical results. Nevertheless, this argument isuseful because it also provides a plausible explanation for thegaussian distributions of the matrix elements. For the scaling,we have provided an alternate argument based on the trace oflocal operators, which turns out to work well especially for γ all .The present work raises a number of questions for furtherstudy. As a new characterization of integrability, the double-peak structure of the A αβ distribution deserves to be betterunderstood. The relative weight of the inner peak is presum-ably connected to the distribution of sizes of the many sub-spaces that the Hilbert space is divided into, due to the manyconservation laws present at integrability. At present, we donot have a quantitative understanding of the exact connec-tion between the subspace distribution and the non-gaussiandistribution of the off-diagonal matrix elements, although thepresence of many subspaces provides a plausible explanationfor the double-peak form. A related question is the type of deviation from the gaussian shape of the A αβ distributionfor systems with a few (nonzero but O ( L ) ) conservationlaws. It would also be interesting to find out whether the two-component versus gaussian (single-component) structures canbe related to differences in real-time relaxation and fluctua-tion behaviors between integrable and non-integrable systems.Also, it is possible that our findings for near-integrable pointsmight have consequences for “pre-thermalization” behaviors[57, 68–71]. ACKNOWLEDGMENTS
We thank J. Dubail, P. Ribeiro, L. Santos, and J. M. St´ephanfor interesting discussions.
Appendix A: Time evolution and off-diagonal matrix elements
In this Appendix we outline some of the connections totime evolution which motivates the study of off-diagonal ma-trix elements. We consider a isolated quantum system withHamiltonian H , with eigenvalues E α and eigenstates | ψ α (cid:105) .Under this Hamiltonian, the time evolution of the initial state | Ψ(0) (cid:105) , that may be the result of a quench at t = 0 from an-other Hamiltonian, is given by | Ψ( t ) (cid:105) = (cid:80) α c α e − i E α t | ψ α (cid:105) ,where c α = (cid:104) ψ α | Ψ(0) (cid:105) are the expansion coefficients in theeigenstate basis. Given an observable A , its expectation valueevolves as (cid:104) A (cid:105) ( t ) = (cid:104) Ψ( t ) | ˆ A | Ψ( t ) (cid:105) = (cid:88) α,β c ∗ α c β A αβ e i( E α − E β ) t . (A1)The long-time average of this quantity is (cid:104) A (cid:105) ( t ) = lim T →∞ T (cid:90) T (cid:104) A (cid:105) ( t ) dt. (A2)For a non-degenerate spectrum, the off-diagonal terms do notcontribute, so that (cid:104) A (cid:105) ( t ) = (cid:80) α | c α | A αα .While the A αα determine the long-time average, these di-agonal matrix elements do not say anything about the tempo-ral fluctuations f A ( t ) ≡ (cid:104) A (cid:105) ( t ) − (cid:104) A (cid:105) ( t ) around the average.A representative value for the magnitude of temporal fluctua-tions is its root-mean-square ( σ t A ) ≡ [ f A ( t )] = lim T →∞ T (cid:90) T [ f A ( t )] dt. (A3)Using Eq. (A1), one finds that ( σ t A ) = (cid:88) α,βα (cid:54) = β | c α | | c β | | A αβ | , (A4)under the assumption that the spectrum is incommensurate,i.e., when there are no degeneracies and E α + E β = E γ + E δ implies that ( α, β ) = ( γ, δ ) or ( α, β ) = ( δ, γ ) .The fluctuation amplitude ( σ t A ) can be considered as a cor-relator of f A ( t ) with itself. Generalizing to correlators at dif-ferent times, we get the autocorrelation function, f A ( t ) f A ( t + τ ) = (cid:88) α,βα (cid:54) = β | c α | | c β | | A αβ | e i( E β − E α ) τ , (A5)which appears in formulations of non-equilibrium fluctuation-dissipation relations [26, 33]. The Fourier transform of thisquantity is s ( ω ) = (cid:88) α,βα (cid:54) = β | c α | | c β | | A αβ | δ ( ω − ( E β − E α )) . (A6)The strength of the fluctuations at frequency E β − E α is equalto | c α | | c β | | A αβ | .Eqs. (A4) and (A6) demonstrate the roles of A αβ in real-time considerations. The quantity γ in our work can be re-garded as a general version of the right hand side of (A4)which is independent of any particular quench protocol. Thequantity S ( ω ) is similarly a smoothed version of the righthand side of (A6), again omitting reference to specific initialstates. Appendix B: Fit to the distribution of A αβ In Sec. IV, we have fitted the sum of two gaussian distri-butions g ( A ) [Eq. (8)] to the actual distribution d ( A ) of theoff-diagonal elements in a small frequency window. The fitparameters in this distribution are a , the mutual weight of thetwo terms, and σ and σ , the standard deviations. For thefits shown in Fig. 3 and for the data plotted in Fig. 4, we im-pose that the fitted distribution g ( A ) has the same variance σ and kurtosis κ as the actual data. This yields the equations σ = aσ + (1 − a ) σ and κ − a (1 − a )( σ − σ ) /σ .By solving these equations for given σ and k = κ − , weobtain expressions for σ and σ in terms of a , σ , = σ (cid:18) ∓ (cid:113) k (1 − a ) /a (cid:19) . (B1)We have imposed σ ≤ σ ≤ σ . The remaining variable a can be obtained in several ways. For Figs. 3 and 4, we haveobtained a by numerically minimizing the integrated squaredifference between the cumulative density function of g ( A ) and that of the data d ( A ) . This method yields an “optimal”value of a , which is substituted into Eq. (B1) in order to obtain σ and σ . However, when the cumulative density distributionof the data behaves erratically due to very few states beinginvolved, this procedure might fail and give an optimal valueof a outside the range [0 , (e.g., the λ = 5 data in Fig. 4). [1] J. M. Deutsch, Phys. Rev. A , 2046 (1991).[2] M. Srednicki, Phys. Rev. E , 888 (1994).[3] M. Rigol, V. Dunjko, and M. Olshanii, Nature , 854 (2008).[4] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore,Rev. Mod. Phys. , 863 (2011).[5] J. Gemmer, M. Michel, and G. Mahler, Quantum Thermody-namics (Springer Berlin Heidelberg, 2009).[6] M. Rigol, Phys. Rev. Lett. , 100403 (2009).[7] M. Rigol and L. F. Santos, Phys. Rev. A , 011604 (2010).[8] G. Biroli, C. Kollath, and A. M. L¨auchli, Phys. Rev. Lett. ,250401 (2010).[9] T. N. Ikeda, Y. Watanabe, and M. Ueda, Phys. Rev. E ,021130 (2011).[10] T. N. Ikeda, Y. Watanabe, and M. Ueda, Phys. Rev. E ,012125 (2013).[11] M. Rigol and M. Srednicki, Phys. Rev. Lett. , 110601(2012).[12] R. Steinigeweg, J. Herbrych, and P. Prelovˇsek, Phys. Rev. E ,012118 (2013).[13] C. Neuenhahn and F. Marquardt, Phys. Rev. E , 060101(2012).[14] W. Beugeling, R. Moessner, and M. Haque, Phys. Rev. E ,042112 (2014).[15] V. K. B. Kota, A. Rela˜no, J. Retamosa, and M. Vyas, J. Stat.Mech. , P10028 (2011).[16] A. Motohashi, Phys. Rev. A , 063631 (2011).[17] M. Rigol, Phys. Rev. A , 053607 (2009).[18] E. Khatami, M. Rigol, A. Rela˜no, and A. M. Garc´ıa-Garc´ıa,Phys. Rev. E , 050102 (2012).[19] C. Gramsch and M. Rigol, Phys. Rev. A , 053615 (2012). [20] K. He, L. F. Santos, T. M. Wright, and M. Rigol, Phys. Rev. A , 063637 (2013).[21] P. R. Zangara, A. D. Dente, E. J. Torres-Herrera, H. M.Pastawski, A. Iucci, and L. F. Santos, Phys. Rev. E , 032913(2013).[22] S. Ziraldo and G. E. Santoro, Phys. Rev. B , 064201 (2013).[23] P. Reimann, Phys. Rev. Lett. , 190403 (2008).[24] A. J. Short, New J. Phys. , 053009 (2011).[25] L. C. Venuti and P. Zanardi, Phys. Rev. E , 012106 (2013).[26] E. Khatami, G. Pupillo, M. Srednicki, and M. Rigol, Phys. Rev.Lett. , 050403 (2013).[27] G. P. Brandino, A. De Luca, R. M. Konik, and G. Mussardo,Phys. Rev. B , 214435 (2012).[28] M. C. Ba˜nuls, J. I. Cirac, and M. B. Hastings, Phys. Rev. Lett. , 050405 (2011).[29] S. Genway, A. F. Ho, and D. K. K. Lee, Phys. Rev. A ,023609 (2012).[30] S. Genway, A. F. Ho, and D. K. K. Lee, Phys. Rev. Lett. ,130408 (2013).[31] M. Fagotti, M. Collura, F. H. L. Essler, and P. Calabrese, Phys.Rev. B , 125101 (2014).[32] M. Fagotti and F. H. L. Essler, Phys. Rev. B , 245107 (2013).[33] M. Srednicki, J. Phys. A , 1163 (1999).[34] S. Ziraldo, A. Silva, and G. E. Santoro, Phys. Rev. Lett. ,247205 (2012).[35] L. F. Santos, F. Borgonovi, and F. M. Izrailev, Phys. Rev. E ,036209 (2012).[36] E. J. Torres-Herrera and L. F. Santos, Phys. Rev. E , 042121(2013); Phys. Rev. A , 043620 (2014). [37] N. Linden, S. Popescu, A. J. Short, and A. Winter, Phys. Rev.E , 061103 (2009).[38] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, Phys. Rev.Lett. , 050405 (2007).[39] T. Barthel and U. Schollw¨ock, Phys. Rev. Lett. , 100601(2008).[40] A. Iucci and M. A. Cazalilla, Phys. Rev. A , 063619 (2009).[41] L. F. Santos and M. Rigol, Phys. Rev. E , 031130 (2010).[42] A. C. Cassidy, C. W. Clark, and M. Rigol, Phys. Rev. Lett. ,140405 (2011).[43] M. A. Cazalilla, A. Iucci, and M.-C. Chung, Phys. Rev. E ,011133 (2012).[44] K. He and M. Rigol, Phys. Rev. A , 063609 (2012).[45] M. Collura, S. Sotiriadis, and P. Calabrese, Phys. Rev. Lett. , 245301 (2013).[46] M. Fagotti and F. H. L. Essler, J. Stat. Mech. , P07012(2013).[47] M. Fagotti, Phys. Rev. B , 165106 (2013).[48] J.-S. Caux and F. H. L. Essler, Phys. Rev. Lett. , 257203(2013).[49] M. Collura, M. Kormos, and P. Calabrese, J. Stat. Mech. ,P01009 (2014).[50] E. Canovi, D. Rossini, R. Fazio, G. E. Santoro, and A. Silva,Phys. Rev. B , 094431 (2011).[51] J. Sirker, N. P. Konstantinidis, F. Andraschko, and N. Sedlmayr,Phys. Rev. A , 042104 (2014).[52] R. Nandkishore and D. A. Huse, arXiv:1404.0686 (2014).[53] J. M. Hickey, S. Genway, and J. P. Garrahan, arXiv:1405.5780(2014).[54] S. Dubey, L. Silvestri, J. Finn, S. Vinjanampathy, and K. Ja-cobs, Phys. Rev. E , 011141 (2012).[55] R. Steinigeweg, A. Khodja, H. Niemeyer, C. Gogolin, andJ. Gemmer, Phys. Rev. Lett. , 130403 (2014). [56] A. Khodja, R. Steinigeweg, and J. Gemmer, arXiv:1408.0187(2014).[57] M. Marcuzzi, J. Marino, A. Gambassi, and A. Silva, Phys. Rev.Lett. , 197203 (2013).[58] H. Niemeyer, K. Michielsen, H. De Raedt, and J. Gemmer,Phys. Rev. E , 012131 (2014).[59] C. Kollath, A. M. L¨auchli, and E. Altman, Phys. Rev. Lett. ,180601 (2007).[60] G. Roux, Phys. Rev. A , 021608 (2009).[61] G. Roux, Phys. Rev. A , 053604 (2010).[62] J. M. Zhang, C. Shen, and W. M. Liu, Phys. Rev. A , 063622(2011).[63] J. M. Zhang, C. Shen, and W. M. Liu, Phys. Rev. A , 013637(2012).[64] S. Sorg, L. Vidmar, L. Pollet, and F. Heidrich-Meisner,arXiv:1405.5404 (2014).[65] M. Cramer, C. M. Dawson, J. Eisert, and T. J. Osborne, Phys.Rev. Lett. , 030602 (2008).[66] M. Cramer and J. Eisert, New J. Phys. , 055020 (2010).[67] Here, Tr( A ) = (cid:0) L + N b L − (cid:1) = 2 D ( L + N b ) N b /L ( L + 1) whichequals DL/ ( L + 1) if N b = L , i.e., at unit filling.[68] M. Moeckel and S. Kehrein, Phys. Rev. Lett. , 175702(2008).[69] M. Eckstein and M. Kollar, Phys. Rev. Lett. , 120404(2008).[70] M. Kollar, F. A. Wolf, and M. Eckstein, Phys. Rev. B ,054304 (2011).[71] J. M. Zhang, F. C. Cui, and J. Hu, Phys. Rev. E85