Full counting statistics in the spin1/2 Heisenberg XXZ chain
Abstract
The spin1/2 Heisenberg chain exhibits a quantum critical regime characterized by quasi longrange 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 sineGordon model.
pacs:
64.70.TgI 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 manyparticle 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 HaldaneShastry model(15). Here we consider the (staggered) subsystem magnetization in the anisotropic onedimensional spin Heisenberg XXZ chain
(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 lowenergy behaviour of the model (1) is described by Luttinger liquid theory or equivalently a free, compact boson(16); (17); (18); (19). The longdistance asymptotics of spinspin correlation functions is of the form
(2) 
where explicit expressions for the amplitudes in (2) are known (20); (45); (21); (23) and is related to the anisotropy parameter by
(3) 
It follows from (2) that throughout the critical regime the dominant correlations are those of the staggered magnetizations in the xyplane. The XXZ chain thus exhibits antiferomagnetic quasilong range order in the XY plane in spin space. Twopoint 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
(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
(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
(6) 
It is easy to see that they have the following properties
(7) 
The last relation allows us to restrict our attention to the interval and can be obtained e.g. from the representation
(8) 
Defining
(9) 
the probability distributions of interest can be expressed as
(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
(11) 
The variances have the following asymptotic expansions for large subsystem sizes
(12) 
For sufficiently large values of we expect the coefficients and to be equal to the corresponding quantities for the entire system, i.e.
(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 nonuniversal quantities they are expected to be susceptible to shortdistance 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
(14) 
If these ratios are universal, the modified generating functions
(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
(16) 
where and are Bose fields with commutation relations , the dots indicate perturbations that are irrelevant in the renormalization group sense and
(17) 
The bosonization formulas for the spin operators are
(18)  
(19) 
where is the lattice spacing and the amplitudes , , are known exactly (23). For large subsystem sizes we thus have
(20) 
Applying the bosonization prescription to our generating functions and ignoring subleading terms we obtain
(21) 
where is the Fock vacuum. The representation (21) reveals that maps onto a simple vertex operator twopoint function in the free boson theory, whereas correspond to expectation values of nonlocal 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 noninteracting spinless fermions by means of a JordanWigner transformation. Using standard techniques(28) we can derive the following determinant representation for the longitudinal generating function
(22) 
Here is the correlation matrix of the free fermion chain obtained by the Jordan Wigner transformation
(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)
(24) 
where
(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
(26) 
, , while keeping fixed.
In this regime the behaviour depends on the parity of the subsystem size
(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
(28)  
where
(29) 
As the leading singularities of the integrand occur when we now make the further approximation
(30) 
This results in
(31) 
The right hand side of (31) is equal to the partition function of a boundary sineGordon 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
(32) 
The function can be computed very efficiently from the solution of the singleparticle Schrödinger equation
(33) 
Denoting by and the solutions to (33) with asymptotics
(34) 
we have
(35) 
where denotes the Wronskian.
Let us now turn to the longitudinal generating function . It has an integral representation
(36) 
This expression needs to be regularized because
(37) 
and ranges from at the isotropic point to infinity when the ferromagnetic point is approached (). The righthand side of (36) can again be related to the partition function of a boundary sineGordon model, but the boundary interaction for is now irrelevant. This suggests that will be independent of and equal to the result at , i.e.
(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 SuzukiTrotter decomposition with imaginary timestep . 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 groundstate energy densities for different values of and compare them to the known exact result
(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
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)
(40) 
Here is related to the function in (35) by
(41) 
where is a nonuniversal 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
(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 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

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.

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.

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

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

The distribution for attractive and moderately repulsive interactions is bimodal, 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
(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 sineGordon 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 xyplane. 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 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
(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
(45) 
We conclude that the fluctuations of have a very simple form: the second moment is
(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
(47) 
Here , , , and are dependent parameters that we fix by considering the dependencies of (for even) and (for odd). Our numerical results suggest that
(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 deltafunctions 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 zdirection in this case, we present results for the subsystem magnetization in the ground state of the Hamiltonian
(49) 
Here a simple determinant formula for is known (39)
(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
(51) 
where
(52) 
Here is the Barnes function and
(53)  
(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
(55) 
We note that at zero magnetic field the generating function is real. This is because all odd cumulants vanish as a consequence of particlehole 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 powerlaw decay
(56) 
where and . An analytic expression for the amplitudes of the leading terms in (56) was conjectured in Ref. (40)
(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 deltafunctions 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.