# Ground state parameters, finite-size scaling, and low-temperature properties of the two-dimensional 𝑆=1/2 XY model

## I Introduction

Studies of effective continuum field theories have resulted in detailed predictions for the low-energy physics of quantum spin systems in two dimensions. In order to make use of these predictions for a given model Hamiltonian, the ground state parameters appearing in the Lagrangian formulation have to be determined. Spin-wave theory  can in some cases give good estimates,  but in general some numerical method has to be employed in order to obtain accurate results. Since the theories can also predict the finite-size scaling behavior of various physical quantities, numerical results for a series of lattice sizes can be used to extract the ground state parameters. Such calculations are also important for testing theoretical predictions.

With todays computers, Lanczos and related exact diagonalization methods can be used for square lattices with up to $`6\times 6`$ spins.  This relatively small maximum size, and the small number of different lattices available, can make the finite-size scaling procedures problematic if sub-leading corrections are significant. It is therefore important to consider also alternative methods that can reach larger lattice sizes.

Here we discuss quantum Monte Carlo (QMC) results for the $`S=1/2`$ XY model, defined by the Hamiltonian

$$H=J\underset{i,j}{}[S_i^xS_j^x+S_i^yS_j^y],$$
(1)   
where $`S_i^x`$ and $`S_i^y`$ are the $`x`$\- and $`y`$-components of a spin-1/2 operator at site $`i`$, and $`i,j`$ denotes a pair of nearest-neighbor sites on a square lattice. Numerical studies of this model have a long history. Exact finite-lattice calculations gave the first indications that the O(2) symmetry is spontaneously broken at $`T=0`$ in the thermodynamic limit.  This was later proved rigorously. QMC simulations have resulted in a quite precise value of the ordered moment  and have also shown  that there is a Kosterliz-Thouless (KT) transition at a temperature $`T_{\mathrm{KT}}0.343J`$.  For $`TT_{\mathrm{KT}}`$ the system acquires quasi-long-range order, i.e., power law decay of the spin-spin correlation function and a non-zero spin stiffness constant. In accordance with the Mermin-Wagner theorem, true long-range order develops only at $`T=0`$.

Our motivation for carrying out calculations for the XY model to even higher precision is to test predictions of effective Lagrangian theories. In particular, the chiral perturbation calculations by Hasenfratz and Niedermayer have resulted in detailed predictions for the finite-size and finite-temperature corrections of various quantities, in some cases beyond leading order. The finite-size and finite-temperature scaling behavior of the O(3) symmetric Heisenberg model has been the topic of numerous studies, and the agreement with the predictions has been confirmed to high precision. The predictions for the O(2) symmetric XY model have not yet been tested exhaustively, however. A recent finite-size scaling study of exact energies for systems with up to $`6\times 6`$ spins has quite convincingly demonstrated agreement with the leading finite-size behavior. Here we consider also several other physical observables. This allows us to carry out a number of independent tests of the consistency of the scaling predictions. With access to larger system sizes, we can also improve considerably on the accuracy of the extrapolated ground state parameters and finite-size corrections. Our data are sufficiently accurate for addressing also sub-leading corrections. In addition to calculations in the ground state, we have also studied systems at finite temperature, on lattices sufficiently large to enable extraction of the leading finite-temperature corrections in the thermodynamic limit.

We have used a numerically exact finite-temperature QMC method based on “stochastic series expansion”, i.e., importance sampling of diagonal matrix elements of the power series expansion of $`\mathrm{exp}(\beta H)`$, where $`\beta =1/T`$ is the inverse temperature. The ground state can be obtained by choosing a sufficiently large $`\beta `$ value. This method has previously been used to study the finite-size scaling behavior and the ground state parameters of the square-lattice Heisenberg model. Recently, a new way of sampling the terms of the expansion was proposed (the “operator-loop” algorithm), which enables calculations to even higher precision. We have used this improved method to study the ground state of the XY model on lattices with $`N=L\times L`$ spins with $`L`$ up to $`16`$. We have carried out calculations at finite temperature for systems with $`L`$ up to $`64`$.

In Sec. II we discuss the physical quantities that we have calculated. The numerical ground state data and the results for finite-size corrections and extrapolations to infinite size are presented in Sec. III. In Sec. IV we discuss calculations at $`T>0`$. Sec. V concludes with a brief summary.

## II Calculated quantities

The QMC algorithm  can be implemented in two different ways; with the spin quantization either along the $`z`$ or $`x`$ axis. The internal energy, $`E=H/N`$, can be easily calculated in both cases. In the $`z`$-representation, the spin stiffness constant $`\rho `$ can also be evaluated. It is defined as the second derivative of the internal energy with respect to a twist $`\varphi `$ under which, for all nearest-neighbor spin pairs in either the $`x`$ or $`y`$ lattice direction, the spin-spin interaction is modified according to

$`S_i^xS_j^x`$ $`+S_i^yS_j^y(S_i^xS_j^x+S_i^yS_j^y)\mathrm{cos}(\varphi )`$ (3)   
$`S_i^xS_j^y\mathrm{sin}(\varphi )+S_i^yS_j^x\mathrm{sin}(\varphi ).`$   
In analogy with the superfluid density of a boson system, one can show that

$$\rho =\frac{^2E(\varphi )}{^2\varphi }=\frac{W_x^2+W_y^2}{2N\beta },$$
(4)   
where $`E(\varphi )`$ is the internal energy per spin and $`W_x`$ and $`W_y`$ are the “winding numbers”, i.e., the net spin currents across the periodic boundaries in the $`x`$ and $`y`$ directions that characterize configurations in simulations carried out in the $`z`$-representation. We also calculate the spin susceptibility, given by

$$\chi =\frac{\beta }{N^2}\left(\underset{i=1}{\overset{N}{}}S_i^z\right)^2.$$
(5)   
In the $`z`$-representation, it is not convenient to calculate spin-spin correlations involving the $`x`$ or $`y`$ spin components. We therefore also use an algorithm implemented in the $`x`$ representation, and there calculate the squared magnetization $`M_x^2`$;

$$M_x^2=\frac{1}{N^2}\left(\underset{i=1}{\overset{N}{}}S_i^x\right)^2.$$
(6)   
Since the O(2) symmetry is not broken on a finite lattice, the magnetization $`m`$ per spin in the thermodynamic limit is given by the infinite-size extrapolated value of

$$m=\sqrt{2M_x^2}.$$
(7)   
We use periodic boundary condition and include each spin pair only once in Eq. (1). For technical details on the simulation algorithm we refer to previously published work.

## III Ground state results

In order to obtain the ground state, we have carried out simulations at inverse temperatures as high as $`\beta =1024`$ for lattices with $`L=416`$. In the $`z`$-representation simulations, we have calculated the ground state energy to within relative statistical errors lower than $`10^5`$. The relative accuracy of the stiffness is on the order of $`10^3`$. Simulations in the $`x`$-representation result in statistical errors for the energy roughly twice as large as in the $`z`$-representation (for simulations of comparable length). The squared magnetization is the least accurate quantity, with statistical errors of roughly $`0.3\%`$ for the largest system size considered in this case ($`L=12`$). In Table I we list the results used in the finite-size analysis presented below. The results for the energy for $`L=4`$ and $`6`$ are in perfect agreement with previous  exact diagonalizations. The magnetization for $`L=4`$ also agrees with the exact diagonalization result. However, for $`L=6`$ there is an $`5\sigma `$ deviation from the result presented in Ref. . We do not know the reason for this disagreement. The QMC method has previously been shown to give perfect agreement with isotropic Heisenberg results for both $`L=4`$ and $`L=6`$, and therefore a failure for (only) the XY-model magnetization for $`L=6`$ would be surprising.

For a model with O(2) symmetry, chiral perturbation theory predicts the size dependence of the energy as 

$$E_0(L)=E_0+\frac{e_3}{L^3}+\frac{e_5}{L^5}+\mathrm{},$$
(8)   
with no $`O(1/L^4)`$ term. The leading-order correction is given by

$$e_3=\gamma c/2,$$
(9)   
where the constant $`\gamma =1.437745`$ and $`c`$ is the spin-wave velocity. In order to obtain a good fit to our $`L=416`$ data, the $`O(1/L^5)`$ term has to be included. The fit then has a chi-squared value per degree of freedom of $`0.7`$. Using instead an $`O(1/L^4)`$ term gives chi-squared $`6`$. We can therefore conclude that our results support the prediction  that $`e_4=0`$. In Figure 1 we show the data for $`L8`$ along with the best fit. The extrapolation $`L\mathrm{}`$ gives $`E_0=0.548824(2)`$. The finite-size correction constants are $`e_3=0.807(2)`$ and $`e_5=1.07(3)`$. Using Eq. (9) we obtain the spin-wave velocity $`c=1.123(2)`$. These results are in agreement with the previous extrapolations using exact diagonalization data for $`L=2,3,4,5,6`$, but the statistical errors are considerably smaller. The energy is also in perfect agreement with a previous Green’s function Monte Carlo calculation, which gave $`E_0=0.54883(1)`$.

We are not aware of any predictions for the size dependence of the spin stiffness of the XY model. For the Heisenberg model, the leading correction is $`O(1/L)`$. In contrast, our XY data can be very well fit with only an $`O(1/L^3)`$ term;

$$\rho (L)=\rho +\frac{r_3}{L^3}+\mathrm{},$$
(10)   
with the infinite-size value $`\rho =0.2696(2)`$ and the cubic correction $`r_3=0.47(1)`$. This fit is shown in Figure 2.

The stiffness result is in perfect agreement with the value obtained using exact diagonalization and finite-size scaling of the ground state energy in higher magnetization sectors. The energy per spin in the sector with magnetization $`j=_iS_i^z`$ should scale as 

$$E_jE_0=\frac{j^2c^2}{2\rho L^4}+O(j^4,1/L^6).$$
(11)   
We have also carried out some QMC calculations in the $`j=1`$ sector. Figure 3 shows results for $`L=8,10`$, and $`12`$, along with the exact result for $`L=6`$. Using also the exact $`L=4`$ result, a fit to Eq. (11), including an $`O(1/L^6)`$ term, gives $`c^2/\rho =4.70(4)`$. This value agrees with the above separate estimates of $`c`$ and $`\rho `$. The consistency between the results, obtained in two different ways, clearly gives very strong support to the Lagrangian theory.

The susceptibility can be obtained from $`c`$ and $`\rho `$ using the standard hydrodynamic relation

$$\chi =\rho /c^2,$$
(12)   
which with our values of $`c`$ and $`\rho `$ gives $`\chi =0.2138(8)`$.

For the magnetization, we assume a finite-size scaling

$$M_x^2(L)=M_x^2+\frac{a_1}{L}+\frac{a_2}{L^2}+\mathrm{}.$$
(13)   
A fit to this form, including the quadratic term, is shown in Figure 4 and gives $`M_x^2=0.0956(6)`$ and $`a_1=0.1212(5)`$. According to Eq. (7), the magnetization is thus $`m=0.437(2)`$, in good agreement with a previous Green’s function Monte Carlo result \[$`m=0.441(5)`$\]. For the Heisenberg model, the linear correction factor is related to ground state parameters according to $`a_1=\alpha M_x^2/(c\chi )`$, where $`\alpha =0.62075`$. The $`a_1`$ obtained here is instead consistent to within a few percent with $`a_1=\alpha M_x^2/(2c\chi )`$. As in the scaling of the energy, the leading size correction is hence proportional to the number of gapless modes in the symmetry-broken system.

## IV Finite-temperature results

We now discuss calculations aimed at extracting the leading finite-temperature corrections to the ground state. Previously, extensive simulations were carried out in order to study the KT transition, and a critical temperature $`T_{\mathrm{KT}}0.343J`$ was found. Here we focus on the susceptibility and the internal energy density at lower temperatures.

Figure 5 shows the susceptibility, defined in Eq. (5), calculated for lattices with $`L=32`$ and $`64`$ for temperatures $`T/J0.05`$. Within statistical errors there are no differences between the data for the two system sizes, and hence the results represent the thermodynamic limit. Chiral perturbation theory predicts a temperature-independent susceptibility for the XY model. Within error bars, our data is temperature independent for $`T/J0.15`$, with the constant value estimated at $`\chi =0.2096(2)`$. This result is lower (by about five standard deviations, or $`2`$%) than the value extracted in Sec III from the ratio $`\chi =\rho /c^2`$. The spin-wave velocity was obtained from the leading finite-size correction to the ground state energy. It is possible that this correction is affected by the presence of sub-leading corrections that are not completely taken into account by only including up to $`O(1/L^5)`$ terms in Eq. (8). Using the data shown in Fig. 5 is a more direct way of extracting the $`T=0`$ susceptibility, and we therefore consider it more reliable. An improved estimate of $`c`$ is then $`c=\rho /\chi =1.134(2)`$.

For a finite lattice, the gap between the $`j=0`$ ground state and the finite-magnetization states implies an exponential decay of $`\chi `$ to zero as $`T0`$. Chiral perturbation theory predicts that the form of this decay for $`TL^2\rho /c^21`$ is given by 

$$\chi =\frac{1}{TL^2}\mathrm{exp}[c^2/(2\rho TL^2)].$$
(14)   
In Figure 6 we show QMC data for $`L=416`$ along with this prediction, where we have used $`c^2/\rho =4.77`$, corresponding to our best estimates of $`c`$ and $`\rho `$. The agreement is not perfect, but satisfactory considering that there are no adjustable parameters and that there should also be corrections to Eq. (14). For $`L=4`$, we have also calculated $`\chi `$ using exact diagonalization down to much lower temperatures. As shown in Figure 7, there is a very good agreement with the predicted form over a sizable low-temperature range. Note, however, that the asymptotic $`T0`$ decay is always purely exponential, without the $`1/T`$ factor in Eq. (14), as only the lowest $`j=\pm 1`$ states contribute.

Chiral perturbation theory predicts that the low-temperature form of the internal energy is given by

$$E(T)=E(0)+\frac{\zeta (3)}{\pi c^2}T^3+O(T^5),$$
(15)   
where $`\zeta (3)=1.20206`$. In Figure 8 we show results for $`T/J0.05`$ calculated for $`L=64`$ lattices. We also show a comparison between $`L=16,32`$, and $`64`$, which suggests that $`L=64`$ gives the thermodynamic limit within error bars. A fit to the $`L=64`$ data with $`O(T^3)`$ and $`O(T^5)`$ terms, and including also the $`T=0`$ energy extracted in Sec. III, gives the cubic correction $`\zeta (3)/\pi c^2=0.284(5)`$. This corresponds to $`c=1.16(1)`$, which is higher than the value $`1.134(2)`$ extracted above using the estimates of $`\rho `$ and $`\chi `$ by about 2.5 error bars. Again, this small deviation may indicate some influence of higher-order corrections in the energy scaling.

## V Summary

In summary, we have presented extensive QMC calculations for the two-dimensional $`S=1/2`$ XY model. We have carried out finite-size and finite-temperature scaling of several physical quantities. The results are consistent with the predictions of effective Lagrangian theory  to within $`12`$%.

The best estimates of the ground state parameters resulting from our calculations are

$`E_0`$ $`=`$ $`0.548824(2),`$ (16)   
$`\rho `$ $`=`$ $`0.2696(2),`$ (17)   
$`m`$ $`=`$ $`0.437(2),`$ (18)   
$`\chi `$ $`=`$ $`0.2096(2),`$ (19)   
$`c`$ $`=`$ $`1.134(2).`$ (20)   
The ground state energy, $`E_0`$, the spin stiffness, $`\rho `$, and the square of the magnetization, $`m^2`$, were all calculated directly in the ground state for systems of linear size $`L=416`$ and extrapolated to infinite size. The susceptibility $`\chi `$ perpendicular to the XY spin-plane was calculated at finite temperature and extrapolated to $`T=0`$. Only the spin-wave velocity $`c`$ was obtained by a more indirect procedure, using the relation $`c^2=\rho /\chi `$. The results are in good agreement with previous exact diagonalization  and QMC work, but the precision is considerably improved over other estimates. Our results for both the energy and the magnetization are in remarkably good agreement with a series expansion calculation.  The magnetization is also in excellent agreement with second order spin-wave theory. The energy obtained in spin-wave theory  is only $`0.2`$% higher than the numerical result obtained here. It would be interesting to carry out high-order spin-wave calculations also for the other quantities discussed here.

## VI Acknowledgments

We would like to thank Jaan Oitmaa for discussions and comments on the manuscript. A. W. S. would like to thank the School of Physics at the University of New South Wales for hospitality and financial support during a visit. Support from the NSF under Grant No. DMR-9712765 is also acknowledged.