Full counting statistics in the spin-1/2 Heisenberg XXZ chain

TL;DR

This paper investigates the probability distributions of subsystem magnetization in the spin-1/2 Heisenberg XXZ chain, revealing universal scaling functions and quantifying quantum fluctuations in its critical ground state.

Contribution

It introduces a method to determine universal scaling functions of magnetization distributions using free fermion techniques and boundary sine-Gordon model relations.

Findings

Identification of scaling behavior in magnetization distributions

Derivation of universal scaling functions

Quantification of quantum fluctuations in the ground state

Abstract

The spin-1/2 Heisenberg chain exhibits a quantum critical regime characterized by quasi long-range magnetic order at zero temperature. We quantify the strength of quantum fluctuations in the ground state by determining the probability distributions of the components of the (staggered) subsystem magnetization. Some of these exhibit scaling and the corresponding universal scaling functions can be determined by free fermion methods and by exploiting a relation with the boundary sine-Gordon model.

Figures

Tables

Table 1. Table 1: Energy densities from iTEBD and exact formula ( 39 ) for several values of Δ Δ \Delta .

$\Delta$	$E_{\mathrm{MPS}}$	$E_{\mathrm{GS}}$	$\delta E(\times {10}^{-8})$
$-0.8$	$-0.256339667$	$-0.256339677$	$0.96$
$-0.6$	$-0.267640618$	$-0.267640628$	$1.03$
$-0.4$	$-0.282089887$	$-0.282089903$	$1.62$
$-0.2$	$-0.299086657$	$-0.299086680$	$2.23$
$0$	$-0.318309858$	$-0.318309886$	$2.84$
$0.2$	$-0.339564266$	$-0.339564304$	$3.84$
$0.4$	$-0.362727187$	$-0.362727227$	$4.06$
$0.6$	$-0.387725863$	$-0.387725910$	$4.71$
$0.8$	$-0.414528779$	$-0.414528832$	$5.35$

Table 2. Table 2: Numerical values for the variance of transverse fluctuations extracted from the two-point function.

$\Delta$	$-0.8$	$-0.6$	$-0.4$	$-0.2$	$0$	$0.2$	$0.4$	$0.6$
$s_x$	$0.101$	$0.081$	$0.064$	$0.049$	$0.036$	$0.024$	$0.013$	$0.005$
DMRG	$0.101$	$0.079$	$0.063$	$0.049$	$0.037$	$0.026$	$0.017$	$0.009$

Equations

$$H=J\sum _{j=1}^L{S}_j^x{S}_{j+1}^x+S_j^y{S}_{j+1}^y+\Delta S_j^z{S}_{j+1}^z.$$

$$⟨\mathrm{GS}|S_{j+n}^x{S}_j^x|\mathrm{GS}⟩$$

$$⟨\mathrm{GS}|S_{j+n}^z{S}_j^z|\mathrm{GS}⟩$$

$$\Delta =-\cos (\pi \eta ).$$

$$S^{\alpha }(\ell )=\sum _{j=1}^{\ell }{S}_j^{\alpha },N^{\alpha }(\ell )=\sum _{j=1}^{\ell }(-1{)}^j{S}_j^{\alpha }.$$

$$P_S^{\alpha }(m,\ell )$$

$$P_N^{\alpha }(m,\ell )$$

$$G_{\ell }^{\alpha }(\theta )\equiv ⟨\mathrm{GS}|e^{i\theta S^{\alpha }(\ell )}|\mathrm{GS}⟩,F_{\ell }^{\alpha }(\theta )\equiv ⟨\mathrm{GS}|e^{i\theta N^{\alpha }(\ell )}|\mathrm{GS}⟩.$$

$$X_{\ell }^{\alpha }(0)$$

$$e^{i\theta S^{\alpha }(\ell )}=\prod _{j=1}^{\ell }\left[\cos (\theta /2)+i\sin (\theta /2){\sigma }_j^{\alpha }\right].$$

$${\tilde{X}}_{\ell }^{\alpha }(r)={\int }_{-\pi }^{\pi }\frac{d\theta }{2\pi }{e}^{-ir\theta }{X}_{\ell }^{\alpha }(\theta ),X=F,G,$$

$$P_{N}^{\alpha}(m,\ell)=\begin{cases}\sum_{r\in\mathbb{Z}}\widetilde{F}^{\alpha}_{\ell}(r)\ \delta(m-r)&\text{if\ }\ell\text{ is even,}\\ \sum_{r\in\mathbb{Z}}\widetilde{F}^{\alpha}_{\ell}\big{(}r+\frac{1}{2}\big{)}\ \delta\big{(}m-r-\frac{1}{2}\big{)}&\text{if\ }\ell\text{ is odd.}\end{cases}$$

$$⟨\mathrm{GS}|S^{\alpha }(\ell )|\mathrm{GS}⟩=0=⟨\mathrm{GS}|N^{\alpha }(\ell )|\mathrm{GS}⟩.$$

$$\langle{\rm GS}|\big{(}S^{\alpha}(\ell)\big{)}^{2}|{\rm GS}\rangle=\ell s_{\alpha}+o(\ell)\ ,\quad\langle{\rm GS}|\big{(}N^{\alpha}(\ell)\big{)}^{2}|{\rm GS}\rangle=\ell n_{\alpha}+o(\ell)\ .$$

$$s_{x}=s_{y}=\lim_{L\to\infty}\frac{1}{L}\langle{\rm GS}|\big{(}S^{x}_{L}\big{)}^{2}|{\rm GS}\rangle\ ,\quad n_{\alpha}=\lim_{L\to\infty}\frac{1}{L}\langle{\rm GS}|\big{(}N^{\alpha}_{L}\big{)}^{2}|{\rm GS}\rangle.$$

$$\frac{\big{\langle}\left(S^{\alpha}(\ell)\right)^{2n}\big{\rangle}}{\big{\langle}\left(S^{\alpha}(\ell)\right)^{2}\big{\rangle}^{n}}\ ,\qquad\frac{\big{\langle}\left(N^{\alpha}(\ell)\right)^{2n}\big{\rangle}}{\big{\langle}\left(N^{\alpha}(\ell)\right)^{2}\big{\rangle}^{n}}\ .$$

$$\langle{\rm GS}|e^{i\theta_{Y}Y^{\alpha}(\ell)}|{\rm GS}\rangle\ ,\quad\theta_{Y}=\frac{\theta}{\sqrt{\langle{\rm GS}|\big{(}Y^{\alpha}(\ell)\big{)}^{2}|{\rm GS}\rangle}}\ ,\quad Y=S,N,$$

$$\mathcal{H}(\Delta )=\frac{v}{2}\int dx\left[K({\partial }_x\theta {)}^2+\frac{1}{K}({\partial }_x\varphi {)}^2\right]+\dots ,$$

$$v=\frac{\pi }{2}\frac{\sqrt{1-{\Delta }^2}}{\arccos \Delta },K=\frac{\pi }{2}\frac{1}{\pi -\arccos \Delta },$$

$$S_j^z$$

$$S_j^x$$

$$S^z(\ell )$$

$$N^z(\ell )$$

$$N^x(\ell )$$

$$G_{\ell }^z(\theta )$$

$$F_{\ell }^z(\theta )$$

$$F_{\ell }^x(\theta )$$

$$F_{\ell }^z(\theta )$$

$${\mathbb{C}}_{nm}=\frac{\sin\big{(}\frac{\pi}{2}(m-n)\big{)}}{\pi(m-n)}\ .$$

$$F^{z}_{\ell}(\theta)=\big{(}\cos(\theta)\big{)}^{\frac{\ell}{2}}\sum_{m=-\infty}^{\infty}(-1)^{m\ell}(2\ell)^{-2\big{(}m+\beta(\theta)\big{)}^{2}}G^{2}\big{(}m+1+b(\theta)\big{)}G^{2}\big{(}1-m-b(\theta)\big{)}\left[1+\frac{{c}_{2}\big{(}m+b(\theta)\big{)}}{\ell^{2}}+\dots\right]$$

$$b(\theta )$$

$$F_{\ell }^z(\theta )\sim e^{-z^2/4}.$$

$$F^{z}_{\ell}(\theta)\sim\sqrt{2}\,G^{2}(1/2)G^{2}(3/2)\,e^{-y^{2}/4}\times\left\{\begin{array}[]{ll}(-1)^{\ell/2}\ell^{-1/2}\,,&\quad\quad\ell{\rm\quad even}\\ &\\ (-1)^{(\ell-1)/2}\,\ell^{-1}[c+\log(2\ell)]\;y/\pi\,,&\quad\quad\ell{\rm\quad odd}\end{array}\right.$$

$$F_{\ell }^x(\theta )$$

Full text

I Introduction

Universality is a key organizing principle for continuous phase transitionscardy ; sachdev . It posits that certain quantities are independent of microscopic details and coincide in different physical systems that belong to the same “universality class”. The latter are determined by properties such as symmetries and dimensionality and are amenable to field theory descriptions. In 1+1 dimensions this permits the exact description of universal properties such as critical exponents and correlation functions at conformally invariant quantum critical points. As emphasized in Ref. LF, , less familiar quantities like the order parameter probability distribution function display universal scaling as well. In quantum theory these probability distributions describe the statistics of measurements on identical systems, which generally give rise to different outcomes. Their analysis provides very detailed information about the physical properties of many-particle systems and has been explored in a variety of areas including condensed matterCM ; CM2 and cold atom physicsCA1 ; CA2 ; CA3 ; CA4 . Theoretical results on full counting statistics in quantum critical systems are relatively scarce. The list of available results includes phase fluctuations in Luttinger liquidsgritsev ; IGD06 ; IGD08 ; lovas , the order parameter statistics in the Ising field theoryLF , the transverse magnetization in the Ising chaincd-07 and the magnetization in the Haldane-Shastry modelstephan . Here we consider the (staggered) subsystem magnetization in the anisotropic one-dimensional spin-$1/2$ Heisenberg XXZ chain

[TABLE]

The XXZ chain is a paradigmatic model for quantum critical behaviour in 1+1 dimensions. It features a critical line parametrized by the exchange anisotropy $-1\leq\Delta\leq 1$. The special values $\Delta=\pm 1$ correspond to the isotropic antiferromagnet and ferromagnet respectively. In the regime $-1<\Delta\leq 1$ the low-energy behaviour of the model (1) is described by Luttinger liquid theory or equivalently a free, compact bosonLutherPeschel ; affleck89 ; GNT ; gia . The long-distance asymptotics of spin-spin correlation functions is of the form

[TABLE]

where explicit expressions for the amplitudes in (2) are known lukyanov ; Affleck98 ; lukyanov2 ; lukyanovterras and $\eta$ is related to the anisotropy parameter $\Delta$ by

[TABLE]

It follows from (2) that throughout the critical regime the dominant correlations are those of the staggered magnetizations in the xy-plane. The XXZ chain thus exhibits antiferomagnetic quasi-long range order in the XY plane in spin space. Two-point functions such s (2) are a standard means for characterizing physical properties and identifying ground state “phases” in quantum critical systemsgia . A key objective of our work is to provide a complementary characterization of ground state properties in the critical XXZ chain by determining the quantum mechanical fluctuations of the subsystem magnetization in the ground state. More precisely we consider the probability distributions of the following observables

[TABLE]

The quantities $S^{\alpha}(\ell)$ and $N^{\alpha}(\ell)$ describe the smooth and staggered components of the $\alpha$-component of the magnetization of the subsystem consisting of sites $1$ to $\ell$, where $\ell\ll L$. We note that whereas $S^{z}(L)$ is a conserved quantity, $S^{z}(\ell)$ is not. The probabilities of the observables (4) taking some value $m$ when the system is prepared in the ground state and a measurement is then performed are

[TABLE]

As we have already mentioned, probability distributions like (5) are experimentally measurable in cold atom experiments. The central objects of our analysis are the generating functions of the moments of the probability distributions (5) are

[TABLE]

It is easy to see that they have the following properties

[TABLE]

The last relation allows us to restrict our attention to the interval $0\leq\theta<2\pi$ and can be obtained e.g. from the representation

[TABLE]

Defining

[TABLE]

the probability distributions of interest can be expressed as

[TABLE]

An analogous equation holds for $P_{S}^{\alpha}(m,\ell)$.

I.1 Moments of the probability distributions

As we are not imposing a magnetic field and spontaneous symmetry breaking of the U(1) symmetry of the Heisenberg Hamiltonian is forbidden in one spatial dimension, translational invariance implies that the averages of $S^{\alpha}(\ell)$ and $N^{\alpha}(\ell)$ vanish

[TABLE]

The variances have the following asymptotic expansions for large sub-system sizes $\ell$

[TABLE]

For sufficiently large values of $\ell$ we expect the coefficients $s_{\alpha}$ and $n_{\alpha}$ to be equal to the corresponding quantities for the entire system, i.e.

[TABLE]

As $S^{z}(L)$ is a conserved quantity and our system is translationally invariant we have $s_{z}=0$. It is instructive to consider the calculation of the variance of the subsystem magnetization by field theory methods. As the variances are non-universal quantities they are expected to be susceptible to short-distance physics, and this is indeed borne out by the explicit calculation summarized in Appendix A.

While the moments themselves depend on microscopic details, certain ratios can be universalprivman ; goldenfeld ; LF . In particular one may expect the following ratios to exhibit universal behaviour

[TABLE]

If these ratios are universal, the modified generating functions

[TABLE]

will be universal functions of the parameter $\theta$. This means in particular that they can be calculated by field theory methods. In practice (15) tells us that the moment generating functions calculated from field theory and computed directly in the lattice model should agree up to an overall rescaling of the parameter $\theta$.

II Field theory description of the XXZ chain

It is well established that the long distance behaviour of local equal time correlation functions in the critical XXZ chain is well described by (perturbed) Luttinger liquid theoryAGSZ89 ; lukyanov ; Affleck98 ; lukyanovterras ; Lyon . In absence of a magnetic field the Hamiltonian can be cast in the form

[TABLE]

where $\phi$ and $\theta$ are Bose fields with commutation relations $[\phi(t,x),\theta(t,y)]=(i/2){\rm sgn}(x-y)$, the dots indicate perturbations that are irrelevant in the renormalization group sense and

[TABLE]

The bosonization formulas for the spin operators are

[TABLE]

where $a_{0}$ is the lattice spacing and the amplitudes $b_{0}$, $c_{1}$, $b_{1}$ are known exactly lukyanovterras . For large subsystem sizes we thus have

[TABLE]

Applying the bosonization prescription to our generating functions and ignoring subleading terms we obtain

[TABLE]

where $|0\rangle$ is the Fock vacuum. The representation (21) reveals that $G^{z}_{\ell}(\theta)$ maps onto a simple vertex operator two-point function in the free boson theory, whereas $F^{\alpha}_{\ell}(\theta)$ correspond to expectation values of non-local operators. The alert reader will have noted that we did not provide a bosonized expression for $G^{x}_{\ell}(\theta)$. The reason is that the field theory calculation of $G^{x}_{\ell}(\theta)$ is easier in a somewhat different setup and we return to this issue in section IV.2.1.

III Generating functions for the staggered subsystem magnetization

We start by considering the probability distributions of the staggered subsystem magnetizations $P^{\alpha}_{N}(m,\ell)$ and the corresponding generating functions $F^{\alpha}_{\ell}(\theta)$. We first present analytic results in certain limits and then compare these to numerical ones.

III.1 The XX point $\Delta=0$

At the XX point the Heisenberg model can be mapped to non-interacting spinless fermions by means of a Jordan-Wigner transformation. Using standard techniquesLSM we can derive the following determinant representation for the longitudinal generating function $F^{z}_{\ell}(\theta)$

[TABLE]

Here $\mathbb{C}_{nm}$ is the correlation matrix of the free fermion chain obtained by the Jordan Wigner transformation

[TABLE]

The matrix $\mathbb{B}_{nm}$ is a Toeplitz matrixremark1 . Its properties have been analyzed in great detail in the context of entanglement entropies in Ref. ce-10, . For large values of the subsystem size $\ell$ one obtains the following asymptotic expansionce-10

[TABLE]

where

[TABLE]

It follows from (24) that $F^{z}_{\ell}(\theta)$ is very small except in the vicinities of $\theta\approx 0,\pi$. To analyze the behaviour in these regions it is useful to define the following scaling limits:

• S1:

$\theta\to 0$, $\ell\to\infty$, while keeping $z=\theta\ell^{1/2}$ fixed.

In this regime the generating function reduces to a simple Gaussian in the scaling variable

[TABLE]

• S2:

$\theta\to\pi$, $\ell\to\infty$, while keeping $y=(\pi-\theta)\ell^{1/2}$ fixed.

In this regime the behaviour depends on the parity of the subsystem size

[TABLE]

where $c=2\log(2)+\gamma_{E}$ and $\gamma_{E}\approx 0.577216$.

We will see in the following that the two limits S1 and S2 are useful for analyzing numerical results for $F^{z}_{\ell}(\theta)$.

III.2 Field theory approach

In the field theory approach we are tasked with evaluating the expressions (21) for $F^{\alpha}_{\ell}(\theta)$. This can be done by following the analysis of Refs gritsev, ; IGD06, ; IGD08, , which considered generating function for phase fluctuations in Luttinger liquids. Expanding in powers of $\theta$ we obtain

[TABLE]

where

[TABLE]

As the leading singularities of the integrand occur when $u_{j}\approx v_{k}$ we now make the further approximation

[TABLE]

This results in

[TABLE]

The right hand side of (31) is equal to the partition function of a boundary sine-Gordon model, for which exact results are available in the literatureFLS1 ; FLS2 ; FS ; BLZ1 ; DT99 ; BLZ2 . Using a result of Ref. BLZ2, for the right hand side of (31) one has

[TABLE]

The function $A^{({\rm vac})}(\lambda)$ can be computed very efficiently from the solution of the single-particle Schrödinger equation

[TABLE]

Denoting by $\Psi^{+}(x,E)$ and $\chi^{+}(x,E)$ the solutions to (33) with asymptotics

[TABLE]

we have

[TABLE]

where $W[f,g]=f\partial_{x}g-g\partial_{x}f$ denotes the Wronskian.

Let us now turn to the longitudinal generating function $F^{z}_{\ell}(\theta)$. It has an integral representation

[TABLE]

This expression needs to be regularized because

[TABLE]

and $K$ ranges from $1/2$ at the isotropic point $\Delta=1$ to infinity when the ferromagnetic point is approached ($\Delta\to-1$). The right-hand side of (36) can again be related to the partition function of a boundary sine-Gordon model, but the boundary interaction for $\Delta<1$ is now irrelevant. This suggests that $F^{z}_{\ell}(\theta)$ will be independent of $\Delta$ and equal to the result at $\Delta=0$, i.e.

[TABLE]

Eqns (32), (35) and (38) provide us with explicit expressions for the generating functions $F^{\alpha}_{\ell}(\theta)$ that now can be compared to numerical results for the lattice model.

III.3 Numerical Method

Our numerical approach is based on the iTEBD algorithmiTEBD ; iTEBD2 . A translationally invariant MPS representation of the ground state of the XXZ Hamiltonian is obtained as follows. We initialize the system in the simple product state $\bigotimes_{j\in\mathbb{Z}}(|\uparrow\rangle_{j}+|\downarrow\rangle_{j})/\sqrt{2}$, which admits an MPS representation with auxiliary dimension $\chi=1$. We then evolve the state in imaginary time by the operator $\exp(-\tau H)$ by means of a second order Suzuki-Trotter decomposition with imaginary time-step $\tau J=10^{-3}$. During imaginary time evolution the MPS loses its canonical form, which we then restore before taking expectation values of operators.

In order to control the convergence of the imaginary time algorithm we keep track of the energy density. In practice we run the algorithm until the energy density becomes stationary (within machine precision). We repeat this procedure with auxiliary dimensions of up to $\chi=256$. In Table 1 we present our best estimates for the ground-state energy densities for different values of $\Delta$ and compare them to the known exact result

[TABLE]

Our numerical results for the energy densities differ from the exact values by $O(10^{-8})$, which is quite satisfactory given that the model is gapless. Once we have obtained the MPS description of the ground state, we can straightforwardly evaluate the generating functions $F^{\alpha}_{\ell}(\theta)$ and $G^{\alpha}_{\ell}(\theta)$ with a computational cost that scales as $O(\ell\chi^{3})$. A useful check on the numerical accuracy of our results can be obtained by considering the noninteracting case $\Delta=0$, where the exact determinant formula (22) for the generating function of the longitudinal staggered magnetization is available. The discrepancy between the iTEBD data and the exact result increases as expected with the subsystem size $\ell$. However up to subsystem sizes of $\ell=200$ the relative error of our iTEBD result is less than $0.1\%$.

III.4 Numerical results for the transverse generating function

$F^{x}_{\ell}(\theta)$

Numerical results for $F^{x}_{\ell}(\theta)$ as a function of $\theta$ for several values of the subsystem size $\ell$ are shown in Fig. 1.

We see that the generating function is very small everywhere except in the vicinities of $\theta=0,\pi$. We also observe that the oscillatory behaviour as a function of $\theta$ becomes more pronounced in the attractive regime $\Delta<0$.

Based on the field theory analysis of section III.2 we expect the $\theta\approx 0$ regime to exhibit scaling with a universal scaling function given by (32), (35)

[TABLE]

Here ${\cal F}^{x}_{0}(z)$ is related to the function $A^{({\rm vac})}(\lambda)$ in (35) by

[TABLE]

where $c$ is a non-universal $\Delta$-dependent constant that arises from the fact that while the ratios (14) are universal, the second moment itself is not, cf. the discussion preceding eqn (15). In practice we determine $c$ by carrying out a best fit of our numerical data to (41). In Fig. 2 we present a comparison of our numerical results for $F^{x}_{\ell}(\theta)$ to the field theory prediction (32), (35). We see that the numerical data exhibits scaling collapse and the agreement with the theoretical scaling function is clearly very good. This holds for all values of $\Delta$ we have considered in the critical regime $-1<\Delta\leq 1$. We again see that in the attractive regime $\Delta<0$ the oscillatory behaviour away from $\theta=0$ becomes more pronounced.

We now turn to the other region in which $F^{x}_{\ell}(\theta)$ is sizeable, namely $\theta\approx\pi$. Interestingly, as shown in Figs 3 and 4, we observe scaling behaviour here as well. There is a strong parity effect in the subsystem size $\ell$ which requires us to consider even and odd $\ell$ separately.

Our numerical data in the vicinity of $\theta=\pi$ is well described by the scaling ansatz

[TABLE]

where $e/o$ refers to even and odd subsystem size $\ell$ respectively. Inspection of Figs 3 and 4 shows that the ansatz is in excellent agreement with the data. We note that at $\theta=\pi$ the numerical data (for $\ell$ even) exhibit a perfect algebraic decay $\sim\ell^{-1/4}$, independent of the value of the interaction $\Delta$. The form (42) suggests that for very large subsystem sizes in the thermodynamic limit the feature at $F^{x}_{\ell}(\theta\approx\pi)$ becomes less and less important compared to $F^{x}_{\ell}(\theta\approx 0)$. At present no analytic results on $F^{x}_{\ell}(\theta\approx\pi)$ are known. It should in principle be possible to calculate $F^{x}_{e/o}(z)$ using field theory methods.

III.5 Probability distribution $P^{x}_{N}(m,\ell)$ of the

transverse, staggered subsystem magnetization

We are now in a position to determine the probability distribution $P^{x}_{N}(m,\ell)$ from the generating function $F^{x}_{\ell}(\theta)$ using eqn (10). In Fig. 5 we show results for $P^{x}_{N}(m,\ell)$ as a function of $m$ for several values of the exchange anisotropy $\Delta$ and subsystem sizes $\ell$. As $P^{x}_{N}(m,\ell)$ is a sum over $\delta$-functions, cf. eqn (10), we plot the corresponding weights $\widetilde{F}_{\ell}^{x}$ at the appropriate values of $m$.

We observe that

There is a strong even/odd effect in $m$. The results for even and odd $m$ follow different smooth curves. The separations between even and odd curves slowly tend to zero as $\ell^{-1/4}$ as the subsystem size $\ell$ is increased. 2. 2.

There is a weaker even/odd effect in the subsystem size $\ell$. This effect remains visible even for the large subsystem sizes we consider here. The magnitude of this effect grows with $\Delta$ and is strongest for $\Delta\to 1$, i.e. when we approach the isotropic antiferromagnet. 3. 3.

The probability distributions are quite broad, implying strong quantum fluctuations in the staggered transverse subsystem magnetization. 4. 4.

The width of $P^{x}_{N}(m,\ell)$ increases as the interaction becomes more attractive and the distribution flattens. 5. 5.

The distribution for attractive and moderately repulsive interactions is bi-modal, while for $\Delta\approx 1$ it displays a single maximum.

These observations can be understood in terms of our scaling analysis of the generating function $F^{x}_{\ell}(\theta)$. The probability distribution $P_{N}(m,\ell)$ is dominated by the behaviour of $F^{x}_{\ell}(\theta)$ in the regions $\theta=0,\pi$, and exploiting the observed scaling behaviour of these contributions we conclude that

[TABLE]

where $\mathcal{\widetilde{F}}^{x}_{0}$ and $\mathcal{\widetilde{F}}^{x}_{e/o}$ are obtained by Fourier transforming the functions ${\cal F}^{x}_{0}(z)$ and ${\cal F}^{x}_{e/o}(z)$ that describe the scaling behaviour of the generating function around $\theta=0$ and $\theta=\pi$ respectively. For large subsystem sizes $\ell$ the even/odd effect in $m$ disappears and we are left with $\ell^{\frac{\eta}{2}-1}\mathcal{\widetilde{F}}^{x}_{0}(m/\ell^{1-\eta/2})$, which can be calculated exactly using the boundary sine-Gordon mapping. The corresponding contribution is shown by a solid red line in Figs 5 and 6. We see that for attractive and moderately strong repulsive interactions there is an enhanced probability to form a large positive or negative staggered moment in the xy-plane. However, this enhancement is not particularly pronounced. The effect is strongest close to the ferromagnet at $\Delta=-1$ as can be seen in Fig. 6, which presents results for $P^{x}_{N}(m,\ell)$ at $\Delta=-0.95$. We observe that for large subsystem sizes $\ell$ the probability distribution becomes fairly flat over most of allowed range of staggered magnetizations $-\ell/2\leq m\leq\ell/2$ except for an enhancement close to the maximal possible values $m\approx\pm\ell/2$.

III.6 Longitudinal generating function $F^{z}_{\ell}(\theta)$

and probability distribution $P^{z}_{N}(m,\ell)$

We now turn to the longitudinal generating function $F^{z}_{\ell}(\theta)$. As we will see, its behaviour is rather different from $F^{x}_{\ell}(\theta)$. According to the field theory approach discussed in section III.2 we expect $F^{z}_{\ell}(\theta)$ to be described by the scaling function

[TABLE]

where $\gamma$ is a $\Delta$-dependent constant that encodes the fact that appropriate ratios of moments are universal, while the second moment itself is not, cf. eqn (15). Numerical results for $F^{z}_{\ell}(\theta)$ are shown in Fig. 7 and are seen to be in excellent agreement with the scaling form (44).

The coefficient $\gamma$ is found to be consistent with

[TABLE]

We conclude that the fluctuations of $N^{z}(\ell)$ have a very simple form: the second moment is

[TABLE]

while all higher cumulants vanish. We now turn to the behaviour of $F^{z}_{\ell}(\theta\approx\pi)$. Guided by the exact result (27) for $\Delta=0$ we have attempted to describe our numerical data by the ansatz

[TABLE]

Here $A$, $B$, $\alpha$, $\beta$ and $c$ are $\Delta$-dependent parameters that we fix by considering the $\ell$-dependencies of $F^{z}_{\ell}(\pi)$ (for $\ell$ even) and $\partial_{\theta}F^{z}_{\ell}(\theta)$ (for $\ell$ odd). Our numerical results suggest that

[TABLE]

In Figs 8 and 9 we compare numerical results for $F^{z}_{\ell}(\theta\approx\pi)$ for several values of $\Delta$ and $\ell$ to the scaling ansatz (47), (48).

The agreement is seen to be quite satisfactory in all cases. Having determined the generating function $F^{z}_{\ell}(\theta)$ we can now use it to obtain the probability distribution of the longitudinal staggered subsystem magnetization $P_{N}^{z}(m,\ell)$ by Fourier transform. Results for several values of $\Delta$ and $\ell$ are presented in Fig. 10.

The probability distribution is again a sum over delta-functions that fix the allowed values of $m$ and we plot the corresponding weights. We observe that in all cases $P_{N}^{z}(m,\ell)$ exhibits a single maximum centred at $m=0$ and is significantly narrower than its transverse counterpart $P_{N}^{x}(m,\ell)$. There is again an even/odd effect in $m$ that increases in magnitude as $\Delta$ approaches $1$, but it is generally weaker than in the transverse case. There also is an even/odd effect in the subsystem size $\ell$.

IV Generating functions for the subsystem magnetization

We now turn to the probability distribution of the (smooth) subsystem magnetization. We first consider the longitudinal generating function $G^{z}_{\ell}(\theta)$, as analytic results are readily available for it.

IV.1 Longitudinal generating function

$G^{z}_{\ell}(\theta)$

For large subsystem sizes the longitudinal generating function can be determined by standard methods: at $\Delta=0$ free fermion techniques apply, while for general values of $\Delta$ Luttinger liquid methods provide detailed predictions.

IV.1.1 The XX point $\Delta=0$

As it is straightforward to take into account a magnetic field along the z-direction in this case, we present results for the subsystem magnetization in the ground state of the Hamiltonian

[TABLE]

Here a simple determinant formula for $G^{z}_{\ell}(\theta)$ is known Colomo93

[TABLE]

The Toeplitz determinant (50) is related to a determinant that has been analyzed in great detail in the context of entanglement entropiesce-10 . Using the results of Ref. ce-10, the large-$\ell$ asymptotics of $G^{z}_{\ell}(\theta)$ can be expressed in the form

[TABLE]

where

[TABLE]

Here ${G(z)}$ is the Barnes $G$-function and

[TABLE]

The leading terms in (51) correspond to $j=0,1$ and have been considered previously in Ref. bss-07, . The constant $c_{1}(\beta)$ has been conjectured in Ref. aiq-11, . At zero magnetic field we have $k_{F}=\pi/2$ and

[TABLE]

We note that at zero magnetic field the generating function is real. This is because all odd cumulants vanish as a consequence of particle-hole symmetry. This ceases to be the case at finite magnetic fields, but odd cumulants still vanish in the large $\ell$ limit (as for the gas CMV-12 ).

IV.1.2 Luttinger liquid description

The large-$\ell$ behaviour of the expression (21) for $G^{z}_{\ell}(\theta)$ has been previously determined by bosonization methods in Refs bss-07, ; aem-07, . In zero magnetic field this gives a power-law decay

[TABLE]

where $\nu(\theta)=\frac{\theta^{2}}{4\pi^{2}\eta}$ and $\eta=1-\frac{1}{\pi}{\rm arccos}(\Delta)$. An analytic expression for the amplitudes of the leading terms in (56) was conjectured in Ref. bss-07,

[TABLE]

and $D_{-1}(\theta)=D_{0}(|\theta|-2\pi)$.

IV.1.3 Comparison to iTEBD results

In Fig. 11 we show the power law decay of $G^{z}_{\ell}(\theta)$ with subsystem size $\ell$ for several fixed values of $\theta$. In all cases the agreement with the Luttinger liquid prediction (56) is seen to be excellent.

A comparison with the exact results at $\Delta=0$ provides a useful accuracy check for our iTEBD data. As expected, the discrepancy grows with increasing subsystem size $\ell$. Moreover, it also depends on $\theta$ and grows as $\theta$ approaches $\pm\pi$. However up to subsystem sizes of $\ell=200$ the relative error of our iTEBD result is less than $0.1\%$.

The probability distribution $P_{S}^{z}(m,\ell)$ is then readily obtained by Fourier transforming $G^{z}_{\ell}(\theta)$. Plotting again the weights of the delta-functions that fix the possible values of $m$ gives the results shown in Fig. 12.

The probability distribution is centred around $m=0$ and is very narrow for all anisotropies $\Delta$. Moreover, there is very little subsystem size dependence for the large values of $\ell$ considered.

IV.2 Transverse generating function $G^{x}_{\ell}(\theta)$

The generating function $G^{x}_{\ell}(\theta)$ cannot be easily analyzed by either free fermion or bosonization methods. As we will see in section IV.2.1 it is however possible to determine it by field theory methods in particular limits. Our numerical results indicate that $G^{x}_{\ell}(\theta)$ decays exponentially in the subsystem size $\ell$ for all values of $\theta$ except $\theta=\pm\pi$, where it displays a power-law decay for even subsystem sizes (and vanishes for odd $\ell$)

[TABLE]

In order to analyze our numerical data for other values of $\theta$ we have carried out fits to the following functional form

[TABLE]

The resulting fits to our numerical data are shown in Figure 13 for some representative values of $\theta$.

The agreement between the fits and the data is seen to be very good. In Figure 14 we show the functions $\xi^{-1}(\theta,\Delta)$ and $\alpha(\theta,\Delta)$ resulting from our fits for several values of the anisotropy parameter $\Delta$. We checked the stability of the results against the “fit window” $[\ell_{min},\ell_{max}]$ of subsystem sizes used by suitably varying $\ell_{min}$ and $\ell_{max}$ in the interval $[50,200]$.

The inverse decay length $\xi^{-1}(\theta,\Delta)$ is seen to take its maximum around $\theta=\frac{\pi}{2}$ and is generally quite small. As the anisotropy $\Delta$ approaches unity $\xi^{-1}(\theta,\Delta)$ is seen to approach zero. This is expected because for $\Delta=1$ the spin rotational symmetry imposes $G^{x}_{\ell}(\theta)=G^{z}_{\ell}(\theta)$ and as we have seen the latter decays as a power law in $\ell$. Similarly, for the correlation length diverges for $\theta\rightarrow\pm\pi$, which indicates power-law behaviour in $\ell$ for these values of $\theta$. The exponent $\alpha(\theta,\Delta)$ of the power-law factor in (59) appears to be a monotonically increasing function of $\theta\in[0,\pi]$. For small $\theta$ it behaves as $\alpha(\theta\ll 1,\Delta)\sim\theta^{2}$.

Finally we present results for the probability distribution $P_{S}^{x}(m,\ell)$ of the smooth, transverse subsystem magnetization in Fig. 15.

We see that the probability distribution has a single maximum at $m=0$ and is generally quite narrow (albeit not as narrow as $P^{z}_{S}(m,\ell)$). Its width increases with diminishing $\Delta$ and is largest when $\Delta\to-1$.

IV.2.1 Perturbed Luttinger liquid approach to the transverse

generating function

In order to make progress, we consider the vicinity of the XXX point $\Delta=1$ and choose our anisotropy axis to lie along the $x$ direction, i.e.

[TABLE]

We now bosonize at the XXX point and then take the anisotropy into account as a perturbation. In the low energy limit the Hamiltonian can be written in the form

[TABLE]

where $\tilde{g}_{a}\propto(\Delta-1)$ and $J^{a}$ and $\bar{J}^{a}$ are the right and left chiral currents respectively. Defining chiral fields by $\Phi(x,t)=\phi_{R}(vt-x)+\phi_{L}(vt+x)=\bar{\varphi}(\bar{z})+\varphi(z)$ where $z=v\tau+ix$, we have the following expressions for the currents

[TABLE]

By virtue of the global $U(1)$ symmetry of $H(\Delta)$ we have only two independent coupling constants and $g_{3}=g_{2}$. The couplings fulfill the Kosterlitz-Thouless RG-equations GNT

[TABLE]

where $L$ and $b$ are hard long and short-distance cutoffs. To proceed it is convenient to define new couplings by

[TABLE]

The combination $g_{1}^{2}-g_{2}^{2}=\mu^{2}$ is an RG invariant and can be obtained by matching to Bethe-Ansatz calculationsgia

[TABLE]

The low-energy projection of the spin operator in the $z$-direction becomes

[TABLE]

where $a_{0}$ is the lattice spacing and the amplitude $a_{1}$ is known exactlylukyanovterras . If we now make the assumption that the staggered piece of the spin operator can be neglected (as was the case for the longitudinal generating function considered above), we have

[TABLE]

We will now determine (67) for large subsystem sizes $r=la_{0}$ by means of RG-improved perturbation theorycollins ; affleck98 ; BarzykinAffleck in the anisotropic current-current interactions (61). The Euclidean action corresponding to (61) is $S=\frac{1}{2}\int\mathrm{d}^{2}z(\partial_{\mu}\Phi)^{2}+S_{\text{int}}$ with

[TABLE]

To second order in perturbation theory in $S_{\rm int}$ we have

[TABLE]

where $T$ is the imaginary time ordering operator and $\left<.\right>_{0}$ is the path integral average with respect to the gaussian action. The perturbative expansion of the transverse generating function $G^{x}_{r}(\theta)$ thus reads

[TABLE]

where ${\cal T}_{1}$ and ${\cal T}_{2,a}$ denote the contributions at first and second order. The first order contribution can be evaluated following Ref. LudwigCardy, , which gives a logarithmic divergence in the short-distance cutoff $b$

[TABLE]

Due to electro-neutrality the various interaction terms in (68) do not mix in second order perturbation theory. The contribution proportional to $g^{2}_{2}$ can again be evaluated following Ref. LudwigCardy, with the result

[TABLE]

The contribution ${\cal T}_{2,2}$ proportional to $(g_{1}+g_{2})^{2}$ is of the form

[TABLE]

The logarithmically divergent parts of (73) can be determined by adapting the results of Refs LeDoussal, ; DL, to our real-space cutoff regularization scheme. This gives

[TABLE]

where $\gamma_{E}$ is the Euler-Mascheroni constant and $\psi(x)$ is the Digamma function. To lighten notations in what follows we define

[TABLE]

The third and final contribution in second order of perturbation theory is proportional to $(g_{1}-g_{2})^{2}$

[TABLE]

While the leading $\log^{2}(r/b)$ contribution can be easily extracted analytically, we resorted to a numerical integration for determining the subleading $\log(r/b)$ term

[TABLE]

Our numerical results for ${\cal T}_{2,3}$ are well-described by the functional form (77) as can be see in Fig. 16.

We now want to RG-improve the perturbative result (70) for $G^{x}_{r}(\theta)$ by solving the Callan-Symanzik equation

[TABLE]

The anomalous dimension is extracted from the perturbative expansion and is given by

[TABLE]

Solving the Callan-Symanzik equation by the method of characteristics we obtain

[TABLE]

The running couplings at scale $\log(r/b)$ are obtained by integrating (63)

[TABLE]

We now expand (80) in powers of the coupling constant and match the result to the perturbative expression (70). This provides us with an expansion of the function $F$ in (80). Putting everything together we arrive at the following expression for the RG improved correlator

[TABLE]

The overall amplitude $\mathcal{A}(\theta)$ is obtained by fitting to iTEBD results for $G^{x}_{\ell}(\theta)$. This leaves us with one free parameter, namely the initial coupling $g_{1}(0)$. We fix this by fitting (82) to iTEBD results for one value of $\theta$. In Fig. 17 we compare (82) obtained in this way to numerical results obtained by iTEBD for $\Delta=0.95$.

The agreement is seen to be quite satisfactory. For larger subsystem sizes we expect $G^{x}_{\ell}(\theta)$ to exhibit exponential decay in the subsystem size $\ell$. This clearly goes beyond RG-improved perturbation theory. It is an interesting problem how to obtain the corresponding correlation length in the framework of the perturbed Luttinger liquid (61).

V Full counting statistics for the Heisenberg ferromagnet

($\Delta=-1$) in the zero magnetization sector $S^{z}=0$

At $\Delta=-1$ the model (1) undergoes a first order quantum phase transition, where the total magnetisation of the ground state discontinuously jumps from $S_{z}=0$, for $\Delta>-1$, to $S_{z}=\pm L$, for $\Delta<-1$. However, if we restrict ourselves to the sector of the Hilbert space with zero total magnetisation, the ground state for $\Delta=-1$ is continuously connected the ground state for $\Delta+1=0^{+}$. The ground state in the $S^{z}=0$ sector at $\Delta=-1$ on a lattice with $L$ sites is

[TABLE]

where $S^{-}=\sum_{j=1}^{L}S^{x}_{j}-iS^{y}_{j}$. Here ${\cal N}$ is a normalization factor and the unitary transformation with $U$ maps the Hamiltonian (1) at $\Delta=-1$ to the isotropic ferromagnet. We note that the state [S^{-}]^{\frac{L}{2}}|\uparrow_{1}\dots\uparrow_{L}\rangle\ admits an exact MPS representation with auxiliary dimension $L+1$ in terms of the vector-valued matrix $\mathbf{\Gamma}_{\alpha\,\beta}=\delta_{\alpha\,\beta}|\!\uparrow\rangle+\delta_{\alpha+1\,\beta}|\!\downarrow\rangle$, and boundary vectors $v^{l}_{\alpha}=\delta_{\alpha\,L+1}$, $v^{r}_{\alpha}=\delta_{\alpha\,1}$.

V.1 Full counting statistics

The ground state (83) has the following useful representation

[TABLE]

which makes it possible to obtain closed-form expressions for the generating functions on the (staggered) subsystem magnetization. A straightforward combinatorial analysis gives the following results for the generating functions for a finite chain of $L$ sites

[TABLE]

In the limit $L\to\infty$ (86) simplify to

[TABLE]

were $P_{\ell}(z)$ are Legendre polynomials. In the vicinity of $\theta=0$ the generating functions (86) exhibit scaling for large subsystem sizes $\ell$

[TABLE]

where $I_{n}(z)$ and $J_{n}(z)$ are (modified) Bessel functions. We note that the universal scaling function in eqn (44) reduces to (87) in the limit $\Delta\to-1$.

V.1.1 Probability distribution functions

From (87) we can extract the following analytic expressions for the probability distribution functions at $\Delta=-1$

[TABLE]

where $K_{n}(z)$ are modified Bessel functions. In Fig. 18 we compare (88) to iTEBD results for finite $\ell$. The agreement is clearly excellent.

VI Summary and Conclusions

We have carried out a detailed study of the probability distributions of the components of the smooth and staggered subsystem magnetizations in the ground state of the critical spin-1/2 Heisenberg XXZ chain. We have shown through a combination of field theory and numerical calculations that appropriate ratios of the moments of these probability distributions are universal. The probability distributions of the longitudinal staggered subsystem magnetization is essentially Gaussian. This is in contrast to the transverse component $N^{x}(\ell)$, which can be thought of as the order parameter of the magnetic quasi long-range order in the XXZ chain. The corresponding probability distribution in the ground state $P^{x}(m,\ell)$ is very broad and for attractive and weakly repulsive interactions it exhibits two slight maxima at fairly large values of $N^{x}(\ell)$. We have shown that the dominant features for large subsystem sizes can be understood in terms of a mapping to the boundary sine-Gordon field theory. The subleading contributions also exhibit scaling, and their calculation in a field theory framework is an interesting open problem. The behaviour of the smooth subsystem magnetization $S^{\alpha}(\ell)$ is rather different from the staggered one. The generating function for the moments of the longitudinal component decays as a power law in subsystem size and can be accurately determined using Luttinger liquid theory. The corresponding probability distribution $P_{S}^{z}(m,\ell)$ is extremely narrow and centred around zero. This is perhaps not surprising as $S^{z}(L)$ is conserved as a result of the U(1) symmetry of the XXZ Hamiltonian. The probability distribution for the transverse component $S^{x}(\ell)$ is narrow and exhibits a single maximum at zero as well. The corresponding generating function of the moments $G^{x}_{\ell}(\theta)$ decays exponentially in $\ell$. We have shown that close to the antiferromagnetic point $\Delta=1$ its behaviour for intermediate values of $\ell$ can be determined by renormalization group improved perturbation theory. This calculation does not account for the exponential decay seen for large values of $\ell$. The description of the large-$\ell$ regime by field theory methods remains an interesting open problem.

Acknowledgements.

We are grateful to Pasquale Calabrese for collaboration in the early stages of this project and for numerous important discussions. We thank Paul Fendley and Austen Lamacraft for enlightening conversations. This work was supported by the EU Horizon 2020 research and innovation programme under Marie Sklodowska-Curie Grant Agreement No. 701221 (MC), the EPSRC under grant EP/N01930X/1 (FHLE) and by the Clarendon Scholarship fund (SG).

Appendix A Variance of the subsystem magnetization

The two-point functions in the XXZ chain have been determined in the framework of perturbed Luttinger liquid theory by Lukyanov and Terraslukyanovterras . After inversion of the spin quantization axes on all odd sites their result is given by (2), where the amplitudes are

[TABLE]

The variance $s_{x}$ (13) is obtained from the two-point function by

[TABLE]

The key identities for calculating $s_{x}$ are

[TABLE]

where $\zeta(x)$ is the Riemann zeta function. These identities show that all terms in the expansion of the two-point functions in fact contribute to $s_{x}$, irrespective of how fast their power law decays are. However, for large $\gamma$ the dominant contribution comes from the $n=1$ terms in the corresponding sums, i.e. from the non-universal short-distance behaviour. This shows that it is useful to take the short-distance behaviour of the correlators into account as precisely as possible. For example, we can decompose $s_{x}$ as

[TABLE]

Using the Lukyanov-Terras result (2) we the obtain the following approximate expression

[TABLE]

The nearest neighbour correlator can be simply obtained from the ground state energy per site and equals

[TABLE]

where we have defined $\nu=\frac{1}{\pi}{\rm arccos}(-\Delta)$. We note that the exact next-nearest-neighbour spin-spin correlators are also available in the literature and can be taken into account in the same way. The amplitude $B$ has an unphysical singularity as $\eta\rightarrow 2/3$. This merely means that the perturbative calculation of Ref. lukyanov, needs to be redone for $\eta\approx 2/3$. As the contribution of the $B$ term to $s_{x}$ becomes important only as we approach the XXX point, we simply drop it in the following. We can compare the field theory results to direct DMRG computations of the variance. In Table 2 we show the results for several values of $\Delta$. We see that in the attractive regime we have good agreement, while in the repulsive regime the agreement is worse.