A Variance of the subsystem magnetization

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

## 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.

64.70.Tg

## I Introduction

Universality is a key organizing principle for continuous phase transitions(1); (2). 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. (3), 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 matter(4); (5) and cold atom physics(6); (7); (8); (9). Theoretical results on full counting statistics in quantum critical systems are relatively scarce. The list of available results includes phase fluctuations in Luttinger liquids(10); (11); (12); (13), the order parameter statistics in the Ising field theory(3), the transverse magnetization in the Ising chain(14) and the magnetization in the Haldane-Shastry model(15). Here we consider the (staggered) subsystem magnetization in the anisotropic one-dimensional spin- Heisenberg XXZ chain

 H=JL∑j=1SxjSxj+1+SyjSyj+1+ΔSzjSzj+1 . (1)

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 . The special values correspond to the isotropic antiferromagnet and ferromagnet respectively. In the regime the low-energy behaviour of the model (1) is described by Luttinger liquid theory or equivalently a free, compact boson(16); (17); (18); (19). The long-distance asymptotics of spin-spin correlation functions is of the form

 ⟨GS|Sxj+nSxj|GS⟩ = Missing or unrecognized delimiter for \right ⟨GS|Szj+nSzj|GS⟩ = −14π2ηn2(1+~Bzn4/η−44−3η2−2η)+(−1)nAz4n1/η(1−Bzn2/η−2)+… (2)

where explicit expressions for the amplitudes in (2) are known (20); (45); (21); (23) and is related to the anisotropy parameter by

 Δ=−cos(πη). (3)

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 systems(19). 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

 Sα(ℓ)=ℓ∑j=1Sαj ,Nα(ℓ)=ℓ∑j=1(−1)jSαj . (4)

The quantities and describe the smooth and staggered components of the -component of the magnetization of the subsystem consisting of sites to , where . We note that whereas is a conserved quantity, is not. The probabilities of the observables (4) taking some value when the system is prepared in the ground state and a measurement is then performed are

 PαS(m,ℓ) = ⟨GS|δ(Sα(ℓ)−m)|GS⟩=∫∞−∞dθ2π e−iθm ⟨GS|eiθSα(ℓ)|GS⟩ , PαN(m,ℓ) = ⟨GS|δ(Nα(ℓ)−m)|GS⟩=∫∞−∞dθ2π e−iθm ⟨GS|eiθNα(ℓ)|GS⟩ . (5)

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

 Gαℓ(θ)≡⟨GS|eiθSα(ℓ)|GS⟩ ,Fαℓ(θ)≡⟨GS|eiθNα(ℓ)|GS⟩ . (6)

It is easy to see that they have the following properties

 Xαℓ(0) = 1 ,Xαℓ(−θ)=(Xαℓ(θ))∗ ,Xαℓ(θ+2π)=(−1)ℓXαℓ(θ) ,X=F,G. (7)

The last relation allows us to restrict our attention to the interval and can be obtained e.g. from the representation

 eiθSα(ℓ)=ℓ∏j=1[cos(θ/2)+isin(θ/2)σαj]. (8)

Defining

 ˜Xαℓ(r)=∫π−πdθ2π e−irθ Xαℓ(θ) ,X=F,G, (9)

the probability distributions of interest can be expressed as

 PαN(m,ℓ)={∑r∈Z˜Fαℓ(r) δ(m−r)if\ ℓ is even,∑r∈Z˜Fαℓ(r+12) δ(m−r−12)if\ ℓ is odd. (10)

An analogous equation holds for .

### 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 and vanish

 ⟨GS|Sα(ℓ)|GS⟩=0=⟨GS|Nα(ℓ)|GS⟩. (11)

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

 ⟨GS|(Sα(ℓ))2|GS⟩=ℓsα+o(ℓ) ,⟨GS|(Nα(ℓ))2|GS⟩=ℓnα+o(ℓ) . (12)

For sufficiently large values of we expect the coefficients and to be equal to the corresponding quantities for the entire system, i.e.

 Missing or unrecognized delimiter for \big (13)

As is a conserved quantity and our system is translationally invariant we have . 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 universal(24); (25); (3). In particular one may expect the following ratios to exhibit universal behaviour

 Missing or unrecognized delimiter for \big (14)

If these ratios are universal, the modified generating functions

 ⟨GS|eiθYYα(ℓ)|GS⟩ ,θY=θ√⟨GS|(Yα(ℓ))2|GS⟩ ,Y=S,N, (15)

will be universal functions of the parameter . 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 .

## 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 theory(26); (20); (45); (23); (27). In absence of a magnetic field the Hamiltonian can be cast in the form

 H(Δ)=v2∫dx[K(∂xθ)2+1K(∂xϕ)2]+… , (16)

where and are Bose fields with commutation relations , the dots indicate perturbations that are irrelevant in the renormalization group sense and

 v=π2√1−Δ2arccosΔ,K=π21π−arccosΔ, (17)

The bosonization formulas for the spin operators are

 Szj ≃ −a0√π∂xϕ(x)+(−1)jc1sin(√4πϕ(x))+…, (18) Sxj ≃ b0(−1)jcos(√πθ(x))+ib1sin(√πθ(x))sin(√4πϕ(x))+…, (19)

where is the lattice spacing and the amplitudes , , are known exactly (23). For large subsystem sizes we thus have

 Sz(ℓ) ≈ −1√π[ϕ(ℓa0)−ϕ(0)]+… , Nz(ℓ) ≈ c1a0∫ℓ0dx sin(√4πϕ(x))+… , Nx(ℓ) ≈ b0a0∫ℓ0dx cos(√πθ(x))+… . (20)

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

 Gzℓ(θ) ≈ ⟨0|e−iθ√πϕ(ℓa0)eiθ√πϕ(0)|0⟩ , Fzℓ(θ) ≈ ⟨0|e−iθc1a0∫ℓ0dx sin(√4πϕ(x))|0⟩ , Fxℓ(θ) ≈ ⟨0|e−iθb0a0∫ℓ0dx cos(√πϕ(x))|0⟩ , (21)

where is the Fock vacuum. The representation (21) reveals that maps onto a simple vertex operator two-point function in the free boson theory, whereas correspond to expectation values of non-local operators. The alert reader will have noted that we did not provide a bosonized expression for . The reason is that the field theory calculation of 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 and the corresponding generating functions . We first present analytic results in certain limits and then compare these to numerical ones.

### iii.1 The XX point Δ=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 techniques(28) we can derive the following determinant representation for the longitudinal generating function

 Fzℓ(θ) = det[B] ,Bnm=eiθ2(−1)nδn,m+2sin(θ/2)(−i)(−1)n+(−1)mCn,m . (22)

Here is the correlation matrix of the free fermion chain obtained by the Jordan Wigner transformation

 Cnm=sin(π2(m−n))π(m−n) . (23)

The matrix is a Toeplitz matrix(29). Its properties have been analyzed in great detail in the context of entanglement entropies in Ref. (30). For large values of the subsystem size one obtains the following asymptotic expansion(30)

 Fzℓ(θ)=(cos(θ))ℓ2∞∑m=−∞(−1)mℓ(2ℓ)−2(m+β(θ))2G2(m+1+b(θ))G2(1−m−b(θ))[1+c2(m+b(θ))ℓ2+…] (24)

where

 b(θ) = 12πiln[cos(θ/2)−sin(θ/2)cos(θ/2)+sin(θ/2)] ,c2(β)=−β2(1+8β2)6 . (25)

It follows from (24) that is very small except in the vicinities of . To analyze the behaviour in these regions it is useful to define the following scaling limits:

• , , while keeping fixed.

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

 Fzℓ(θ)∼e−z2/4 . (26)
• , , while keeping fixed.

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

 Fzℓ(θ)∼√2G2(1/2)G2(3/2)e−y2/4×⎧⎪⎨⎪⎩(−1)ℓ/2ℓ−1/2,ℓeven(−1)(ℓ−1)/2ℓ−1[c+log(2ℓ)]y/π,ℓodd (27)

where and .

We will see in the following that the two limits S1 and S2 are useful for analyzing numerical results for .

### iii.2 Field theory approach

In the field theory approach we are tasked with evaluating the expressions (21) for . This can be done by following the analysis of Refs (10); (11); (12), which considered generating function for phase fluctuations in Luttinger liquids. Expanding in powers of we obtain

 Fxℓ(θ) ≈ ∞∑n=0(iθb02a0)2n1(2n)!∫ℓ0dx1…∫ℓ0dx2n ⟨2cos(√πθ(x1))…2cos(√πθ(x2n))⟩ (28) = ∞∑n=01(n)!(iθℓb04πa0)2n Z(0)2n(L,ℓ/L,K) ,

where

 Z(0)2n(L,z,K) = ∫2π0du12π…dun2π∫2π0dv12π…dvn2π∏nj,k=1G(uj−vk)∏nj

As the leading singularities of the integrand occur when we now make the further approximation

 Z(0)2n(L,ℓ/L,K)≈Z(0)2n(ℓ,1,K) . (30)

This results in

 Fxℓ(θ)≈∞∑n=01(n!)2(iθℓb02a0)2nZ(0)2n(ℓ,1,K). (31)

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 literature(31); (32); (33); (34); (35); (36). Using a result of Ref. (36) for the right hand side of (31) one has

 Missing or unrecognized delimiter for \big (32)

The function can be computed very efficiently from the solution of the single-particle Schrödinger equation

 −∂2xΨ(x)+[x4K−2−14x2]Ψ(x)=EΨ(x) . (33)

Denoting by and the solutions to (33) with asymptotics

 Ψ+(x,E) ∼ √πx2Kfor x→0 , χ+(x,E) ∼ x−2K−12e−x4K/4Kfor x→∞ , (34)

we have

 A(vac)(λ)=12W[χ+,Ψ+]∣∣∣E=ρλ2 ,ρ=(8K)2−1/2KΓ2(1−1/4K), (35)

where denotes the Wronskian.

Let us now turn to the longitudinal generating function . It has an integral representation

 Fzℓ(θ) ≈ ∞∑n=0(iθc12a0)2n1(2n)!∫ℓ0dx1…∫ℓ0dx2n ⟨2sin(√4πϕ(x1))…2sin(√4πϕ(x2n))⟩ . (36)

This expression needs to be regularized because

 ⟨ei√4πϕ(x) e−i√4πϕ(x)⟩=(Lπa0sin(πxL))−2K , (37)

and ranges from at the isotropic point to infinity when the ferromagnetic point is approached (). 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 is now irrelevant. This suggests that will be independent of and equal to the result at , i.e.

 Fzℓ(θ)≈e−z2/4 ,z=θℓ1/2. (38)

Eqns (32), (35) and (38) provide us with explicit expressions for the generating functions that now can be compared to numerical results for the lattice model.

### iii.3 Numerical Method

Our numerical approach is based on the iTEBD algorithm(37); (38). 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 , which admits an MPS representation with auxiliary dimension . We then evolve the state in imaginary time by the operator by means of a second order Suzuki-Trotter decomposition with imaginary time-step . 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 . In Table 1 we present our best estimates for the ground-state energy densities for different values of and compare them to the known exact result

 EGS=14∫∞−∞dλ12πcosh(λ1−η)4(Δ2−1)cosh(2λ)−Δ+Δ4. (39)

Our numerical results for the energy densities differ from the exact values by , 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 and with a computational cost that scales as . A useful check on the numerical accuracy of our results can be obtained by considering the noninteracting case , 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 . However up to subsystem sizes of the relative error of our iTEBD result is less than .

### iii.4 Numerical results for the transverse generating function Fxℓ(θ)

Numerical results for as a function of for several values of the subsystem size are shown in Fig. 1.

We see that the generating function is very small everywhere except in the vicinities of . We also observe that the oscillatory behaviour as a function of becomes more pronounced in the attractive regime .

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

 Fxℓ(θ≈0)=Fx0(z) ,z=θℓ1−η/2 . (40)

Here is related to the function in (35) by

 Fx0(z)=A(vac)(cz) , (41)

where is a non-universal -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 by carrying out a best fit of our numerical data to (41). In Fig. 2 we present a comparison of our numerical results for 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 we have considered in the critical regime . We again see that in the attractive regime the oscillatory behaviour away from becomes more pronounced.

We now turn to the other region in which is sizeable, namely . 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 which requires us to consider even and odd separately.

Our numerical data in the vicinity of is well described by the scaling ansatz

 Fxℓ(θ≈π)≃(−1)⌊ℓ/2⌋ℓ−1/4Fxe/o(z) ,z=(θ−π)ℓ1−η/2 , (42)

where refers to even and odd subsystem size respectively. Inspection of Figs 3 and 4 shows that the ansatz is in excellent agreement with the data. We note that at the numerical data (for even) exhibit a perfect algebraic decay , independent of the value of the interaction . The form (42) suggests that for very large subsystem sizes in the thermodynamic limit the feature at becomes less and less important compared to . At present no analytic results on are known. It should in principle be possible to calculate using field theory methods.

### iii.5 Probability distribution PxN(m,ℓ) of the transverse, staggered subsystem magnetization

We are now in a position to determine the probability distribution from the generating function using eqn (10). In Fig. 5 we show results for as a function of for several values of the exchange anisotropy and subsystem sizes . As is a sum over -functions, cf. eqn (10), we plot the corresponding weights at the appropriate values of .

We observe that

1. There is a strong even/odd effect in . The results for even and odd follow different smooth curves. The separations between even and odd curves slowly tend to zero as as the subsystem size is increased.

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

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

4. The width of increases as the interaction becomes more attractive and the distribution flattens.

5. The distribution for attractive and moderately repulsive interactions is bi-modal, while for it displays a single maximum.

These observations can be understood in terms of our scaling analysis of the generating function . The probability distribution is dominated by the behaviour of in the regions , and exploiting the observed scaling behaviour of these contributions we conclude that

 PxN(m,ℓ)≃ℓη2−1˜Fx0(m/ℓ1−η/2)+(−1)⌊ℓ/2⌋+⌊m⌋ℓη2−5/4˜Fxe/o(m/ℓ1−η/2), (43)

where and are obtained by Fourier transforming the functions and that describe the scaling behaviour of the generating function around and respectively. For large subsystem sizes the even/odd effect in disappears and we are left with , 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 as can be seen in Fig. 6, which presents results for at . We observe that for large subsystem sizes the probability distribution becomes fairly flat over most of allowed range of staggered magnetizations except for an enhancement close to the maximal possible values .

### iii.6 Longitudinal generating function Fzℓ(θ) and probability distribution PzN(m,ℓ)

We now turn to the longitudinal generating function . As we will see, its behaviour is rather different from . According to the field theory approach discussed in section III.2 we expect to be described by the scaling function

 Fzℓ(θ≈0)=e−γz2/4 ,z=θℓ1/2 , (44)

where is a -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 are shown in Fig. 7 and are seen to be in excellent agreement with the scaling form (44).

The coefficient is found to be consistent with

 γ=12−2η. (45)

We conclude that the fluctuations of have a very simple form: the second moment is

 ⟨GS|(Nz(ℓ))2|GS⟩=ℓ8−8η+o(ℓ) , (46)

while all higher cumulants vanish. We now turn to the behaviour of . Guided by the exact result (27) for we have attempted to describe our numerical data by the ansatz

 Fzℓ(θ≈π)={(−1)ℓ/2A ℓαe−γz2/4ℓ even(−1)(ℓ−1)/2Bc+log(2ℓ)ℓβ ze−γz2/4ℓ odd ,z=(π−θ)ℓ1/2. (47)

Here , , , and are -dependent parameters that we fix by considering the -dependencies of (for even) and (for odd). Our numerical results suggest that

 α=14η ,β=12η . (48)

In Figs 8 and 9 we compare numerical results for for several values of and to the scaling ansatz (47), (48).

The agreement is seen to be quite satisfactory in all cases. Having determined the generating function we can now use it to obtain the probability distribution of the longitudinal staggered subsystem magnetization by Fourier transform. Results for several values of and are presented in Fig. 10.

The probability distribution is again a sum over delta-functions that fix the allowed values of and we plot the corresponding weights. We observe that in all cases exhibits a single maximum centred at and is significantly narrower than its transverse counterpart . There is again an even/odd effect in that increases in magnitude as approaches , but it is generally weaker than in the transverse case. There also is an even/odd effect in the subsystem size .

## 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 , as analytic results are readily available for it.

### iv.1 Longitudinal generating function Gzℓ(θ)

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

#### The XX point Δ=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

 H=JL∑j=1(SxjSxj+1+SyjSyj+1)−hL∑j=1Szj ,0

Here a simple determinant formula for is known (39)

 Gzℓ(θ) = eiθℓ2det[I+(e−iθ−1)C] , Cnm = sin(kF(m−n))π(m−n) ,kF=arccos(h/J) . (50)

The Toeplitz determinant (50) is related to a determinant that has been analyzed in great detail in the context of entanglement entropies(30). Using the results of Ref. (30) the large- asymptotics of can be expressed in the form

 Gzℓ(θ)=eiℓθ/2∞∑j=−∞ρℓ(j+θ2π) (51)

where

 ρℓ(β)=e−2iβkFℓ(2ℓsin(kF))−2β2G2(1+β)G2(1−β)[1+c1(β)ℓ+c2(β)ℓ2+…]. (52)

Here is the Barnes -function and

 c1(β) = 2iβ3cot(kF), (53) c2(β) = β26(−1+7β2+12β4−3β2(5+4β2)csc2kF) . (54)

The leading terms in (51) correspond to and have been considered previously in Ref. (40). The constant has been conjectured in Ref. (41). At zero magnetic field we have and

 c1(β)=0 ,c2(β)=−16β2(1+8β2). (55)

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 limit (as for the gas (42)).

#### Luttinger liquid description

The large- behaviour of the expression (21) for has been previously determined by bosonization methods in Refs (40); (43). In zero magnetic field this gives a power-law decay

 Gzℓ(θ)≃∞∑j=−∞e−iπℓjDj(θ) ℓ−ν(θ+2πj) , (56)

where and . An analytic expression for the amplitudes of the leading terms in (56) was conjectured in Ref. (40)

 D0(θ) = ⎡⎢ ⎢⎣Γ(η2−2η)2√πΓ(12−2η)⎤⎥ ⎥⎦θ2/(4ηπ2)exp[−∫∞0dtt(sinh2θ2πtsinhtcosh(1−η)tsinhηt−θ2e−2t4ηπ2)], (57)

and .

#### Comparison to iTEBD results

In Fig. 11 we show the power law decay of with subsystem size for several fixed values of . In all cases the agreement with the Luttinger liquid prediction (56) is seen to be excellent.

A comparison with the exact results at provides a useful accuracy check for our iTEBD data. As expected, the discrepancy grows with increasing subsystem size . Moreover, it also depends on and grows as approaches . However up to subsystem sizes of the relative error of our iTEBD result is less than .

The probability distribution is then readily obtained by Fourier transforming . Plotting again the weights of the delta-functions that fix the possible values of gives the results shown in Fig. 12.

The probability distribution is centred around and is very narrow for all anisotropies . Moreover, there is very little subsystem size dependence for the large values of considered.