A Variance of the subsystem magnetization

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


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.


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


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


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


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


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


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


It is easy to see that they have the following properties


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




the probability distributions of interest can be expressed as


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


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


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


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


If these ratios are universal, the modified generating functions


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


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


The bosonization formulas for the spin operators are


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


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


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

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


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


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)




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

  • , , while keeping fixed.

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


    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




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


This results in


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


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


Denoting by and the solutions to (33) with asymptotics


we have


where denotes the Wronskian.

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


This expression needs to be regularized because


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.


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

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

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

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

Figure 1: Staggered transverse generating function , for representative values of and .

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)


Here is related to the function in (35) by


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.

Figure 2: Staggered transverse generating function for several values of . The numerical results (symbols) are seen to exhibit scaling collapse in the variable and are well described by the universal scaling function (32), (35) calculated from the boundary sine-Gordon model (red line).

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.

Figure 3: Scaling behaviour of the staggered transverse generating function for , even subsystem sizes and several values of exchange anisotropy .
Figure 4: Same as Fig. 3 for odd .

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


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

Figure 5: Probability distribution functions for , and . As is a sum over -functions we plot the corresponding weights at the appropriate values of . The solid red lines show the field theory result, which becomes exact in the large- limit.

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


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 .

Figure 6: Probability distribution functions for . As the system become more and more ”ferromagnetic” the probability distribution for the staggered subsystem magnetisation tends to become broader and flat.

iii.6 Longitudinal generating function and probability distribution

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


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

Figure 7: Staggered longitudinal generating function for several values of . The numerical data exhibit scaling collapse that is in excellent agreement with the universal scaling function (red line).

The coefficient is found to be consistent with


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


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


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


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

Figure 8: Scaling behaviour of the staggered longitudinal generating function for even subsystem sizes and several values of . For the point the full red line represents the analytical scaling function .
Figure 9: Scaling behaviour of the staggered longitudinal generating function for odd subsystem sizes and several values of . For the case, the red line represent the analytical scaling function in Eq. (27). The constant equals for and has been set to for and for respectively.

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.

Figure 10: Probability distribution functions for , and . In the noninteracting case the full lines are the exact results obtained using the determinant formula. In the interacting cases, the red dashed lines are the leading smooth contribution coming from the scaling behaviour of the generating function in the vicinity of . Notice in particular that, in the ferromagnetic regime, the sub-leading staggered corrections are almost invisible for the sizes considered here. Otherwise, as the antiferromagnetic regime is approached, the sub-leading parity effects become more significant.

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

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

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


Here a simple determinant formula for is known (39)


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




Here is the Barnes -function and


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


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


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


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.

Figure 11: as a function of subsystem size for different values of the anisotropy and . iTEBD results (symbols) are compared to the Luttinger liquid prediction (56) (full lines), apart from XX case () where the exact formula (50) has been used.

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.

Figure 12: Probability distribution functions for , and . In the noninteracting case the full lines are the exact results obtained integrating the determinant formula.

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.

iv.2 Transverse generating function