Closure of the entanglement gap at quantum criticality: The case of the Quantum Spherical Model

The study of entanglement spectra is a powerful tool to detect or elucidate universal behaviour in quantum many-body systems. We investigate the scaling of the entanglement (or Schmidt) gap $\delta\xi$, 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 that the Schmidt gap vanishes at the critical point, although only logarithmically. For a system on a torus and the half-system bipartition, the entanglement gap vanishes as $\pi^2/\ln(L)$, with $L$ the linear system size. The entanglement gap is nonzero in the paramagnetic phase and exhibits a faster decay in the ordered phase. The rescaled gap $\delta\xi\ln(L)$ exhibits a crossing for different system sizes at the transition, 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 to a non-trivial structure of the eigenvector, whereas in the ordered phase it is flat. We also show that the vanishing of the entanglement gap at criticality can be qualitatively but not quantitatively captured by neglecting the structure of the zero-mode eigenvector.

In this work we investigate the ES in critical twodimensional quantum many-body systems. We focus on * swald@pks.mpg.de the lowest laying entanglement gap δξ defined as where ξ 0 and ξ 1 are the lowest and the first excited ES level, respectively. The behaviour of the entanglement gap at quantum critical points has not been thoroughly addressed, except for one-dimensional systems [5-7, 13, 29, 30, 34, 45, 46]. Several exact results suggest that at one-dimensional quantum critical points δξ vanishes. For instance, in CFT systems δξ decays logarithmically as ∝ 1 ln(L) with the subsystem's length L [40]. Similar scaling is found in corner transfer matrix calculations [45] (see also [8] for a review). Higher-dimensions are far less explored. Interestingly, it has been argued that the closing of the entanglement gap does not necessarily signal critical behaviour [23]. Similar conclusions have been reached by considering the ES of a bipartition in momentum space [47]. Still, the ES can be useful to distinguish different phases of matter. This is the case for systems that exhibit order by breaking of a continuous symmetry [31]. It has been suggested that deep in the ordered phase the lower part of the ES contains the fingerprints of symmetry breaking, being reminiscent of the so-called Anderson tower-of-states [48][49][50]. This has been verified by analytical calculations in the quantum rotor model [31], numerical simulations in the two-dimensional Bose-Hubbard model in the superfluid phase [33] (see also [39]), and also in two-dimensional Heisenberg models on the square [36] and on the kagome lattice [38]. A signature of the tower-of-states scenario is that the gaps in the lower part of the ES decay as a power-law with the subsystem volume, with multiplicative logarithmic cor- We also define the ratios ω x(y) = x(y) L. We mostly consider the the case with ωy = 1.
rections [31]. Higher ES levels are expected to exhibit a much slower decay [31,33,37]. The behaviour of the entanglement gap upon approaching the critical point has not been investigated thoroughly.
Here we address this issue in the quantum spherical model [51][52][53][54][55] (QSM). The QSM is a paradigmatic manybody system in which the effects of strongly interacting degrees of freedom may be studied at a considerably low cost, as the model can be mapped to free bosons subject to an additional external constraint. Despite its simplicity it exhibits several salient features of realistic quantum many-body systems. For instance, its classical version served as a testing ground for the theory of critical phenomena and finite size scaling [56]. In two dimensions the QSM exhibits a standard paramagnetic (disordered) phase and a ferromagnetic (ordered) one, which are separated by a second order quantum phase transition. The universality class of the transition is that of the threedimensional classical O(N ) vector model [57] in the large N limit [52,53,58]. Surprisingly, entanglement properties of the QSM are rather unexplored, although there is recent interest [59][60][61]. We should stress that although the results that we are going to derive for the ES cannot be considered general, they certainly represent an interesting case study, and can be useful to understand the generic behaviour of ES in quantum many-body systems.
Here we consider a two-dimensional lattice of linear size L with periodic boundary conditions in both directions. The typical bipartitions that we use are reported in Fig. 1. Figure 1 (a) shows a bipartition with a straight boundary between A and its complement, with A spanning the full lattice along theŷ direction. This is not the case in Fig. 1 (b), where the boundary has a corner. The effect of corners in the scaling of the entanglement entropies is nontrivial, and it has been studied intensely in the last decade [4,[62][63][64][65][66][67][68][69]. Since the QSM is mappable to free bosons, entanglement-related observables can be calculated from the two-point correlations functions [8].
Here we show that δξ (cf. (2)) is nonzero in the paramagnetic 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 quantum critical point, in the case of straight boundary the entanglement gap vanishes as π 2 ln(L). However, we show that logarithmic corrections are present, which make it difficult to robustly verify the finite-size scaling of δξ. We also show that the behaviour of the entanglement gap is reflected in the zero-mode eigenvector of the spin-spin correlation matrix. As the transition is approached from the paramagnetic side, the eigenvector flattens, meaning that all its components become equal. This reflects the presence of a zero mode. Exactly at criticality, the eigenvector is not flat in the thermodynamic limit, due to the presence of strong fluctuations, whereas it is flat in the ordered phase. Interestingly, we show that by neglecting the structure of the eigenvector at the critical point, i.e., by approximating the eigenvector with the flat vector, we obtain that δξ = A ln(L), which accounts for the vanishing of the entanglement gap, although it is not quantitatively accurate. We clarify how the behaviour as A ln(L) arises from some interesting multiplicative logarithmic corrections in the expectation values of the QSM correlators with the flat vector. Interestingly, the constant A depends only on low-energy properties of the model and on the geometry of the bipartition.
The manuscript is organised as follows. In section II we introduce the QSM and its phase diagram. In section III we define the quantities of interest. In section IV we 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 subsections IV B and IV C we derive the finite-size scaling of the spin and momentum correlation functions, respectively. In section V we investigate the critical behaviour of δξ. Our prediction is discussed in section V A, and it is compared against numerical results in section V B. We describe the behaviour of δξ across the phase diagram of the QSM in subsection V B 1, whereas we address the vanishing of δξ and its finite-size scaling in subsections V B 2 and V B 3, respectively. In section VI we discuss how the entanglement gap is related to the zero-mode eigenvector of the correlator, which we introduce in subsection VI A. In subsection VI B we show that by assuming that the eigenvector is flat at criticality one can qualitatively explain the vanishing of the entanglement gap. We conclude in section VII. In Appendix A we report the derivation of the finite-size scaling of the correlation functions in the QSM. In Appendix B we derive the expectation values of the correlators with the flat vector.

II. QUANTUM SPHERICAL MODEL
The QSM [52][53][54] on a two dimensional cubic lattice of linear size L and volume V = L 2 is defined by the Hamiltonian Here, n = (n x , n y ) ∈ [1, . . . , L] 2 denotes a generic lattice site, and ⟨n, m⟩ a lattice bond joining two nearestneighbour sites. J > 0 is the ferromagnetic exchange constant and we choose J = 1 in the remainder of the paper. The canonically conjugated variables s n and p n satisfy the standard bosonic commutation relations We refer to p n as momentum variable, and to the parameter g as quantum coupling as the model reduces to the famous classical spherical model [70,71] in the limit g → 0. The Lagrange multiplier µ is called spherical parameter and fixes the spherical constraint, i.e. n ⟨s 2 n ⟩ = V.
This means that all allowed configurations of the QSM are located around the sphere in configuration space that is defined by Eq. (5). Critical properties of the QSM are determined through the self-consistent behaviour of µ [53]. The two dimensional QSM does not exhibit a finite temperature phase transition [70,71], although it possesses a ground-state transition, i.e., at T = 0 [52][53][54].
We now briefly review how to diagonalise the Hamiltonian (3) and describe its critical behaviour. First, we exploit the translational invariance of the model by performing a Fourier transform as Here the sum over k = (k x , k y ) runs in the first Brillouin zone k i = 2π L j, with j ∈ [−L 2, L 2] integer. The Hamiltonian (3) in Fourier space reads with the single-particle dispersion relation In order to fully diagonalise (7) we introduce bosonic ladder operators b k and b † k obeying standard bosonic commutation relations viz.
with the parameter α 2 k = g 2Λ −1 k . In terms of these ladder operators, the Hamiltonian (7) is diagonal and reads Entanglement-related properties of Gaussian systems such as the QSM stem from the two-point correlation functions ⟨s n s m ⟩ and ⟨p n p m ⟩. In equilibrium at zero temperature T = 0, the eigenmodes k of the system are occupied according to (11), we can thus immediately derive the twopoint correlation functions [54] S nm = ⟨s n s m ⟩ = Importantly, from (12) and (13) one obtains the relation which allows to relate the critical behaviour of the spin correlator to that of the momentum correlator. From (12), one can rewrite the spherical constraint (5) as This equation is also called gap equation [72] and implies that only the average number of bosons is fixed. From the finite-size expressions (12) (13) and (16), the thermodynamic limit L → ∞ is obtained in the usual way by replacing A crucial observation is that the correlator (12) and the spherical parameter (16) exhibit a singularity for k = 0, due to the zero mode. We anticipate that this will play an important role in the behaviour of the entanglement gap. This contribution of the zero mode to the entanglement entropy was previously investigated focusing on the harmonic chain [73]. We now summarise the zero-temperature critical behaviour of the QSM. In two dimensions the model exhibits a second order phase transition at a critical value g c . For g < g c the ground-state of (3) exhibits magnetic order. At g > g c the ground-state is paramagnetic. The behaviour of the QSM is determined by the scaling of the spherical parameter µ. In the thermodynamic limit, in the paramagnetic phase one has that µ is finite and nonzero. On the other hand, one has µ = 0 at the critical point, and in the ordered phase. The value of g c can be determined analytically. In the thermodynamic limit the spherical constraint (16) is rewritten as with the complete elliptic integral [74] K(x) = The critical coupling g c follows by imposing the condition µ = 0. This yields The different phases of the model correspond to different finite-size scaling behaviours of µ. In the paramagnetic phase one has µ = O(1) in the limit L → ∞. At the critical point one can show that µ = O(1 L 2 ), whereas in the ordered phase µ = O(1 L 4 ) (see section IV). These behaviours are numerically illustrated in Fig. 2. The universality class of the ground state transition [53] is that of the large-N vector model in three dimensions, as expected from general renormalisation group arguments. Critical properties of the large-N vector model have been characterised analytically [56] and finite-size corrections have also been investigated [75][76][77][78][79].

III. ENTANGLEMENT SPECTRA AND ENTANGLEMENT GAPS
Here we are interested in the ground-state entanglement spectrum of the QSM, focussing on the two bipartitions depicted in Fig. 1. The lattice, with periodic boundary conditions, is divided into two regions A andĀ. Region A is of size A = x × y and we define the corresponding aspect ratios ω x = x L and ω y = y L, with 0 ≤ ω x,y ≤ 1. In Fig. 1 (a) the subsystem A spans the full lattice along thê y direction implying that the boundary between the two subsystems A andĀ is straight. This case corresponds to ω y = 1. In Fig. 1 (b), the boundary presents a corner and is thus not straight. The presence of corners has striking consequences for entanglement entropies, giving rise to sub-leading universal logarithmic corrections [62][63][64][65][66]80]. The effects of corners in the scaling of the ES have not been investigated yet. For the case of a straight boundary with periodic boundary conditions the momentum k y is a good quantum number for the correlation matrices (12) and (13), and for the ES. This will be exploited in section V to reduce the computation of the ES of the QSM to that of an effective one-dimensional model. This dimensional reduction has been employed to study symmetry-resolved entanglement entropies [81]. This rather simple observation will also allow to obtain analytically the scaling of the entanglement gap at the critical point, by exploiting corner transfer matrix results [5][6][7]45].
We now review the calculation of entanglement-related quantities in the QSM. Since the QSM is essentially mappable to a free bosonic model (see section II), its entanglement properties are derived from the two-point correlation functions (12) and (13) (see Ref. [8] for a review). The crucial ingredient is the correlation matrix C restricted to the subsystem A, viz.
with S A and P A being the correlation matrices defined in (12) and (13), restricted to the subsystem A. Since in the remainder we mostly consider the restricted correlation 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 subsystem A is a quadratic operator and is written as [8] Here H A is the so-called entanglement Hamiltonian, k are the single-particle ES levels, b k are free-bosonic operators and Z ensures the normalisation of the reduced density matrix Trρ A = 1. The spectrum {e k } k=1,..., A of the correlation matrix C A is simply related to that of H A viz.
The normalisation factor Z is obtained as The ES, i.e., the spectrum of the entanglement Hamiltonian H A , is obtained by filling the single-particle levels k in all the possible ways. To construct the ES, it is convenient to introduce the bosonic occupation numbers α k = 0, 1, . . . in the levels k . The generic ES level ξ({α k }) is written as The eigenvalues e k satisfy the constraint e k > 1 4, implying that k > 0. Clearly, the lowest ES level ξ 0 corresponds to the vacuum state with α k = 0 for all k. Let us order the k as 1 ≤ 2 ≤ ⋅ ⋅ ⋅ ≤ A . The first excited ES level is obtained by populating the smallest single particle level 1 . Thus, the lowest entanglement gap δξ (Schmidt gap) is defined as Here we focus on δξ, although one can define higher gaps [82].

IV. FINITE-SIZE CRITICAL CORRELATORS IN THE QSM
As explained in section III, entanglement-related observables, and also the entanglement gap, in the QSM are entirely encoded in the two-point correlation functions (12) and (13). In the following sections we derive the finitesize behaviour of these two-point correlation functions.
In section IV A we discuss the gap equation (16). In sections IV B and IV C we the focus on the spin and momentum correlators respectively. For the classical spherical model similar results were obtained [75,83].

A. Spherical parameter
Here we derive the finite-size scaling of the spherical parameter µ at the quantum phase transition. The result is not new [56] but it is a useful initiation for the discussion of the correlators. To treat the sum over k in (12) we observe that the following identity holds where the prime in the sum indicates that the l = 0 contribution is removed, and I ν are modified Bessel functions of the first kind [74]. To derive (27), we introduce an auxiliary integration [72] over t to represent the term (µ + ω k ) −1 2 , then we employ Poisson's summation formula. Further details are reported in Appendix A. The first term in the brackets in (27) does not depend explicitly 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 for large L. In this limit the integral in (27) can be treated by using a saddle point approximation.
In order to use (16), we decompose the diagonal cor-relator S nn as with the thermodynamic contribution corresponding to the term I 0 (t 2 ) 2 in (27). The remaining terms in (27) are collected in S (L) nn . 1 After expanding the square in (27), we observe that S (L) nn is written as In order to extract the large L behaviour of (30) we employ a standard saddle point approximation. The calculation is straightforward and details are reported in Appendix A. A striking simplification occurs at the critical point and in the ordered phase, where µ → 0. One can verify numerically that at the thermodynamical critical point µ ∝ 1 L 2 . This is expected because µ ∝ m 2 = 1 ξ 2 corr , with m the mass of the theory and ξ corr the correlation length, and at the critical point ξ corr ∝ L. In the limit µ → 0, one obtains the surprisingly elegant result (see Appendix A) Interestingly, in (31) the first term is of one-dimensional nature, and it is obtained by isolating the terms with either l = 0 or l ′ = 0 in the sum in (30). In the second term in (31) the scaling as µ ∝ 1 L 2 gives rise to a nontrivial 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 converge quickly. Double sums as in (31) appear often in lattice calculations, and have been investigated in the past [75,76,83]. In some cases they can be expressed in terms of generalised Riemann zeta functions [84]. Using Eqs. (29) and (31) in the gap equation (16) at criticality yields The integral in (32) has to be considered carefully due to a ∝ 1 L contribution in the µ → 0 limit which can be extracted as [61] dk where the dots denote subleading terms in 1 L. The second term in (33) is the singular term that determines the critical behaviour of three-dimensional QSM at the thermal phase transition [61]. This is not surprising because the universality class of the quantum phase transition in two dimensions is the same [52,53]. Based on the expected finite-size scaling µ ∝ 1 L 2 it is convenient to define where the constant γ 2 is to be determined and the factor 2 is for later convenience. We substitute the Ansatz (34) in the gap equation (32), and use the spherical constraint in the thermodynamic limit (16) at criticality, where µ = 0. This yields where the first term is (33) and the other two are obtained from (30). Eq. (35) can be solved numerically to obtain the universal constant γ 2 ≃ 1.51196. Note that Eq. (35) has also been found in the context of the large N limit of the three dimensional N −vector model [56,77]. The behaviour of µ in the different regions of the phase diagram of the QSM and the accuracy of (34) are verified in Fig. 2 where we show the numerical solution of Eq. (16). In the paramagnetic region for g > g c one has µ = O(1). At the critical point and in the ferromagnetic phase µ → 0 in the limit L → ∞. The dashed-dotted line is the analytic result (34) with γ 2 obtained from (35).
Below the critical point we expect µ ∝ 1 L 4 [56], which is confirmed by the fit (dashed line).

B. Spin-spin correlation function Snm
We now discuss the finite-size scaling of the spin-spin correlation function (12) at the quantum critical point. We only discuss the final result, reporting the details of the derivation in Appendix A. First, one can again decompose the correlator as with the thermodynamic contribution As in Eq. (29) there is an implicit dependence on L via µ. The finite-size part has the surprisingly simple form Here we defined The prime in the sum means that the term (l, l ′ ) = (0, 0) has been removed. Again, Eq. (38) holds in the limit L → ∞ and µ → 0. The general expression, which is valid also in the paramagnetic phase, is reported in Appendix 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 due to translation invariance. Moreover, one has that S nm is periodic along the two directions, i.e., it is invariant under n y −m y → n y −m y ±L and n x −m x → n x −m x ±L. This is enforced by the infinite sums over l, l ′ . For a bipartition with straight boundary between the two subsystems ( Fig. 1 (a)) the invariance under n y −m y → n y −m y ±L remains true also for the correlator restricted to A. Finally, nm exhibits an interesting singularity structure. For ω y = 1 the denominator in Eq. (38) is singular, whereas it is regular for ω y < 1. Specifically, the terms with l = 0 and l ′ = ±1 in (38) exhibit a singularity in the limit n x − m x → 0 and n y − m y → ±L. On the other hand, terms with l ′ > 1 or l > 1 in (38) are not singular.
The same singularity appears if ω x = 1 and ω y < 1. We anticipate that these singularities will give rise to multiplicative logarithmic corrections in the expectation value of the correlators that we will show in section VI.

C. Momentum correlation function Pnm
The same finite-size analysis as in section IV B can be carried out for the momentum correlator P nm (cf. (13)). Following the decomposition with the finite-size part P (L) nm has the same structure as (38), and it reads This expression is obtained from the spin-spin correlator, cf. Eq (38), by using (15). As for (38), the finite-size term (42) is singular if subsystem A spans the full lattice in one of the two directions, i.e., if ω x = 1 or ω y = 1. For ω y = 1 the singularity occurs for l = 0 and l ′ = ±1 in the limit n x − m x → 0 and n y − m y → ±L. Note that the first term in Eq. (42) exhibits a stronger singularity than the second one.

V. CRITICAL BEHAVIOUR OF THE ENTANGLEMENT GAP
We now discuss the critical behaviour of the entanglement gap δξ. In subsection V A, by using a dimensional reduction, we provide an exact result for the case of a smooth boundary between the subsystems. In subsection V B we discuss numerical results. We first discuss the behaviour of the entanglement gap across the phase diagram of the QSM in subsection V B 1. In subsection V B 2 we show that at the critical point the entanglement gap vanishes logarithmically with the system size. Finally, in subsection 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 = 1 (see Fig. 1 a). Periodic boundary conditions along theŷ direction imply that the momentum k y is a good quantum number for the correlation matrix C A (cf. (21)) restricted to subsystem A. Moreover, translation invariance implies that by performing a Fourier transform along theŷ direction the Hamiltonian (3) can be written as the sum of L decoupled quadratic one-dimensional systems [8]. This dimensional reduction is effective for any free system, and has been recently employed to study the so-called symmetry-resolved entanglement entropies [81]. The fact that k y is a good quantum number implies that the correlation matrix C A has a block structure with each block corresponding to a different k y viz.
It is straightforward to diagonalise a given block with fixed k y by imposing that the eigenvectors of C A are also eigenvectors of the momentum alongŷ with the given eigenvalue k y . Since we are interested only in the largest eigenvalue e 1 of C A a further simplification occurs. As the critical behaviour is associated with the formation of a uniform magnetization, it is natural to expect that e 1 is in the sector with k y = 0. This can be readily checked numerically. Thus, in the following we restrict the calculation to k y = 0. By imposing that the eigenvectors of C A are "flat" alongŷ, i.e., they do not depend on y, the problem is reduced to the diagonalisation of the reduced correlation matrix where we defined the reduced spin and momentum correlators as Eqs. (45) and (46) depend only on the coordinates n x −m x along thex direction, and subsystem A is the interval of length x . Here α kx 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 massive harmonic chain with frequency Ω = √ 2µ [8]. The full ES of the massive harmonic chain for the bipartition in two semi-infinite chains has been calculated by using the corner transfer matrix approach [8]. The reduced density matrix ρ A , up to a trivial renormalisation, is written as with the corner transfer matrix Hamiltonian where β j are bosonic ladder operators. Here K(x) is the complete elliptic integral of the first kind (see Eq. (19)). The parameter κ is given in terms of Ω as [81] Eq. (47) holds if A is the half-infinite line. In this limit, as it is clear from Eq. (48), the single-particle ES levels are equally spaced [8] with spacing . To determine the finite-size scaling of the entanglement gap δξ we use the fact that for L → ∞ at criticality µ ∝ 1 L 2 (see Eq. (34)). By substituting (34) in the corner transfer matrix results (48) and (49), we obtain that in the large L limit δξ decays logarithmically with L as Note the dependence on the universal constant γ 2 . To derive (50), one can also observe that close to the critical point, on the paramagnetic side, Eq. (48) gives δξ = π 2 ln(ξ corr ). Eq. (50) then follows from standard scaling arguments. A similar decay of the entanglement gap as in (50) is obtained for critical one-dimensional systems [8], both fermionic and bosonic ones. An important remark is that the corner transfer matrix calculation is valid for the bipartition in two semi-infinite systems, which implies that there is only one boundary between the two subsystems, in contrast with the bipartitions Fig. 1, which contain two boundaries because we are using periodic boundary conditions. Despite that, as it will be clear in section V B, Eq. (50) gives the leading behaviour for large L of δξ. We anticipate that a logarithmic subleading term as O(ln −2 (L)), which is missing in Eq. (50), is present. From Eqs. (22) and (50) one obtains that the eigenvalue e 1 of C A is given as Importantly, the missing O(ln −2 (L)) term in (50) will give a O(ln(L)) contribution in (51).

B. Numerical results
In this section we discuss numerical results confirming the validity of the logarithmic scaling of the entanglement gap at criticality. We provide numerical evidence that the prefactor of the logarithmic decay obeys the standard finite-size scaling behaviour. For instance, it exhibits a crossing for different system sizes at the critical point. However, logarithmic corrections are present, and a precise finite-size scaling analysis is very challenging.

Overview
Before discussing the scaling of δξ at the critical point, it is useful to focus on its behaviour across the phase diagram of the QSM, see Fig. 3. The figure shows δξ as a function of g for several system sizes L. The entanglement spectrum is calculated for the bipartition with straight boundary, i.e., ω y = 1 and ω x = 1 2 (see Fig. 1 (a)). In Fig. 3 the solid line is δξ as obtained by using the value of the spherical constraint µ in the thermodynamic limit L → ∞ (cf. (32)). This yields µ = O(1) in the paramagnetic phase and µ = 0 in the ferromagnetic phase and at criticality (g ≤ g c ). The thermodynamic entanglement gap is obtained by substituting the thermodynamic value of µ in the finite-size expressions for the correlators (cf. (12) and (13)) and taking the limit L → ∞ after. This procedure gives the correct thermodynamic behaviour of δξ, at least away from the critical point. Although we use the finite-size expressions for the correlators, we observe that δξ converges quickly to its thermodynamic value. This is expected because the behaviour of the QSM is determined by the scaling of µ.
In the ordered phase and at the critical point the spin  Fig. 1 (a) with ωx = 1 2 and ωy = 1.
e1 is plotted versus linear size L. In the ordered phase (diamonds) we observe a fast increase with L, whereas in the paramagnetic phase e1 = O(1). Note the logarithmic divergence as e1 ∝ ln 2 (L) at the critical point at gc. The dashed dotted line is a fit to e1 = 1 π 4 ln 2 (8L γ2) + A0 + A1 ln(8L γ2), with A0, A1 fitting parameters.
correlator (12) diverges due to the zero mode. Thus, we regularise the zero-mode by fixing µ = 10 −6 . As it is clear from Fig. 3, this analysis, although it is not rigorous, suggests that δξ = 0 in the ordered phase, whereas µ is finite and nonzero in the paramagnetic phase.
Let us now discuss the finite-size behaviour of δξ. In the paramagnetic phase, i.e. g > g c , the approach to the thermodynamic limit is exponential, which is expected because the model is massive. For g < g c , i.e., in the ferromagnetic phase, the data suggest a vanishing gap. The scaling of the entanglement gap in magnetically ordered phases has been investigated extensively [31,33,36,38,39]. For instance, in Ref. [31] it was predicted that in the presence of continuous symmetry breaking in generic dimension d, δξ should decay as δξ ∝ (L d−1 ln(L)) −1 .
In d = 1 one recovers the logarithmic decay as 1 ln(L), reflecting the absence of symmetry breaking. In d > 1 Eq. (52) yields a "fast" power-law decay with a multiplicative logarithmic correction. An important remark is that Eq. (52) applies to the gaps in the lower part of the entanglement spectrum, i.e., the part which is related to the Anderson tower of states. Gaps in the higher part of the entanglement spectrum are expected to vanish logarithmically [31].  Fig. 4. The line is a fit to A0 + A1 ln(8L γ2), with A0, A1 fitting parameters. The fit gives A1 ≃ 0.041. The inset shows e1 obtained by using µ = γ 2 2 (2L) and fixing γ2 = 8. e1 is plotted versus ln 2 (L). The line is a fit to A ′ 0 + 1 π 4 ln 2 (L).
of δξ we, equivalently, consider the scaling of the largest eigenvalue e 1 of C A . We show our numerical results for e 1 in Fig. 4 as a function of L (note the logarithmic scale on the x-axis). To highlight the different scaling as compared to other regions of the phase diagram, we report also data in the paramagnetic phase (square symbols) and in the ferromagnetic phase (diamonds). Within the ordered phase e 1 increases faster than logarihmically. In the paramagnetic region e 1 exhibits a mild increase for small L, saturating at L → ∞. This is a consequence of the finite correlation length in the paramagnetic phase. A dramatically different behaviour is visible at criticality (circles), for which we report data up to L ∼ 40000. 2 Interestingly, for moderately large L the behaviour of δξ is compatible with a logarithmic increase, although Eq. (51) suggests a ln 2 (L) scaling. This should be attributed to the presence of a sub-leading logarithmic term ln(L) (cf. (51)). A fit to A 2 ln 2 (8L γ 2 )+A 0 +A 1 ln(8L γ 2 ) (dashed-dotted line) gives A 2 ≃ 0.01, which is in good agreement with the prediction 1 π 4 . One also obtains A 1 ≃ 0.04 and A 0 ≃ 0.16. Note that A 0 ≃ 1 6, as predicted by (51).
To further corroborate our results, in Fig. 5 we show e 1 − 1 π 4 ln(8L γ 2 ) versus L using a logarithmic scale on the x-axis. The data are the same as in Fig. 4. The continuous line is a fit to with A 0 and A 1 fitting constants. The logarithmic behaviour is perfect. Note that this logarithmic term is 2 Note that since ωy = 1, we can use dimensional reduction to attain large system sizes (see section V A). not predicted by (51). Its origin could be attributed to the fact that the corner transfer matrix result is obtained for the semi-infinite system, i.e., the biparititon with one boundary. It is interesting to investigate the dependence on γ 2 of the constant A 1 in (53). In the inset in Fig. 5 we show e 1 obtained by fixing µ = γ 2 2 (2L 2 ) with γ 2 = 8 in (12) and (13). In the inset e 1 is plotted versus ln 2 (L). The dashed-dotted line is a fit to 1 π 4 ln 2 (L) + A ′ 0 . The perfect linear behaviour suggests that the subleading logarithmic term is absent or its prefactor is small. A fit to 1 π 4 ln 2 (8L γ 2 ) + A ′ 0 + A ′ 1 ln(L) gives A ′ 1 ≈ 0.0007. It would be interesting to investigate this behaviour more systematically. One possible scenario is that the prefactor of the logarithmic term is of the form A ′ 1 = ln(γ 2 8).

Finite-size scaling analysis
Having established the logarithmic vanishing of δξ at the critical point, it is natural to investigate its behaviour in the vicinity of the quantum phase transition. A natural idea is that δξ obeys standard finite-size scaling [85] δξ where the dots stand for scaling corrections, f (x) is a scaling function and ν is the exponent that governs the divergence of the correlation length at the critical point. For the QSM one has ν = 1 [53] . The scaling function f (x) is determined by the universality class of the QSM, and, in principle, can be calculated. Under the assumption that the f (x) is analytic, one can expand (54) near g c to obtain δξ ln(L) = f (0) + (g − g c )L 1 ν + . . . , From the analysis in section V A one should expect f (0) = π 2 ≃ 9.8. Eq. (55) implies that the data for δξ for different system sizes should exhibit a crossing at g c . This crossing method for the entanglement gap has been used to detect a quantum phase transition in a system of coupled one-dimensional models [35]. However, since δξ has logarithmic corrections, one should expect strong limitations, as we are going to show. The scaling Ansatz (54) implies that by plotting the rescaled gap δξ ln(L) as a function of the scaling variable (g − g c )L 1 ν one should observe a data collapse for different system sizes, provided that scaling corrections can be neglected.
Our finite-size data for δξ as a function of g for several system sizes L are shown in Fig. 6 focussing on the vicinity g ≈ g c . We only show data for moderately large system sizes L ≲ 200. Clearly, the data exhibit a crossing at g ≈ 9.6, which is close to the critical point g c ≃ 9.67826. This is quite remarkable because logarithmic corrections are present. In fact, we observe that even including larger system sizes, it is challenging to obtain a more precise estimate of g c . In Fig. 7 we perform a data collapse analysis plotting the rescaled entanglement gap δξ ln(L) versus the scaling variable (g−g c )L 1 ν . Since we expect that the scaling behaviour is determined by the QSM universality class, we fix ν = 1. Due to the logarithmic scaling corrections, the data collapse is poor. From section V A one should expect f (0) = π 2 . On the other hand, the data up to L ≲ 10 4 suggest f (0) ≈ 7, which is quite far from the expected value f (0) ≃ 9.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 of the entanglement gap can be used to estimate the position of the quantum critical point, although extracting the critical exponent ν and the scaling function requires knowledge of the logarithmic corrections.

VI. ENTANGLEMENT GAP AND THE ZERO-MODE EIGENVECTOR
In this section we discuss how the vanishing of the entanglement gap is reflected in the eigenstate of the corre- Eigenvector's components are rescaled by A 1 2 . On the xaxis i is a label. In the ordered phase for g < gc the eigenvector becomes flat in the thermodynamic limit, in contrast with the behaviour at the critical point at gc, and in the paramagnetic phase.
lation matrix that corresponds to the zero mode. Moreover, we show that assuming a flat structure of the zeromode eigenvector at criticality allows one to capture qualitatively the logarithmic vanishing of the entanglement gap. Within this approximation the vanishing of δξ is related to some interesting multiplicative logarithmic corrections in the correlators. Finally, the result suggests that the presence of corners in the bipartition affects the vanishing of the gap.

A. The zero-mode eigenvector
Let us consider the eigenvector ψ 0 ⟩ corresponding to the largest eigenvalue of the spin-spin correlator S A . This eigenvector is closely related to that of C A corresponding to e 1 , which gives the smallest single-particle ES level. Its behaviour is summarised in Fig. 8, showing the components of the eigenvector for different system sizes and in different regions of the phase diagram. We consider the bipartition with straight boundary ω y = 1 and ω x = 1 2 (see Fig. 1 a). Upon increasing L all the components decay to zero. Thus, it is convenient to rescale by A 1 2 = x y (see Fig. 1). We define the flat vector 1⟩ in region A as It is clear from Fig. 8 in the thermodynamic limit in the ordered phase one has that ψ 0 ⟩ → 1⟩, up to an irrelevant global phase. The structure of ψ 0 ⟩ for g > g c can be understood as follows. Deep in the paramagnetic phase the correlation length is small. In the limit g → ∞ spin-spin correlators become ultra-local, viz.
with ε vanishing for g → ∞. In the case ω y = 1, it is straightforward to determine the eigenvector of (57) corresponding to the largest eigenvalue in the sector with k y = 0. Due to ω y = 1, the eigenvector is "flat" alonĝ y, and has a non-trivial dependence only on the x coordinate. The components of the eigenvector are given as The dotted line in Fig. 8 shows the eigenvector ψ 0 ⟩ for g = 10 and the data are in perfect agreement with (58). Upon approaching the quantum critical point, the zero-mode eigenvector flattens, reflecting that the system develops ferromagnetic order. To understand that, let us consider the spin correlator (12) in the thermodynamic limit. Upon increasing L, as µ → 0, the correlator develops a singularity for k = 0 which encodes the critical behaviour of the QSM. In the limit of large L one can isolate the contribution of the zero mode as [73] where c is a constant. Here the first term is obtained by setting µ = 0 and by replacing the sum in (12) with an integral and the second term is the contribution of the zero mode k = 0. The second contribution in (59) does not depend on n and m, and is divergent in the limit µ → 0. In this limit one has that the flat vector becomes an exact eigenvector of S nm with an eigenvalue that is proportional to L. However, the decomposition in (59) is not justified because the limit µ → 0 and the limit L → ∞ cannot be taken independently, because µ ∝ 1 L 2 . Figure 8 shows that at the critical point the rescaled components of ψ 0 ⟩ collapse on the same curve. The structure of the eigenvector is not flat. On the other hand, in the ordered phase, where µ ∝ 1 L 4 (see Fig. 2) upon increasing L the eigenvector becomes flat. This suggests that the decomposition (59) holds if µ decays sufficiently fast for large L.

B. An interesting logarithmic correction
In this section we investigate the scaling of the entanglement gap assuming that the eigenvector ψ 0 ⟩ is flat also at the critical point, and that the decomposition in Eq. (57) holds. A similar analysis for the massive harmonic chain was presented in Ref. [73]. Here we assume that S nm can be decomposed as and we assume that S ′ nm is negligible. The product P ⋅ S  12)) over the flat vector 1⟩. Symbols are numerically exact data for the bipartition with several values of ωx and ωy (see Fig. 1). The dashed dotted line is the analytic result s0L. Note that s0 is obtained by summing (67) and (68). is thus decomposed as where we suppress the indices n, m to lighten the notation. Consistently with (60), we are going to neglect the second term in (61). The matrix P ⋅ S is not hermitian, whereas P and S are hermitian. This means that one has to introduce right and left eigenvectors. We define two vectors u R and u L as It is now straightforward to check that u R and u L are the right and left eigenvectors of P ⋅ S, respectively. The eigenvalue is given as Eq. (64) implies that the problem of calculating the eigenvalue e 1 of C A (cf. (21)) is reduced to the simpler problem of calculating the flat-vector expectation values in (64).
In the following we are going to calculate Note that (65) has the same form as the spin susceptibility. To obtain (65) and (66), we use the expansion of the spin and momentum correlators discussed in section IV B and section IV C. Importantly, both the thermodynamic and the finite-size contributions in (36) and (40) have to be taken into account. We start discussing the expectation value ⟨1 S 1⟩ and first consider the contribution of the thermodynamic part of the correlator in (37). From (37) we can perform the sums over n and m, and after using the explicit form of the spherical parameter (34), taking the limit L → ∞, we obtain for a bipartition with generic ω x and ω y Note that this expectation value grows linearly with L.
The constant γ 2 is defined in (35). The integral in (67) depends on the universal low-energy behaviour of the QSM, 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 linear contribution in L in (65). Indeed, it is straightforward to take the limit L → ∞ in (38) to obtain Note that the integral in (68) is finite, although the denominator in (68) is singular for l = 0 and l ′ = ±1 (see section IV B). It is straightforward to integrate the contributions (67) and (68) numerically. We conclude that the expectation value (65) grows linearly with L in the limit L → ∞. The accuracy of (67) and (68) is numerically verified in Fig. 9. The symbols are exact numerical data for (65), whereas the dashed-dotted lines are the theoretical predictions obtained by summing (67) and (68). We now show that, surprisingly, the expectation value (66) decays as ln(L) L, i.e., it exhibits a multiplicative logarithmic correction. The derivation is quite cumbersome, although it requires standard methods such as Poisson's summation formula and the Euler-Maclaurin formula. The details are reported in Appendix B. Here we 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 simplicity we consider the bipartition with ω x = 1 p and ω y = 1 q, with p, q ∈ N. Clearly, for ω y < 1 the boundary between the two subsystems is not straight, i.e., it has a corner (see Fig. 1 b). One obtains The function η p ′ ,q ′ (k x , k y ) reads as The dots in the square brackets denote terms of higher powers of 1 (k x +p ′ p) that can be derived systematically by using the Euler-Maclaurin formula (see Appendix B). The function ψ ′ (x) is the first derivative of the digamma function ψ(x) with respect to its argument [74]. As anticipated above, the behaviour as ln(L) L is clearly visible in (70). As for (67) and (68), it is clear that η p ′ ,q ′ is determined by the low-energy part of the dispersion of the QSM.
Let us now consider the finite-size contribution (42). Interestingly, as it is clear from (42), the finite-size correlator is smooth for ω y < 1 and ω x < 1, whereas it exhibits a singularity if either ω y = 1 or ω x = 1, i.e., if the boundary between A and its complement is straight. Similar to (69), the singular contribution is Interestingly, the minus sign in (71) suggests that the presence of corners increases the prefactor of the logarithmic correction. Finally, after combining Eqs. (67), (68) and (69), (71) with (64), one obtains that e 1 ∝ ln(L). The prefactor of the logarithmic growth depends on the low-energy properties of the QSM. As anticipated, by approximating the zero-mode eigenvector with the flat vector one obtains that δξ decays logarithmically upon increasing L. However, from (22) one obtains that δξ ∝ 1 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 twodimensional critical QSM. Our main result is that in the QSM there is a relationship between critical behaviour and vanishing of the entanglement gap.
There are several intriguing directions for future research. First, it would be interesting to study the behaviour of the entanglement gap below the transition, i.e., in the ordered phase. Furthermore, an interesting question is how the scenario outlined in this work survives beyond the large N limit. This, however, is a very demanding task because entanglement-related observables cannot be calculated efficiently at finite N . Still, the flatvector approximation discussed in section VI could be generalized, at least perturbatively in 1 N . It would be interesting to check whether the logarithmic correction that is responsible of the vanishing of the entanglement gap persists at finite N . Another natural direction is to understand if the vanishing of the entanglement gap at the critical point is an artifact of the large N limit. The question is whether at finite N a spurious transition appears, as observed in Ref. [23].
It would be also interesting to study the negativity spectrum [86][87][88][89] at the quantum phase transition, and in particular the effect of the zero mode. A very interesting direction is to understand how the fluctuations of the number of particles between the two subsystems is reflected in the entanglement spectrum and the entanglement gap. Very recently, the symmetry resolved entanglement entropies emerged as ideal tools to do that [30,41,81,[90][91][92][93][94][95][96][97][98][99][100][101][102][103][104][105][106]. However, an important remark is that in the QSM the number of bosons is not conserved, and the symmetry-resolved entanglement entropies are not well defined. The particle number conservation is only enforced on average via the gap equation (2). Still, it should be possible to generalize the QSM to investigate this issue, e.g. by studying spin-anisotropy in the QSM [54]. It would be also important to understand how our results can be generalized to long-range spherical models. Finally, it would be interesting to consider higher-dimensional fermionic models. An interesting question is whether the area-law violation [107][108][109][110][111][112][113] affects the scaling of the entanglement gap. In this appendix we derive the large L behaviour of the correlation function S nm in the QSM. Specifically, we provide exact expressions for the leading and the first subleading terms in powers of 1 L. The correlator to evaluate is defined as (cf. Eq. (12)) The correlation depends only on the distance d = n − m, reflecting translation invariance. Eq. (A1) can be rewritten as +ikj dj (A2) We now apply Poisson's summation formula which, for a periodic function G(q) = G(q + 2π), is stated as The application of (A3) to (A2) yields .
(A4) Here I n (t) is the modified Bessel function of the first kind [74]. It is convenient to isolate the terms with l x = l y = 0 in (A4), viz.
The first term on the right-hand side gives the thermodynamic contribution to the correlator S nm , i.e., in the limit L → ∞, whereas the other terms are finite-size corrections. The prime in the sum is to stress that the terms with l x = l y = 0 is removed. We now derive the large L behaviour of (A5). Upon expanding (A5), it is clear that we have to derive the asymptotic behaviour of integrals of the type (A6) Without loss of generality we can restrict ourselves to the case with l, l ′ > 0. The generalization to arbitrary l, l ′ is straightforward by using the symmetry of the Bessel function I −n = I n . It is convenient to change variables in (A6) to z 2 = t 2 (lL + x), viz.
(A9) Finally, a standard calculation yields Here we defined and the function η(t) as The main ingredient to derive (A10) is the asymptotic behaviour 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 the correlators, it is useful to consider the limit µ → 0, because µ vanishes at criticality. Specifically, we consider the limit L → ∞ with µ ∝ 1 L 2 . In this limit we obtain the expression where we fixed g = g c . Finally, we now obtain that in the large L limit and for µ → 0 the correlator S nm is given as where K l,l ′ is defined in (A14) and the prime in the sum is to stress that the term with l = l ′ = 0 has been removed. In (A15) one can recognize the two contributions in (37) and (38). Note that the finite-size term (second term in (A15)) is O(1 L), whereas the thermodynamic one (first term in (A15)) is O(1). In (A15) we neglect higher order corrections in powers of 1 L. The large L expansion for the momentum correlator P nm can be obtained from (A15) by using (15).
Appendix B: Derivation of the flat-vector expectation value ⟨1 P 1⟩ In this appendix we derive the large L behaviour of the expectation value of the momentum correlator with the flat vector ⟨1 P 1⟩. We consider the leading, i.e, the thermodynamic, as well as the first subleading contribution.
The main goal is to show that the expectation value exhibits multiplicative logarithmic corrections. Two types of contributions are present. One originating from the thermodynamic limit of the correlator, whereas the second one is due to the first subleading term. The latter is present only for a straight boundary between the two subsystems, and it vanishes if the bipartition has corners.

Thermodynamic contribution
Here derive the thermodynamic contribution, which is given as ⟨1 P (th) 1⟩, cf. (40). Here 1⟩ is the flat vector restricted to region A, i.e, The thermodynamic part of the momentum correlator reads After performing the sum over n and m in (B2), and after changing variables to k ′ x = Lω x k x π and k ′ y = Lω y k y π, we obtain Lωx k x sin 2 π Lωy k y In order to extract the large L behaviour of (B3) it is useful to split the integration domains [0, Lω x 2] and [0, Lω y 2] and to write We now restrict ourselves to the case with ω x = 1 p and ω y = 1 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 We now focus on the behaviour at the quantum phase transition. We set g = g c , µ = γ 2 2 (2L 2 ), and we ex-pand (B5) in the limit L → ∞, using the periodicity of the sine function. This yields Importantly, as a result of the large L limit, Eq. (B6) depends only on the low-energy part of the dispersion of the QSM, although it contains non-universal information.
To proceed we determine the large L behaviour of the sum over l x , l y in (B6), i.e., of the function η p ′ ,q ′ (k x , k y ) defined as The asymptotic behaviour of η p,q in the limit L → ∞ can be obtained by using the Euler-Mclaurin formula. Given a function f (x) this is stated as Here the dots denote terms with higher derivatives of f (x) calculated at the integration boundaries x 1 and x 2 , that can be derived to arbitrary order. To proceed, we first isolate the term with either l x = 0 or l y = 0 in (B7). The remaining sum after fixing l x = 0 or l y = 0 can be treated with (B8). We define this contribution to the large L behaviour of η p ′ ,q ′ as η 0 , which is given as In the derivation of (B9) we neglected the boundary terms in (B8) because they are subleading. We are now left with the sums over l x ∈ [1, L (2p)] and l y ∈ [1, L (2q)] in (B7). These can be evaluated again by using (B8). We first apply (B8) to the sum over l x and obtain two contributions. The first one is obtained after evaluating the integral in (B8) at x 2 = L (2p). After expanding the result for L → ∞, we find the contribution η 1 given as 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 Here ψ ′ (z) is the first derivative of the digamma function ψ(z) with respect to its argument [74]. The second contribution is obtained by evaluating the integral in (B8) at x 1 = 1. The remaining sum over l y cannot be evaluated analytically. However, one can again treat the sum over l y with the Euler-Mclaurin formula (B8). After neglecting the boundary terms in (B8), which are subleading for large L, and after evaluating the integral in (B8) at x 2 = L (2q), we obtain the contribution η 2 as To obtain the full contribution of the sum over l x in (B7) we now have to consider the effect of the boundary terms in (B8). Before doing that we check the accuracy of (B11) and (B12) by defining J is obtained by neglecting the terms with either l x = 0 or l y = 0 in (B7), which were treated in (B9), and by approximating the sum over 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 − η 1 − η 2 )L versus 1 L. For large L the data show a linear behaviour attaining a finite value in the limit L → ∞. This shows that the leading order term ∝ ln(L) L of J is fully captured by η 1 + η 2 , the remaining contribution being ∝ 1 L, which we neglect.
Having discussed the contribution which derives from approximating the sum over l x in (B7) with the integral in (B8), we finally focus on the effect of the boundary terms in (B8). Let us consider the first boundary term (first term in the second row in (B8)). The contribution as ln(L) L is obtained by fixing l x = 1, other contributions are subleading. After performing the sum over l y , one obtains the first boundary contribution η b1 as In a similar way the second boundary term (last term in (B8)) gives η b2 = 2 3 √ g c π 3 q (1 + k x + p ′ p) 3
In Fig. 12 we check the accuracy of (B16), showing the function η p ′ ,q ′ for fixed values of q = 1, which corresponds to a straight partition between the two subsystems, and p = 1 2. For all values of p ′ , q ′ and k x , k y that we consider η p ′ ,q ′ is well described by (B16).

Finite-size contribution
In this section we derive the leading behaviour in the large L limit of ⟨1 P (L) 1⟩. Interestingly, we show that in the presence of a straight boundary between the two subsystems the expectation value behaves as ⟨1 P (L) 1⟩ ∝ ln(L) L. On the other hand, in the presence of corners, Large L behaviour of the function η p ′ ,q ′ (kx, ky) defined in (B7). Here we consider a bipartition with ωx = 1 p and ωy = 1 q (see Fig. 1). We fix q = 1 considering q ′ = 0 and p ′ = 0, 1 (empty and full symbols, respectively). Symbols are numerical results. Dashed-dotted lines are the asymptotic behaviours in (B16).
the multiplicative logarithmic correction is absent. The finite-size correlator to calculate reads as P (L) nm = − Crucially, if ω x < 1 and ω y < 1, the denominators in (B17) are never singular. This implies that the logarithmic correction is not present, which can be straightforwardly checked numerically. Let us now consider the situation with ω x < 1 and ω y = 1. The other case with ω x = 1 and ω y < 1 can be treated similarly. A singularity appears in the limit L → ∞ for l = 0 and l ′ = ±1. We numerically observe that only the first term in (B17) gives rise to a singular behaviour. Thus, we neglect the second term and fix l = 0, obtaining ⟨1 P (L) 1⟩ ≃ − (B19) Again, the singular behaviour occurs for x ≈ 0 and y ≈ −lL, with l ′ = ±1. In this limit we can neglect the exponential in (B20) because it is regular. Thus, we obtain To proceed, we consider the case l = 1 and it is clear that the contribution from l = −1 is the same. We can restrict the sum over x in (B20) to x > 0 because of the symmetry x → −x. We also restrict to y < 0 because the singularity in (B20) occurs for y < 0. We now have ⟨1 P (L) 1⟩ ≃ Now the strategy is to treat the sum (B21) by using the Euler-Mclaurin formula (B8). For instance, one can first apply (B8) to the sum over x and obtain that the leading term in the large L limit is obtained by evaluating the integral in (B8) at ω x L. One can also verify that the boundary terms in (B8) can be neglected. A straightfor- where the contribution of l = −1 in (B20) has been taken into account. The validity of (B22) is numerically confirmed in Fig. 13.