The meson from lattice QCD.
K. Jansen, C. Michael, C. Urbach
DESY, Zeuthen, Platanenallee 6, D-15738 Zeuthen, Germany
Theoretical Physics Division, Dept. of Mathematical Sciences,
University of Liverpool, Liverpool L69 7ZL, UK
Institut für Elementarteilchenphysik, Fachbereich Physik,
Humbolt Universität zu Berlin, D-12489, Berlin, Germany
We study the flavour singlet pseudoscalar mesons from first principles using lattice QCD. With flavours of light quark, this is the so-called meson and we discuss the phenomenological status of this. Using maximally twisted-mass lattice QCD, we extract the mass of the meson at two values of the lattice spacing for lighter quarks than previously discussed in the literature. We are able to estimate the mass value in the limit of light quarks with their physical masses.
There is considerable interest in understanding hadronic decays involving and in the final state. The phenomenological study of hadronic processes involving flavour singlet pseudoscalar mesons makes assumptions about their composition. Here we address the issue of the nature of the and from QCD directly, making use of lattice techniques.
Lattice QCD directly provides a bridge between the underlying quark description and the non-perturbative hadrons observed in experiment. The amplitudes to create a given meson from the vacuum with a particular operator made from quark fields are measurable, an example being the determination of . It also allows a quantitative study of the disconnected quark contributions that arise in the flavour singlet sector. The lattice approach provides other information such as that obtained by varying the number of quark flavours and their masses.
The disconnected diagram responsible for giving the flavour-singlet pseudoscalar meson a mass (in the chiral limit) is closely related to the fluctuations in the topological charge. This link will be discussed shortly in this paper and in more detail in a later publication .
From the chiral perturbation theory description (for a review, see ref. ), one expects the mixing of and to be most simply described in a quark model basis. In the flavour singlet sector, for pseudoscalar mesons, we then have contributions to the mass squared matrix with quark model content and (which we label as and respectively):
Here corresponds to the mass of the flavour non-singlet eigenstate and is the contribution to the mass coming from connected fermion diagrams while corresponds to the contribution from disconnected fermion diagrams. Thus is the pion mass. Because of mixing, as will be discussed, does not correspond to any specific meson, but its mass can be estimated assuming that the non-singlet pseudoscalar mass-squared is linear in the quark mass. So , that is , leading to GeV. For a discussion from lattice results of small corrections to this assumption, see ref. . Chiral perturbation theory also gives corrections to linearity coming from loop corrections.
The mixing between the and flavour singlet channels must produce the experimental and . One approximation used historically is that of SU(3) flavour symmetry for which the two physical states will be a flavour octet and a flavour singlet . The is then identified as primarily the flavour octet while the is the flavour singlet.
Using as input and and requiring that the output masses ( and ) are correctly reproduced, the three mixing parameters cannot be fully determined. It is usual to express the resulting one parameter freedom in terms of a mixing angle, here defined by
In our lattice study, we have two degenerate light quarks in the sea but no explicit strange quark. Since we do not consider partially quenched QCD here, we only study the top left corner of the mixing matrix. The flavour singlet pseudoscalar meson (called in this case) will then have mass-squared .
From the phenomenological analysis of the full mixing matrix, one can then estimate the mass of the meson. One such phenomenological analysis, motivated by lattice input (basically the magnitude of for strange quarks), gave  values GeV. This assignment corresponds to a mixing angle close to and to a mass of the meson of 0.776 GeV. This mass value is plausible since it is close to the (mass-squared weighted) average mass of the and . Moreover, it is somewhat lighter than the meson, which has an additional contribution to its mass from strange quark loops (). However, here we plan to determine the meson mass more directly, avoiding some of the assumptions made above.
Note that the mass of the meson comes from two components: which decreases with decreasing quark mass and which, in the phenomenological study, increases with decreasing quark mass. It is thus important to study the meson at light quark masses to explore this. For instance, with the parameters of the phenomenological model described above, the meson made of light quarks of mass equal to the strange quark mass would have mass-squared giving mass 0.862 GeV. Thus we can expect a rather flat dependence of the mass on the mass of the underlying quarks (here a degenerate pair of valence quarks with sea quarks of the same nature).
Lattice QCD is able to access flavour-singlet states, but at a cost which has limited these explorations in practice. Since the correlator of a flavour singlet meson created at and annihilated at will have two contributions: a connected contribution where both quark and antiquark propagate from to and a disconnected contribution where there is a quark loop at and another at . This latter contribution needs to be measured at many values of and , so it is optimal to use a stochastic method to achieve this. One then evaluates the disconnected contributions (loops) at all , and from combining them pair-wise, one can evaluate the disconnected contribution to the meson correlator. The stochastic method introduces noise, which can be reduced by appropriate variance reduction techniques. The goal is to reduce the noise to a level lower than the inherent noise coming from the underlying gauge-field variation. We shall show that the formalism known as twisted mass lattice QCD allows a very effective variance reduction method to be applied.
We present details of our lattice methods and explain the stochastic method used to evaluate the flavour singlet pseudoscalar meson correlators (both connected and disconnected). Because the signal from the disconnected correlator is relatively noisy, we discuss several strategies for reducing this error. From a combination of methods, we are able to present quite precise results, which enable us to discuss the continuum and chiral limit of the flavour singlet pseudoscalar mass. We compare with previous determinations and summarise.
2 Twisted mass lattice QCD action
In the gauge sector we employ the so-called tree-level Symanzik improved gauge action (tlSym) , viz.
where is the untwisted bare quark mass, is the bare twisted quark mass, is the third Pauli matrix acting in flavour space and
is the mass-less Wilson-Dirac operator. and are the forward and backward covariant difference operators, respectively.
The mass term is tuned to maximal twist (as described in ref. ) at our lightest parameter. This guarantees improvement . This was shown to work excellently in the quenched approximation [8, 9], and there are good indications that this holds true also with dynamical fermions [10, 11]. Moreover, this formalism has been found to give a very effective way to reach light quark masses [7, 12, 13]. As such, it provides an attractive route to evaluate the properties of the flavour-singlet pseudoscalar mesons.
There is one complication, however, namely that twisted mass lattice formalism breaks flavour and parity symmetries at finite value of the lattice spacing . These symmetries are restored in the continuum limit and the theory is well defined (so providing a valid regularisation) at finite lattice spacing. For our present purpose, this flavour-breaking causes the and mesons to have a mass splitting (of order ). Moreover the has contributions from disconnected diagrams [14, 7]. The analysis of the flavour-singlet pseudoscalar meson is not significantly affected by this and it can be studied in the same way as was used previously in Wilson-based lattice formalisms [3, 15]. In the twisted mass formalism, the meson can mix in principle with the neutral meson. This is an order mixing which could induce an order contribution to the mass observed at finite lattice spacing. We expect this effect to be negligible, both because the meson is heavy and because we are working at maximal twist. We check this by varying the lattice spacing . For a recent review see ref. .
For later convenience we introduce the following notation
for the fermion matrix of the two degenerate quark flavours (here labelled and ) separately.
The results presented in this paper are based on gauge configurations as produced by the European Twisted Mass collaboration (ETMC). The details of the ensembles are described in ref. . In particular, we concentrated for this paper on the ensembles labelled and and . We have compiled the details for those ensembles in table 1. The -ensembles correspond to a value of the lattice spacing of about and the -ensembles to . The spatial lattice size is of about for all ensembles used here, apart from , which has identical parameters to but .
3 Neutral particles in twisted mass QCD
In quenched or partially-quenched lattice QCD, the disconnected contribution to the flavour singlet meson does not combine properly with the connected contribution to give a physical state. To avoid this problem, it is mandatory to study flavour singlet states in full QCD - with sea quarks having the same properties as valence quarks. Then the spectrum of flavour singlet states is well defined and can be extracted from the -dependence of the full correlator. Here we focus on the case where there are two degenerate light quarks (called ) which is a consistent theory in which to study the flavour singlet mesons.
One promising way to study light quarks is using twisted mass QCD at maximal twist. In the case of twisted mass fermions, the bilinear operator appropriate to create the state is which on transformation into the twisted basis (used in lattice evaluation) will become , where acts in the flavour space. This amounts to evaluating, in the lattice basis, the difference of the disconnected loop between and quarks. As already reported [17, 18], this enables a very efficient variance reduction technique to be used to evaluate the disconnected diagram relevant to the .
The key observation is the following relation for the inverse and :
Consider now the disconnected loop where is some -matrix (here ) and/or colour-matrix and the sum is over space. The conventional method involves solving with stochastic volume sources . The desired result is then obtained from
This has signal/noise which has a volume dependence which is much more favourable than the conventional method with signal/noise (here ). For ensemble , we find the zero-momentum disconnected loop at a given time-value has a standard deviation of for the variation with gauge configuration and time-slice (here called the intrinsic variation and obtained by extrapolating to an infinite number of stochastic samples) whereas the stochastic noise on this loop has a standard deviation of from 24 samples of volume source (conventional method) but only from the above method (with 12 samples). The relevant standard deviation for any analysis is the folding of the intrinsic variation () with the stochastic error. So with a stochastic error of the net standard deviation is 19.5 whereas with the conventional method it would be 89. Thus there is a significant improvement. So 12 inversions give the disconnected correlator from all to all with no significant increase in errors from the stochastic evaluation.
We also use a non-local source/sink for the meson, constructed using ”fuzzed” links of length as described in refs. [21, 18], which can be evaluated by replacing by the corresponding fuzzed gauge links. Thus twisted mass lattice QCD allows a very efficient way to evaluate flavour-singlet pseudoscalar disconnected loops, and combining them, the required disconnected correlator .
Evaluating the connected correlator for the neutral pseudoscalar meson as described elsewhere , then we have the full information needed to explore the spectrum in the channel. We have a matrix of such correlators available (local or non-local at sink/source). We could also explore pseudoscalar correlators obtained by creating the state with which on twisting becomes . This does not allow our very efficient variance reduction to be applied, so the disconnected contributions are evaluated using a hopping parameter variance reduction [22, 18]. Even so, they are sufficiently noisy that we have been unable to use them to constrain the fits or to determine the decay constant.
We measure the auto-correlation versus trajectory number for our correlators. The largest auto-correlation is found for the disconnected contribution. To explore this more fully, we investigated the auto-correlation for the time-slice sum (i.e. zero momentum) of the quark loop at a given time versus trajectory, finding comparable auto-correlation times to that of the plaquette . This is not unexpected, since the disconnected contribution is related to the topological charge density and thus both it and the plaquette are sensitive to the vacuum structure encoded in the gauge configurations. For further discussion see ref. . In order to cope with this measured auto-correlation, we block the data into sufficiently large blocks (more than trajectories) to remove any effect on our final error estimates. Note that the gauge configurations we use for our measurements are usually well separated in units of trajectories, such that even after blocking trajectories we are left with a sufficiently large number of (then independent) measurements.
To explore the signal, we plot the effective mass versus for ensemble in fig. 1, where we have used the variational basis from -values and to optimise the ground state contribution. The zero momentum results (open symbols) show a plateau, although statistical errors are large, as will be discussed later. The results for the two lattice spatial volumes available in this case agree well.
In order to explore options to reduce the errors, we also evaluate the full correlator for momentum (in units of ). Since we evaluate the disconnected correlator for momentum 1 in each of the 3 spatial directions, we obtain an improved estimate of (for example with relative error reduced to 11% at compared to error of 14% for zero-momentum). Since the symmetry classification of mesonic states is less sharp at non-zero momentum, the interpretation of the lattice spectrum is thus less straightforward. In practice, we find that the small reduction in error on does not translate into a significant reduction in the error in determining the underlying mass of the state.
We make a fit to the combined matrix of zero-momentum correlators for a -range of 3-10 and with two states. The results from an uncorrelated fit to blocked data are presented in table 2. These fits are stable to varying the -range.
Before analysing these results further, we now discuss why the errors are so large for the disconnected correlator, despite the fact that we measure all and , we use many gauge configurations and our stochastic errors are small. As described above, we have investigated the autocorrelation (versus trajectory) of the disconnected contribution and find no statistically significant evidence of any autocorrelation beyond 40 trajectories. We block our data into larger blocks (typically covering trajectories) to avoid any increase in error from this source.
The origin of the large error is that the signal for the disconnected part of the correlator comes from only a small part of the total data sample. Here we concentrate on the zero-momentum method to illustrate this. We study where labels the gauge configuration and we explore the distribution versus at fixed . For instance (for ensemble from -values for 888 gauge configurations) with , of the gauge configurations contribute of the observed signal (i.e. the sum over all gauge configurations). Thus the statistical impact of the data set is smaller than expected since parts of the data have big fluctuations (in a fermionic loop related to topological charge density).
Another way to illustrate this error problem is that from the first gauge configurations (covering trajectories) we find at , while for the second gauge configurations we find . Positivity requires . So the error estimate for certain quantities from even such a big ensemble as trajectories can be underestimated compared to having trajectories, while other quantities do not show such a problem on the same set of gauge configurations.
So even more configurations, than we have here, would be needed to get reliable and small errors in the case of disconnected contributions. This same conclusion has been obtained before, most strongly in a study involving staggered fermions .
4 Error reduction strategies
Because the statistical error on the mass is still large, we now discuss various strategies to reduce it while retaining the same set of gauge configurations. Some of these strategies do indeed reduce the statistical error, but at the expense of introducing a systematic error that has to be discussed.
4.1 Excited state removal
We now consider replacing the connected neutral pseudoscalar correlator by just the ground state contribution. This has been considered previously . The basic idea is that the disconnected contributions () are only big for the lightest flavour-singlet pseudoscalar (the ) and not for excited states. In other words, and are split but and are almost degenerate. This might happen because the topological charge fluctuations in the vacuum, which give the flavour singlet states a different mass, are more strongly coupled to the ground state than to its excited states. Note that this is an assumption which needs to be checked by looking at the lattice results.
In order to have no significant excited state contributions in the appropriate flavour-singlet correlator , one should remove them from . Then neither nor will have excited state contributions, and a one state fit to should be possible down to quite low -values, as indeed we shall find. The extension of this argument to twisted mass lattice QCD is not trivial, since there are also disconnected contributions to the neutral pion itself. Near the continuum limit, however, the argument goes through as before.
In detail, we use fits of the form to the zero momentum data for the connected neutral correlator in the -range 10-23 to determine the mass and the coupling . Note that this mass corresponds to a neutral pion only if it is interpreted as a (partially-quenched or mixed-action) Osterwalder-Seiler state ; we shall call it the ‘connected pion’. We then use our best fit parameters of and to construct the connected pion correlator, and we assign it with an error using the bootstrap method. This constructed connected pion correlator is finally used together with the disconnected contribution to build the full correlator, which we use even below . As shown in fig. 2, the effective mass is now essentially flat. Hence our result is consistent with dominance by a single state down to surprisingly small -values, where the data have small statistical errors. In contrast, fig. 1 shows the effective mass obtained without this excited state removal assumption, and it is less flat and has larger errors, although the plateau value is consistent.
With the excited state removal assumption, we can then make one state fits to the resulting matrix of correlators with -range 2-10 and the results are reported in table 3. This shows that we achieve a significant reduction in the error on the mass. This reduction comes at the cost of a possible systematic error coming from the assumption we have made.
4.2 Point to point correlators
We discussed above the possibility of using correlators at smaller to reduce errors. Another approach, which we now discuss, is to reduce errors at all -values by focussing on a different quantity. The method we describe has been used recently  in a study of the disconnected contributions to the meson.
It is mandatory to use an all-to-all method to estimate the disconnected loops at all and . The disconnected correlator of interest reads
assuming that has no vacuum expectation value, as is the case for the flavour-singlet pseudoscalar meson because of the combined flavour-parity symmetry of TMQCD. The approach we have considered above is to sum over the spatial coordinates and . This has the advantage that the zero momentum correlator is studied which has the simplest theoretical interpretation: as a sum of exponentials for each state. Now is actually peaked when and is small at large . Moreover we find that the absolute error is largely independent of , thus the relative error is smallest when . We illustrate the dependence of on in fig. 3. The zero momentum contribution is a sum over this -distribution with weight approximately and this enhances the larger region which has larger relative statistical errors. Basically the zero-momentum sum picks up (inherent gauge-time variation) noise by summing over large where the signal is unimportant but the noise is still significant.
For example, we find for ensemble with that the relative error on when () is 50% of the relative error on summed over both and . Thus the peak signal to noise is twice as big as that summed over relative spatial separation. It might be thought that taking would give an even smaller relative error, since 6 spatial directions are averaged (for spatial separations on axis). Although we find less error reduction than for , we do make use of a further small error reduction and we select (on axis).
We now discuss how to extract the conventional physics (i.e. mass values) from correlators at fixed (rather than summed over all ). Consider a correlator corresponding to a meson of mass . In the region of where the ground state contribution dominates, we expect the correlator to be dependent explicitly on the meson mass when point operators are used at source and sink. Indeed one way to extract this behaviour is by evaluating the 4-dimensional lattice Fourier transform:
to model the data, where, to respect the periodic boundary conditions (for a meson), and . In eq. 6, we take which corresponds to the most local lattice realisation of the derivative in the Klein Gordon action. Other implementations of the derivative would shift the mass by corrections of order as discussed long ago . We have checked that this expression correctly reproduces the -dependence for the charged pion correlator in the -region where the ground state dominates, using the mass value previously determined from fitting the zero momentum correlator. Thus, though less familiar, one can use the expression of eq. 6 to extract the mass value from the -dependence at fixed , in place of the usual expression used at zero momentum.
To exploit this error reduction for the , we combine the disconnected contribution (with as discussed above) with the neutral pion connected contribution also evaluated with . This we evaluate using point sources. However, it is optimal, as discussed in the preceding section, to use only the ground state contribution to the connected neutral pion correlator, as we now discuss.
For this ground-state connected component, we use eq. 6 with as input the mass and coupling determined from the zero-momentum fit to the neutral (connected) pion. The dependence of the resulting correlator on (at fixed ) can then be parametrised again by eq. 6. One way to illustrate this, as shown in fig. 2, is by plotting an ‘effective mass’ (here defined as the mass parameter that solves the dependence given by eq. 6 for two adjacent -values). This shows that we again have a rather constant ‘effective mass’ even from low -values, consistent with a description by one state only. Moreover, the errors on the ‘effective mass’ are smaller than from the fixed momentum methods used above.
To determine the mass by this method, we then fit the resulting combined correlator to a two parameter expression (free mass and coupling) according to eq. 6. Since this approach is less familiar, we illustrate in fig. 4 the fit versus and also in fig. 3 the -dependence of the disconnected contribution compared to the measured data. In detail, we use fits to the zero momentum data for the connected neutral correlator for the -range 10-23 to determine the connected pion mass and the coupling , as above. From the bootstrap ensemble of values of and , we determine, at each -value, the point-to-point connected pion correlator at fixed with its associated error. This is then used together with the measured point-to-point disconnected contribution to construct the full correlator. The result is then fitted with one state using eq. 6 and a -range 2-10 to obtain the mass value given in table 4.
Using different -ranges gives the same mass value, within errors. Compared to using the zero-momentum approach described previously, the statistical error on the mass is much smaller. Here we have used local source and sink operators (since in this case the fixed behaviour is given without extra parameters). What is especially helpful, however, is that the resultant data is well fitted by only one contribution for .
We have presented this new and powerful method for one ensemble () but we also apply it (using ) to the other cases considered and the results are collected in table 4.
This point-to-point approach provides a useful reduction in the statistical error, though at the potential cost of introducing a source of systematic error in relating the quantity measured to that actually required. For local meson creation and destruction operators, the -dependence at fixed is given by the same parameters as the conventional zero-momentum case, up to possible non rotation-invariant contributions of order which we expect to be insignificant at large . For spatially smeared (or fuzzed) operators this invariance is not present and we would need more parameters to describe the point-to-point correlators. We also find strong evidence of ground state dominance, which enables us to use the more precise data from smaller -values.
5 Summary of mass
In order to compare the error-reduced results from different lattice spacings, we plot them using the -values given in table 1 to create dimensionless quantities. We see in fig. 5(a) that the results are consistent with each other and are also consistent with the results presented above using the standard fixed-momentum analysis technique. This increases the statistical impact of our results.
We also plot previous results [27, 15, 28] that are in the quark mass and lattice spacing region we are exploring in fig. 5(b). Different groups using different lattice formulations obtain consistent results. Results from our two lattice spacings are consistent with each other and we are able to evaluate the flavour singlet mass closer to the continuum limit (and in an order improved formalism since we work at maximal twist) than hitherto.
An extrapolation to the chiral limit of our results gives a value of the mass of assuming a linear dependence. From our lattice results, we actually see no statistically significant evidence for any slope, and if we assume a constant behaviour, we would obtain . Using fm obtained  from the ETMC evaluation of , we obtain an mass of 0.865(65) GeV where the error reflects that from the linear extrapolation. This updates the phenomenological estimate of 0.776 GeV discussed above, which was in turn mainly motivated by lattice results at heavier quark masses than those we now have available. We have additional systematic errors which are not fully under control: from the chiral extrapolation, the excited state removal assumption, the continuum extrapolation and any possible volume dependence. We do not see any sign that these systematic errors are required to our current statistical precision, but, as a conservative estimate, we consider these systematic errors to be comparable to the statistical error. So we quote an mass of 0.865(65)(65) GeV.
One of the main conclusions of our study is that we find strong evidence that the flavour singlet mass remains finite as the quark mass is reduced to the chiral limit. This has implications for the topological charge: it implies that the topological charge susceptibility must decrease as in the chiral limit . We will discuss our results for topological charge more fully elsewhere .
We thank Jaume Carbonell and Mariane Brinet for assistance with making point source propagators available for testing point to point propagation. We thank all members of ETMC for useful comments, help and a fruitful collaboration. One of the authors (CM) acknowledges a useful discussion with Paul Rakow. The computer time for this project was made available to us by the John von Neumann-Institute for Computing on the JUMP and Jubl systems in Jülich and apeNEXT system in Zeuthen, by UKQCD on the QCDOC machine at Edinburgh, by INFN on the apeNEXT systems in Rome, by BSC on MareNostrum in Barcelona (www.bsc.es), by the NW Development cluster at Liverpool and by the Leibniz Computer centre in Munich on the Altix system. We thank these computer centres and their staff for all technical advice and help. On QCDOC we have made use of Chroma  and BAGEL  software and we thank members of UKQCD for assistance.
This work has been supported in part by the DFG Sonderforschungsbereich/ Transregio SFB/TR9-03 and the EU Integrated Infrastructure Initiative Hadron Physics (I3HP) under contract RII3-CT-2004-506078. We also thank the DEISA Consortium (co-funded by the EU, FP6 project 508830), for support within the DEISA Extreme Computing Initiative (www.deisa.org).
-  ETMC, in preparation (2007).
-  T. Feldmann, Int. J. Mod. Phys. A15, 159 (2000), [hep-ph/9907491].
-  UKQCD, C. McNeile and C. Michael, Phys. Lett. B491, 123 (2000), [hep-lat/0006020].
-  P. Weisz, Nucl. Phys. B212, 1 (1983).
-  ALPHA, R. Frezzotti, P. A. Grassi, S. Sint and P. Weisz, JHEP 08, 058 (2001), [hep-lat/0101001].
-  R. Frezzotti and G. C. Rossi, JHEP 08, 007 (2004), [hep-lat/0306014].
-  ETMC, P. Boucaud et al., Phys. Lett. B650, 304 (2007), [hep-lat/0701012].
-  L F , K. Jansen, A. Shindler, C. Urbach and I. Wetzorke, Phys. Lett. B586, 432 (2004), [hep-lat/0312013].
-  L F , K. Jansen, M. Papinutto, A. Shindler, C. Urbach and I. Wetzorke, JHEP 09, 071 (2005), [hep-lat/0507010].
-  ETMC, C. Urbach, PoS LAT2007, 022 (2006), [arXiv:0710.1517].
-  ETMC, P. Dimopoulos, R. Frezzotti, G. Herdoiza, C. Urbach and U. Wenger, PoS LAT2007, 102 (2007), [arXiv:0710.2498].
-  ETMC, B. Blossier et al., JHEP 04, 020 (2008), [arXiv:0709.4574].
-  ETMC, C. Alexandrou et al., arXiv:0803.3190 [hep-lat].
-  L F , K. Jansen et al., Phys. Lett. B624, 334 (2005), [hep-lat/0507032].
-  UKQCD, C. R. Allton et al., Phys. Rev. D70, 014501 (2004), [hep-lat/0403007].
-  A. Shindler, arXiv:0707.4093 [hep-lat].
-  ETMC, C. Michael and C. Urbach, PoS LAT2007, 122 (2007), [arXiv:0709.4564].
-  ETMC, P. Boucaud et al., arXiv:0803.0224 [hep-lat].
-  UKQCD, M. Foster and C. Michael, Phys. Rev. D59, 074503 (1999), [hep-lat/9810021].
-  UKQCD, C. McNeile and C. Michael, Phys. Rev. D73, 074506 (2006), [hep-lat/0603007].
-  UKQCD, P. Lacock, A. McKerrell, C. Michael, I. M. Stopher and P. W. Stephenson, Phys. Rev. D51, 6403 (1995), [hep-lat/9412079].
-  UKQCD, C. McNeile and C. Michael, Phys. Rev. D63, 114503 (2001), [hep-lat/0010019].
-  E. B. Gregory, A. C. Irving, C. M. Richards and C. McNeile, Phys. Rev. D77, 065019 (2008), [arXiv:0709.4224].
-  H. Neff, N. Eicker, T. Lippert, J. W. Negele and K. Schilling, Phys. Rev. D64, 114509 (2001), [hep-lat/0106016].
-  C. DeTar and L. Levkova, PoS LAT2007, 116 (2007), [arXiv:0710.1322 [hep-lat]].
-  C. Michael and I. Teasdale, Nucl. Phys. B215, 433 (1983).
-  CP-PACS, V. I. Lesk et al., Phys. Rev. D67, 074503 (2003), [hep-lat/0211040].
-  K. Hashimoto and T. Izubuchi, arXiv:0803.0186 [hep-lat].
-  H. Leutwyler and A. Smilga, Phys. Rev. D46, 5607 (1992).
-  SciDAC, R. G. Edwards and B. Joo, Nucl. Phys. Proc. Suppl. 140, 832 (2005), [hep-lat/0409003].
-  P. Boyle, http://www.ph.ed.ac.uk/~paboyle/bagel/Bagel.html.