Constraint Effective Potential of the Staggered Magnetization in an Antiferromagnet
We employ an improved estimator to calculate the constraint effective potential of the staggered magnetization in the spin quantum Heisenberg model using a loop-cluster algorithm. The first and second moment of the probability distribution of the staggered magnetization are in excellent agreement with the predictions of the systematic low-energy magnon effective field theory. We also compare the Monte Carlo data with the universal shape of the constraint effective potential of the staggered magnetization and study its approach to the convex effective potential in the infinite volume limit. In this way the higher-order low-energy parameter is determined from a fit to the numerical data.
Systematic low-energy effective field theory is a very powerful tool for investigating the low-energy physics of Goldstone bosons. In particle physics this tool has been used to describe the low-energy physics of pions . Pions arise as the pseudo-Nambu-Goldstone bosons of the spontaneously broken chiral symmetry of QCD. Chiral perturbation theory provides a systematic momentum expansion of the low-energy physics of pions . Based on symmetry considerations, the effective theory makes detailed quantitative predictions, which depend on the values of some a priori unknown low-energy parameters, such as the pion decay constant and the chiral order parameter. Systematic effective field theories have also been developed for quantum field theories and classical spin systems with a spontaneously broken symmetry . Effective field theories with an symmetry have been used to describe the low-energy physics of magnons — the Goldstone bosons of quantum antiferromagnetism [4, 5, 6, 7]. In this case, the relevant low-energy parameters are the spin stiffness , the spinwave velocity , and the staggered magnetization density . In particular, the effective theory makes detailed predictions for the finite-size and finite-temperature effects of systems in large volumes at low temperatures. These effects have been worked out by Hasenfratz and Niedermayer even at the 2-loop level . By comparing the effective field theory predictions with Monte Carlo data obtained with an efficient loop-cluster algorithm , the low-energy parameters , , and have been determined with high accuracy for the spin quantum Heisenberg model, both on the square [10, 11] and on the honeycomb lattice , as well as for the - model with 2- and 4-spin interactions on the square lattice . Similar methods have been applied to the 4-d -symmetric quantum field theory describing the Higgs sector of the standard model . By comparing with effective field theory predictions for the susceptibility of the order parameter, both the vacuum expectation value of the Higgs field and the Higgs boson mass have been determined from Monte Carlo data for the 4-d lattice model [15, 16] obtained with the Wolff cluster algorithm . The systematic low-energy effective field theory also predicts the shape of the constraint effective potential of the order parameter, which has been worked out in great detail by Göckeler and Leutwyler [18, 19]. Their predictions for the shape of the constraint effective potential have been tested against Monte Carlo simulations of classical 3-d and 4-d lattice models .
In this paper, we develop an improved estimator for the probability distribution of the staggered magnetization in order to extract the constraint effective potential from loop-cluster simulations of the spin quantum Heisenberg model on the square lattice. The first moments of the probability distribution are in excellent agreement with the predictions of the magnon effective theory. The Monte Carlo data approach the universal shape of the effective potential in the infinite volume limit, which is known to be a convex function . Our study confirms that the effective field theory provides quantitative predictions for a wide class of low-energy observables, which are exact, order by order in a systematic momentum expansion.
The rest of this paper is organized as follows. Section 2 summarizes the predictions of the low-energy magnon effective field theory that are relevant for our study, and section 3 discusses the determination of the low-energy parameters. In section 4 we describe the improved estimator that allows us to obtain very accurate Monte Carlo data for the probability distribution of the staggered magnetization. In section 5 the results of numerical simulations obtained with a loop-cluster algorithm are compared with the effective field theory predictions. Finally, section 6 contains our conclusions.
2 Effective Field Theory Predictions
In this section we summarize the predictions of Göckeler and Leutwyler [18, 19] which are derived from a 3-d -symmetric effective field theory. Although these authors had applications to the 3-d classical Heisenberg model in mind, their results also apply to the antiferromagnetic -d quantum Heisenberg model. This microscopic model is defined by the Hamiltonian
where and are nearest-neighbor sites on a square lattice with spacing . The spin operators obey the standard commutation relations
Note that we work in natural units in which . Furthermore, represents an external staggered magnetic field and
is the staggered magnetization order parameter. In the infinite volume limit at zero temperature and with , develops a non-zero vacuum expectation value signaling the spontaneous breakdown of the spin symmetry down to its subgroup.
The corresponding low-energy effective field theory is formulated in terms of the staggered magnetization order parameter field
where is a point in Euclidean space-time. The leading terms in the effective action for the staggered magnetization field take the form
where is the spin stiffness, is the spinwave velocity, and is the staggered magnetization density. The corresponding partition function is given by
Antiferromagnetic magnons have a “relativistic” dispersion relation with the spinwave velocity playing the role of the velocity of light. In fact, by introducing , the effective action can be written in the manifestly Euclidean rotation-invariant form
which ensures Lorentz-invariance after analytic continuation from Euclidean to real time. It should be pointed out, however, that Euclidean rotation invariance is just an accidental symmetry of the leading terms of the effective action. Since the underlying quantum Heisenberg model itself does not have this symmetry, some of the higher-order four-derivative terms in the effective action break Euclidean rotation-invariance down to the discrete rotation subgroup of the square lattice.
It is interesting to note that the ferromagnetic quantum Heisenberg model (which differs from the antiferromagnet only by the sign of the Hamiltonian) has very different symmetry properties at low energies. Unlike quantum antiferromagnets, quantum ferromagnets have a conserved order parameter — the uniform magnetization. Consequently, ferromagnetic magnons have a nonrelativistic dispersion relation and the corresponding effective action contains an additional Wess-Zumino term , which breaks rotation invariance between space and Euclidean time already at leading order. The resulting effective field theory for ferromagnetic magnons has been studied in detail in .
Here we concentrate entirely on antiferromagnets which are described by the effective action of eq.(2.5). We consider the system in a periodic cubic space-time volume with the inverse temperature fixed at and with . The space-time average of the staggered magnetization is given by
In contrast to [18, 19], we have included a factor in the definition of because the quantum spins of the underlying Heisenberg model have while the effective field is normalized to 1. Due to the symmetry, the probability distribution
of the mean staggered magnetization vector only depends on the magnitude . The distribution is normalized by
The constraint effective potential represents the free energy density of configurations constrained to a fixed mean staggered magnetization and is determined by
with the normalization derived in  given by
Here is the staggered magnetization per spin, while and are shape-coefficients of the space-time box. For the exactly cubic space-time volume considered here and . In the infinite volume and zero temperature limit the constraint effective potential approaches the infinite volume effective potential which is known to be a convex function of . In a finite volume, on the other hand, is not convex. An extensive variant of the intensive quantity is
Here and are -independent functions of the rescaled variable
The leading order contribution to the constraint effective potential is given by an inverse Laplace transform
of the function
Here the are shape-coefficients of the finite space-time box described in detail in appendix B of . As was pointed out in , the function is entirely kinematical and thus, unlike , the quantity is universal, i.e. completely independent of the low-energy parameters. Consequently, is the same for all 3-d systems with an symmetry spontaneously broken down to , including the 3-d classical and the -d quantum Heisenberg model. The correction to the leading contribution is given by
Here is a low-energy constant related to the higher-order terms
in the effective action, which is given by
It should be noted that the value of has no impact on eq.(2.12) because it affects the normalization factor only at higher orders in .
Just as a non-zero quark mass in QCD explicitly breaks chiral symmetry and gives the pion its mass, a non-zero staggered magnetic field explicitly breaks the spin symmetry and turns the magnons into pseudo-Nambu-Goldstone bosons with a non-zero mass determined at leading order by
The constant also appears in the -dependence of the field expectation value
It should be noted that eq.(2.22) was derived in the -regime of chiral perturbation theory in which while , . In particular, in eq.(2.22) one cannot make (and thus ) arbitrarily small, because one would otherwise enter the -regime in which .
Besides the constraint effective potential, Göckeler and Leutwyler have also derived analytic predictions for the first and second moment of the probability distribution , including the 2-loop level. The resulting expressions are
For the cubic box considered here, the additional shape-coefficient is given by .
Hasenfratz and Niedermayer have used the effective theory to derive the finite-size and finite-temperature effects of the staggered susceptibility
from a 2-loop calculation in the -regime of magnon chiral perturbation theory . Similarly, the uniform susceptibility takes the form
Here determines the shape of an approximately cubic space-time box of size , with . The functions , , and are known shape-coefficients [3, 8]. For an exactly cubical space-time volume (i.e. for ) the result of eq.(2.24) agrees with eq.(2) since
The factor 3 arises due to the three components of the staggered magnetization vector.
3 Determination of the Low-Energy Parameters , , and
The susceptibilities and have been calculated numerically for the antiferromagnetic spin quantum Heisenberg model on the square lattice using the very efficient loop-cluster algorithm . By comparing the Monte Carlo data with the effective theory predictions of eqs.(2.24) and (2.25), the low-energy parameters , , and have been determined with high precision. At very low temperatures, one enters the cylindrical regime of space-time volumes with in which the -expansion of chiral perturbation theory applies. In this case, the staggered magnetization acts as a quantum mechanical rotor, again resulting in characteristic finite-volume effects . By simulating cylindrical space-time volumes using the continuous-time variant of the loop-cluster algorithm, and again comparing with the corresponding predictions of the low-energy effective theory, the values of the low-energy parameters, previously obtained from the cubical space-time regime, have been verified independently . Using the value of the exact mass gap of the 2-d model , a result of Chakravarty, Halperin, and Nelson for the finite-temperature correlation length  was extended by Hasenfratz and Niedermayer  who obtained
This expression is valid in space-time volumes with a slab geometry, i.e. for . While the data in the cubical and cylindrical regimes determine the ratio with higher precision than and individually, Monte Carlo data for the very long correlation length in the slab regime determine very precisely . The combined analysis of all numerical data in cubical, cylindrical, and slab geometries resulted in and . In a recent study using a zero-temperature valence-bond projector method, Sandvik and Evertz obtained the very accurate result . Although the discrepancy between the two results for is only in the permille range, it is statistically significant. In particular, it is important to clarify the discrepancy because our present very accurate study is sensitive to such small effects. For this purpose we have generated new data in the cubical regime for volumes ranging from to 80. The largest volumes are substantially bigger than those of the original study . In order to have an independent handle on the spinwave velocity , in addition to the uniform susceptibility , which is given by the temporal winding number , we have also measured the spatial winding numbers . The condition
determines an exactly cubical space-time box with . By varying until eq.(3.2) is satisfied, we have determined in excellent agreement with the result of . By fitting the new data for and in the cubical regime to eqs.(2.24) and (2.25), we have obtained and . The value for is again in excellent agreement with the previous results obtained in the cubical, cylindrical, and slab regimes. The new result for agrees with the one of the original study in the cubical regime  and is consistent with the result of , which is about two permille lower than the result of . We attribute this small but statistically significant discrepancy to an underestimation of the systematic errors of in the cylindrical regime data of , related to the termination of the Seeley expansion described in . Until this issue is completely clarified, we discard the cylindrical regime data and instead include the result of . The best estimate of the low-energy parameters obtained in cubic and slab geometries as well as at zero temperature is then given by
4 Improved Estimator for the Distribution of the Staggered Magnetization
The loop-cluster algorithm [9, 10, 11] is a very efficient numerical tool that allows us to perform high-accuracy numerical simulations of the quantum Heisenberg model.The cluster algorithm connects the spin variables to closed loop-clusters, which are completely independent of one another. All spins in a given cluster are then flipped simultaneously with 50 percent probability. A given spin configuration containing clusters is just one member of a sub-ensemble of equally probable configurations. An improved estimator substantially increases the statistics by analytically averaging a given observable over all configurations in the sub-ensemble. For the quantum Heisenberg model, improved estimators have been constructed previously for the staggered and uniform susceptibilities and as well as for the energy density . The improved estimator for the distribution of the staggered magnetization to be constructed here is similar to the improved estimator for the topological charge distribution in the meron-cluster algorithm for the 2-d classical model , which has been combined with a re-weighting technique .
In the loop-cluster algorithm, every cluster contributes additively to the total 3-component of the staggered magnetization. While it is straightforward to implement the improved estimator in continuous Euclidean time, it is most easily explained in the discrete-time variant of the loop-cluster algorithm [9, 10]. In that case, the cluster size (i.e. the number of lattice points in a given cluster) determines the 3-component of the staggered magnetization of the cluster , which is proportional to . Under cluster flip the staggered magnetization of a cluster changes sign. The distribution of the staggered magnetization is recorded as a histogram which is built iteratively using one cluster after the other. The initial histogram (with , where is the number of space-time lattice points) is constructed from the first cluster as
The two entries of the initial histogram correspond to the two possible orientations of the first cluster, each arising with probability . In the -th iteration step (with ) a new histogram is built from the previous one as
After steps, all clusters have been incorporated and the final histogram is given by . Examples of histograms obtained for two individual spin configurations are shown in figure 1. The example in the left panel contains one cluster that is bigger than all the other clusters together. Hence, the region around is not sampled. In addition, there are two relatively large clusters that give rise to the multiple peaks in the distribution. In the example shown in the right panel, on the other hand, there are two clusters of similar size, such that the region around is also sampled.
The final probability distribution of the 3-component of the staggered magnetization
is the average of the histograms for all configurations in the Markov chain generated by the cluster algorithm. By construction, this distribution is properly normalized as
It should be noted that the numerical effort to build the improved estimator is proportional to the number of lattice points and, in addition, proportional to the number of clusters. Since the number of clusters is proportional to the volume, the evaluation of the improved estimator requires a computational effort proportional to , and thus becomes rather time-consuming for large volumes. This is in contrast to the improved estimators for the susceptibilities and which only require a computational effort proportional to . Of course, the improved estimator increases the statistics by a factor of which is exponential in the volume. Hence, investing a polynomial effort should still be justified. However, one should not forget that the configurations in a sub-ensemble are not statistically independent. As we will see, the improved estimator works very well and by far outperforms calculations done without it.
The mean value of the 3-component of the staggered magnetization corresponding to a given value of is
The factor 2 in the denominator arises because we are dealing with quantum spins . Hence, one can identify the probability distribution of the 3-component of the mean staggered magnetization as
Due to the time-discretization chosen in  (which has four Trotter steps), all cluster sizes are multiples of 8. Consequently, the non-zero entries of the histogram correspond to values of which are also multiples of 8. This implies that
Altogether, we thus obtain
with restricted to multiples of 8. By construction, in the Euclidean time continuum limit the resulting probability distribution is normalized as
We have simulated the spin quantum Heisenberg model on a square lattice with between and at inverse temperatures using the loop-cluster algorithm in its discrete-time variant [9, 10]. We have worked at a sufficiently small lattice spacing in Euclidean time, such that the systematic discretization error is negligible compared to the statistical errors. The probability distribution of the 3-component of the mean staggered magnetization has been obtained using the improved estimator described above. A typical distribution is shown in figure 2.
As we will see below, the information about the vicinity of the minimum of the constraint effective potential is contained in the region of where changes rapidly. Figures 2 and 3 compare Monte Carlo data obtained with and without the improved estimator, investing the same amount of computer time in both cases. The error reduction of the improved estimator is very substantial.
By employing the re-weighting technique used in [29, 28], one could concentrate the statistics in the relevant region in which changes rapidly. This should further improve the efficiency of our numerical method.
Due to the symmetry, for a fixed magnitude of the mean staggered magnetization vector , its 3-component has a flat distribution given by
Here is a step function which is equal to 1 for and zero otherwise. The probability distribution of the 3-component and the one of the magnitude are related by
Hence, given the Monte Carlo data for , we can extract the probability distribution of the magnitude of the staggered magnetization as
This indeed ensures the correct normalization of eq.(2.10) because
Using eq.(4.13), we have determined the probability distributions from obtained using the improved estimator. Some results for are shown in figure 4.
As the volume increases the mean value of decreases and the width of the distribution becomes narrower. It should be noted that the distribution is not symmetric around its maximum. In the infinite volume limit, the distribution turns into a -function centered at .
We also like to compute the first and second moment and of the distribution . Of course, this is trivial once we have computed using the improved estimator which requires a computational effort proportional to . If one only wanted to compute the two moments but not itself, one may wonder how to achieve this in the most efficient manner. As we noted before, the combination is proportional to the staggered susceptibility which can be obtained using an improved estimator requiring a computational effort proportional to . How can one determine itself? Since is the magnitude of the mean staggered magnetization vector , of which only the 3-component is easily accessible in a quantum Monte Carlo simulation, this seems not entirely straightforward. However, using eq.(4.13) one obtains
Hence, by just measuring in a standard unimproved manner, one can determine with a computational effort proportional to . The crucial question is how the statistical error achieved in this way compares with the one obtained by extracting from (whose construction requires a computational effort proportional to ). In order to investigate which of the two methods is more efficient, we have first determined from using the improved estimator described above. Then we have invested the same amount of computer time in an unimproved measurement of . For and the statistical error of obtained in this manner is a factor of 1.5 larger than when one uses the improved estimator. Hence, despite its computational effort proportional to , thanks to the average over the large number of configurations in the sub-ensemble, the improved estimator is slightly better than the unimproved method, even if one is only interested in the first moment and not in the distribution itself.
In order to compare our Monte Carlo data with the expectation value of eq.(2.22), we also like to switch on a staggered magnetic field . The corresponding probability distribution
can be sampled using the loop-cluster algorithm with an additional Metropolis accept-reject step applied to each cluster flip, which takes into account the contribution to the Boltzmann weight. Using an ordinary unimproved estimator, one then simply measures .
5 Comparison of Monte Carlo Simulations and Effective Theory Predictions
We have determined the first and second moment and of the distributions , which are compared with the effective field theory predictions of eq.(2) in table 1. The errors of the theoretical predictions result from the uncertainties in the low-energy parameters of eq.(3.3). For the first moment the agreement is very good for . The absolute value of the second moment is very small and its statistical error is relatively large. Still, there are systematic discrepancies between the Monte Carlo data and the effective theory predictions of eq.(2). This discrepancy is well accounted for by additional corrections. Such corrections involve new low-energy parameters multiplying higher-order terms in the effective action. Their evaluation would require a 3-loop calculation which has not been worked out in the effective theory. Parameterizing the 3-loop terms with unknown coefficients and , i.e.
one obtains a good fit to the Monte Carlo data for and . One may conclude that precise calculations of the two moments allow the determination of some combination of sub-leading low-energy parameters.
Using , the probability distributions of figure 3 are readily converted into the corresponding constraint effective potentials shown in figure 5.
With increasing volume the constraint effective potential approaches the convex shape of the infinite volume effective potential. Using the rescaled variable , one can also consider the extensive quantity which is shown in figure 6.
Expanding , we have used the Monte Carlo data for between and to determine and . Some values of the function extracted from the numerical data are compared with the analytic result of eq.(2.16) in figure 7. It should be pointed out that the observed agreement does not rely on any adjustable parameters. Even the normalization constant of eq.(2.12), which fixes an additive constant in the constraint effective potential, is predicted by the effective theory.
As quantified in table 2, in the interval , i.e. around the minimum of the constraint effective potential, the theoretical values and the simulation data for agree within error bars.
|0.10||- 3.782(62)||- 3.818|
|0.15||- 4.271(39)||- 4.302|
|0.20||- 4.578(22)||- 4.589|
|0.25||- 4.646(18)||- 4.653|
|0.30||- 4.458(20)||- 4.462|
|0.35||- 3.983(36)||- 3.986|
In order to determine the low-energy parameter , we have considered the expectation value as a function of the external staggered magnetic field given in eq.(2.22). Since these results do not require the probability distribution , they could be obtained on larger volumes up to . The results summarized in table 3 give a good fit for
Using the theoretical prediction for , figure 8 compares values of the function determined from the Monte Carlo data with the analytic result of eq.(2.18), for the two values of at the edge of the corresponding error band. In the interval around the minimum of the constraint effective potential, the Monte Carlo data are consistent with the theoretical predictions of Göckeler and Leutwyler. We attribute the deviations outside this interval to effects of higher order. Indeed, as one sees in figure 6, the finite volume effects become larger when moves away from the minimum of the constraint effective potential.
We have constructed an improved estimator for the probability distribution of the staggered magnetization in the quantum Heisenberg model. Using the improved estimator in a loop-cluster algorithm simulation, we have determined the first and second moment and of the distribution of the magnitude of the mean staggered magnetization vector , as well as the constraint effective potential (obtained from ) for different space-time volumes. The Monte Carlo data are in excellent quantitative agreement with analytic predictions which Göckeler and Leutwyler derived from a systematic low-energy effective field theory. This demonstrates that the magnon effective theory indeed provides correct predictions, order by order in a systematic low-energy expansion. Thanks to the very efficient loop-cluster algorithm, the Heisenberg model is an excellent testing ground for the effective field theory method. Its quantitative success provides encouragement to also test the recently constructed systematic effective field theories for holes [30, 31] and electrons  doped into an antiferromagnet against numerical simulations. First results confirming the effective theory have already been obtained for the - model on the honeycomb lattice . The results obtained in the study presented here should also be encouraging for lattice QCD simulations, were the numerical problem is much harder. Eventually, one may expect agreement between lattice QCD and chiral perturbation theory at the same level of accuracy as achieved in the condensed matter problem investigated in this paper.
We have benefited from correspondence and discussions with B. B. Beard, M. Göckeler, P. Hasenfratz, K. Jansen, F. Niedermayer, and H. Leutwyler. C. P. H. would like to thank the members of the Institute for Theoretical Physics at Bern University for their hospitality during a visit where this project was initiated. The work of C. P. H. is supported by CONACYT Grant No. 50744-F and by Grant Proyecto Cuerpo-Academico-56-UCOL. This work is supported in part by funds provided by the Schweizerischer Nationalfonds (SNF). The “Center for Research and Education in Fundamental Physics” at Bern University is supported by the “Innovations- und Kooperationsprojekt C-13” of the Schweizerische Universitätskonferenz (SUK/CRUS).
-  S. Weinberg, Physica 96 A (1979) 327.
-  J. Gasser and H. Leutwyler, Nucl. Phys. B250 (1985) 465.
-  P. Hasenfratz and H. Leutwyler, Nucl. Phys. B343 (1990) 241.
-  S. Chakravarty, B. I. Halperin, and D. R. Nelson, Phys. Rev. B39 (1989) 2344.
-  H. Neuberger and T. Ziman, Phys. Rev. B39 (1989) 2608.
-  D. S. Fisher, Phys. Rev. B39 (1989) 11783.
-  P. Hasenfratz and F. Niedermayer, Phys. Lett. B268 (1991) 231.
-  P. Hasenfratz and F. Niedermayer, Z. Phys. B92 (1993) 91.
-  H. G. Evertz, G. Lana, and M. Marcu, Phys. Rev. Lett. 70 (1993) 875.
-  U.-J. Wiese and H.-P. Ying, Z. Phys. B93 (1994) 147.
-  B. B. Beard and U.-J. Wiese, Phys. Rev. Lett. 77 (1996) 5130.
-  F. J. Jiang, F. Kämpfer, M. Nyfeler, and U.-J. Wiese, Phys. Rev. B78 (2008) 214406.
-  F. J. Jiang, M. Nyfeler, S. Chandrasekharan, and U.-J. Wiese, J. Stat. Mech. (2008) P02009.
-  A. Hasenfratz, K. Jansen, J. Jersák, C. B. Lang, H. Leutwyler, and T. Neuhaus, Z. Phys. C46 (1990) 257.
-  M. Göckeler, K. Jansen, and T. Neuhaus, Phys. Lett. B273 (1991) 450.
-  M. Göckeler, H. A. Kastrup, T. Neuhaus, and F. Zimmermann, Nucl. Phys. B404 (1993) 517.
-  U. Wolff, Phys. Rev. Lett. 62 (1989) 361.
-  M. Göckeler and H. Leutwyler, Nucl. Phys. B350 (1991) 228.
-  M. Göckeler and H. Leutwyler, Phys. Lett. B253 (1991) 193.
-  I. Dimitrovic, J. Nager, K. Jansen, and T. Neuhaus, Phys. Lett. B268 (1991) 408.
-  L. O’Raifeartaigh, A. Wipf, and H. Yoneyama, Nucl. Phys. B271 (1986) 653.
-  H. Leutwyler, Phys. Rev. D49 (1994) 3033.
-  C. P. Hofmann, Phys. Rev. B60 (1999) 388; Phys. Rev. B60 (1999) 406; Phys. Rev. B65 (2002) 094430; AIP Conference Proceedings vol.623 (2002) 305.
-  P. Hasenfratz, M. Maggiore, and F. Niedermayer, Phys. Lett. B245 (1990) 522.
-  P. Hasenfratz and F. Niedermayer, Phys. Lett. B245 (1990) 529.
-  B. B. Beard, R. J. Birgeneau, M. Greven, and U.-J. Wiese, Phys. Rev. Lett. 80 (1998) 1742.
-  A. W. Sandvik and H. G. Evertz, arXiv:0807.0682 [cond-mat.str-el].
-  W. Bietenholz, A. Pochinsky, and U.-J. Wiese, Phys. rev. Lett. 75 (1995) 4524.
-  U.-J. Wiese, Nucl. Phys. B318 (1989) 153.
-  F. Kämpfer, M. Moser, and U.-J. Wiese, Nucl. Phys. B729 (2005) 317.
-  C. Brügger, F. Kämpfer, M. Moser, M. Pepe, and U.-J. Wiese, Phys. Rev. B74 (2006) 224432.
-  C. Brügger, C. P. Hofmann, F. Kämpfer, M. Moser, M. Pepe, and U.-J. Wiese, Phys. Rev. B75 (2007) 214405.