Universal terms in the overlap of the ground state of the spin-1/2 XXZ chain with the N\'eel state

TL;DR

This paper investigates universal finite-size scaling terms in the overlap between the N9el state and the ground state of the spin-1/2 XXZ chain, connecting boundary conformal field theory predictions with exact asymptotic calculations.

Contribution

It provides an exact asymptotic analysis of the overlap's order one term using a finite-size determinant formula, confirming boundary CFT predictions and extracting finite-size correction exponents.

Findings

Confirmed boundary conformal field theory predictions for the overlap.

Derived exact asymptotic expressions for the order one term.

Determined the exponent of the leading finite-size correction.

Abstract

We analyze universal terms that appear in the large system size scaling of the overlap between the N\'eel state and the ground state of the spin-1/2 XXZ chain in the antiferromagnetic regime. In a critical theory, the order one term of the asymptotics of such an overlap may be expressed in terms of $g$-factors, known in the context of boundary conformal field theory. In particular, for the XXZ model in its gapless phase, this term provides access to the Luttinger parameter. In its gapped phase, on the other hand, the order one term simply reflects the symmetry broken nature of the phase. In order to study the large system size scaling of this overlap analytically and to compute the order one term exactly, we use a recently derived finite-size determinant formula and perform an asymptotic expansion. Our analysis confirms the predictions of boundary conformal field theory and enables us to determine the exponent of the leading finite-size correction.

Equations

$$H_L=\sum _{j=1}^L\left({\sigma }_j^x{\sigma }_{j+1}^x+{\sigma }_j^y{\sigma }_{j+1}^y+\Delta {\sigma }_j^z{\sigma }_{j+1}^z\right),\Delta =\cos \gamma ,0\le \gamma <\pi .$$

$$K=\frac{\pi }{2(\pi -\gamma )}.$$

$${\mathcal{O}}_L={\left|⟨N⟩{\Psi }_L\right|}^2.$$

$$\mathcal{O}_{L}\;=\;\sqrt{K}e^{-\alpha L}\Big{[}1\;+\;AL^{-\delta}\;+\;\ldots\Big{]}\,.$$

$$\mathcal{O}_{L}\;=\;\frac{1}{\sqrt{2}}e^{-\alpha L}\Big{[}1\;+\;\mathcal{O}\big{(}\ln^{-1}(L)\big{)}\Big{]}\,.$$

$$\mathcal{O}_{L}\;=\;\frac{1}{2}e^{-\alpha L}\Big{[}1\;+\;\mathcal{O}\big{(}e^{-cL}\big{)}\Big{]}\,.$$

$$e^{-\text{uptau}{H}_L}|a⟩\underset{\text{uptau}\to \infty }{\sim }{e}^{-\text{uptau}{E}_0}|{\Psi }_L⟩⟨{\Psi }_L⟩a.$$

$${\mathcal{O}}_L=\underset{\text{uptau}\to \infty }{\lim }\frac{⟨N⟩e^{-\text{uptau}{H}_L}N⟨N⟩e^{-\text{uptau}{H}_L}N}{⟨N⟩e^{-2\text{uptau}{H}_L}N}.$$

$${\mathcal{O}}_L=\underset{\text{uptau}\to \infty }{\lim }\frac{{\left[Z_{\mathrm{cyl}}(L,\text{uptau})\right]}^2}{Z_{\mathrm{cyl}}(L,2\text{uptau})}.$$

$$H_{\text{uptau}}=E\text{uptau}+e+\frac{\pi v}{\text{uptau}}\left[\frac{(S^z{)}^2}{2K}+\sum _{k=1}^{\infty }k\left(a_k^{†}{a}_k+\frac{1}{2}\right)\right].$$

$$Z_{\rm cyl}(L,\uptau)\;=\;{\rm Tr}\big{\{}e^{-LH_{\uptau}}\big{\}}\,,$$

$$Z_{\mathrm{cyl}}(L,\text{uptau})=\frac{e^{-EL\text{uptau}-eL+\pi vL/(24\text{uptau})}}{{\prod }_{n=1}^{\infty }\left(1-e^{-\pi vLn/\text{uptau}}\right)}\sum _{n\in \mathbb{Z}}{e}^{-\pi vL{n}^2/(2K\text{uptau})}.$$

$$Z_{\mathrm{cyl}}(L,\text{uptau})=\sqrt{K}\frac{e^{-EL\text{uptau}-eL+\pi \text{uptau}/(24vL)}}{{\prod }_{n=1}^{\infty }\left(1-e^{-\pi \text{uptau}n/(vL)}\right)}\sum _{n\in \mathbb{Z}}{e}^{-\pi K\text{uptau}{n}^2/(2vL)}.$$

$${\mathcal{O}}_L=\underset{\text{uptau}\to \infty }{\lim }\frac{[Z(L,\text{uptau}){]}^2}{Z(L,2\text{uptau})}=\sqrt{K}{e}^{-eL}.$$

$${\left[\frac{\sinh ({\lambda }_j+\frac{i\gamma }{2})}{\sinh ({\lambda }_j-\frac{i\gamma }{2})}\right]}^L$$

$${\left[\frac{{\lambda }_j+\frac{i}{2}}{{\lambda }_j-\frac{i}{2}}\right]}^L$$

$${\left[\frac{\sin ({\lambda }_j+\frac{i\eta }{2})}{\sin ({\lambda }_j-\frac{i\eta }{2})}\right]}^L$$

$$\lambda_{1}\;\simeq\;\left\{\begin{array}[]{cl}\displaystyle\frac{\gamma}{L}&\text{(planar)}\\[12.91663pt] \displaystyle\frac{1}{L}&\text{(isotropic)}\\[12.91663pt] \displaystyle\frac{\pi}{2K\mspace{-1.5mu}(k)}\frac{\eta}{L}&\text{(axial)}\end{array}\right.\,,\qquad\lambda_{L/4}\;\simeq\;\left\{\begin{array}[]{cl}\displaystyle\frac{\gamma}{\pi}\ln\!\Big{(}\frac{2L}{\pi}\Big{)}&\text{(planar)}\\[12.91663pt] \displaystyle\frac{1}{\pi}\ln\!\Big{(}\frac{2L}{\pi}\Big{)}&\text{(isotropic)}\\[12.91663pt] \displaystyle\frac{\pi}{2}\,-\,\frac{\pi}{2kK\mspace{-1.5mu}(k)}\frac{\eta}{L}&\text{(axial)}\end{array}\right.\,.$$

$$\rho_{L}(\lambda)\;=\;\frac{1}{L}\sum_{j=1}^{L/2}\delta(\lambda-\lambda_{j})\;=\;\frac{1}{L}\sum_{j=1}^{L/4}\big{[}\delta(\lambda-\lambda_{j})+\delta(\lambda+\lambda_{j})\big{]}\,,$$

$$\rho(\lambda)\;=\;\lim_{L\to\infty}\rho_{L}(\lambda)\;=\;\left\{\begin{array}[]{lll}\displaystyle\frac{1}{2\gamma\cosh\!\big{(}\frac{\pi\lambda}{\gamma}\big{)}}\,,&\lambda\in\mathbb{R}&\text{(planar)}\\[12.91663pt] \displaystyle\frac{1}{2\cosh(\pi\lambda)}\,,&\lambda\in\mathbb{R}&\text{(isotropic)}\\[12.91663pt] \displaystyle\frac{1}{2\pi}\sum_{n\in\mathbb{Z}}\frac{\cos(2n\lambda)}{\cosh(\eta n)}\,,&\lambda\in[-\tfrac{\pi}{2},\tfrac{\pi}{2}]&\text{(axial)}\end{array}\right.\,.$$

$${\mathcal{O}}_L={\mathcal{P}}_L\cdot {\mathcal{R}}_L\text{with}{\mathcal{P}}_L=\prod _{j=1}^{L/4}p({\lambda }_j),{\mathcal{R}}_L=\left|\frac{{\det }_{L/4}(G^{+})}{{\det }_{L/4}(G^{-})}\right|,$$

$$p(\lambda)\;=\;\left\{\begin{array}[]{ll}\displaystyle\frac{\tanh(\lambda+\frac{i\gamma}{2})\tanh(\lambda-\frac{i\gamma}{2})}{4\sinh^{2}(2\lambda)}&\text{(planar)}\\[12.91663pt] \displaystyle\frac{\lambda^{2}+\frac{1}{4}}{16\lambda^{2}}&\text{(isotropic)}\\[12.91663pt] \displaystyle\frac{\tan(\lambda+\frac{i\eta}{2})\tan(\lambda-\frac{i\eta}{2})}{4\sin^{2}(2\lambda)}&\text{(axial)}\end{array}\right.\,.$$

$$G_{jk}^{\pm }$$

$${\tilde{\rho }}_L(\lambda )$$

$$K_{\beta}(\lambda)\;=\;\left\{\begin{array}[]{ll}\displaystyle\frac{1}{2\pi}\,\frac{\sin(2\beta\gamma)}{\sinh^{2}(\lambda)+\sin^{2}(\beta\gamma)}&\text{(planar)}\\[12.91663pt] \displaystyle\frac{1}{2\pi}\,\frac{2\beta}{\lambda^{2}+\beta^{2}}&\text{(isotropic)}\\[12.91663pt] \displaystyle\frac{1}{2\pi}\,\frac{\sinh(2\beta\eta)}{\sin^{2}(\lambda)+\sinh^{2}(\beta\eta)}&\text{(axial)}\end{array}\right.\,,\qquad\beta\;=\;1/2,\,1\,.$$

$$\ln {\mathcal{O}}_L+\alpha L=\ln {\mathcal{P}}_L+\alpha L+\ln {\mathcal{R}}_L=\sum _{j=1}^{L/4}\ln [p({\lambda }_j)]+\alpha L+\ln \left|\frac{{\det }_{L/4}(G^{+})}{{\det }_{L/4}(G^{-})}\right|.$$

$$P(L)=\ln ({\mathcal{P}}_L)+\alpha L=\sum _{j=1}^{L/4}\ln [p({\lambda }_j)]+\alpha L,$$

$$P(L)\simeq P_0+P_1(L)+\dots ,$$

$$\mathfrak{a}(\lambda )={\left[\frac{\sinh (\lambda -\frac{i\gamma }{2})}{\sinh (\lambda +\frac{i\gamma }{2})}\right]}^L\prod _{k=1}^M\frac{\sinh (\lambda -{\lambda }_k+i\gamma )}{\sinh (\lambda -{\lambda }_k-i\gamma )}.$$

$$\sum _{j=1}^{M}f({\lambda }_j)={\oint }_{\mathcal{C}}\frac{\mathrm{d}\omega }{2\pi i}\frac{{\mathfrak{a}}^{\prime }(\omega )}{1+\mathfrak{a}(\omega )}f(\omega ).$$

$$\ln [\mathfrak{a}(\lambda )]=i\gamma L+L\ln \left[\frac{\sinh (\lambda -\frac{i\gamma }{2})}{\sinh (\lambda +\frac{i\gamma }{2})}\right]+i{\oint }_{\mathcal{C}}{K}_1(\lambda -\omega )\ln \left[1+\mathfrak{a}(\omega )\right]\mathrm{d}\omega ,$$

$$\mathfrak{b}(x)={\mathfrak{a}}^{-1}(x-i\gamma /2),\bar{\mathfrak{b}}(x)=\mathfrak{a}(x+i\gamma /2),x\in \mathbb{R}.$$

$$\ln [\mathfrak{b}(x)]$$

Full text

1 Introduction

Recent years have witnessed renewed interest in the study of quantum phase transitions and quantum critical points in low-dimensional systems at zero temperature. Our theoretical understanding of such points relies mainly on two different but complementary approaches. One is field theory, which captures the important low-energy degrees of freedom that emerge at criticality. Another one is a more direct approach, starting from a concrete microscopic model —typically a system with $L$ locally coupled elementary degrees of freedom (e.g. spins or particles)— and studying the thermodynamic limit (TDL) $L\to\infty$.

The case of one dimension (1d) is special as very powerful techniques are available on both sides. Most of the critical points in 1d can be described by conformal field theories (CFT) [1, 2] which are exactly solvable in the sense that it is in principle possible to determine all correlation functions. On the side of the microscopic approach, powerful numerical techniques such as the density matrix renormalization group (DMRG) algorithm [3, 4, 5, 6] or exact diagonalization [7, 8] are available. An important class of microscopic systems are the so-called Bethe ansatz integrable models. For such models an exact determination of all energy eigenstates is possible [9, 10]. Also some correlation functions can be accessed analytically, and the results can be compared with CFT.

Both sets of techniques have different advantages. One obvious limitation of CFT is that it can, by definition, only access universal numbers that describe the corresponding universality class. So, it is insensitive to the exact details of the underlying microscopic model. Integrability methods based on Bethe ansatz, however, allow to access both, universal and non-universal numbers. In that sense, they provide more information. On the other hand, CFT can describe the low-energy physics of non-integrable models. Said differently, the low-energy effective field theory can still be exactly solvable even though the underlying microscopic model is not. Connecting the two sets of methods is of paramount importance and has been a subject of tremendous effort over the last few years [11, 12, 13, 14, 15, 16, 17].

One key model where such a connection can be successfully performed is the anisotropic spin-1/2 Heisenberg chain (XXZ model) whose Hamiltonian reads

[TABLE]

The operators $\sigma_{j}^{\alpha}$ act only non-trivially on lattice site $j$ as Pauli matrices $\sigma^{\alpha}$, $\alpha=x,y,z$, and satisfy periodic boundary conditions, $\sigma_{L+1}^{\alpha}=\sigma_{1}^{\alpha}$. For simplicity, we assume throughout the paper that the system size $L$ is a multiple of four. The parameter $\Delta$ is called the anisotropy of the model. Here, it is restricted to the interval $-1<\Delta\leq 1$, where the XXZ model becomes gapless in the TDL. Later, in Sec. 3, we will also consider the gapped phase $\Delta>1$, where we parameterize the anisotropy by $\Delta=\cosh\eta$, $\eta>0$.

The low-energy behavior of the gapless model (1) can be described by a free compact boson CFT, also known as Tomonaga-Luttinger liquid [18, 19, 20], or Gaussian free field in the mathematical literature [21]. Much effort has also been devoted to the calculation of correlation functions and their large-distance asymptotics, which is a challenging problem [11, 12, 13, 14, 15]. This has allowed to confirm many predictions of the Tomonaga-Luttinger liquid theory in this system. In particular, one can extract the Luttinger parameter $K$ from the low-energy spectrum or, similarly, from the exponents governing the algebraic decay of bulk correlation functions. It depends on the anisotropy $\Delta=\cos\gamma$ as follows,

[TABLE]

The knowledge of $K$ uniquely determines the underlying CFT.

In this paper, we demonstrate an alternative way to access the Luttinger parameter, namely, by making a connection to boundary CFT, a subject which has been largely developed by Cardy [22, 23, 24]. We consider the overlap $\mathcal{O}_{L}$ of the ground state $\ket{\Psi_{L}}$ of the XXZ model (1) with the Néel state $\ket{N}=\ket{\uparrow\downarrow}^{\otimes L/2}$,

[TABLE]

As was shown in Refs. [25] and [26], its large system size scaling can be analyzed by means of boundary CFT. The overlap may be interpreted as a ratio of partition functions of 2d classical systems (see Sec. 2) with boundary conditions that can be argued to renormalize to conformal invariant boundary conditions[23]. The order one contribution to the overlap is determined by the “universal ground state degeneracy”, or “the boundary $g$-factor” [27]. It can be viewed as an analogue of the central charge in boundary critical phenomena. In particular, it decreases under RG flow.

Crucially, the overlap with the Néel state can also be evaluated in closed form by means of certain determinant formulas[28, 29, 30, 31, 32] whose derivation is based on algebraic Bethe ansatz. This works, in principle, for all eigenstates of the XXZ model, not only for its ground state. Overlaps of XXZ eigenstates with the Néel state were subject of intense study, as they proved instrumental in the understanding of thermalization after quantum quenches [30, 33, 34, 35]. We refer the interested reader to the review article [36]. In the following, however, we are interested in the overlap with the ground state only and, in particular, in its large system size scaling. For this analysis, we will use the determinant formula that is analytically most tractable, and which was derived in Ref. [31].

To be more specific, our main result is the following formula, valid in the ‘planar regime’ $-1<\Delta<1$ and obtained by a careful asymptotic analysis of the exact finite-size determinant formula (see Sec. 3),

[TABLE]

Here, $K$ is the Luttinger parameter of Eq. (2). The exponential factor $e^{-\alpha L}$ is called the ‘exponential decay’ of the overlap with ‘decay rate’ $\alpha$. Furthermore, we shall refer to the prefactor (here $\sqrt{K}$) as the ‘order one term’, since it determines the $\mathcal{O}(1)$ contribution of $\ln(\mathcal{O}_{L})$, as well as to $A$ and $\delta$ as the ‘amplitude’ and the ‘exponent’ of the leading finite-size correction. The dots in Eq. (4) represent higher subleading terms. All coefficients are non-trivial functions of $\Delta$. The order one term, which is in exact agreement with the prediction of boundary CFT (see Eq. (14) in Sec. 2), as well as the exponent $\delta=\min\{1,4K-2\}$ are universal whereas the decay rate $\alpha$, which has the dimension of an inverse length, and the amplitude $A$, which has the dimension length to the $\delta$, are not. We nevertheless derive a fully explicit integral formula for $\alpha$ and compute the amplitude $A$ (for several values of $\Delta$) to high numerical accuracy (see Sec. 3.3.2). In this sense, we are able to extract even the leading finite-size correction $A/L^{\delta}$ to the overlap, which is not accessible by CFT per se (even though the exponent $\delta$ may be obtained by a perturbed CFT treatment, see Refs. [37, 38, 39] for such kind of studies).

We further provide results for the ‘isotropic point’ $\Delta=1$ as well as for the ‘axial regime’ $\Delta>1$ (see Sec. 3). The isotropic point can be accessed by taking the limit $\gamma\to 0$ (i.e. $K\to 1/2$ $\Rightarrow$ “$\delta\to 0$”), in which Eq. (4) turns into

[TABLE]

The main difference to the planar regime is that the leading finite-size correction is no longer algebraic ($\sim 1/L^{\delta}$) but logarithmic ($\sim 1/\ln L$).

The result for the overlap in the axial regime reads

[TABLE]

This formula looks a bit simpler than the result (4) of the planar regime, as expected for a gapped phase. Note that the finite-size corrections in the axial regime are exponentially small, and that the order one term is independent of the anisotropy $\Delta>1$. The value of the latter can be intuitively understood by considering the limit $\Delta\to\infty$, where the ground state is given by $\ket{\Psi}_{L}=(\ket{\uparrow\downarrow}^{\otimes L/2}+\ket{\downarrow\uparrow}^{\otimes L/2})/\sqrt{2}$. Hence, the overlap is $\mathcal{O}_{L}=1/2$ (implying that $\alpha=0$ and $c=\infty$). For $\Delta>1$ but finite, $\ket{\Psi}_{L}$ is of course more complicated, but the order one term still reflects the symmetry broken nature of the ground state.

2 Overlap in boundary conformal field theory

We focus in this section on the XXZ model in its gapless phase $-1<\Delta\leq 1$, described by the Hamiltonian (1). The ground state of the XXZ model can be seen as the result of an infinite imaginary-time evolution starting from an arbitrary state $\ket{a}$, provided it has nonzero overlap with the ground state,

[TABLE]

Choosing the Néel state $\ket{a}=\ket{N}$, for convenience, and using the previous equation, the overlap can be recast as

[TABLE]

This can be interpreted as a ratio of partition functions of two-dimensional classical systems,

[TABLE]

Here, $Z_{\rm cyl}(L,\uptau)$ denotes the partition function of a cylinder of circumference $L$ and height $\uptau$ (see Fig. 1). The boundary conditions at the top and at the bottom are set by the Néel state. Now, we make the crucial assumption that the Néel boundary condition renormalizes to a conformal invariant boundary condition [22, 24]. This argument is very general and only requires to determine the appropriate conformal boundary condition in the continuum limit, if it exists. It was already exploited in Ref. [26], in order to compute the ratio (9) by means of boundary CFT. It can be applied to other models corresponding to different universality classes as well [25, 40, 41]. Generalizations to chains with open boundary conditions and with universal power-law corrections to the partition function [42] are also possible [26, 43, 44, 45].

The conformal boundary condition can be determined as follows. First, we map spin configurations onto discrete height configurations $h_{j}$, where the height picks up an additive factor, say, $+1$ when encountering an up spin and $-1$ when encountering a down spin, $h_{j}=\sum_{i=1}^{j}\sigma_{i}^{z}$. Thus, for the Néel state, we obtain $h_{j}=(1+(-1)^{j+1})/2$. In the continuum limit, this discrete ‘height field’ $h_{j}$ is expected to renormalize to a free compact bosonic field, through a procedure known as bosonization [20]. Since the discrete height configuration of the Néel state is flat on average, it is natural to expect that the field be constant at the boundary that is set by the Néel state. Such a conformal invariant boundary condition is called Dirichlet boundary condition.

Let us now come back to the computation of the overlap (9). We are left with partition functions of a ‘free compact boson CFT’ on a cylinder with Dirichlet boundary condition, which are well-known [2]. We find it instructive to reproduce the derivation of Ref. [46], as it allows to make contact with other works which are similar in spirit to ours. Instead of looking at a periodic system which is let evolve in imaginary time $\uptau$, it is equivalent to interpret the CFT partition function as a finite temperature average of a system of width $\uptau$ with free spins at the boundary and a periodic time direction, where the time is now represented by the system size $L$. This is illustrated in Fig. 1.

The low-energy physics of an open XXZ chain of length $\uptau$ is well-described by the bosonized form of the XXZ Hamiltonian [46],

[TABLE]

The coefficients $E$ and $v$ can be interpreted as the extensive part of the ground state energy and as a velocity, respectively. The order one contribution $e$ originates from the boundaries solely. The total magnetization $S^{z}$ of the chain is an integer, and the $a_{k},a_{k}^{\dagger}$ are bosonic operators that satisfy the Heisenberg commutation relations $[a_{k},a_{l}^{\dagger}]=\delta_{kl}$ and $[a_{k}^{\dagger},a_{l}^{\dagger}]=[a_{k},a_{l}]=0$. It is also understood that $\sum_{k=1}^{\infty}k=\zeta(-1)=-1/12$ (Zeta regularization). The partition function we are looking for can be written as

[TABLE]

where $L$ now plays the role of an inverse temperature. The trace has to be performed over all bosonic modes and all magnetization sectors. Inserting the Hamiltonian (10) into this definition yields

[TABLE]

This result can also be expressed in terms of Dedekind eta and Jacobi theta functions. Several works have looked at the large-$L$ scaling of this partition function and have compared with Bethe ansatz results for open chains at low temperatures, see e.g. Refs. [47] and [29]. However, we are interested in performing a different limit, namely $\uptau\to\infty$, which is not apparent from the above formula. To access this limit, it is convenient to express the partition function (12) in a different form. Using the Poisson re-summation formula (or, equivalently, the known modular properties of eta and theta functions), we eventually obtain

[TABLE]

Note the appearance of a prefactor $\sqrt{K}$. In the CFT literature, such factors are sometimes called $g$-factors, whose importance has been pointed out in Ref. [27]. In general, they depend on the CFT as well as on the type of conformal boundary condition, which is Dirichlet here. Using Eq. (13), formula (9) for the overlap simplifies to

[TABLE]

Of course, this derivation takes only the contributions from low-energy degrees of freedom into account. In particular, the Hamiltonian (10) is not exact any more at high energies. Treating high-energy effects would require to add multi-boson interactions and to put a cut-off on the large but finite number of bosons, which would make exact computations much harder. It is expected, but by no means obvious, that the prefactor $\sqrt{K}$ will remain unaffected by their presence, whereas the coefficient $e$ of the exponential decay will change. Lattice effects will also be responsible for additional but subleading corrections to Eq. (14).

These observations make it desirable to explicitly check the validity of the boundary CFT formula (14) for the case of the spin-1/2 XXZ chain. In Ref. [26] this was done exactly at the free fermion point $\Delta=0$, i.e. $K=1$, where $\mathcal{O}_{L}=2^{-L/2}$ [48], as well as numerically for other values of $\Delta$ by means of exact diagonalization techniques for small systems up to $L=28$. Agreement within a few percent was found. However, since there is now an exact Bethe ansatz-based formula for the overlap available [31], which is valid for any system size $L$, an analytical treatment of the large system size scaling of $\mathcal{O}_{L}$ is possible. In the next section, we extract the leading terms of its asymptotic expansion. This enables us to analytically confirm the presence of the universal $\sqrt{K}$ term (analytically in a certain sense, as discussed below). We also explicitly determine the coefficient of the exponential decay, which does not coincide, as expected, with the order one contribution $e$ to the ground state energy of an open XXZ chain. Furthermore, we compute the leading finite-size correction, which shows an universal algebraic decay.

3 Asymptotics from the exact finite-size determinant formula

The starting point of our asymptotic analysis is a recently derived Bethe ansatz-based determinant formula for the overlap of the Néel state with an XXZ eigenstate [31], valid for any system size $L$. Before presenting the details of this analysis, we first outline in Sec. 3.1 some important properties of the XXZ model and, in particular, of its ground state. In Sec. 3.2 we provide the explicit form of the determinant formula in terms of so-called Bethe roots. Subsequently, in Secs. 3.3 and 3.4, we analyze the large system size scaling of the overlap. We calculate the thermodynamic limit (TDL) $L\to\infty$ of properly defined quantities, e.g. the decay rate, the order one term, and (the exponent of) the leading finite-size correction, which altogether characterize the asymptotic behavior of the overlap as predicted by Eqs. (4)–(6).

3.1 Properties of the XXZ model and its ground state

The low-energy spectrum of the XXZ model (1) with anisotropy $-1<\Delta\leq 1$ becomes gapless in the TDL, whereas it remains gapped for $\Delta>1$. Eigenstates can be explicitly determined within the framework of coordinate Bethe ansatz [9, 10]. They depend on a set of parameters (Bethe roots) which satisfy, due to periodic boundary conditions, a set of so-called Bethe equations. Depending on the value of $\Delta$, this set of equations looks slightly different. We distinguish three cases:

• •

the ‘planar regime’ $-1<\Delta<1$ with parameterization $\Delta=\cos\gamma$, $0<\gamma<\pi$,

• •

the ‘isotropic’ (also ‘phase-separating’) point $\Delta=1$, and

• •

the ‘axial regime’ $\Delta>1$ with parameterization $\Delta=\cosh\eta$, $\eta>0$.

The corresponding Bethe equations read as follows,

[TABLE]

Although the Bethe equations of these three cases look similar, the structures of their solutions show perceptible differences. In the axial regime, for instance, real parts of Bethe roots are bounded in the TDL, whereas they are unbounded in the two other cases (see e.g. Eq. (16) below). In the planar regime, the dependence of a Bethe state on Bethe roots shows a periodicity of $\pi$ in imaginary direction. At the isotropic point, i.e. $\gamma\to 0$ in $\Delta=\cos\gamma$, this periodicity gets lifted after rescaling all Bethe roots by $\gamma$.

The XXZ ground state is determined by $L/2$ many different Bethe roots $\{\lambda_{j}\}_{j=1}^{L/2}$ (no magnetic field) that are symmetrically distributed on the real line with distances between them as small as possible. We therefore only need to know the positive roots $\{\lambda_{j}\}_{j=1}^{L/4}$ where, by definition, $0<\lambda_{1}<\ldots<\lambda_{L/4}$ (remember that $L$ is assumed to be a multiple of four). For later purposes, we show in B that the smallest and largest positive ground state Bethe roots scale in the large system size limit as

[TABLE]

In the expressions of the axial regime ($\Delta>1$), $K(k)$ is the elliptic integral of the first kind with modulus $k$, which is uniquely determined by the anisotropy $\Delta=\cosh(\eta)$ via $\eta=\pi K(k)/K(\sqrt{1-k^{2}})$. Here, $K(k)$ should not be confused with the Luttinger parameter $K$ of Eq. (2), defined for $-1<\Delta\leq 1$ only. We explicitly see from this result that the positive ground state Bethe roots in the planar regime and at isotropic point are unbounded in the TDL, $\lambda_{L/4}\sim\ln(2L/\pi)$, while they are bounded by $\pi/2$ in the axial regime, approaching this bound in the TDL.

The discrete distribution of ground state Bethe roots (‘ground state distribution’),

[TABLE]

is a symmetric function, which is normalized to $1/2$ by definition. It becomes smooth in the TDL,

[TABLE]

The explicit form of $\rho$ has been known for a long time (as solution of a linear integral equation (LIE) [49, 50, 51]), though a rigorous derivation of the LIE (at half filling) was lacking. Recently, in Ref. [52], a proof of this LIE was found.

In Fig. 2 we show some TDL ground state distributions for several values of the anisotropy parameter $\Delta$. The symmetry $\rho_{L}(-\lambda)=\rho_{L}(\lambda)$ is a necessary condition for the overlap with the Néel state not to vanish. If an XXZ Bethe state does not fulfil the “condition of parity invariance” [31], i.e. $\{\lambda_{j}\}_{j=1}^{L/2}=\{-\lambda_{j}\}_{j=1}^{L/2}$, then the exact finite-size determinant formula yields exactly zero, as was shown in Ref. [32].

Despite all these differences (see e.g. the ground state distributions (18) of the planar, isotropic, and axial cases), many Bethe ansatz-based formulas look similar for the three cases (see e.g. Eqs. (15) and Ref. [13]). They can be transformed into each other by replacing all hyperbolic sine functions (planar regime) by sine functions (axial regime) or by its argument (isotropic point), and so forth. This also holds for the exact finite-size overlap formula whose explicit expression is presented in the next subsection.

3.2 Finite-size determinant formula

We exploit the determinant formula of Ref. [31] for the overlap of an XXZ Bethe state with the Néel state by inserting the positive ground state Bethe roots $\{\lambda_{j}\}_{j=1}^{L/4}$. In all three cases (planar, isotropic, axial), the exact overlap formula has the following structure (‘prefactor’ $\mathcal{P}_{L}$ times ‘ratio of determinants’ $\mathcal{R}_{L}$),

[TABLE]

It can be shown [51, 52] that $\lim\nolimits_{L\to\infty}\tilde{\rho}_{L}(\lambda)=\rho(\lambda)$ with the TDL ground state distribution $\rho$ of Eq. (18).

We are interested in the large system size scaling of this overlap $\mathcal{O}_{L}$, which we expect to be, in the three cases ‘planar’, ‘isotropic’, and ‘axial’, of the form (4), (5), and (6), respectively. To show this, we first of all observe that the prefactor $\mathcal{P}_{L}$ alone accounts for the exponential decay $e^{-\alpha L}$ of $\mathcal{O}_{L}$. This will become manifest in our analysis in Sec. 3.4, where we show that $\mathcal{R}_{L}$ converges to a constant and, hence, contributes at most to the order one term. Thus, the decay rate can be expressed as $\alpha=-\lim_{L\to\infty}\ln(\mathcal{P}_{L})/L$. In A we explicitly compute it for all $\Delta\geq-1$.

For later convenience, we consider the logarithm of the overlap $\mathcal{O}_{L}$ plus $L$ times the decay rate $\alpha$,

[TABLE]

In the following, we calculate the order one contributions as well as the leading finite-size corrections of $\ln(\mathcal{P}_{L})+\alpha L$ and $\ln(\mathcal{R}_{L})$ explicitly, proving that Eqs. (4)–(6) are indeed the correct asymptotic expansions of $\mathcal{O}_{L}$ for large system size $L$. In Sec. 3.3 we study the different contributions of the ‘prefactor term’ $\ln(\mathcal{P}_{L})+\alpha L$ while Sec. 3.4 is devoted to the analysis of the $\ln(\mathcal{R}_{L})$-term. In Sec. 3.5 we summarize the results by comparing the contributions of the prefactor term with those of the ratio of determinants. We further identify the source of the leading finite-size correction of the overlap in the three different cases ‘planar’, ‘isotropic’, and ‘axial’.

3.3 Asymptotic analysis of the prefactor

Let us define the ‘prefactor term’ $P(L)$ by

[TABLE]

where $p$ is the prefactor function of Eq. (19b), and $\alpha$ is the decay rate as computed in A. We are interested in an asymptotic large-$L$ expansion of $P(L)$,

[TABLE]

where we refer to the zeroth order term $P_{0}$ as the order one contribution of $P(L)$ and to the first order term $P_{1}$ as its leading finite-size correction. The dots denote higher orders, i.e. subleading finite-size corrections.

In Sec. 3.3.1 we derive an explicit expression for $P_{0}$ in all three cases (planar, isotropic, axial), which is the most difficult part and which will be the main result of our asymptotic analysis. Furthermore, in Sec. 3.3.2 we extract the leading finite-size correction $P_{1}(L)$, which turns out to be qualitatively different in the three cases ‘planar’, ‘isotropic’, and ‘axial’. In the planar regime, for instance, we will see that the prefactor term scales as $P(L)\simeq\ln(K)+A/L^{\delta}$ with an algebraically decaying finite-size correction $P_{1}(L)=A/L^{\delta}$ whereas, at the isotropic point ($K=1/2$), it scales as $P(L)\simeq-\ln(2)+\mathcal{O}(\ln^{-2}(L))$. In the axial regime, on the contrary, the correction to the zeroth order term $P_{0}=-2\ln(2)$ is exponentially small, $P(L)\simeq-2\ln(2)+\mathcal{O}(e^{-cL})$.

3.3.1 Order one term.

One idea to compute the order one contribution $P_{0}$ of the prefactor term (22) could be to use the Euler-Maclaurin formula. By means of this formula, a sum over ground state Bethe roots as in Eq. (21) can be approximated for large system size $L$ by an integral (times $L$) plus order one and subleading corrections (see e.g. Ref. [53] and Eq. (62) in B). In our case, the function under the integral, $\ln[p(\lambda)]$, shows logarithmic singularities at the boundaries of the integration interval. In B we demonstrate that an application of the Euler-Maclaurin formula to such a function encounters technical difficulties. Nevertheless, by making use of the expansion of the log-gamma function [54, 55], an analytic treatment of the singularity at the lower boundary is possible in all three cases (planar, isotropic, axial). In contrast, the singularity at the upper boundary can be treated only for the axial case. For the planar and isotropic cases, unfortunately, the singularity of $\ln[p(\lambda)]$ at infinity leads to insurmountable difficulties. This is why this approach works well in the axial regime, while it does not work in the planar regime and at the isotropic point.

Another approach to compute sums over Bethe roots is based on the well-established technique of non-linear integral equations (NLIE). This technique was intended for and successfully applied to the computation of finite-temperature thermodynamics [56, 57]. It was shown that it is also applicable to finite systems (of arbitrary length) at zero temperature [58]. More specifically, the NLIE can be formulated for any excited (low-lying) state of the finite system. “Zero temperature” then just means to consider the ground state. This approach works well in all three cases (planar, isotropic, axial). The only difference is that the driving terms and kernels of the NLIE look slightly different [59]. The structure of the NLIE, however, as well as the logic of the derivation of the order one term $P_{0}$ are exactly the same. Note that our analysis of the NLIE does not only provide an exact expression for the order one contribution $P_{0}$, but also enables us to determine the scaling of the leading finite-size correction $P_{1}(L)$ (see Sec. 3.3.2).

In the following, we focus on the second approach in the regime $0\leq\Delta<1$. The other part of the planar regime, $-1<\Delta<0$, as well as the isotropic point $\Delta=1$ and the axial regime $\Delta>1$ can be treated similarly. The NLIE approach is based on the following formulas [60, 58]. We consider the ‘auxiliary function’

[TABLE]

Obviously, due to the Bethe equations (15a), $\mathfrak{a}(\lambda_{j})=-1$. Furthermore, we consider a contour $\mathcal{C}$ in the complex plane that encircles the Bethe roots $\{\lambda_{j}\}_{j=1}^{M}$ in a counter-clockwise manner and that does not encircle hole solutions of the Bethe equations (if present, e.g. for excited states of the XXZ chain). However, since we are interested in the ground state of the XXZ model, characterized by roots only, we do not need to be concerned too much about hole solutions. For a function $f$ that is analytic inside the contour $\mathcal{C}$ the sum of $f(\lambda_{j})$ over all Bethe roots $\lambda_{j}$ can be then expressed by the following contour integral,

[TABLE]

The auxiliary function $\mathfrak{a}$ satisfies the following non-linear integral equation (NLIE) [58, 60],

[TABLE]

where $K_{1}$ is the kernel function of Eq. (19e). In principle, the full $L$-dependence of a sum over Bethe roots is completely determined by the last two equations. In practice, however, it is a difficult task (if not impossible) to solve the NLIE analytically, in order to get an explicit expression for $\mathfrak{a}(\lambda)$, even in the limit of large system size $L\to\infty$.

In the case of the ground state, which consists of $L/2$ symmetrically distributed real Bethe roots, it is convenient to choose the contour $\mathcal{C}$ such that the entire real axis is enclosed. This way, $\mathcal{C}$ becomes independent of the detailed positions of the ground state Bethe roots and, in particular, independent of $L$. The $L$-dependence of the sum in Eq. (24) is fully captured by the ‘driving term’ of the NLIE (25). Still, an analytical solution of the NLIE is out of reach. Due to the complex contour $\mathcal{C}$, even a numerical evaluation is problematic. This second obstacle can be cleared by expanding the contour $\mathcal{C}$ to an infinitly long rectangle with edges parallel to the real axis (at levels $\pm i\gamma/2$) and by defining two new auxiliary functions,

[TABLE]

Using this specific contour and plugging these definitions into Eq. (25), a set of two coupled non-linear integral equations (NLIEs) for the functions $\mathfrak{b}$ and $\bar{\mathfrak{b}}$ can be derived [58]. Since the contributions from the vertical parts of the rectangular contour can be neglected, it reads

[TABLE]

Now, we would like to rewrite the right-hand side of Eq. (24) in terms of the auxiliary functions $\mathfrak{b}$ and $\bar{\mathfrak{b}}$. First of all, we note that the driving term $d$ of these new NLIEs is related to the TDL ground state distribution $\rho$ via $id^{\prime}(x+i\gamma/2)=2\pi\rho(x)$. Provided that $f$ is analytic in the whole strip $\{\lambda\in\mathbb{C}\;|\,-\frac{\gamma}{2}<{\rm Im}(\lambda)<\frac{\gamma}{2}\,\}$ the sum over ground state Bethe roots can be explicitly written as

[TABLE]

Here, $\delta$ is the delta distribution defined by $(\delta\ast g)(x)=g(x)$. The derivation of Eq. (28) is straightforward. First, we deform the contour $\mathcal{C}$ in Eq. (24) to the rectangle described above Eqs. (26) and subsequently insert the definitions (26) themselves. Since on the lower contour (lower edge of the rectangle) the inverse of $\mathfrak{a}$ is used in the definition of $\mathfrak{b}$, also the derivative of $\ln(\mathfrak{b})$ occurs. In order to express everything by derivatives of $\ln(1+\mathfrak{b})$ and $\ln(1+\bar{\mathfrak{b}})$ solely, we plug in the derivative of the NLIE (27a) with respect to $x$. After an expedient shift of the integration contour in one of the four convolution terms, we finally arrive at Eq. (28). In the very last step we explicitly used the analyticity condition of $f$ mentioned above Eq. (28).

We further note that ${\rm Re}(\mathfrak{b})={\rm Re}(\bar{\mathfrak{b}})$ and ${\rm Im}(\mathfrak{b})=-{\rm Im}(\bar{\mathfrak{b}})$, which reflects the symmetry of the ground state Bethe roots and basically means that the set (27) of two NLIEs can be reduced to one NLIE for one unknown function $\mathfrak{b}$. However, solving the whole set (27) is computationally not much more costly, especially if we take into account that the real parts of $\mathfrak{b}$ and $\bar{\mathfrak{b}}$ are symmetric functions and that the imaginary parts are antisymmetric. If we require that $f(x\pm i\gamma/2)$ as a function of real argument $x$ has the same symmetry as the auxiliary functions $\mathfrak{b}$ and $\bar{\mathfrak{b}}$, i.e. ${\rm Re}\,f(x\pm i\gamma/2)={\rm Re}\,f(-x\pm i\gamma/2)$ and ${\rm Im}\,f(x\pm i\gamma/2)=-{\rm Im}\,f(-x\pm i\gamma/2)$, the sum over positive ground state Bethe roots can be simplified to

[TABLE]

Since we sum over positive ground state Bethe roots only, we integrate from zero to infinity in the first term, and an additional factor of $1/2$ appears in the second term.

Regarding an evaluation of Eqs. (27) and (28) or (29), there are three big advantages as compared to Eqs. (24) and (25). Firstly, both functions $\mathfrak{b}$ and $\bar{\mathfrak{b}}$ have to be evaluated only for purely real arguments. Secondly, they are very small close to the origin ($\sim x^{L}$), as can be most easily seen from their definitions, i.e. Eqs. (26) together with Eq. (23). Thirdly, the integrations are convolutions along the real axis. To evaluate them we can apply efficient numerical techniques such as the Fast Fourier Transform algorithm. We solve the NLIEs (27) for a fixed system size $L$ by iteration, which converges after a few steps (a couple of hundreds or less) and provides the functions $\mathfrak{b}$ and $\bar{\mathfrak{b}}$ to high numerical accuracy. Since $L$ enters only through the prefactor of the driving term $d$, we are able to compute solutions for system sizes up to $L=10^{10}$ or even larger. Subsequently, after an integration by parts, we plug the solution $\mathfrak{b}$ into Eq. (29) and compute the integral. In this sense, we can calculate sums over ground state Bethe roots exactly, i.e. to high numerical accuracy (up to $10^{-14}$ if we want), for any system size $L$ and, in particular, for very large $L$. This enables us, for instance, to calculate the order one term $P_{0}$ of the prefactor term $P(L)$ as well as its leading finite size correction $P_{1}(L)$, both to high numerical accuracy (for details see the following discussion as well as Sec. 3.3.2 and, in particular, Fig. 3 therein).

We are interested in the special case

[TABLE]

We observe that the function $f$ satisfies the symmetry requirements (see statement right before Eq. (29)), but it does not satisfy the analyticity condition (see right before Eq. (28)). The function $f^{\prime}$ has a pole on the real axis at zero (with residue $-2$), and therefore $f$ is not analytical at zero and shows a branch cut similar to that of the function $\ln(x)$. This leads to an additional explicit term which can be easily computed by considering the contour integration in Eq. (24) directly. We integrate the right-hand side by parts and expand the contour $\mathcal{C}$ to the rectangle that eventually leads to the $\mathfrak{b}$-$\bar{\mathfrak{b}}$-formulation (27) of the NLIE. This order of steps is most convenient. Though, the other way around leads, of course, to the same result. Then, from this pole at zero, we encounter an explicit term which has to be subtracted from the right-hand side of Eq. (29). Therefore, Eq. (29) will be modified by the following extra term,

[TABLE]

where $D_{\epsilon}(0)$ is the disc around zero with small radius $\epsilon$, and where we used that $\mathfrak{a}(0)=1$ which is due to the symmetry of ground state Bethe roots, $\{\lambda_{j}\}_{j=1}^{L/2}=\{-\lambda_{j}\}_{j=1}^{L/2}$. Taking this term into account and using that the first integral in Eq. (29) is the decay rate $\alpha$ (see A), we eventually obtain

[TABLE]

This formula gives us access to the prefactor term $P(L)$ for any system size $L$, either numerically to very high accuracy (as discussed before) or, in a certain sense, analytically in the TDL (as discussed in the following).

The NLIEs (27) can be analyzed in the TDL, similarly to the low-temperature analysis of Ref. [61]. One observation from our numerical investigation of the NLIEs (27) is that the functions $\ln^{\prime}(1+\mathfrak{b})$ and $\ln^{\prime}(1+\bar{\mathfrak{b}})$ change significantly only in a small region around $\pm\mathcal{L}$ with $\mathcal{L}=\frac{\gamma}{\pi}\ln(2L)$. Furthermore, they vanish for $x\to 0$ (rapidly) and for $x\to\infty$ (slowly), i.e. $\ln^{\prime}(1+\mathfrak{b})$ and $\ln^{\prime}(1+\bar{\mathfrak{b}})$ show ‘bumps’ around $\pm\mathcal{L}$. Thus, it is convenient to define ‘scaling functions’ by

[TABLE]

Using these definitions in the NLIEs (27) and neglecting contributions from convolution integrals which connect bumps of $\ln^{\prime}(1+\mathfrak{b})$ or $\ln^{\prime}(1+\bar{\mathfrak{b}})$ around $\mathcal{L}$ with those around $-\mathcal{L}$, we obtain the following NLIEs for the scaling functions,

[TABLE]

where the integration kernels are the same as in Eq. (27d). Note that the driving terms of these NLIEs are independent of $L$. We interpret the scaling functions (i.e. solutions of these ‘TDL-NLIEs’) as the TDL of the auxiliary functions $\mathfrak{b}$ and $\bar{\mathfrak{b}}$. Taking into account that the function $\ln(p)$ scales for large arguments as

[TABLE]

we derive from Eq. (32) the following expression,

[TABLE]

Here, we used $\int_{-\infty}^{\infty}[\delta(x)-F(x)]{\rm d}x=K$ with the Luttinger parameter $K=\pi/[2(\pi-\gamma)]$ and, for the second term in Eq. (36), that $\lim_{x\to\pm\infty}\ln[1+\mathfrak{b}(x)]=\ln(2)$. To derive the third term on the right-hand side of Eq. (36), we used an integration by parts, the definition (33), and that $\lim_{x\to\pm\infty}\ln^{\prime}[\sinh^{2}(2x)]=4$. The result (36), with $\mathfrak{b}_{+}$ as the (unique) solution of Eqs. (34), is exact. Unfortunately, we were unable to simplify it further. However, by solving the TDL-NLIEs (34) numerically to high accuracy (essentially machine precision), we could show that the integral in Eq. (36) is given by

[TABLE]

This eventually leads to the final result of this subsection,

[TABLE]

The limit $\gamma\to 0$ (isotropic point, $K=1/2$) can be easily taken and yields $P_{0}=-\ln(2)$.

In the axial regime, the analysis of the prefactor term $P(L)$ is similar but less involved. This is due to the fact that Bethe roots are bounded and integrations (of the convolutions and of all other integrals) run from $-\pi/2$ to $\pi/2$. Hence, the integral in Eq. (32) is subleading, rather than of order one as in the planar regime, and does not contribute to $P_{0}$. Instead, there is a second additional explicit term, which is similar to that of Eq. (31) and stems from the pole of $\ln^{\prime}(p)$ at $\pi/2$ (also with residue $-2$),

[TABLE]

Here, we used that in the axial regime the hyperbolic sine functions in the definition (23) of the auxiliary function $\mathfrak{a}$ are replaced by usual sine functions. Hence, $\mathfrak{a}(\pi/2)=1$, again due to the symmetry of ground state Bethe roots $\{\lambda_{j}\}_{j=1}^{L/2}=\{-\lambda_{j}\}_{j=1}^{L/2}$. Therefore, we obtain for the order one term $P_{0}$ of the prefactor term $P(L)$ in the axial regime the following explicit, rigorously derived expression,

[TABLE]

which is independent of the anisotropy parameter $\Delta>1$.

3.3.2 Leading finite-size corrections.

The analysis of the leading finite-size correction $P_{1}(L)$ of the prefactor term $P(L)$ proceeds as follows. There are a priori two main contributions to $P_{1}(L)$, which we discuss separately. The first comes from the terms that connect bumps of the auxiliary functions around $\mathcal{L}$ with bumps around $-\mathcal{L}$, which were neglected in the derivation of the TDL-NLIEs (34). Shifting the integration variables of the convolution integrals by $\pm\mathcal{L}$ and inserting the definition (33) of the scaling functions into the NLIEs (27), one observes that the leading correction to the TDL-NLIEs (34) stems from convolution terms with $\ln(1+\mathfrak{b}_{-})$, $\ln(1+\bar{\mathfrak{b}}_{-})$ (not from those with $\ln(1+\mathfrak{b}_{+})$, $\ln(1+\bar{\mathfrak{b}}_{+})$), see e.g. Eqs. (6)–(8) of Ref. [61]. It is therefore of the form

[TABLE]

where $\beta$ is the exponent of the leading exponential decay of the kernel function $F$. The proportionality factor is a function of $x$ and, to leading order, independent of $\mathcal{L}$. So, it is unimportant to us. In summary, the leading corrections to the scaled auxiliary functions are of order $\exp\{-2\mathcal{L}\beta\}$. The second contribution is an integral proportional to $F(\mathcal{L})\simeq\exp\{-\mathcal{L}\beta\}$, which was neglected in going from Eq. (32) to its scaled version (36). Hence, the leading finite-size corrections of the scaled auxiliary functions are subleading with respect to this integral and, therefore, $P_{1}(L)$ is proportional to $\exp\{-\mathcal{L}\beta\}$ with $\mathcal{L}=\gamma\ln(2L)/\pi$. The exponent $\beta$ can be easily obtained by investigating the poles of the Fourier transform of $F$ in the complex plane, see Eq. (19e). The result is simply $\beta=\text{min}\{\frac{\pi}{\gamma},\frac{2\pi}{\pi-\gamma}\}$.

Thus, by inspecting the decay of the kernel function $F$ of the NLIEs (27) for large arguments, we have shown that the leading finite-size correction $P_{1}(L)$ of the prefactor term $P(L)$ decays algebraically as $A/L^{\delta}$ with amplitude $A=A(\Delta)$ and exponent111Note that a priori this amplitude and this exponent can be different from the amplitude and the exponent of the leading finite-size correction of $\ln(\mathcal{O}_{L})$ (in the planar regime), cf. Eq. (4). As we will discuss later in Sec. 3.5, they are in fact the same. We therefore denote them by the same letters here.

[TABLE]

This formula is one of the main results of this subsection. It is valid in the entire planar regime $-1<\Delta=\cos(\gamma)<1$.

In Fig. 3 we show the exponent $\delta=\delta(\Delta)$ as a function of the anisotropy $\Delta$ in the regime $0<\Delta<1$ (black solid line). It is compared to ‘finite-size exponents’

[TABLE]

which we computed by means of the NLIEs (27) for several values of $L$ up to $L=10^{10}$ (with a properly chosen $\Delta L\ll L$). We further checked for several values of $\Delta$ with $-1<\Delta<0$ (by solving Bethe equations) that the exponents $\delta_{L}$ are very close to one already for $L\gtrsim 100$ (results not shown in the figure). We also computed ‘finite-size amplitudes’ (see right panel of Fig. 3),

[TABLE]

again to very high accuracy by solving the NLIEs (27).

Furthermore, we can take the isotropic limit of the kernel function $F$. For large arguments $x$ it decays algebraically, $F(x)\sim 1/x^{2}$, instead of exponentially as in the planar regime. Following the line of arguments of the planar case, we infer that the leading finite-size correction $P_{1}(L)$ at the isotropic point is logarithmic, $P_{1}(L)\simeq\mathcal{O}(\ln^{-2}(L))$. We verified this by solving the corresponding non-linear integral equations up to $L=10^{15}$. Note that solving the Bethe equations (15b) up to $L\approx 8000$ is by far not sufficient to reach the asymptotic regime (see lower curve of the right panel of Fig. 4 in Sec. 3.4 and, in particular, the inset).

In the axial regime, on the other hand, the leading finite-size correction $P_{1}(L)$ is exponentially small, $P_{1}(L)\sim e^{-cL}$. Thus, it is sufficient to solve the Bethe equations (15c) for system sizes up to $L=320$ and to estimate the value of $c$ by an extrapolation (see e.g. Sec. 3.4 and, in particular, the inset of Fig. 5).

3.4 Asymptotic analysis of the ratio of determinants

This section deals with the asymptotics of the ratio of determinants $\mathcal{R}_{L}$, Eq. (19a), for large system size $L$. The formula for the ratio of determinants reads (see Sec. 3.2)

[TABLE]

where the kernel functions $K_{1}^{\pm}$ and the function $\tilde{\rho}_{L}$ are given by Eqs. (19e) and (19d), respectively. It is valid for all system sizes $L$ (with the only restriction that $L$ has to be a multiple of four). We first demonstrate how $\mathcal{R}_{L}$ converges to a ratio of Fredholm determinants in the TDL, which we evaluate analytically using a theorem of Ref. [62]. We also would like to draw the interested reader’s attention to Ref. [63] for a review on the asymptotic analysis of such types of determinants. Subsequently, the finite-size corrections to this TDL value are studied numerically by solving the Bethe ansatz equations (15) for finite but large system sizes up to $L=8192$.

Let us first focus on the planar and isotropic cases, which are the most complicated. In the TDL, the distribution of ground state Bethe roots becomes a smooth function $\rho$ on the real line, Eq. (18), which is the limit of the sequence of functions $\tilde{\rho}_{L}$ [52]. Neglecting the finite-size corrections, which we shall analyze at the end of this subsection, Eq. (45) may be recast as a limit of ratio of two Fredholm determinants (denoted in the following by ‘Det’ to distinguish it from the finite-size determinant, which is denoted by ‘det’),

[TABLE]

Here, $K_{1}(\lambda)$ is the kernel defined in Eq. (19e). The discussion below applies to any real-valued kernel under reasonable smoothness assumptions. In analogy to regular matrices, we used the notation $T$ for Toeplitz (resp. $H$ for Hankel) for the part that depends only on the difference (resp. sum) of the arguments $\lambda$ and $\mu$. The symbol $\hat{K}_{1}$ stands for the cosine transform of $K_{1}$, $\hat{K}_{1}(q)=2\int_{0}^{\infty}K_{1}(\lambda)\cos(q\lambda)\>{\rm d}\lambda$. Said differently, $T_{\Lambda}[\hat{K}_{1}]$ and $H_{\Lambda}[\hat{K}_{1}]$ act as shown in Eqs. (46b) and (46c), where $K_{1}$ is the inverse cosine transform of $\hat{K}_{1}$. Such types of determinants have been widely studied in the mathematical literature. In particular, the limit (46a) follows from Theorem 5 of Ref. [62]. We obtain

[TABLE]

Computing the integral $\hat{K}_{1}(0)=2\int_{0}^{\infty}K_{1}(\lambda)\,{\rm d}\lambda$ in Eq. (47) is straightforward. We find

[TABLE]

where $K$ is the Luttinger parameter of Eq. (2). We note in passing that the large $\Lambda$ limit of each determinant in (46a) may also be determined separately using the aforementioned theorem. We find (for any smooth integral operator $V$)

[TABLE]

with $\hat{S}(q)=\ln[1+\hat{V}(q)]$ and the inverse cosine transform $S(\lambda)=\frac{1}{\pi}\int_{0}^{\infty}\hat{S}(q)\cos(q\lambda)\,{\rm d}q$.

In the axial case, the algebra is only slightly different. Formula (46a) still holds but the kernel now acts on $L^{2}([0,\pi/2])$ instead of $L^{2}([0,\Lambda])$. Another difference is that both determinants are finite. We obtain

[TABLE]

Note that even though the kernel itself depends on $\Delta=\cosh(\eta)$, see Eq. (19e), the final result in the axial regime does not. To summarize, we have found

[TABLE]

In all three cases, the order one contribution is half of the order one contribution of $\ln(\mathcal{P}_{L})$ with opposite sign (see Sec. 3.3.1).

We also performed numerical checks of these formulas by computing the exact ratio (46a) by means of the Gaussian quadrature method, demonstrated in Ref. [64] and successfully applied, for instance, in Ref. [65] in the context of low-temperature, large-distance asymptotics of two-point functions of the XXZ model. In the planar regime, the kernels are regular and decay very fast at infinity. We find that the ratio of Fredholm determinants can be evaluated to almost arbitrary accuracy. For example, choosing $\Lambda=30$ and a few hundred points for the Gaussian quadrature, typically allows to confirm Eq. (50) to more than $30$ digits of numerical accuracy. The situation becomes less favorable when approaching the isotropic limit from below. In fact, exactly at the isotropic point the kernel decays algebraically, and we were able to confirm the result $\lim_{L\to\infty}\mathcal{R}_{L}=\sqrt{2}$ only up to relative error of the order $10^{-7}$. Finally, in the axial case, the numerical evaluation of the ratio of Fredholm determinants becomes tremendously accurate, presumably due to the fact that the kernels now act on a finite interval. This holds even close to the isotropic point. Choosing, for instance, $\Delta=1.01$ and $300$ quadrature points yields about $50$ digits of accuracy.

Another way to approach the ratio of determinants $\mathcal{R}_{L}$ is to compute it directly by solving the Bethe equations (15) and inserting the solution of the ground state Bethe roots into the exact finite-size formula (45). This procedure is limited by the computational time needed to solve the Bethe equations to a given accuracy. We solved Bethe equations of the three different cases up to $L=8192$. This way, we could numerically confirm that $\mathcal{R}_{L}$ scales indeed to the TDL result (50). Furthermore, at the same time, we got access to the finite-size corrections of $\ln(\mathcal{R}_{L})$. In the planar regime, we could numerically verify (see left panel of Fig. 4; regime $-1<\Delta<0$ not shown in the figure) that the finite-size correction behaves as $\sim L^{-2\delta}$, with $\delta=\text{min}\{1,4K-2\}$ being the exponent of the leading finite-size correction of the prefactor term $P(L)$.

At the isotropic point $\Delta=1$, the correction appears to be logarithmic, namely of the form $\sim 1/\ln(L)$ (see upper curve in the right panel of Fig. 4). In the axial regime, the finite size correction of $\ln(\mathcal{R}_{L})$ is exponentially small ($\sim e^{-cL}$, see inset of Fig. 5).

3.5 Results

We summarize the results of the previous subsections. In the planar regime $-1<\Delta<1$, we found

[TABLE]

Here, the decay rate $\alpha$ is determined by Eq. (55) and $K$ is the Luttinger parameter of Eq. (2). The exponent $\delta$, Eq. (42), as well as the amplitude $A$ describe the algebraic decay of the leading finite-size correction. Fig. 3 illustrates the behavior of the two parameters $\delta$ and $A$. The leading finite-size correction of $\ln(\mathcal{R}_{L})$ is completely dominated by the leading correction of $\ln(\mathcal{P}_{L})$ (see left panel of Fig. 4). By exponentiating the asymptotic expression for $\ln(\mathcal{O}_{L})$, we finally obtain the result (4).

The asymptotic formula (51) still holds in the isotropic limit, with the only difference that the leading finite-size correction is no longer algebraic but logarithmic. Since the correction of the prefactor term is of order $\mathcal{O}(\ln^{-2}(L))$, the leading correction of the overlap $\ln(\mathcal{O}_{L})$ stems here from the ratio of determinants, which is  $\sim 1/\ln(L)$ (see left panel of Fig. 4).

The result in the axial regime $\Delta>1$ (gapped phase) reads

[TABLE]

where $\alpha$ is determined by Eq. (58). The leading finite-size corretion is exponentially small, with the same $c=c(\Delta)$ for the corrections of $\ln(\mathcal{P}_{L})$ and $\ln(\mathcal{R}_{L})$. We further observe that $\lim_{\Delta\to 1}c(\Delta)=0$ and $\lim_{\Delta\to\infty}c(\Delta)=\infty$ (see inset of Fig. 5). The former is expected since the correction becomes logarithmic in the isotropic limit. The latter is also expected since there is no finite-size correction to the order one term in the Ising limit (see below Eq. (6)). Again, exponentiating the asymptotic expression for $\ln(\mathcal{O}_{L})$ yields the expected result, Eq. (6).

In Fig. 5 we show the order one term (black solid and dashed lines) in all three cases (planar, isotropic, axial) together with some finite-size calculations based on solutions of Bethe equations up to $L=8192$. Note that the order one term as function of anisotropy is not continuous at $\Delta=1$, a clear signature of the phase transition.

4 Conclusion

In this paper, we presented an extensive asymptotic analysis of the overlap between the Néel state and the ground state of the XXZ spin chain. We found that such an overlap always decays exponentially, which is expected for a non-trivial ground state of a quantum many body system. Nevertheless, we computed the corresponding decay rate exactly (see A). Subleading corrections to this behavior are most interesting, as they reveal important physical properties of the underlying quantum system at low temperature. In the gapped case (axial regime $\Delta>1$), for instance, the order one correction is exactly $1/2$ and independent of the anisotropy $\Delta$, reflecting the symmetry broken nature of the ground state.

We put particular emphasis on the gapless phase of the model (planar and isotropic cases), where our results are perhaps most significant. In Sec. 2 we argued (following Refs. [25, 26]) that the order one correction is universal and exactly given by $\sqrt{K}$ with the Luttinger parameter $K$. The derivation relies on a standard, but crucial CFT assumption. In an imaginary time picture, the Néel state plays the role of a boundary condition which renormalizes to a conformal invariant boundary condition [23], which is Dirichlet in our case. We checked this assumption in Sec. 3 by means of a microscopic approach, taking full advantage of the integrability of the XXZ spin chain. This approach is technically considerably more involved, but nevertheless doable using a recently derived exact finite-size formula for the overlap [31]. As we have shown, the order one term in the asymptotic expansion of the overlap is in perfect agreement with the prediction of CFT.

The evidence gathered in Sec. 3 for the correctness of the CFT formula in the XXZ spin chain is overwhelming. We checked the validity of identity (4) to machine precision (about $16$ digits) for several values of the anisotropy $\Delta$ by solving the NLIE equations (27) iteratively and subsequently computing the integral (32). In this sense, the order one term is exactly given by the TDL formulas (34) and (36), but its explicit expression (38) could not be rigorously derived yet (in a mathematical strict sence). Aside from this issue, the derivation of the CFT term presented here was fully analytical. Furthermore, we derived the asymptotic formula (6) of the axial (gapped) regime, in Sec. 3 by means of the NLIE approach as well as in B by means of the Euler-Maclaurin formula (EMF). Again, from a mathematically rigorous perspective, there is a bottleneck to prove formula (6) by means of the EMF. It lies in the TDL ground state root distribution (18) whose applicability to sums of functions that we consider in this paper, see Eqs. (20) and (19b), is not proven yet. We refer to Ref. [52] for a proof in case of slightly more regular functions. We hope that further elaboration on the ideas of Ref. [52] will close this gap (in the line of arguments of the EMF approach).

Another important feature of the order one term of the overlap is that it is universal, in the sense that a different microscopic model that renormalizes to the same CFT would have the same order one term. This can be beautifully illustrated by considering the simpler Haldane-Shastry chain[66, 67] which scales to a Luttinger liquid with Luttinger parameter $K=1/2$ (the corresponding point in XXZ is the isotropic point). The overlap between the ground state and any product state may be expressed in square root determinant form [68, 69], leading to the exact Néel overlap

[TABLE]

with decay rate $\alpha=\ln(2)-1/2\simeq 0.19314718$. As expected, the universal order one term is $\sqrt{K}=\sqrt{1/2}$, but the non-universal decay rate is different from the XXZ decay rate (57), in agreement with the discussion above.

We further extracted the form of the leading finite-size correction of the overlap, which takes the form of a power law with an exponent that we determined exactly to be $\delta={\rm min}\{1,4K-2\}$. At the isotropic point ($\gamma=0$, $K=1/2$), this turns into an $1/\ln L$ correction. The exponent and other subleading corrections can be derived, in principle, also from a perturbed CFT analysis. To do so, one needs to identify the least irrelevant bulk/boundary operators authorized by the symmetries, and then to compute the generated corrections. Such corrections depend on the scaling dimensions of the operators considered. The leading exponent, for instance, can be explained by considering perturbations by the stress-tensor (which gives $L^{-1}$) combined with a cosine term (which gives $L^{-4K+2}$). Such calculations typically require a little bit of hindsight from the underlying microscopic model. This analysis was not necessary here. The determination of the corrections follows from a standard reasoning based on the structure of the kernels that appear in the non-linear integral equations. This approach can also be generalized to analyze higher subleading terms. We note in passing that the knowledge of these corrections might prove useful in the study of certain basis-dependent Rényi-Shannon entropies [25, 26, 70], of which our overlap is a particular case. These entropies have also been studied using CFT, but the lack of tractable finite-size formulas makes the determination of subleading corrections even more crucial.

There are several other potential directions in which the present work may be generalized. For simplicity, we focused our attention on the ground state solely. However, the determinant formula we relied on holds for any eigenstate, and an asymptotic analysis for low-lying energy states should be also possible. In this case, however, boundary CFT does not predict any extra contributions. Still, this could be verified using the techniques presented in this paper. Another possible generalization concerns the asymptotic analysis of the partition function of the six-vertex model on a torus or cylinder, the latter being related [28] to the overlap we have studied. It would be interesting to try to push our analysis one step further and to derive the full CFT partition function (13) as an order one correction to the partition function of the discrete model in the TDL. A similar CFT analysis may be applied also to open chains [43, 26], where boundary contributions for low-lying excitations display a richer structure [44, 45]. This includes order one as well as power-law corrections originating from the Cardy-Peschel formula [42]. In that case, tractable formulas for the overlap with the Néel state are unfortunately not available yet. Another intriguing idea would be to go back to the ground state of the periodic XXZ chain and to realize a different conformal invariant boundary condition, e.g. a Neumann boundary condition. This can be done by computing the overlap of the ground state with a normalized equal-weight superposition of all product states [71]. We hope that the present paper stimulates future works in such directions.

Acknowledgements

We are grateful to F. Göhmann, Y. Ikhlef, M. Karbach, A. Klümper, K. K. Kozlowski, P. Pearce, and M. Rajabpour for useful discussions. MB acknowledges financial support by the DFG research group FOR-2316 funded by the DFG.

Appendix A Computation of the decay rate

The decay rate $\alpha=\alpha(\Delta)$ can be obtained by taking the TDL of the logarithm of the overlap $\mathcal{O}_{L}$ (or equivalently of the prefactor $\mathcal{P}_{L}$) divided by the system size $L$,

[TABLE]

Hence, for any value of $\Delta\geq-1$, it is determined by the integral of the logarithm of the prefactor function $p$, Eq. (19b), times the TDL ground state distribution $\rho(\lambda)$, Eq. (18). This relation has been often used in the past, though a rigorous proof is still lacking for our purpose. Note, however, that such type of results can be rigorously proved [52] for arbitrary anisotropy $\Delta>-1$ and for a wide class of functions $f$ slightly more regular than the one in Eq. (54).

Let us first consider the planar regime. The positive ground state Bethe roots are unbounded in the TDL (see e.g. Eq. (16) or Fig. 2). Therefore, we have to integrate in Eq. (54) over the whole positive real axis. Thus, the decay rate reads ($\Delta=\cos\gamma$)

[TABLE]

In particular, for $\Delta=\Delta_{m}=\cos\gamma_{m}$ with $\gamma_{m}=\pi/(2m+1)$, $m\in\mathbb{N}_{0}$, the integral in Eq. (55b) can be further simplified. By doing some manipulations (i.e. adding a factor $e^{ikx}$ under the intergal, taking the $x$-derivative, using the periodicity of the integrand in imaginary direction, choosing as a proper closed contour an infinitely long rectangle of height $i\pi$, applying Cauchy’s residue theorem, integrating with respect to $x$, identifying the correct integration constant by looking at the asymptotics for large imaginary $x$, and finally sending $x$ to zero), we find the following explicit expression,

[TABLE]

with the digamma function $\psi$ and Euler’s constant $\gamma_{E}$. Simple examples of Eq. (56) are $m=0,1,2$: $\alpha(-1)=\ln(2)$, $\alpha(1/2)=\ln(3\sqrt{3}/4)$, and $\alpha(\cos(\pi/5))=\ln[4/(1+\sqrt{5})]$. Further explicit expressions are $\alpha(0)=\ln(\sqrt{2})$ and $\alpha(1/\sqrt{2})=3\ln(2)/4-G/\pi$ with Catalan’s constant $G$.

The decay rate at the isotropic point $\Delta=1$ can be computed by taking the limit $\gamma\to 0$ in Eq. (55b). The result reads

[TABLE]

Let us now consider the axial regime. The positive ground state Bethe roots are bounded, $0<\lambda_{j}<\pi/2$, and $\rho$ is supported on the interval $\lambda\in[-\pi/2,\pi/2]$ (see e.g. Eq. (16) and Fig. 2). Thus, in Eq. (54), we now have to integrate from zero to $\pi/2$ rather than from zero to infinity as in the planar and isotropic cases. The decay rate can be simplified to ($\Delta=\cosh\eta$)

[TABLE]

It vanishes in the Ising limit $\Delta\to\infty$, $\lim_{\Delta\to\infty}\alpha(\Delta)=0$, as expected (see text below Eq. (6) in Sec. 1).

The curve $\alpha(\Delta)$ is shown in Fig. 6, together with the explicit values $\alpha(-1)$, $\alpha(0)$, and $\alpha(1)$. It is a smooth function in the whole regime $\Delta>-1$, even at the point $\Delta=1$.

Appendix B Approximating sums over ground state Bethe roots by means of the Euler-Maclaurin formula

In this section we study how sums over the positive ground state Bethe roots $\{\lambda_{j}\}_{j=1}^{L/4}$ with $0<\lambda_{1}<\ldots<\lambda_{L/4}$ can be approximated by means of the Euler-Maclaurin formula (see e.g. Ref. [53]). In particular, we are interested in sums like $\sum_{j=1}^{L/4}f(\lambda_{j})$ where $f$ is a function with logarithmic singularities at the endpoints, as e.g. in Eq. (21).

The two endpoints are given by the smallest and largest positive ground state Bethe roots, $\lambda_{1}$ and $\lambda_{L/4}$, respectively. Therefore, we have to analyze how these two Bethe roots scale in the TDL. To this end, we compare integrals of the discrete ground state distribution $\rho_{L}$, Eq. (17), with integrals of the corresponding TDL distribution $\rho$, Eq. (18), both over the same interval. In the planar regime, for instance, the primitive functions of $\rho_{L}$ and $\rho$ respectively read

[TABLE]

where $\Theta$ is the Heavyside step function, i.e. $\Theta(x)=0$ for $x<0$, $\Theta(0)=1/2$, and $\Theta(x)=1$ for $x>0$. Therefore, $\sigma_{L}(\lambda_{1})-\sigma_{L}(0)=1/(2L)$. Identifying this with $\sigma(\lambda_{1})-\sigma(0)$, we obtain $\lambda_{1}=\frac{2\gamma}{\pi}{\rm artanh}[\tan(\frac{\pi}{2L})]$ which implies $\lambda_{1}\simeq\gamma/L$ in the limit of large system size $L$. Similarly, in order to compute the large system size scaling of $\lambda_{L/4}$, we identify $\sigma_{L}(\infty)-\sigma_{L}(\lambda_{L/4})=1/(2L)$ with

[TABLE]

where we used in the second step that $e^{-\frac{\pi}{\gamma}\lambda_{L/4}}$ is small. Thus, we obtain $\lambda_{L/4}\simeq\frac{\gamma}{\pi}\ln(\frac{2L}{\pi})$. The isotropic limit $\gamma\to 0$ (after scaling all Bethe roots by $\gamma$) can be performed easily, yielding $\lambda_{1}\simeq 1/L$ and $\lambda_{L/4}\simeq\frac{1}{\pi}\ln(\frac{2L}{\pi})$. The analysis of the scaling in the axial regime is not much more involved and proceeds similarly. Without presenting the details of this calculation, we refer to the final result, Eq. (16).

Let us now consider a sum $\sum_{j=1}^{L/4}f(\lambda_{j})$ over the positive ground state Bethe roots $\{\lambda_{j}\}_{j=1}^{L/4}$. We are interested in the large system size scaling of this sum. In particular, we would like to compute the extensive part, which is proportional to $L$, as well as the order one contribution by means of the Euler-Maclaurin formula (EMF). Denoting the smallest and largest positive ground state Bethe roots by $\epsilon_{L}=\lambda_{1}$ and $\Lambda_{L}=\lambda_{L/4}$, respectively, a slightly modified version of the EMF reads [53]

[TABLE]

Here, $B_{2n}$ are the Bernoulli numbers and g(\lambda)\big{|}^{\Lambda}_{\epsilon}=g(\Lambda)-g(\epsilon). The main modification as compared to the original EMF (see e.g. Ref. [55]) is that we replaced the measure ${\rm d}j$ under the integral by $L\rho(\lambda){\rm d}\lambda$ as well as the derivative $\frac{\rm d}{{\rm d}j}$ by $\frac{1}{L\rho(\lambda)}\frac{\rm d}{{\rm d}\lambda}$.

In the following analysis, we use for $\rho$ the TDL ground state distribution of Eq. (18). This is of course a delicate point [52] since we completely neglect finite-size corrections to the root distribution $\rho$. The above replacements are fine for functions that are regular in the interval $[\epsilon_{L},\Lambda_{L}]$ (see Ref. [52] for more details), in particular at the endpoints $\lim_{L\to\infty}\epsilon_{L}=0$ and $\lim_{L\to\infty}\Lambda_{L}=\Lambda_{\infty}$. Here, regular means that $f$ is bounded and almost constant in the intervals $[0,\epsilon_{L}]$ and $[\Lambda_{L},\Lambda_{\infty}]$. Then, the sum in the second line of Eq. (62) is subleading. Further, the integral in Eq. (62) can be simplified in the following way, using that $\int_{0}^{\epsilon_{L}}\rho(\lambda)\,{\rm d}\lambda=\int_{\Lambda_{L}}^{\Lambda_{\infty}}\rho(\lambda)\,{\rm d}\lambda=1/(2L)$,

[TABLE]

Hence, the exptensive part of the left-hand side of Eq. (62) is given by the first term, $L\int_{0}^{\Lambda_{\infty}}f(\lambda)\rho(\lambda){\rm d}\lambda$, and there is in total no order one contribution. We checked this for several examples of functions $f$ by solving Bethe equations and computing the sum of the left-hand side of Eq. (62) explicitly.

On the contrary, for a function $f$ that has logarithmic singularities at the endpoints, e.g. $f=\ln(p)$, the analysis is more involved. Let us first focus on the lower boundary, since the calculation is basically the same in all three cases ‘planar’, ‘isotropic’, and ‘axial’. Let us assume that the asymptotics of $f$ for small positive arguments is given by $f(\lambda)\simeq-2\ln(\lambda)+\mathcal{O}(\lambda)$, which is fulfilled for $f=\ln(p)$. Then, we find

[TABLE]

In the very last step, we used the expansion of the log-gamma function (see e.g. Ref. [54])

[TABLE]

with argument $x=L\rho(\epsilon_{L})\epsilon_{L}\simeq 1/2$. Therefore, the contribution to the order one term of the sum $\sum_{j=1}^{L/4}\ln[p(\lambda_{j})]$ from the lower endpoint is in total  $-\ln(2)$. The same line of arguments holds for the upper endpoint in the axial regime, $\Lambda_{L}\simeq\pi/2-\pi\eta/(2kK\mspace{-1.0mu}(k)L)$, see Eq. (16), yielding as a second contribution again  $-\ln(2)$. This confirms the order one term (40), which we also rigorously derived in Sec. 3.3.1 by a completely independent approach (NLIE). The upper endpoint in the two other cases (planar and isotropic), $\Lambda_{\infty}=\infty$, causes serious problems. The main obstacle is that the finite-size corrections of $\rho$, which are neglected in Eq. (62), become important for functions that scale for large arguments $\lambda$ as  $\sim\ln(\lambda)$. Again, we have checked this by solving Bethe equations for finite, but large system sizes up to $L=8192$. Computing the finite-size corrections of $\rho$ is a difficult task.

In conclusion, the EMF approach to compute the sum $\sum_{j=1}^{L/4}\ln[p(\lambda_{j})]$ works well in the axial regime but, unfortunately, leads to unsurmountable difficulties in the planar regime and at the isotropic point.