Closure of the entanglement gap at quantum criticality: The case of the Quantum Spherical Model
CClosure of the entanglement gap at quantum criticality: The case of the QuantumSpherical Model
Sascha Wald, ∗ Ra´ul Arias,
2, 3 and Vincenzo Alba Max-Planck-Institut f¨ur Physik Komplexer Systeme, N¨othnitzer Straße 38, D-01187, Dresden, Germany SISSA and INFN, via Bonomea 265, 34136 Trieste, Italy Instituto de F´ısica La Plata - CONICET and Departamento de F´ısica,Universidad Nacional de La Plata C.C. 67, 1900, La Plata, Argentina Institute for Theoretical Physics, Universiteit van Amsterdam,Science Park 904, Postbus 94485, 1098 XH Amsterdam, The Netherlands
The study of entanglement spectra is a powerful tool to detect or elucidate universal behaviour inquantum many-body systems. We investigate the scaling of the entanglement (or Schmidt) gap δξ ,i.e., the lowest laying gap of the entanglement spectrum, at a two-dimensional quantum critical point.We focus on the paradigmatic quantum spherical model, which exhibits a second-order transition,and is mappable to free bosons with an additional external constraint. We analytically show thatthe Schmidt gap vanishes at the critical point, although only logarithmically. For a system on atorus and the half-system bipartition, the entanglement gap vanishes as π / ln ( L ) , with L the linearsystem size. The entanglement gap is nonzero in the paramagnetic phase and exhibits a faster decayin the ordered phase. The rescaled gap δξ ln ( L ) exhibits a crossing for different system sizes at thetransition, although logarithmic corrections prevent a precise verification of the finite-size scaling.Interestingly, the change of the entanglement gap across the phase diagram is reflected in the zero-mode eigenvector of the spin-spin correlator. At the transition quantum fluctuations give rise toa non-trivial structure of the eigenvector, whereas in the ordered phase it is flat. We also showthat the vanishing of the entanglement gap at criticality can be qualitatively but not quantitativelycaptured by neglecting the structure of the zero-mode eigenvector. I. INTRODUCTION
In the last two decades the study of quantum entangle-ment has revolutionised our understanding of quantummany-body systems [1–4]. The main ingredient to ad-dress entanglement-related questions in a quantum sys-tem S is the reduced density matrix ρ A of a subsystem A ⊂ S . Given the ground-state ∣ Ψ ⟩ of S and a spatialbipartition of S = A ∪ ¯ A (see e.g. Fig. 1), ρ A is defined as ρ A = Tr ¯ A ∣ Ψ ⟩⟨ Ψ ∣ . (1)The entanglement spectrum (ES) { ξ i = − ln ( λ i ) ∣ λ i ∈ spec ( ρ A )} has been the subject of intense investigation.Pioneering studies [5–8] were fueled by the rapid successof the density matrix renormalisation group [9, 10] tosimulate one-dimensional quantum many-body systems.The interest in the ES was revived after it was discov-ered that for fractional quantum Hall states the lowerpart of the ES contains universal information about theedge modes and the conformal field theory (CFT) de-scribing them [11]. This sparked intense theoretical ac-tivity to clarify the nature of the ES in fractional quan-tum Hall systems [12–23], topologically ordered phasesof matter [24–26], frustrated and magnetically orderedsystems [23, 27–39], CFT systems [40–43], and systemswith impurities [44].In this work we investigate the ES in critical two-dimensional quantum many-body systems. We focus on ∗ [email protected] the lowest laying entanglement gap δξ defined as δξ = ξ − ξ , (2)where ξ and ξ are the lowest and the first excitedES level, respectively. The behaviour of the entangle-ment gap at quantum critical points has not been thor-oughly addressed, except for one-dimensional systems [5–7, 13, 29, 30, 34, 45, 46]. Several exact results suggestthat at one-dimensional quantum critical points δξ van-ishes. For instance, in CFT systems δξ decays logarith-mically as ∝ / ln ( L ) with the subsystem’s length L [40].Similar scaling is found in corner transfer matrix calcu-lations [45] (see also [8] for a review). Higher-dimensionsare far less explored. Interestingly, it has been arguedthat the closing of the entanglement gap does not neces-sarily signal critical behaviour [23]. Similar conclusionshave been reached by considering the ES of a bipartitionin momentum space [47]. Still, the ES can be useful todistinguish different phases of matter. This is the case forsystems that exhibit order by breaking of a continuoussymmetry [31]. It has been suggested that deep in theordered phase the lower part of the ES contains the fin-gerprints of symmetry breaking, being reminiscent of theso-called Anderson tower-of-states [48–50]. This has beenverified by analytical calculations in the quantum rotormodel [31], numerical simulations in the two-dimensionalBose-Hubbard model in the superfluid phase [33] (see also[39]), and also in two-dimensional Heisenberg models onthe square [36] and on the kagome lattice [38]. A sig-nature of the tower-of-states scenario is that the gaps inthe lower part of the ES decay as a power-law with thesubsystem volume, with multiplicative logarithmic cor- a r X i v : . [ c ond - m a t . s t a t - m ec h ] S e p FIG. 1. Bipartition of the two dimensional lattice as A ∪ ¯ A .Periodic boundary conditions in both directions are used. (a)A bipartition with straight boundary between A and ¯ A . A contains ∣ A ∣ = L × (cid:96) x sites and spans the full lattice along theˆ y direction. (b) Bipartition with a corner. Now ∣ A ∣ = (cid:96) x (cid:96) y .We also define the ratios ω x ( y ) = (cid:96) x ( y ) / L . We mostly considerthe the case with ω y = rections [31]. Higher ES levels are expected to exhibit amuch slower decay [31, 33, 37]. The behaviour of the en-tanglement gap upon approaching the critical point hasnot been investigated thoroughly.Here we address this issue in the quantum sphericalmodel [51–55] (QSM). The QSM is a paradigmatic many-body system in which the effects of strongly interactingdegrees of freedom may be studied at a considerably lowcost, as the model can be mapped to free bosons subjectto an additional external constraint. Despite its simplic-ity it exhibits several salient features of realistic quantummany-body systems. For instance, its classical versionserved as a testing ground for the theory of critical phe-nomena and finite size scaling [56]. In two dimensionsthe QSM exhibits a standard paramagnetic (disordered)phase and a ferromagnetic (ordered) one, which are sep-arated by a second order quantum phase transition. Theuniversality class of the transition is that of the three-dimensional classical O ( N ) vector model [57] in the large N limit [52, 53, 58]. Surprisingly, entanglement proper-ties of the QSM are rather unexplored, although there isrecent interest [59–61]. We should stress that althoughthe results that we are going to derive for the ES cannotbe considered general, they certainly represent an inter-esting case study, and can be useful to understand thegeneric behaviour of ES in quantum many-body systems.Here we consider a two-dimensional lattice of linearsize L with periodic boundary conditions in both direc-tions. The typical bipartitions that we use are reportedin Fig. 1. Figure 1 (a) shows a bipartition with a straightboundary between A and its complement, with A span-ning the full lattice along the ˆ y direction. This is notthe case in Fig. 1 (b), where the boundary has a corner.The effect of corners in the scaling of the entanglemententropies is nontrivial, and it has been studied intenselyin the last decade [4, 62–69]. Since the QSM is mappableto free bosons, entanglement-related observables can becalculated from the two-point correlations functions [8].Here we show that δξ (cf. (2)) is nonzero in the param-agnetic phase, whereas it vanishes in the ordered phase,as expected [31]. This is compatible with the numeri- cal results in [33] (see also [36, 38]). At the quantumcritical point, in the case of straight boundary the en-tanglement gap vanishes as π / ln ( L ) . However, we showthat logarithmic corrections are present, which make itdifficult to robustly verify the finite-size scaling of δξ . Wealso show that the behaviour of the entanglement gap isreflected in the zero-mode eigenvector of the spin-spincorrelation matrix. As the transition is approached fromthe paramagnetic side, the eigenvector flattens, meaningthat all its components become equal. This reflects thepresence of a zero mode. Exactly at criticality, the eigen-vector is not flat in the thermodynamic limit, due to thepresence of strong fluctuations, whereas it is flat in theordered phase. Interestingly, we show that by neglect-ing the structure of the eigenvector at the critical point,i.e., by approximating the eigenvector with the flat vec-tor, we obtain that δξ = A /√ ln ( L ) , which accounts forthe vanishing of the entanglement gap, although it is notquantitatively accurate. We clarify how the behaviouras A /√ ln ( L ) arises from some interesting multiplicativelogarithmic corrections in the expectation values of theQSM correlators with the flat vector. Interestingly, theconstant A depends only on low-energy properties of themodel and on the geometry of the bipartition.The manuscript is organised as follows. In section IIwe introduce the QSM and its phase diagram. In sec-tion III we define the quantities of interest. In section IVwe discuss the finite-size scaling in the QSM. Specifically,in subsection IV A we focus on the so-called gap equation,which ensures the external constraint in the QSM. In sub-sections IV B and IV C we derive the finite-size scaling ofthe spin and momentum correlation functions, respec-tively. In section V we investigate the critical behaviourof δξ . Our prediction is discussed in section V A, andit is compared against numerical results in section V B.We describe the behaviour of δξ across the phase dia-gram of the QSM in subsection V B 1, whereas we ad-dress the vanishing of δξ and its finite-size scaling insubsections V B 2 and V B 3, respectively. In section VIwe discuss how the entanglement gap is related to thezero-mode eigenvector of the correlator, which we intro-duce in subsection VI A. In subsection VI B we show thatby assuming that the eigenvector is flat at criticality onecan qualitatively explain the vanishing of the entangle-ment gap. We conclude in section VII. In Appendix Awe report the derivation of the finite-size scaling of thecorrelation functions in the QSM. In Appendix B we de-rive the expectation values of the correlators with the flatvector. II. QUANTUM SPHERICAL MODEL
The QSM [52–54] on a two dimensional cubic lattice oflinear size L and volume V = L is defined by the Hamil-tonian H = g ∑ n p − J ∑ ⟨ n , m ⟩ s n s m + ( µ − ) ∑ n s n . (3)Here, n = ( n x , n y ) ∈ [ , . . . , L ] denotes a generic lat-tice site, and ⟨ n , m ⟩ a lattice bond joining two nearest-neighbour sites. J > J = s n and p n satisfy the standard bosonic commutation relations [ p n , p m ] = [ s n , s m ] = , [ s n , p m ] = i δ nm . (4)We refer to p n as momentum variable, and to the pa-rameter g as quantum coupling as the model reduces tothe famous classical spherical model [70, 71] in the limit g →
0. The Lagrange multiplier µ is called spherical pa-rameter and fixes the spherical constraint, i.e. ∑ n ⟨ s n ⟩ = V. (5)This means that all allowed configurations of the QSMare located around the sphere in configuration space thatis defined by Eq. (5). Critical properties of the QSMare determined through the self-consistent behaviour of µ [53]. The two dimensional QSM does not exhibit afinite temperature phase transition [70, 71], although itpossesses a ground-state transition, i.e., at T = p n = √ V ∑ k e − i nk π k , s n = √ V ∑ k e i nk q k . (6)Here the sum over k = ( k x , k y ) runs in the first Bril-louin zone k i = π / L j , with j ∈ [− L / , L / ] integer. TheHamiltonian (3) in Fourier space reads H = ∑ k g π k π − k + Λ k q k q − k (7)with the single-particle dispersion relationΛ k = √ µ + ω k with ω k = − cos k x − cos k y (8)In order to fully diagonalise (7) we introduce bosonic lad-der operators b k and b † k obeying standard bosonic com-mutation relations viz. q k = α k b k + b † − k √ , π k = i α k b † k − b − k √ α k = √ g / − k . In terms of theseladder operators, the Hamiltonian (7) is diagonal andreads H = ∑ k E k ( b † k b k + / ) , with E k = √ g Λ k . (10)Entanglement-related properties of Gaussian systemssuch as the QSM stem from the two-point correlationfunctions ⟨ s n s m ⟩ and ⟨ p n p m ⟩ . In equilibrium at zero temperature T =
0, the eigenmodes k of the system areoccupied according to ⟨ b k b k ′ ⟩ = ⟨ b † k b † k ′ ⟩ = ⟨ b † k ′ b k ⟩ = , ⟨ b k ′ b † k ⟩ = δ kk ′ . (11)From Eq. (11), we can thus immediately derive the two-point correlation functions [54] S nm = ⟨ s n s m ⟩ = V ∑ k e i ( n − m )⋅ k α k , (12) P nm = ⟨ p n p m ⟩ = V ∑ k e − i ( n − m )⋅ k α − k , (13) K nm = ⟨ s n p m ⟩ = i2 δ nm . (14)Importantly, from (12) and (13) one obtains the relation P nm = g ∫ dµ S nm , (15)which allows to relate the critical behaviour of thespin correlator to that of the momentum correlator.From (12), one can rewrite the spherical constraint (5)as ∑ n S nn = V ⇒ g = V ∑ k E k . (16)This equation is also called gap equation [72] and impliesthat only the average number of bosons is fixed. Fromthe finite-size expressions (12) (13) and (16), the ther-modynamic limit L → ∞ is obtained in the usual way byreplacing2 πk j L → k ′ j , L ∑ k x ,k y → ∏ j = x,y ∫ π − π dk ′ j π . (17)A crucial observation is that the correlator (12) and thespherical parameter (16) exhibit a singularity for k = g c . For g < g c the ground-state of (3) exhibits magneticorder. At g > g c the ground-state is paramagnetic. Thebehaviour of the QSM is determined by the scaling ofthe spherical parameter µ . In the thermodynamic limit,in the paramagnetic phase one has that µ is finite andnonzero. On the other hand, one has µ = g c can bedetermined analytically. In the thermodynamic limit thespherical constraint (16) is rewritten as √ g = π K ( − µ + √ − √ µ ( µ + ) µ + )[( + µ ) ( µ + √ µ ( µ + ) + ) − ] (18)with the complete elliptic integral [74] K ( x ) = ∫ π dθ √ − x sin ( θ ) . (19)The critical coupling g c follows by imposing the condition µ =
0. This yields g c = π K − ( / − /√ ) ≃ . . (20)The different phases of the model correspond to differentfinite-size scaling behaviours of µ . In the paramagneticphase one has µ = O( ) in the limit L → ∞ . At thecritical point one can show that µ = O( / L ) , whereas inthe ordered phase µ = O( / L ) (see section IV). Thesebehaviours are numerically illustrated in Fig. 2. The uni-versality class of the ground state transition [53] is thatof the large- N vector model in three dimensions, as ex-pected from general renormalisation group arguments.Critical properties of the large- N vector model have beencharacterised analytically [56] and finite-size correctionshave also been investigated [75–79]. III. ENTANGLEMENT SPECTRA ANDENTANGLEMENT GAPS
Here we are interested in the ground-state entanglementspectrum of the QSM, focussing on the two bipartitionsdepicted in Fig. 1. The lattice, with periodic boundaryconditions, is divided into two regions A and ¯ A . Region A is of size ∣ A ∣ = (cid:96) x × (cid:96) y and we define the corresponding as-pect ratios ω x = (cid:96) x / L and ω y = (cid:96) y / L , with 0 ≤ ω x,y ≤
1. InFig. 1 (a) the subsystem A spans the full lattice along theˆ y direction implying that the boundary between the twosubsystems A and ¯ A is straight. This case corresponds to ω y =
1. In Fig. 1 (b), the boundary presents a corner andis thus not straight. The presence of corners has strikingconsequences for entanglement entropies, giving rise tosub-leading universal logarithmic corrections [62–66, 80].The effects of corners in the scaling of the ES have notbeen investigated yet.For the case of a straight boundary with periodicboundary conditions the momentum k y is a good quan-tum number for the correlation matrices (12) and (13),and for the ES. This will be exploited in section V toreduce the computation of the ES of the QSM to that ofan effective one-dimensional model. This dimensional re-duction has been employed to study symmetry-resolvedentanglement entropies [81]. This rather simple observa-tion will also allow to obtain analytically the scaling ofthe entanglement gap at the critical point, by exploitingcorner transfer matrix results [5–7, 45].We now review the calculation of entanglement-relatedquantities in the QSM. Since the QSM is essentially map-pable to a free bosonic model (see section II), its entan-glement properties are derived from the two-point cor-relation functions (12) and (13) (see Ref. [8] for a re-view). The crucial ingredient is the correlation matrix C restricted to the subsystem A , viz. C A = S A ⋅ P A , (21)with S A and P A being the correlation matrices definedin (12) and (13), restricted to the subsystem A . Since inthe remainder we mostly consider the restricted correla-tion matrices S A and P A , we will often omit the subscript A to lighten the notation.For free bosons the reduced density matrix of subsys-tem A is a quadratic operator and is written as [8] ρ A = Z − e −H A , H A = ∑ k (cid:15) k b † k b k . (22)Here H A is the so-called entanglement Hamiltonian, (cid:15) k are the single-particle ES levels, b k are free-bosonic op-erators and Z ensures the normalisation of the reduceddensity matrix Tr ρ A =
1. The spectrum { e k } k = ,..., ∣ A ∣ ofthe correlation matrix C A is simply related to that of H A viz. √ e k =
12 coth ( (cid:15) k ) . (23)The normalisation factor Z is obtained as Z = ∣ A ∣ ∏ j = (√ e j + ) . (24)The ES, i.e., the spectrum of the entanglement Hamil-tonian H A , is obtained by filling the single-particle lev-els (cid:15) k in all the possible ways. To construct the ES, itis convenient to introduce the bosonic occupation num-bers α k = , , . . . in the levels (cid:15) k . The generic ES level ξ ({ α k }) is written as ξ ({ α k }) = ln Z + ∣ A ∣ ∑ j = α j (cid:15) j . (25)The eigenvalues e k satisfy the constraint e k > /
4, im-plying that (cid:15) k >
0. Clearly, the lowest ES level ξ cor-responds to the vacuum state with α k = k . Letus order the (cid:15) k as (cid:15) ≤ (cid:15) ≤ ⋅ ⋅ ⋅ ≤ (cid:15) ∣ A ∣ . The first excitedES level is obtained by populating the smallest singleparticle level (cid:15) . Thus, the lowest entanglement gap δξ (Schmidt gap) is defined as δξ ≡ ξ − ξ = (cid:15) . (26)Here we focus on δξ , although one can define highergaps [82]. IV. FINITE-SIZE CRITICAL CORRELATORSIN THE QSM
As explained in section III, entanglement-related observ-ables, and also the entanglement gap, in the QSM are en-tirely encoded in the two-point correlation functions (12)
10 100 L µ g=3 (ordered)g=g c g ≈ (analytical)1/L (fit) FIG. 2. Spherical parameter µ as a function of linear size L at the quantum critical point at g c (circles), in the orderedphase (squares), and in the paramagnetic phase (diamonds).Note the different scaling with L in the different phases andat the critical point. The dashed-dotted line is the analyticbehaviour γ /( L ) . The dashed line is a fit. and (13). In the following sections we derive the finite-size behaviour of these two-point correlation functions.In section IV A we discuss the gap equation (16). In sec-tions IV B and IV C we the focus on the spin and momen-tum correlators respectively. For the classical sphericalmodel similar results were obtained [75, 83]. A. Spherical parameter
Here we derive the finite-size scaling of the spherical pa-rameter µ at the quantum phase transition. The result isnot new [56] but it is a useful initiation for the discussionof the correlators. To treat the sum over k in (12) weobserve that the following identity holds1 L ∑ k √ µ + ω k = ∫ ∞ dt √ π e −( µ + ) t [ I ( t )+ ∞ ∑ ′ l =−∞ I lL ( t )] , (27)where the prime in the sum indicates that the l = I ν are modified Bessel func-tions of the first kind [74]. To derive (27), we introducean auxiliary integration [72] over t to represent the term ( µ + ω k ) − / , then we employ Poisson’s summation for-mula. Further details are reported in Appendix A. Thefirst term in the brackets in (27) does not depend ex-plicitly on L , and gives the thermodynamic contribution.However, there is an implicit dependence on L through µ .The second term is the genuine finite-size contribution.We are interested in the leading finite-size behaviour forlarge L . In this limit the integral in (27) can be treatedby using a saddle point approximation.In order to use (16), we decompose the diagonal cor- relator S nn as S nn = S ( th ) nn + S ( L ) nn , (28)with the thermodynamic contribution S ( th ) nn = π ∫ d k α k (29)corresponding to the term I ( t ) in (27). The remainingterms in (27) are collected in S ( L ) nn . After expanding thesquare in (27), we observe that S ( L ) nn is written as S ( L ) nn = √ g c √ π ∫ ∞ dte −( µ + ) t ∞ ∑ l,l ′ =−∞ I lL ( t ) I l ′ L ( t ) . (30)In order to extract the large L behaviour of (30) we em-ploy a standard saddle point approximation. The cal-culation is straightforward and details are reported inAppendix A.A striking simplification occurs at the critical pointand in the ordered phase, where µ →
0. One can verifynumerically that at the thermodynamical critical point µ ∝ / L . This is expected because µ ∝ m = / ξ ,with m the mass of the theory and ξ corr the correlationlength, and at the critical point ξ corr ∝ L . In the limit µ →
0, one obtains the surprisingly elegant result (seeAppendix A) S ( L ) nn = − √ g c πL [ ln ( − e −√ µL )− ∞ ∑ l,l ′ = e − L √ µ ( l + l ′ ) √ l + l ′ ] . (31)Interestingly, in (31) the first term is of one-dimensionalnature, and it is obtained by isolating the terms witheither l = l ′ = µ ∝ / L gives rise to a non-trivial behaviour of the correlator as it cancels the factor L in the exponential. It also implies that terms with large l, l ′ are exponentially suppressed, and the sums convergequickly. Double sums as in (31) appear often in latticecalculations, and have been investigated in the past [75,76, 83]. In some cases they can be expressed in terms ofgeneralised Riemann zeta functions [84].Using Eqs. (29) and (31) in the gap equation (16) atcriticality yields1 = √ g c √ π ∫ π − π d k √ µ + ω k + √ g c πL ∞ ∑ l,l ′ = e − L √ µ ( l + l ′ ) √ l + l ′ − √ g c πL ln ( − e −√ µL ) . (32)The integral in (32) has to be considered carefully dueto a ∝ / L contribution in the µ → ∫ d k √ µ + ω k = ∫ d k √ ω k − π √ µ + . . . , (33) A similar decomposition as (28) holds for the generic spin-spincorrelator S nm (see section IV B). where the dots denote subleading terms in 1 / L . The sec-ond term in (33) is the singular term that determines thecritical behaviour of three-dimensional QSM at the ther-mal phase transition [61]. This is not surprising becausethe universality class of the quantum phase transitionin two dimensions is the same [52, 53]. Based on theexpected finite-size scaling µ ∝ / L it is convenient todefine µ = γ L , (34)where the constant γ is to be determined and the factor 2is for later convenience. We substitute the Ansatz (34) inthe gap equation (32), and use the spherical constraint inthe thermodynamic limit (16) at criticality, where µ = γ − ∞ ∑ l,l ′ = e − γ √( l + l ′ ) √ l + l ′ + ln ( − e − γ ) = , (35)where the first term is (33) and the other two are ob-tained from (30). Eq. (35) can be solved numerically toobtain the universal constant γ ≃ . N limit of the three dimensional N − vector model [56, 77].The behaviour of µ in the different regions of the phasediagram of the QSM and the accuracy of (34) are ver-ified in Fig. 2 where we show the numerical solution ofEq. (16). In the paramagnetic region for g > g c one has µ = O( ) . At the critical point and in the ferromagneticphase µ → L → ∞ . The dashed-dotted lineis the analytic result (34) with γ obtained from (35).Below the critical point we expect µ ∝ / L [56], whichis confirmed by the fit (dashed line). B. Spin-spin correlation function S nm We now discuss the finite-size scaling of the spin-spin cor-relation function (12) at the quantum critical point. Weonly discuss the final result, reporting the details of thederivation in Appendix A. First, one can again decom-pose the correlator as S nm = S ( th ) nm + S ( L ) nm , (36)with the thermodynamic contribution S ( th ) nm = √ g c √ ( π ) ∫ d k e i k ( n − m ) √ µ + ω k . (37)As in Eq. (29) there is an implicit dependence on L via µ . The finite-size part has the surprisingly simple form S ( L ) nm = √ g c π ∞ ∑ ′ l,l ′ =−∞ e −√ µF ll ′ ( n , m ) F ll ′ ( n , m ) . (38)Here we defined F ll ′ ( n , m ) = √( lL + n x − m x ) + ( l ′ L + n y − m y ) . (39) The prime in the sum means that the term ( l, l ′ ) = ( , ) has been removed. Again, Eq. (38) holds in the limit L → ∞ and µ →
0. The general expression, which isvalid also in the paramagnetic phase, is reported in Ap-pendix A. From Eq. (38), it is clear that the correlators S nm depend only on n x − m x and n y − m y , as expected dueto translation invariance. Moreover, one has that S nm isperiodic along the two directions, i.e., it is invariant un-der n y − m y → n y − m y ± L and n x − m x → n x − m x ± L . Thisis enforced by the infinite sums over l, l ′ . For a biparti-tion with straight boundary between the two subsystems(Fig. 1 (a)) the invariance under n y − m y → n y − m y ± L re-mains true also for the correlator restricted to A . Finally, S ( L ) nm exhibits an interesting singularity structure. For ω y = ω y <
1. Specifically, the terms with l = l ′ = ± n x − m x → n y − m y → ± L . On the other hand,terms with ∣ l ′ ∣ > ∣ l ∣ > ω x = ω y <
1. Weanticipate that these singularities will give rise to multi-plicative logarithmic corrections in the expectation valueof the correlators that we will show in section VI.
C. Momentum correlation function P nm The same finite-size analysis as in section IV B can becarried out for the momentum correlator P nm (cf. (13)).Following the decomposition P nm = P ( th ) nm + P ( L ) nm , (40)with P ( th ) nm = √ g c π ∫ π − π d k e i k ( n − m ) √ µ + ω k , (41)the finite-size part P ( L ) nm has the same structure as (38),and it reads P ( L ) nm =− π √ g c ∞ ∑ ′ l,l ′ =−∞ e −√ µF ll ′ ( n , m ) F ll ′ ( n , m ) [ F ll ′ ( n , m ) + √ µ ] . (42)This expression is obtained from the spin-spin correlator,cf. Eq (38), by using (15). As for (38), the finite-sizeterm (42) is singular if subsystem A spans the full latticein one of the two directions, i.e., if ω x = ω y =
1. For ω y = l = l ′ = ± n x − m x → n y − m y → ± L . Note that the firstterm in Eq. (42) exhibits a stronger singularity than thesecond one. V. CRITICAL BEHAVIOUR OF THEENTANGLEMENT GAP
We now discuss the critical behaviour of the entangle-ment gap δξ . In subsection V A, by using a dimensionalreduction, we provide an exact result for the case of asmooth boundary between the subsystems. In subsec-tion V B we discuss numerical results. We first discuss thebehaviour of the entanglement gap across the phase dia-gram of the QSM in subsection V B 1. In subsection V B 2we show that at the critical point the entanglement gapvanishes logarithmically with the system size. Finally, insubsection V B 3 we investigate the finite-size scaling δξ near criticality. A. Exact result via dimensional reduction
Let us focus on the bipartition with ω y = y directionimply that the momentum k y is a good quantum num-ber for the correlation matrix C A (cf. (21)) restricted tosubsystem A . Moreover, translation invariance impliesthat by performing a Fourier transform along the ˆ y di-rection the Hamiltonian (3) can be written as the sum of L decoupled quadratic one-dimensional systems [8]. Thisdimensional reduction is effective for any free system,and has been recently employed to study the so-calledsymmetry-resolved entanglement entropies [81]. The factthat k y is a good quantum number implies that the cor-relation matrix C A has a block structure with each blockcorresponding to a different k y viz. C A = ⊕ k y C ( k y ) A , k y = πL j, j = , , . . . L − . (43)It is straightforward to diagonalise a given block withfixed k y by imposing that the eigenvectors of C A are alsoeigenvectors of the momentum along ˆ y with the giveneigenvalue k y . Since we are interested only in the largesteigenvalue e of C A a further simplification occurs. Asthe critical behaviour is associated with the formation ofa uniform magnetization, it is natural to expect that e is in the sector with k y =
0. This can be readily checkednumerically. Thus, in the following we restrict the cal-culation to k y =
0. By imposing that the eigenvectors of C A are “flat” along ˆ y , i.e., they do not depend on y , theproblem is reduced to the diagonalisation of the reducedcorrelation matrix C ( k y = ) A = S ( k y = ) A ⋅ P ( k y = ) A , (44)where we defined the reduced spin and momentum cor-relators as S ( k y = ) A ( n x − m x ) = L ∑ k x e i ( n x − m x ) k x α k x , (45) P ( k y = ) A ( n x − m x ) = L ∑ k x e − i ( n x − m x ) k x α − k x . (46) Eqs. (45) and (46) depend only on the coordinates n x − m x along the ˆ x direction, and subsystem A is the interval oflength (cid:96) x . Here α k x corresponds to α k in Eq. (9) with k y =
0. The correlators (45) and (46) and hence (44)are formally the same as those of the so-called massiveharmonic chain with frequency Ω = √ µ [8]. The fullES of the massive harmonic chain for the bipartition intwo semi-infinite chains has been calculated by using thecorner transfer matrix approach [8]. The reduced densitymatrix ρ A , up to a trivial renormalisation, is written as ρ A ∼ e −H ctm , (47)with the corner transfer matrix Hamiltonian H ctm = ∞ ∑ j = (cid:15) ( j + ) β † j β j , (cid:15) = πK (√ − κ ) K ( κ ) , (48)where β j are bosonic ladder operators. Here K ( x ) is thecomplete elliptic integral of the first kind (see Eq. (19)).The parameter κ is given in terms of Ω as [81] κ = ( + Ω − Ω √ + Ω ) . (49)Eq. (47) holds if A is the half-infinite line. In this limit,as it is clear from Eq. (48), the single-particle ES lev-els are equally spaced [8] with spacing (cid:15) . To determinethe finite-size scaling of the entanglement gap δξ we usethe fact that for L → ∞ at criticality µ ∝ / L (seeEq. (34)). By substituting (34) in the corner transfermatrix results (48) and (49), we obtain that in the large L limit δξ decays logarithmically with L as δξ = π ln ( Lγ ) + O( ln − ( L )) , (50)Note the dependence on the universal constant γ . Toderive (50), one can also observe that close to the crit-ical point, on the paramagnetic side, Eq. (48) gives δξ = π / ln ( ξ corr ) . Eq. (50) then follows from standardscaling arguments. A similar decay of the entanglementgap as in (50) is obtained for critical one-dimensionalsystems [8], both fermionic and bosonic ones. An im-portant remark is that the corner transfer matrix cal-culation is valid for the bipartition in two semi-infinitesystems, which implies that there is only one boundarybetween the two subsystems, in contrast with the bipar-titions Fig. 1, which contain two boundaries because weare using periodic boundary conditions. Despite that, asit will be clear in section V B, Eq. (50) gives the leadingbehaviour for large L of δξ . We anticipate that a loga-rithmic subleading term as O( ln − ( L )) , which is missingin Eq. (50), is present. From Eqs. (22) and (50) oneobtains that the eigenvalue e of C A is given as e = + π ln ( Lγ ) + O( ln − ( L )) . (51)Importantly, the missing O( ln − ( L )) term in (50) willgive a O( ln ( L )) contribution in (51). g δξ L=20L=30L=40L=50L=80L=120L=180L= ∞ δξ g ≈ FIG. 3. Entanglement gap δξ as a function of g and linearsize L : Overview across the phase diagram. The results arefor the bipartition in Fig. 1 (a) with (cid:96) x = L /
2. The verticalline marks the critical point at g c . The continuous line isthe result in the thermodynamic limit. Inset: Scaling of theentanglement gap in the ordered phase at g < g c . B. Numerical results
In this section we discuss numerical results confirmingthe validity of the logarithmic scaling of the entangle-ment gap at criticality. We provide numerical evidencethat the prefactor of the logarithmic decay obeys thestandard finite-size scaling behaviour. For instance, itexhibits a crossing for different system sizes at the criti-cal point. However, logarithmic corrections are present,and a precise finite-size scaling analysis is very challeng-ing.
1. Overview
Before discussing the scaling of δξ at the critical point,it is useful to focus on its behaviour across the phase di-agram of the QSM, see Fig. 3. The figure shows δξ asa function of g for several system sizes L . The entan-glement spectrum is calculated for the bipartition withstraight boundary, i.e., ω y = ω x = / δξ as obtained by usingthe value of the spherical constraint µ in the thermody-namic limit L → ∞ (cf. (32)). This yields µ = O( ) inthe paramagnetic phase and µ = g ≤ g c ). The thermodynamicentanglement gap is obtained by substituting the ther-modynamic value of µ in the finite-size expressions forthe correlators (cf. (12) and (13)) and taking the limit L → ∞ after. This procedure gives the correct thermo-dynamic behaviour of δξ , at least away from the criticalpoint. Although we use the finite-size expressions for thecorrelators, we observe that δξ converges quickly to itsthermodynamic value. This is expected because the be-haviour of the QSM is determined by the scaling of µ .In the ordered phase and at the critical point the spin
10 100 1000 10000 L e g ≈ c ≈ ≈ π ln (8 x/ γ )+A +A ln(8 x/ γ ) FIG. 4. Largest eigenvalue e of the correlation matrix. Dataare for the bipartition in Fig. 1 (a) with ω x = / ω y = e is plotted versus linear size L . In the ordered phase (di-amonds) we observe a fast increase with L , whereas in theparamagnetic phase e = O( ) . Note the logarithmic diver-gence as e ∝ ln ( L ) at the critical point at g c . The dasheddotted line is a fit to e = / π ln ( L / γ )+ A + A ln ( L / γ ) ,with A , A fitting parameters. correlator (12) diverges due to the zero mode. Thus, weregularise the zero-mode by fixing µ = − . As it is clearfrom Fig. 3, this analysis, although it is not rigorous, sug-gests that δξ = µ is finiteand nonzero in the paramagnetic phase.Let us now discuss the finite-size behaviour of δξ . Inthe paramagnetic phase, i.e. g > g c , the approach tothe thermodynamic limit is exponential, which is ex-pected because the model is massive. For g < g c , i.e.,in the ferromagnetic phase, the data suggest a vanish-ing gap. The scaling of the entanglement gap in mag-netically ordered phases has been investigated exten-sively [31, 33, 36, 38, 39]. For instance, in Ref. [31] itwas predicted that in the presence of continuous sym-metry breaking in generic dimension d , δξ should decayas δξ ∝ ( L d − ln ( L )) − . (52)In d = / ln ( L ) ,reflecting the absence of symmetry breaking. In d >
2. Vanishing of the entanglement gap at the quantumcritical point
We now focus on the scaling of the entanglement gap atthe quantum critical point g c ≃ .
100 1000 10000 L e - / π l n ( L / γ ) g=g c ≈ +A ln(L) 20 30 40 50 60 70 80 ln (L) e FIG. 5. Largest eigenvalue e of the correlation matrix: Sub-leading logarithmic correction. Plot of e − / π ln ( L / γ ) versus ln ( L ) . The data are the same as in Fig. 4. The lineis a fit to A + A ln ( L / γ ) , with A , A fitting parameters.The fit gives A ≃ . e obtained by us-ing µ = γ /( L ) and fixing γ = e is plotted versus ln ( L ) .The line is a fit to A ′ + / π ln ( L ) . of δξ we, equivalently, consider the scaling of the largesteigenvalue e of C A . We show our numerical results for e in Fig. 4 as a function of L (note the logarithmic scaleon the x -axis). To highlight the different scaling as com-pared to other regions of the phase diagram, we reportalso data in the paramagnetic phase (square symbols)and in the ferromagnetic phase (diamonds). Within theordered phase e increases faster than logarihmically. Inthe paramagnetic region e exhibits a mild increase forsmall L , saturating at L → ∞ . This is a consequence ofthe finite correlation length in the paramagnetic phase.A dramatically different behaviour is visible at critical-ity (circles), for which we report data up to L ∼ Interestingly, for moderately large L the behaviour of δξ is compatible with a logarithmic increase, althoughEq. (51) suggests a ln ( L ) scaling. This should be at-tributed to the presence of a sub-leading logarithmic termln ( L ) (cf. (51)). A fit to A ln ( L / γ )+ A + A ln ( L / γ ) (dashed-dotted line) gives A ≃ .
01, which is in goodagreement with the prediction 1 / π . One also obtains A ≃ .
04 and A ≃ .
16. Note that A ≃ /
6, as pre-dicted by (51).To further corroborate our results, in Fig. 5 we show e − / π ln ( L / γ ) versus L using a logarithmic scale onthe x -axis. The data are the same as in Fig. 4. Thecontinuous line is a fit to e − π ln ( Lγ ) = A + A ln ( Lγ ) (53)with A and A fitting constants. The logarithmic be-haviour is perfect. Note that this logarithmic term is Note that since ω y =
1, we can use dimensional reduction toattain large system sizes (see section V A). g δ ξ l n ( L ) L=20L=30L=40L=50L=60L=80L=100L=120L=140L=180
FIG. 6. Finite-size scaling of the rescaled entanglement gap δξ ln ( L ) plotted as function of g . Here L is the system size.The vertical line marks the critical point. not predicted by (51). Its origin could be attributed tothe fact that the corner transfer matrix result is obtainedfor the semi-infinite system, i.e., the biparititon with oneboundary. It is interesting to investigate the dependenceon γ of the constant A in (53). In the inset in Fig. 5we show e obtained by fixing µ = γ /( L ) with γ = e is plotted versus ln ( L ) .The dashed-dotted line is a fit to 1 / π ln ( L ) + A ′ . Theperfect linear behaviour suggests that the subleading log-arithmic term is absent or its prefactor is small. A fit to1 / π ln ( L / γ ) + A ′ + A ′ ln ( L ) gives A ′ ≈ . A ′ = ln ( γ / ) .
3. Finite-size scaling analysis
Having established the logarithmic vanishing of δξ at thecritical point, it is natural to investigate its behaviour inthe vicinity of the quantum phase transition. A naturalidea is that δξ obeys standard finite-size scaling [85] δξ ln ( L ) = f (( g − g c ) L / ν ) + . . . , (54)where the dots stand for scaling corrections, f ( x ) is ascaling function and ν is the exponent that governs thedivergence of the correlation length at the critical point.For the QSM one has ν = f ( x ) is determined by the universality class of the QSM,and, in principle, can be calculated. Under the assump-tion that the f ( x ) is analytic, one can expand (54) near g c to obtain δξ ln ( L ) = f ( ) + ( g − g c ) L / ν + . . . , (55)From the analysis in section V A one should expect f ( ) = π ≃ .
8. Eq. (55) implies that the data for δξ fordifferent system sizes should exhibit a crossing at g c . Thiscrossing method for the entanglement gap has been used0 -200 0 200 400 (g-g c )L ν δ ξ l n ( L ) L=50L=100L=500L=1000L=2000L=6000L= ∞ ∞ L= 30000L= 10000L= 1000
FIG. 7. Scaling behaviour of the rescaled entanglement gap. δξ ln ( L ) plotted against ( g − g c ) L / ν . Here g c ≃ . ν = to detect a quantum phase transition in a system of cou-pled one-dimensional models [35]. However, since δξ haslogarithmic corrections, one should expect strong limita-tions, as we are going to show. The scaling Ansatz (54)implies that by plotting the rescaled gap δξ ln ( L ) as afunction of the scaling variable ( g − g c ) L / ν one shouldobserve a data collapse for different system sizes, pro-vided that scaling corrections can be neglected.Our finite-size data for δξ as a function of g for sev-eral system sizes L are shown in Fig. 6 focussing on thevicinity g ≈ g c . We only show data for moderately largesystem sizes L ≲ g ≈ .
6, which is close to the critical point g c ≃ . g c . In Fig. 7 we perform a data collapse analy-sis plotting the rescaled entanglement gap δξ ln ( L ) versusthe scaling variable ( g − g c ) L / ν . Since we expect that thescaling behaviour is determined by the QSM universalityclass, we fix ν =
1. Due to the logarithmic scaling cor-rections, the data collapse is poor. From section V A oneshould expect f ( ) = π . On the other hand, the data upto L ≲ suggest f ( ) ≈
7, which is quite far from theexpected value f ( ) ≃ .
8. As it is shown in the inset,a very slow drift towards the asymptotic value is visible,compatible with the presence of logarithmic corrections.In conclusion, our analysis suggests that the scaling ofthe entanglement gap can be used to estimate the posi-tion of the quantum critical point, although extractingthe critical exponent ν and the scaling function requiresknowledge of the logarithmic corrections. VI. ENTANGLEMENT GAP AND THEZERO-MODE EIGENVECTOR
In this section we discuss how the vanishing of the en-tanglement gap is reflected in the eigenstate of the corre- i/|A| -1.4-1.2-1-0.8-0.6-0.4-0.20 | A | / ψ L=100L=500L=1000 g=g c g
Let us consider the eigenvector ∣ ψ ⟩ corresponding to thelargest eigenvalue of the spin-spin correlator S A . Thiseigenvector is closely related to that of C A correspond-ing to e , which gives the smallest single-particle ESlevel. Its behaviour is summarised in Fig. 8, showing thecomponents of the eigenvector for different system sizesand in different regions of the phase diagram. We con-sider the bipartition with straight boundary ω y = ω x = / L all the com-ponents decay to zero. Thus, it is convenient to rescaleby ∣ A ∣ / = √ (cid:96) x (cid:96) y (see Fig. 1). We define the flat vector ∣ ⟩ in region A as ∣ ⟩ = √∣ A ∣ ( , , . . . , ) T . (56)It is clear from Fig. 8 in the thermodynamic limit in theordered phase one has that ∣ ψ ⟩ → ∣ ⟩ , up to an irrelevantglobal phase.The structure of ∣ ψ ⟩ for g > g c can be understood asfollows. Deep in the paramagnetic phase the correlationlength is small. In the limit g → ∞ spin-spin correlators1become ultra-local, viz. S nm = δ nm + ε ( δ ∣ n x − m x ∣ , + δ ∣ n y − m y ∣ , ) , (57)with ε vanishing for g → ∞ . In the case ω y =
1, it isstraightforward to determine the eigenvector of (57) cor-responding to the largest eigenvalue in the sector with k y =
0. Due to ω y =
1, the eigenvector is “flat” alongˆ y , and has a non-trivial dependence only on the x co-ordinate. The components of the eigenvector are givenas ψ n x ,n y = ∣ A ∣ / sin ( πn x (cid:96) x ) . (58)The dotted line in Fig. 8 shows the eigenvector ∣ ψ ⟩ for g =
10 and the data are in perfect agreement with (58).Upon approaching the quantum critical point, thezero-mode eigenvector flattens, reflecting that the sys-tem develops ferromagnetic order. To understand that,let us consider the spin correlator (12) in the thermody-namic limit. Upon increasing L , as µ →
0, the correlatordevelops a singularity for k = L one canisolate the contribution of the zero mode as [73] S nm = S ( th ) nm + c √ µ + . . . , (59)where c is a constant. Here the first term is obtainedby setting µ = k =
0. The second contribution in (59)does not depend on n and m , and is divergent in thelimit µ →
0. In this limit one has that the flat vectorbecomes an exact eigenvector of S nm with an eigenvaluethat is proportional to L . However, the decompositionin (59) is not justified because the limit µ → L → ∞ cannot be taken independently, because µ ∝ / L . Figure 8 shows that at the critical point therescaled components of ∣ ψ ⟩ collapse on the same curve.The structure of the eigenvector is not flat. On the otherhand, in the ordered phase, where µ ∝ / L (see Fig. 2)upon increasing L the eigenvector becomes flat. Thissuggests that the decomposition (59) holds if µ decayssufficiently fast for large L . B. An interesting logarithmic correction
In this section we investigate the scaling of the entan-glement gap assuming that the eigenvector ∣ ψ ⟩ is flatalso at the critical point, and that the decomposition inEq. (57) holds. A similar analysis for the massive har-monic chain was presented in Ref. [73]. Here we assumethat S nm can be decomposed as S nm = s L ∣ ⟩⟨ ∣ + S ′ nm (60)and we assume that S ′ nm is negligible. The product P ⋅ S L 〈 | S | 〉 ω y =1 ω x =1/2 (numerical) ω y =1/2 ω x =1/2 (numerical)analytical FIG. 9. Expectation value ⟨ ∣ S ∣ ⟩ of the correlation matrix S (cf. (12)) over the flat vector ∣ ⟩ . Symbols are numericallyexact data for the bipartition with several values of ω x and ω y (see Fig. 1). The dashed dotted line is the analytic result s L . Note that s is obtained by summing (67) and (68).
10 100 1000 L 〈 | P | 〉 L ω y =1 ω x =1/2 ω y =1 ω x =1/4theory FIG. 10. Rescaled expectation value ⟨ ∣ P ∣ ⟩ L over the flatvector ∣ ⟩ . The data are for the straight bipartition with ω y = ω x = / , /
4. Note the logarithmic scale on the x axis. The dashed-dotted line is the analytical result. is thus decomposed as P ⋅ S = s L P ∣ ⟩⟨ ∣ + P ⋅ S ′ , (61)where we suppress the indices n , m to lighten the nota-tion. Consistently with (60), we are going to neglect thesecond term in (61). The matrix P ⋅ S is not hermitian,whereas P and S are hermitian. This means that one hasto introduce right and left eigenvectors. We define twovectors u R and u L as u R = P ∣ ⟩ (62) u L = ∣ ⟩ . (63)It is now straightforward to check that u R and u L arethe right and left eigenvectors of P ⋅ S , respectively. Theeigenvalue is given as e = ⟨ ∣ S ∣ ⟩⟨ ∣ P ∣ ⟩ . (64)2Eq. (64) implies that the problem of calculating the eigen-value e of C A (cf. (21)) is reduced to the simpler problemof calculating the flat-vector expectation values in (64).In the following we are going to calculate ⟨ ∣ S ∣ ⟩ = ∣ A ∣ ∑ n , m ∈ A S nm , (65) ⟨ ∣ P ∣ ⟩ = ∣ A ∣ ∑ n , m ∈ A P nm . (66)Note that (65) has the same form as the spin susceptibil-ity. To obtain (65) and (66), we use the expansion of thespin and momentum correlators discussed in section IV Band section IV C. Importantly, both the thermodynamicand the finite-size contributions in (36) and (40) have tobe taken into account.We start discussing the expectation value ⟨ ∣ S ∣ ⟩ andfirst consider the contribution of the thermodynamic partof the correlator in (37). From (37) we can perform thesums over n and m , and after using the explicit form ofthe spherical parameter (34), taking the limit L → ∞ , weobtain for a bipartition with generic ω x and ω y ⟨ ∣ S ( th ) ∣ ⟩ = √ g c Lπ ω x ω y ∬ ∞−∞ dk x dk y sin ( k x ω x ) sin ( k y ω y ) k x k y ( γ + k x + k y ) . (67)Note that this expectation value grows linearly with L .The constant γ is defined in (35). The integral in (67)depends on the universal low-energy behaviour of theQSM, i.e., at k x , k y →
0, although it is not fully universal.We now show that the finite-size term (38) yields a linearcontribution in L in (65). Indeed, it is straightforward totake the limit L → ∞ in (38) to obtain ⟨ ∣ S ( L ) ∣ ⟩ = √ g c L πω x ω y ∞ ∑ ′ l,l ′ =−∞ ∬ ω x dxdx ′ ∬ ω y dydy ′ e − γ √( l + x − x ′ ) +( l ′ + y − y ′ ) √( l + x − x ′ ) + ( l ′ + y − y ′ ) . (68)Note that the integral in (68) is finite, although the de-nominator in (68) is singular for l = l ′ = ± L in the limit L → ∞ . The accuracy of (67) and (68) is numericallyverified in Fig. 9. The symbols are exact numerical datafor (65), whereas the dashed-dotted lines are the theoret-ical predictions obtained by summing (67) and (68).We now show that, surprisingly, the expectationvalue (66) decays as ln ( L )/ L , i.e., it exhibits a multi-plicative logarithmic correction. The derivation is quitecumbersome, although it requires standard methods suchas Poisson’s summation formula and the Euler-Maclaurinformula. The details are reported in Appendix B. Herewe solely discuss the final result. Similar to (65) one can treat separately the thermodynamic contribution of (66)(cf. (41)) and the finite-size one (cf. (42)). For simplicitywe consider the bipartition with ω x = / p and ω y = / q ,with p, q ∈ N . Clearly, for ω y < ⟨ ∣ P ( th ) ∣ ⟩ = p − ∑ p ′ = q − ∑ q ′ = ∫ / p dk x ∫ / q dk y sin ( π ( k x + p ′ / p )) sin ( π ( k y + q ′ / q )) η p ′ ,q ′ ( k x , k y ) . (69)The function η p ′ ,q ′ ( k x , k y ) reads as η p ′ ,q ′ ( k x , k y ) = π √ g c [ q ( k x + p ′ / p ) + p ( k y + q ′ / q ) + pψ ′ ( + k y + q ′ / q ) + q + k x + p ′ / p + q ( + k x + p ′ / p ) q ( + k x + p ′ / p ) + . . . ] ln ( L ) L . (70)The dots in the square brackets denote terms of higherpowers of 1 /( k x + p ′ / p ) that can be derived systematicallyby using the Euler-Maclaurin formula (see Appendix B).The function ψ ′ ( x ) is the first derivative of the digammafunction ψ ( x ) with respect to its argument [74]. As antic-ipated above, the behaviour as ln ( L )/ L is clearly visiblein (70). As for (67) and (68), it is clear that η p ′ ,q ′ is de-termined by the low-energy part of the dispersion of theQSM.Let us now consider the finite-size contribution (42).Interestingly, as it is clear from (42), the finite-size corre-lator is smooth for ω y < ω x <
1, whereas it exhibitsa singularity if either ω y = ω x =
1, i.e., if the bound-ary between A and its complement is straight. Similarto (69), the singular contribution is ⟨ ∣ P ( L ) ∣ ⟩ = − √ g c π ln ( L ) L . (71)Interestingly, the minus sign in (71) suggests that thepresence of corners increases the prefactor of the logarith-mic correction. Finally, after combining Eqs. (67), (68)and (69), (71) with (64), one obtains that e ∝ ln ( L ) .The prefactor of the logarithmic growth depends on thelow-energy properties of the QSM. As anticipated, byapproximating the zero-mode eigenvector with the flatvector one obtains that δξ decays logarithmically uponincreasing L . However, from (22) one obtains that δξ ∝ /√ ln ( L ) , instead of the correct behaviour as 1 / ln ( L ) established in section V A. VII. CONCLUSIONS
We investigated the entanglement gap δξ in the two-dimensional critical QSM. Our main result is that in theQSM there is a relationship between critical behaviourand vanishing of the entanglement gap.3There are several intriguing directions for future re-search. First, it would be interesting to study the be-haviour of the entanglement gap below the transition,i.e., in the ordered phase. Furthermore, an interestingquestion is how the scenario outlined in this work survivesbeyond the large N limit. This, however, is a very de-manding task because entanglement-related observablescannot be calculated efficiently at finite N . Still, the flat-vector approximation discussed in section VI could begeneralized, at least perturbatively in 1 / N . It would beinteresting to check whether the logarithmic correctionthat is responsible of the vanishing of the entanglementgap persists at finite N . Another natural direction isto understand if the vanishing of the entanglement gapat the critical point is an artifact of the large N limit.The question is whether at finite N a spurious transitionappears, as observed in Ref. [23].It would be also interesting to study the negativityspectrum [86–89] at the quantum phase transition, andin particular the effect of the zero mode. A very in-teresting direction is to understand how the fluctua-tions of the number of particles between the two sub-systems is reflected in the entanglement spectrum andthe entanglement gap. Very recently, the symmetry re-solved entanglement entropies emerged as ideal tools todo that [30, 41, 81, 90–106]. However, an important re-mark is that in the QSM the number of bosons is notconserved, and the symmetry-resolved entanglement en-tropies are not well defined. The particle number con-servation is only enforced on average via the gap equa-tion (2). Still, it should be possible to generalize the QSMto investigate this issue, e.g. by studying spin-anisotropyin the QSM [54]. It would be also important to under-stand how our results can be generalized to long-rangespherical models. Finally, it would be interesting to con-sider higher-dimensional fermionic models. An interest-ing question is whether the area-law violation [107–113]affects the scaling of the entanglement gap. ACKNOWLEDGMENTS
V.A. would like to thank Paola Ruggiero for drawingto our attention Ref. 81 and for discussions. We alsothank Pasquale Calabrese for useful comments on themanuscript. V.A. acknowledges support from the Eu-ropean Research Council under ERC Advanced grant743032 DYNAMINT.
Appendix A: Critical behaviour of the spincorrelator
In this appendix we derive the large L behaviour of thecorrelation function S nm in the QSM. Specifically, weprovide exact expressions for the leading and the firstsubleading terms in powers of 1 / L . The correlator to evaluate is defined as (cf. Eq. (12)) S nm = √ g √ V ∑ k e i k ( n − m ) √ µ + ω k . (A1)The correlation depends only on the distance d = n − m , reflecting translation invariance. Eq. (A1) can berewritten as S nm = √ g π ∫ ∞ dte −( µ + ) t ∏ j = x,y L ∑ k j e − cos ( k j ) t + ik j d j (A2)We now apply Poisson’s summation formula which, for aperiodic function G ( q ) = G ( q + π ) , is stated as1 L ( L − )/ ∑ n =−( L − )/ G ( πnL ) = ∞ ∑ l =−∞ ∫ π − π dq π G ( q ) e iqlL . (A3)The application of (A3) to (A2) yields S nm = √ g √ π ∫ ∞ dte −( µ + ) t ∏ j = x,y ( ∞ ∑ l j =−∞ I l j L + d j ( t )) . (A4)Here I n ( t ) is the modified Bessel function of the firstkind [74]. It is convenient to isolate the terms with l x = l y = ∏ j = x,y ∞ ∑ l j =−∞ I l j L + d j ( t ) =∏ j = x,y ( I d j ( t ) + ∞ ∑ ′ l j =−∞ I l j L + d j ( t )) . (A5)The first term on the right-hand side gives the thermo-dynamic contribution to the correlator S nm , i.e., in thelimit L → ∞ , whereas the other terms are finite-size cor-rections. The prime in the sum is to stress that the termswith l x = l y = L behaviour of (A5). Upon expanding (A5), it is clear thatwe have to derive the asymptotic behaviour of integralsof the type K l,l ′ ( x, x ′ ) = √ g √ π ∫ ∞ e −( µ + ) t I lL + x ( t ) I l ′ L + x ′ ( t ) . (A6)Without loss of generality we can restrict ourselves tothe case with l, l ′ >
0. The generalization to arbitrary l, l ′ is straightforward by using the symmetry of the Besselfunction I − n = I n . It is convenient to change variablesin (A6) to z = t /( lL + x ) , viz. K l,l ′ ( x, x ′ ) = √ g √ Ll + x √ π ∫ ∞ dze −( µ + )( lL + x ) z I lL + x ( z ( lL + x )) I l ′ L + x ′ ( r ( l ′ L + x ′ ) z ) , (A7)where we introduced the ratio r as r = lL + xl ′ L + x ′ . (A8)4We can now perform a saddle point analysis for large Ll + x . For large L , the integral K l,l ′ is determined bythe saddle point t ∗ = ( ( µ + )( + r ) + √ r + ( µ ( µ + ) + ) r + µ ( µ + )( µ + ) r ) . (A9)Finally, a standard calculation yields K l,l ′ = √ g √ lL + x ( ( l ′ L + x ′ )) / √ rπ × e −( l ′ L + x ′ )( r ( + µ ) t − rη ( t )− η ( t r )) g ′ ( t )√ f ( t ) RRRRRRRRRRR t → t ∗ . (A10)Here we defined g ′ ( t ) = ( t + ) / ( r t + ) / (A11) f ( t ) = − r t − t √ r t + + ( µ + ) r − r ( t − ) t √ t + , (A12)and the function η ( t ) as η ( t ) = ( + t ) + ln ( t + ( + t ) ) . (A13)The main ingredient to derive (A10) is the asymptoticbehaviour of the Bessel function I z ( z ) for z → ∞ [74]together with the standard saddle point analysis [114].Since we are interested in the critical behaviour of thecorrelators, it is useful to consider the limit µ →
0, be-cause µ vanishes at criticality. Specifically, we considerthe limit L → ∞ with µ ∝ / L . In this limit we obtainthe expression K l,l ′ ( x, x ′ ) = √ g c π e −√ µ √( lL + x ) +( l ′ L + x ′ ) √( lL + x ) + ( l ′ L + x ′ ) , (A14)where we fixed g = g c . Finally, we now obtain that in thelarge L limit and for µ → S nm is given as S nm = √ g c √ π ∫ ∞ dte −( µ + ) t I n x − m x ( t ) I n y − m y ( t )+ ∞ ∑ ′ l,l ′ =−∞ K l,l ′ ( n x − m x , n y − m y ) , (A15)where K l,l ′ is defined in (A14) and the prime in the sumis to stress that the term with l = l ′ = O( / L ) , whereas the thermodynamic one(first term in (A15)) is O( ) . In (A15) we neglect higherorder corrections in powers of 1 / L . The large L expan-sion for the momentum correlator P nm can be obtainedfrom (A15) by using (15). Appendix B: Derivation of the flat-vectorexpectation value ⟨ ∣ P ∣ ⟩ In this appendix we derive the large L behaviour of theexpectation value of the momentum correlator with theflat vector ⟨ ∣ P ∣ ⟩ . We consider the leading, i.e, the ther-modynamic, as well as the first subleading contribution.The main goal is to show that the expectation value ex-hibits multiplicative logarithmic corrections. Two typesof contributions are present. One originating from thethermodynamic limit of the correlator, whereas the sec-ond one is due to the first subleading term. The latteris present only for a straight boundary between the twosubsystems, and it vanishes if the bipartition has corners.
1. Thermodynamic contribution
Here derive the thermodynamic contribution, which isgiven as ⟨ ∣ P ( th ) ∣ ⟩ , cf. (40). Here ∣ ⟩ is the flat vectorrestricted to region A , i.e, ∣ ⟩ = √∣ A ∣ ( , , . . . , ) , ∣ A ∣ = (cid:96) x (cid:96) y . (B1)The thermodynamic part of the momentum correlatorreads P ( th ) nm = √ gπ ∫ π − π d k e i k ( n − m ) √ µ + ω k . (B2)After performing the sum over n and m in (B2), andafter changing variables to k ′ x = Lω x k x / π and k ′ y = Lω y k y / π , we obtain ⟨ ∣ P ( th ) ∣ ⟩ = √ √ gL ω x ω y ∫ Lω x / dk x ∫ Lω y / dk y sin ( πk x ) sin ( πk y ) sin ( πLω x k x ) sin ( πLω y k y ) × [ µ + − cos ( πLω x k x ) − cos ( πLω y k y )] . (B3)In order to extract the large L behaviour of (B3) it is useful to split the integration domains [ , Lω x / ] and5 [ , Lω y / ] and to write ⟨ ∣ P ( th ) ∣ ⟩ = √ √ gL ω x ω y L / − ∑ l x ,l y = ∫ ( l x + ) ω x l x ω x dk x ∫ ( l y + ) ω y l y ω y dk y sin ( πk x ) sin ( πk y ) sin ( πLω x k x ) sin ( πLω y k y ) [ µ + − cos ( πLω x k x ) − cos ( πLω y k y )] . (B4)We now restrict ourselves to the case with ω x = / p and ω y = / q , with p, q positive integers. After a simple shift of the integration variables as k x → k x − l x ω x and k y → k y − l y ω y , one obtains ⟨ ∣ P ( th ) ∣ ⟩ = √ p q √ gL p − ∑ p ′ = q − ∑ q ′ = L /( p )− ∑ l x = L /( q )− ∑ l y = ∫ / p dk x ∫ / q dk y sin ( π ( k x + l x + p ′ / p )) sin ( π ( k y + l y + q ′ / q )) sin ( pπL ( k x + l x + p ′ / p )) sin ( qπL ( k y + l y + q ′ / q ))× [ µ + − cos ( pπL ( k x + l x + p ′ / p )) − cos ( qπL ( k y + l y + q ′ / q ))] . (B5)We now focus on the behaviour at the quantum phasetransition. We set g = g c , µ = γ /( L ) , and we ex- pand (B5) in the limit L → ∞ , using the periodicity ofthe sine function. This yields ⟨ ∣ P ( th ) ∣ ⟩ = √ g c π L p − ∑ p ′ = q − ∑ q ′ = L /( p )− ∑ l x = L /( q )− ∑ l y = ∫ / p dk x ∫ / q dk y sin ( π ( k x + p ′ / p )) sin ( π ( k y + q ′ / q ))( k x + l x + p ′ / p ) ( k y + l y + q ′ / q ) ×[ γ π + p ( k x + l x + p ′ / p ) + q ( k y + l y + q ′ / q ) ] . (B6)Importantly, as a result of the large L limit, Eq. (B6)depends only on the low-energy part of the dispersion ofthe QSM, although it contains non-universal information. To proceed we determine the large L behaviour of thesum over l x , l y in (B6), i.e., of the function η p ′ ,q ′ ( k x , k y ) defined as η p ′ ,q ′ ( k x , k y ) = √ g c π L L /( p )− ∑ l x = L /( q )− ∑ l y = √ γ π + p ( k x + l x + p ′ / p ) + q ( k y + l y + q ′ / q ) ( k x + l x + p ′ / p ) ( k y + l y + q ′ / q ) . (B7)The asymptotic behaviour of η p,q in the limit L → ∞ canbe obtained by using the Euler-Mclaurin formula. Givena function f ( x ) this is stated as x ∑ x = x f ( x ) = ∫ x x f ( x ) dx + f ( x ) + f ( x ) + f ′ ( x ) − f ′ ( x ) + . . . (B8)Here the dots denote terms with higher derivatives of f ( x ) calculated at the integration boundaries x and x , that can be derived to arbitrary order. To proceed, wefirst isolate the term with either l x = l y = l x = l y = L behaviour of η p ′ ,q ′ as η , which is given as η = √ g c π [ q ( k x + p ′ / p ) + p ( k y + q ′ / q ) ] ln ( L ) L . (B9)In the derivation of (B9) we neglected the boundaryterms in (B8) because they are subleading. We are6now left with the sums over l x ∈ [ , L /( p )] and l y ∈[ , L /( q )] in (B7). These can be evaluated again by us-ing (B8). We first apply (B8) to the sum over l x andobtain two contributions. The first one is obtained afterevaluating the integral in (B8) at x = L /( p ) . After ex-panding the result for L → ∞ , we find the contribution η given as η = L /( q ) ∑ l y = p √ g c π ( k y + l y + q ′ / q ) ln ( L ) L . (B10)Note the term ln ( L )/ L in (B10). The sum over l y in (B10) can be performed exactly to obtain in the large L limit η = √ g c π pψ ′ ( + k y + q ′ / q ) ln ( L ) L . (B11)Here ψ ′ ( z ) is the first derivative of the digamma function ψ ( z ) with respect to its argument [74]. The second con-tribution is obtained by evaluating the integral in (B8) at x =
1. The remaining sum over l y cannot be evaluatedanalytically. However, one can again treat the sum over l y with the Euler-Mclaurin formula (B8). After neglect-ing the boundary terms in (B8), which are subleadingfor large L , and after evaluating the integral in (B8) at x = L /( q ) , we obtain the contribution η as η = √ g c π q + k x + p ′ / p ln ( L ) L . (B12)To obtain the full contribution of the sum over l x in (B7)we now have to consider the effect of the boundary termsin (B8). Before doing that we check the accuracy of (B11)and (B12) by defining J = √ g c π L ∫ L /( p ) dl xL /( q )− ∑ l y = √ γ π + p ( k x + l x + p ′ / p ) + q ( k y + l y + q ′ / q ) ( k x + l x + p ′ / p ) ( k y + l y + q ′ / q ) . (B13) J is obtained by neglecting the terms with either l x = l y = l x in (B8) with an integral(see first term in (B8)), treating the sum over l y exactly.In Fig. 11 we show (J − η − η ) L versus 1 / L . For large L the data show a linear behaviour attaining a finite valuein the limit L → ∞ . This shows that the leading orderterm ∝ ln ( L )/ L of J is fully captured by η + η , theremaining contribution being ∝ / L , which we neglect.Having discussed the contribution which derives fromapproximating the sum over l x in (B7) with the integralin (B8), we finally focus on the effect of the boundaryterms in (B8). Let us consider the first boundary term(first term in the second row in (B8)). The contribution ( J - η - η ) L p=q=1 p’=q’=0linear fit FIG. 11. Check of the asymptotic behaviour in the large L limit of J (cf. (B13)). The circles are numerical data for ( J − η − η ) L , with η and η as defined in (B11) and (B12).The dashed-dotted line is a linear fit. Data are for p = q = p ′ = q ′ = J − η − η ∝ / L for L → ∞ . as ln ( L )/ L is obtained by fixing l x =
1, other contribu-tions are subleading. After performing the sum over l y ,one obtains the first boundary contribution η b as η b = √ g c π q ( + k x + p ′ / p ) ln ( L ) L . (B14)In a similar way the second boundary term (last termin (B8)) gives η b = √ g c π q ( + k x + p ′ / p ) ln ( L ) L . (B15)Note that boundary terms in (B8) are expected to besmall. Specifically, the k -th term is suppressed as 1 /( k + ) !. The final result for η ( k x , k y , p, p ′ , q, q ′ ) is obtained byputting together Eqs. (B9), (B11), (B12), (B14), (B15)to obtain η p ′ ,q ′ ( k x , k y ) = η + η + η + η b + η b . (B16)In Fig. 12 we check the accuracy of (B16), showing thefunction η p ′ ,q ′ for fixed values of q =
1, which correspondsto a straight partition between the two subsystems, and p = /
2. For all values of p ′ , q ′ and k x , k y that we consider η p ′ ,q ′ is well described by (B16).
2. Finite-size contribution
In this section we derive the leading behaviour in thelarge L limit of ⟨ ∣ P ( L ) ∣ ⟩ . Interestingly, we show thatin the presence of a straight boundary between the twosubsystems the expectation value behaves as ⟨ ∣ P ( L ) ∣ ⟩ ∝ ln ( L )/ L . On the other hand, in the presence of corners,7
10 100 L η p ’ , q ’ ( k x , k y ) p’=0, q’=0 k x =0.05 k y =1p’=0, q’=0 k x =0.1 k y =0.5p’=1, q’=0, k x =0.05 k y =0.1p’=1, q’=0, k x =0.05 k y =0.05analytical FIG. 12. Large L behaviour of the function η p ′ ,q ′ ( k x , k y ) defined in (B7). Here we consider a bipartition with ω x = / p and ω y = / q (see Fig. 1). We fix q = q ′ = p ′ = , the multiplicative logarithmic correction is absent. Thefinite-size correlator to calculate reads as P ( L ) nm = − √ g c π ∞ ∑ ′ l,l ′ =−∞ e −√ µ √( lL + n x − m x ) +( l ′ L + n y − m y ) × [ [( lL + n x − m x ) + ( l ′ L + n y − m y ) ] / + √ µ ( lL + n x − m x ) + ( l ′ L + n y − m y ) ] . (B17)Crucially, if ω x < ω y <
1, the denominators in (B17)are never singular. This implies that the logarithmic cor-rection is not present, which can be straightforwardlychecked numerically. Let us now consider the situationwith ω x < ω y =
1. The other case with ω x = ω y < L → ∞ for l = l ′ = ±
1. We numericallyobserve that only the first term in (B17) gives rise to asingular behaviour. Thus, we neglect the second termand fix l =
0, obtaining ⟨ ∣ P ( L ) ∣ ⟩ ≃ − √ g c πL ω x ∞ ∑ ′ l ′ =−∞ Lω x ∑ n x ,m x = L − ∑ n y ,m y = e −√ µ √( n x − m x ) +( l ′ L + n y − m y ) (( n x − m x ) + ( l ′ L + n y − m y ) ) / . (B18)Only the differences n x − m x and n y − m y appear in (B18). Thus, it is convenient to change variables to x = n x − m x and y = n y − m y , to obtain ⟨ ∣ P ( L ) ∣ ⟩ ≃ − √ g c πL ω x ∞ ∑ ′ l ′ =−∞ Lω x ∑ x =− Lω x L − ∑ y =−( L − ) ( Lω x + − ∣ x ∣)( L − ∣ y ∣) e −√ µ √ x +( l ′ L + y ) ( x + ( l ′ L + y ) ) / . (B19)Again, the singular behaviour occurs for x ≈ y ≈ − lL , with l ′ = ±
1. In this limit we can neglect the ex-ponential in (B20) because it is regular. Thus, we obtain ⟨ ∣ P ( L ) ∣ ⟩ ≃ − √ g c πL ω x ∞ ∑ ′ l ′ =−∞ Lω x ∑ x =− Lω x L − ∑ y =−( L − ) ( Lω x + − ∣ x ∣)( L − ∣ y ∣)( x + ( l ′ L + y ) ) / . (B20)To proceed, we consider the case l = l = − x in (B20) to x > x → − x . We also restrict to y < y <
0. We now have ⟨ ∣ P ( L ) ∣ ⟩ ≃ √ g c πL ω x Lω x ∑ x = L − ∑ y = ( Lω x + − x )( y − L )( x + ( L − y ) ) / . (B21) Now the strategy is to treat the sum (B21) by using theEuler-Mclaurin formula (B8). For instance, one can firstapply (B8) to the sum over x and obtain that the leadingterm in the large L limit is obtained by evaluating theintegral in (B8) at ω x L . One can also verify that theboundary terms in (B8) can be neglected. A straightfor-8
100 1000 L - 〈 | P ( L ) | 〉 L exact1/(g c1/2 π ) ln(L) FIG. 13. Large L behaviour of ⟨ ∣ P ( L ) ∣ ⟩ (cf. (B18) forits definition). The symbols are numerical data obtained byusing (B17). The dashed-dotted line is the analytical re-sult (B22). All the results are for the bipartition with ω x = / ω y = ward calculation gives the final result ⟨ ∣ P ( L ) ∣ ⟩ ≃ − √ g c π ln ( L ) L , (B22)where the contribution of l = − [1] J. Eisert, M. Cramer, and M. B. Plenio, Rev. Mod.Phys. , 277 (2010).[2] L. Amico, R. Fazio, A. Osterloh, and V. Vedral, Rev.Mod. Phys. , 517 (2008).[3] P. Calabrese, J. Cardy, and B. Doyon, Journal ofPhysics A: Mathematical and Theoretical , 500301(2009).[4] N. Laflorencie, Physics Reports , 1 (2016).[5] I. Peschel, M. Kaulke, and . Legeza, Annalen der Physik , 153 (1999).[6] M.-C. Chung and I. Peschel, Phys. Rev. B , 4191(2000).[7] I. Peschel, Journal of Statistical Mechanics: Theory andExperiment , P06004 (2004).[8] I. Peschel and V. Eisler, Journal of Physics A: Mathe-matical and Theoretical , 504003 (2009).[9] S. R. White, Phys. Rev. Lett. , 2863 (1992).[10] U. Schollw¨ock, Annals of Physics , 96 (2011).[11] H. Li and F. D. M. Haldane, Phys. Rev. Lett. ,010504 (2008).[12] R. Thomale, D. P. Arovas, and B. A. Bernevig, Phys.Rev. Lett. , 116805 (2010).[13] A. M. L¨auchli, E. J. Bergholtz, J. Suorsa, andM. Haque, Phys. Rev. Lett. , 156404 (2010).[14] M. Haque, O. Zozulya, and K. Schoutens, Phys. Rev.Lett. , 060401 (2007).[15] R. Thomale, A. Sterdyniak, N. Regnault, and B. A.Bernevig, Phys. Rev. Lett. , 180502 (2010).[16] M. Hermanns, A. Chandran, N. Regnault, and B. A.Bernevig, Phys. Rev. B , 121309(R) (2011).[17] A. Chandran, M. Hermanns, N. Regnault, and B. A.Bernevig, Phys. Rev. B , 205136 (2011).[18] X.-L. Qi, H. Katsura, and A. W. W. Ludwig, Phys.Rev. Lett. , 196402 (2012).[19] Z. Liu, E. J. Bergholtz, H. Fan, and A. M. L¨auchli,Phys. Rev. B , 045119 (2012).[20] A. Sterdyniak, A. Chandran, N. Regnault, B. A.Bernevig, and P. Bonderson, Phys. Rev. B , 125308 (2012).[21] J. Dubail, N. Read, and E. H. Rezayi, Phys. Rev. B , 115321 (2012).[22] J. Dubail, N. Read, and E. H. Rezayi, Phys. Rev. B , 245310 (2012).[23] A. Chandran, V. Khemani, and S. L. Sondhi, Phys.Rev. Lett. , 060501 (2014).[24] F. Pollmann, A. M. Turner, E. Berg, and M. Oshikawa,Phys. Rev. B , 064439 (2010).[25] A. M. Turner, F. Pollmann, and E. Berg, Phys. Rev. B , 075102 (2011).[26] B. Bauer, L. Cincio, B. Keller, M. Dolfi, G. Vidal,S. Trebst, and A. Ludwig, Nature Communications (2014), 10.1038/ncomms6137.[27] D. Poilblanc, Phys. Rev. Lett. , 077202 (2010).[28] J. I. Cirac, D. Poilblanc, N. Schuch, and F. Verstraete,Phys. Rev. B , 245134 (2011).[29] G. De Chiara, L. Lepori, M. Lewenstein, and A. San-pera, Phys. Rev. Lett. , 237208 (2012).[30] V. Alba, M. Haque, and A. M. L¨auchli, Phys. Rev.Lett. , 227201 (2012).[31] M. A. Metlitski and T. Grover, “Entanglement entropyof systems with spontaneously broken continuous sym-metry,” (2011), arXiv:1112.5166 [cond-mat.str-el].[32] V. Alba, M. Haque, and A. M. L¨auchli, Journal ofStatistical Mechanics: Theory and Experiment ,P08011 (2012).[33] V. Alba, M. Haque, and A. M. L¨auchli, Phys. Rev.Lett. , 260403 (2013).[34] L. Lepori, G. De Chiara, and A. Sanpera, Phys. Rev.B , 235107 (2013).[35] A. J. A. James and R. M. Konik, Phys. Rev. B ,241103 (2013).[36] F. Kolley, S. Depenbrock, I. P. McCulloch,U. Schollw¨ock, and V. Alba, Phys. Rev. B ,144426 (2013).[37] L. Rademaker, Phys. Rev. B , 144419 (2015).[38] F. Kolley, S. Depenbrock, I. P. McCulloch, U. Schollw¨ock, and V. Alba, Phys. Rev. B ,104418 (2015).[39] I. Fr´erot and T. Roscilde, Phys. Rev. Lett. , 190401(2016).[40] P. Calabrese and A. Lefevre, Phys. Rev. A , 032329(2008).[41] A. M. L¨auchli, (2013), arXiv:1303.0741 [cond-mat.stat-mech].[42] V. Alba, P. Calabrese, and E. Tonni, Journal of PhysicsA: Mathematical and Theoretical , 024001 (2017).[43] J. Cardy, “The entanglement gap in cfts, talk at the kitpconference ”closing the entanglement gap: Quantuminformation, quantum matter, and quantum fields”.”(2015).[44] A. Bayat, H. Johannesson, S. Bose, and P. Sodano, Na-ture Communications (2014), 10.1038/ncomms4784.[45] T. T. Truong and I. Peschel, Zeitschrift f¨ur Physik BCondensed Matter , 119 (1989).[46] G. D. Giulio and E. Tonni, Journal of Statistical Me-chanics: Theory and Experiment , 033102 (2020).[47] R. Lundgren, J. Blair, M. Greiter, A. L¨auchli, G. A.Fiete, and R. Thomale, Phys. Rev. Lett. , 256404(2014).[48] C. Lhuillier and G. Misguich, in High Magnetic Fields (Springer Berlin Heidelberg, 2002) pp. 161–190.[49] A. J. Beekman, L. Rademaker, and J. van Wezel, Sci-Post Phys. Lect. Notes , 11 (2019).[50] A. Wietek, M. Schuler, and A. M. Luchli, “Studyingcontinuous symmetry breaking using energy level spec-troscopy,” (2017), arXiv:1704.08622 [cond-mat.str-el].[51] G. Obermair, in