{'fraction_non_alphanumeric': 0.03407540514458753, 'fraction_numerical': 0.01627233702705401, 'mean_word_length': 4.47894257064722, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A study of the eta and eta' mesons with improved staggered fermions

We report on a high statistics lattice QCD calculation of the mass of the eta and eta' mesons using ASQTAD improved staggered fermions. The calculation used two ensembles with different lattice spacings and pion masses. We also report results for eta-eta' mixing. The results are in satisfactory agreement with other lattice calculations using other fermion formulations and with experiment, given the the unphysical quark masses used. We see no evidence of abnormal behaviour at the lattice spacings studied.


I. INTRODUCTION
The explanation for the large mass of the η ′ meson relative to the masses of the light pseudo-scalar mesons involves nonperturbative physics, such as topology and the U (1) anomaly [1][2][3], so these masses present an interesting challenge for lattice QCD. There are many phenomenological applications that require knowledge of the properties of the η and η ′ mesons, such as their decay constants. For example, Di Donato et al. [4] have recently noted that a theoretical description of η-η ′ mixing is required for understanding the CP asymmetries of charmed and bottom mesons with η or η ′ in the final state. Fleischer et al [5] have recently studied the dependence of η-η ′ mixing on the decay B 0 s,d → J/ψη (′) , and which could be used to look for physics beyond the standard model. There is currently an active experimental program [6] into the properties of the η or η ′ mesons at experiments such as WASA [7], KLOE-2 [8], and MAMI [9]. Lattice QCD can in principle compute the η-η ′ mixing from first principles, so it should be able to help interpret the experiments.
The improved staggered fermion formulation, known as ASQTAD fermions [10,11], has resulted in some very accurate lattice QCD calculations which have been validated against experiment with a high degree of precision [12,13]. The results from the improved staggered fermions physics program have been recently been reviewed by the MILC collaboration [14].
Despite the success of the phenomenology from staggered fermions, concerns have been expressed over the theoretical basis of the formalism. The central issue arises from the so-called 'rooting' of the staggered determinant [15,16]. This is invoked in order to correct * Current address: Bergische Universität Wuppertal, Gaussstr. 20, D-42119 Wuppertal, Germany the number of flavors/tastes of sea quarks which would otherwise arise. These aspects of the staggered fermions formalism have been reviewed by Dürr [17], Sharpe [18], Creutz [19] and Kronfeld [20]. The physics of the η ′ meson is sensitive to the topology of QCD. Since some of the concerns about rooting the staggered formalism are related to the topology [16,19], the η ′ could be a place where incorrect results are obtained. The properties of the η ′ meson are also thought to be related to the eigenvalue spectrum of the Dirac operator. The eigenvalue spectrum of a type of improved staggered fermions called HISQ has recently been shown to agree with continuum expectations [21][22][23][24], which is suggestive that the physics of the η ′ meson is correct with improved staggered fermions. The recent paper by Donald et al [23] also contains a rebuttal of some previously raised theoretical concerns. A calculation of the properties of flavor singlet mesons is a further and crucial test of the validity of the improved staggered formalism.
The very successful physics program that uses improved staggered fermions has largely been independent of the properties of the η and η ′ mesons [14]. The MILC collaboration has recently used the properties of the η ′ in their study of the topological susceptibility [25]. The η and η ′ mesons could indirectly affect other lattice calculations. For example, the η and η ′ mesons are common decay products of hadrons that decay via the strong force, such as the a 0 meson [26,27] or the proposed exotic hybrid mesons [28]. However the majority of lattice QCD calculations performed by the MILC collaboration are probably independent of the properties of the η and η ′ mesons.
In [29] we reported on methods to compute disconnected diagrams with improved staggered fermions. In that study we found that we needed much higher statistics than those of a typical ensemble generated by the MILC collaboration (at that time). This new work is based on generating many more configurations in order to improve the statistics. We have recently reported [30] results for the masses of the 2 ++ , 0 ++ and 0 −+ glueballs obtained as a byproduct of this calculation.
Lattice QCD calculations of flavor singlet quantities are not as well developed as those of flavor nonsinglet quantities, because they are computationally more expensive. There have been lattice QCD calculations of the η ′ mesons with 2+1 flavors of sea quarks, by the JLQCD collaboration [31], RBC/UKQCD collaboration [32], Hadron Spectrum Collaboration [33], and TWQCD and JLQCD collaborations [34]. A preliminary study of the form factors for the semileptonic decay of D s → η s has been carried out [35]. There is also a preliminary calculation of the masses of η and η ′ mesons, with 2+1+1 flavors of sea quarks, from the ETM collaboration [36].
There have been a few lattice QCD calculations of the flavor singlet pseudoscalar meson with n f =2 sea quarks. However, the strange quarks in the sea play an essential role in lattice QCD calculations of the mass of the η ′ meson. For example Jansen et al. [37] found the ground state of the flavor singlet pseudoscalar meson in n f = 2 to be 0.865(65)(65) GeV compared with the experimental value of the η mass of 0.548 GeV which lattice QCD calculations with 2+1 sea quarks should reproduce.

II. DETAILS OF THE LATTICE CALCULATIONS
We used the improved staggered fermion action called ASQTAD [10,11] with the tadpole improved Symanzik gauge action. We generated two ensembles at two lattice spacings, using QCDOC machines [38]. The basic parameters of the lattice QCD calculations are in Table I. We generated N (total) cfg configurations, but base our final calculation on a subset of N (meas) cfg configurations. The RHMC algorithm was used to generate the configurations [39,40]. The details of the tuning of the RHMC algorithm are in [41].
The use of new measurement techniques has improved the precision of many lattice QCD measurements. However, for many applications of lattice QCD the use of higher statistics is still an essential requirement. In particular, this is so for glueball studies [42], or for anything which includes the calculation of disconnected diagrams.
Our aim was to generate 30,000 trajectories. For comparison the ETM collaboration has typically used ensemble sizes of 5000 or 10000 trajectories [37]. The length of simulation required also depends on the autocorrelation time. There have been some recent high statistics lattice QCD calculations using an anisotropic lattice [43].
In our methods paper [29] on the η ′ we found that the correlators did not have a Gaussian distribution. A similar observation has been made by the ETM collaboration [37] in their calculation with n f = 2 flavors of sea quarks. In this calculation we again observe a long tail in the distribution of the disconnected correlators. In [29] we pointed out that the long-tailed distribution is an inherent challenge to measuring disconnected correlators. Under the assumption (well-justified by numerical observation) that the pseudoscalar fermion loop operators fluctuate within a Gaussian distribution about zero, the product of two loop operators, measurements disconnected correlator d must fall in a long-tailed distribution, specified by where K 0 is a modified Bessel function of the second kind, and A, B and C are constants. The distribution of disconnected correlator measurements will always have its peak at zero, with the signal coming from the asymmetry of the long tails. This behavior may, however, be masked by any binning of measurements. We show some histograms of disconnected correlator measurements in Fig. 1.
The UKQCD/RBC collaboration has, however, not observed any "large statistical excursions" in the distribution of the disconnected correlators in an unquenched calculation with 2+1 flavors of sea quarks [32].
It might therefore be tempting to suggest that the long tail in the distribution of the disconnected correlators is caused by the lack of 'lattice' chiral symmetry in the staggered and twisted mass fermion operators, compared to the domain wall formalism. However, we explicitly checked that we used sufficient numbers of noise sources so that our errors were dominated by the gauge noise. The wall sources used by the RBC/UKQCD collaboration [32] can be thought of as using a single noise source but with partitioning in time. It might be that the distribution of the errors on disconnected correlators in the latter case are not dominated by the gauge noise. This issue is discussed by Wilcox [44] and there is a recent study by Alexandrou et al. [45]. Recently Endres et al. [46] have discussed heavy-tailed correlator distributions in condensed matter simulations and lattice QCD calculations [47], and started to develop new techniques to extract physics in such situations.
We note that the width of the disconnected correlator distribution increases with fuzzing of the sources. The taste-singlet pseudoscalar operator in the staggered formulation has the quark and anti-quark separated by four links. It is possible that the introduction of the gauge links in the covariant shift or in the fuzzing increases the width of the Gaussian distribution of the loop operators and hence the width of distribution of the disconnected correlator described by (1). In other formulations, the distributions may possibly be narrower and less troublesome.
It can be difficult to make a good choice of strange quark mass to use in a lattice QCD calculation, although the choice is getting easier with increased experience of 2+1 calculations. In the original calculations by the MILC collaboration, the estimate of the strange quark mass turned out to be off by as much as 25% [14]. The MILC collaboration then corrected for the mismatch of the strange quark mass by using their extensive partially quenched data sets and later unquenched runs. Their updated estimates for the strange quark mass (in lattice units) were 0.035(7) and 0.0261(7) on the coarse and fine ensembles respectively [14]. However we have used the new value of m s that MILC used on the coarse [48], so as to interpolate with their original value of the strange quark mass. Indeed, our coarse ensemble has already been used in an analysis of the strangeness content of the nucleon, by the MILC collaboration [49].
We first discuss the interpolating operators and twopoint functions for flavor singlet pseudoscalar mesons using Wilson fermions with n f degenerate quarks. A longer discussion can be found in our original paper [29]. Twopoint functions are constructed using the interpolating operators Using standard Wick contractions, the two-point function decomposes into a connected (C) and disconnected (D) part In the staggered formalism the native four tastes of sea quarks are reduced to n degenerate flavors through the n 4 -rooting of the fermion matrix. The four valence tastes manifest themselves in valence loops, which come in too large by a factor of four. A connected correlator is essentially a single valence loop and requires normalizing by a factor of 1 4 to get the single-flavor contribution. A disconnected correlator contains two valence loops (along with an arbitrary number of sea quark loops). Computing the single-flavor disconnected contribution therefore requires normalizing by two factors of 1 4 . This factor is sometimes referred to as 'valence rooting' [18]. We include these factors implicitly in our definition of D and C.
In the staggered formulation we use the Kluberg-Stern notation [50]. The Goldstone pion operators are (γ 5 ×γ 5 ), and the flavor singlet pion operators are (γ 5 × 1), where the first index refers to spin and the second index to taste.
The signal-to-noise ratio for the flavor singlet pseudoscalar mesons rapidly falls with time, so it is important to use spatially smeared operators to project onto the ground and first excited state as soon as possible. We used the fuzzing technique [51] with a fuzzing length of 4 (this must be even to respect the symmetries of the staggered action). We computed a variational correlator matrix C using four basis states -light and strange quark pseudoscalar correlators combined with local and fuzzed sinks and sources: In Figs. 2 and 3 we show the effective mass plots for the local and fuzzed smeared operators for the (γ 5 × 1) pion operator. For the main η-η ′ analysis, we use the factorizing fit model [52] shown in Eq. (6), where i, j run from 1 to 4. Here N is finite and all states are treated as stable. In final fits we used no Bayesian constraints. Physics can also be extracted from the measured variational matrix using the generalized eigenvalue method [53][54][55].
where v n (t, t 0 ) are the generalized eigenvectors. Typically a small time t 0 is chosen. The energies are extracted from the generalized eigenvalues We use the effective mass from Eq. (8) as a check on the masses from the fits to Eq. (6). We used the recently updated value of r 0 =0.4661(38) fm [56] to determine the lattice spacing from our measurements of r 0 /a [30]. The inverse lattice spacing determined from r 0 is 1.61 and 2.14 GeV for the coarse and fine ensembles respectively. With our estimates of the lattice spacing we find the value of the Goldstone pion mass on the coarse and fine ensemble to be 295 and 349 MeV respectively.
The mass of the connected strange-strange pseudoscalar meson m ηs = 0.6858(40) GeV, has recently been determined by the HPQCD collaboration [56]. By using m 2 ηs ∝ m s we find that the mass of the strange quark is mistuned by -20% on the coarse ensemble and +5% on the fine ensemble.
The masses in table II can be used to check the reduction of taste breaking as the lattice spacing is reduced [26]. Defining we find d fine /d coarse = 0.39 (1). The MILC collaboration obtained [26] d fine /d coarse = 0.38(3) at similar parameters. Note that MILC use r 1 in equation 9, but the conversion between r 1 and r 0 drops out in the ratio.

III. ANALYSIS OF DISCONNECTED LOOPS
In [29] we reported on methods to compute disconnected diagrams with improved staggered fermions. We found that the technique proposed by Kilcup and Venkataraman [57] was the most efficient one of those we tested for computing the disconnected diagrams of (γ 5 × 1) and (1 × 1) operators. We therefore used that with 32 Gaussian noise vectors in this calculation.
The integrated autocorrelation time for the light and strange (γ 5 × 1) loops are 42 (7)  An important test of the staggered formalism is to check the long time behavior of the D(t)/C(t) ratio of the correlators in Eq. (4). In quenched QCD, the D(t)/C(t) ratio rises linearly with time but, in unquenched QCD, the ratio should tend to 1 for large times. This test was performed for clover fermions in [58]. In our earlier study [29], the statistics were not large enough to determine the large time behavior of the ratio in unquenched QCD. In Fig. 4 we plot the D(t)/C(t) ratio for the coarse and fine ensembles respectively, as well as a quenched ensemble with the same lattice spacing as the coarse ensemble.
With the higher statistics afforded by the two new ensembles there is now a clearer difference between the quenched and unquenched D(t)/C(t) ratios. The D(t)/C(t) ratio for the two unquenched seem to be reaching a plateau close to 1. As explained in Sec. II, the staggered ratio depends on the 'valence rooting' factors of 1/4, so this is a good test of the staggered formalism. We also include the result from a quenched lattice QCD calculation [29] with a lattice spacing close to the coarse ensemble.

B. Fit results for the masses
We plot the fit results against the effective masses from the variational method (8) in Figs. 5 and 6. Here, we set ∆ = t − t 0 = 2 in (8). The parameters from the final factorising fits are in given Table III. In arriving at these, we varied t min and N exp looking for a sweet spot of stability with respect to these, for χ 2 /dof < ∼ 1, and for error bars small enough to distinguish three states. Plateaux in t min were small due to the competition between excited states at small t and the onset of large noise/signal at larger t. An SVD cut of around 10 −5 was necessary to get fits with reasonable confidence levels.
Note that the statistical errors on the final mass parameters deduced from these correlated fits are generally somewhat less than the individual variational effective masses calculated at each value of t, even though the same number of configurations were used. The global factorizing fits include all the data from the different flavor and fuzzing channels simultaneously, and from a large range of t-values. They also take correlations into account. The individual effective mass estimates from (8), on the other hand, are each determined by relatively small samples of this data set -uncorrelated measurements from two time-slices yield independent estimates for each eigenvalue.
We found it necessary to restrict t max to moderate values (16 and 20) since, at larger values, the correlators can fluctuate below zero. We do not attribute this large-t behavior to unphysical effects. Rather, the unique statistical properties of disconnected correlators mean that a small number of configurations at the tail of the distribution can cause large fluctuations of the mean, as discussed in [29].
We associate the lowest state with η and the first excited state with the η ′ .

IV. COMPARISON WITH EXPERIMENT AND OTHER LATTICE CALCULATIONS
Given that we only have two masses at two different lattice spacings, we cannot make a controlled continuum  We also include data from the UKQCD/RBC [32], Hadron Spectrum Collaboration [33]. The state listed as η(3rd) is the next excited state in our fits. or chiral extrapolation. However, we do wish to check that the resulting parameters are physically reasonable.
In lattice QCD calculations with two flavors of sea quarks, it was found that the mass of the η 2 light flavor singlet pseudoscalar meson had very little mass dependence [37].
In Fig. 7 we plot our results for the masses of the η and η ′ mesons as a function of the square of the pion mass. We also include the recent results from the UKQCD/RBC collaboration [32] and the hadron spectrum collaboration [33]. Note that while the lattice spacing is larger for the coarse ensemble than the fine ensembles, the light quark mass is smaller for the coarse ensemble than the fine.
The mass from the (γ 5 , 1) pion operator is different from the Goldstone (γ 5 , γ 5 ) pion operator, because of taste breaking. In Fig. 8 we plot, as a function of the square of the lattice spacing, the masses of the η and η ′ mesons along with the mass splitting between the (γ 5 , 1) and (γ 5 , γ 5 ) pions.
In a quantum mechanics analysis of the mixing of the η and η ′ [59,60], the starting point is two bare flavor singlet pseudoscalar mesons. As interactions are introduced, the two flavor singlet mesons mix and a mixing angle is introduced. This mixing angle is a parameter in the various decays of the η and η ′ mesons.
In unquenched lattice QCD calculations it is not clear how to start from unmixed states. It is possible to estimate the elements of the mixing matrix in quenched and partially quenched QCD [61,62]. We also note that the properties of the connected strange-strange pseudoscalar meson determined by the HPQCD collaboration [56] are also useful for developing a partially quenched mixing analysis [61].
In QCD, the mixing between the η and η ′ mesons is understood at the amplitude level. In phenomenological analysis of decays, the mixing is parametrized in terms of decay constants and 1 or 2 angles. Many of the analyses of the decays also include additional strong interaction physics, such as form factors [63].
Di Donato et al. [4] and [5] Fleischer et al. have recently discussed the importance of the mixing of η and η ′ mesons on various decays. Shore and collaborators [2] have derived equations relating the masses and decay constants of the η and η ′ to the topological susceptibility.
The phenomenological analysis of η-η ′ mixing is usually described in terms of the axial decay constants [1]. To study the decay constants of flavor singlet meson operators with staggered fermions requires the (γ 4 γ 5 , 1) pion operators. These were not included in our calculation. Since there is no Ward identity for the (γ 4 γ 5 , 1) operators, renormalization factors would need to be computed [64]. The connection between η − η ′ mixing with axial currents or pseudoscalar currents is discussed in section 3.8 of Feldmann's review [1]. To first order the mixing angle is the same with axial or pseudoscalar currents.
There are two popular bases in which to express η/η ′ states. One basis is the SU (3) basis: The η and η ′ mesons states are then linear combinations of the basis states: The η and η ′ mesons can also be described using the quark flavor basis: The two mixing angles are related to each other: In lattice QCD calculations it is not possible to start from pure basis states such as |η q and |η s . Instead the mixing is described in terms of amplitudes. For example we can define an amplitude a η8η via Other amplitudes, for different interpolating operators and mesons, are similarly defined. For example, the amplitudes a ηqη and a ηsη are parameters in the fit model in Eq. (6), with numerical indices replaced by the names of the meson or interpolating operator. With two possible interpolating operators and two states, there are in principle four possible decay constants. The discussion of mixing in quantum mechanics suggests that these four amplitudes can be parameterized in terms of two decay constants and a mixing angle.
It has been found that, with SU (3)-breaking, there cannot be one single mixing angle in the SU (3) basis [65]. For example from a review of the literature, Feldmann [1] notes a spread of 10 degrees in the mixing angles extracted from experiment. In the SU (3) flavor basis the mixing of decay constants is expressed by The decay constants can also be computed in the quark basis: To extract the mixing angles from the matrix in Eq. (17), the following combinations can be used: There are arguments [1,66] that suggest φ q ≈ φ s . Feldmann [1] reviews the determination of the mixing angles from various processes and he finds | φ q − φ s |< 5 degrees. The UKQCD/RBC colaborations use perturbation theory to identify conditions where there is only one angle [32].
In Fig. 9 we show the fit amplitudes plotted in couplets (a η , a η ′ ) demonstrating the mixing angles in the quark flavor basis as in Eq. (17).
Because of SU (3) symmetry breaking, f q is not equal to f s in Eq. (17), hence a rotation of the interpolating basis from the quark to the flavor basis will no longer give a matrix of the same form as Eq. (17). The use of mixing angles to parametrize the mixing matrix is similar to the polar decomposition of a matrix [67], where a matrix can be written as the product of a Hermitian matrix and a unitary matrix. In this application the Hermitian matrix is actually diagonal.
The RBC/UKQCD and Hadron Spectrum Collaboration extract a single angle by computing the following angle so that tan φ est is the geometric mean of tan φ est q and tan φ est s . The relation φ q ≈ φ s is an assumption, under which the four possible parameters are reduced to three, hence it is interesting to test it using lattice QCD. We have determined the mixing angles using different methods and present the results in Table IV. The angles φ est q , φ est s and φ est were obtained from the fits to Eq. (6), using Eqs. (18), (19) and (20). Local operators were used for this. We also found similar results using fuzzed operators.
The angle φ fit , in Table IV, was determined by simultaneously fitting the smeared and local matrix of amplitudes in Eq. (17) using a single angle φ fit = φ q = φ s , with two decay constants for the local matrix, and two decay constants for the fuzzed correlators.
One cannot confidently conclude from the results in Table IV, that a single mixing angle does not describe the data. As an exercise we ignored the statistical errors on the amplitudes from the β = 7.095 ensemble and rotated the matrix in Eq. 17 by an angle α in the interpolating operator basis and computed φ est q and φ est s . Around 9 degrees from the quark basis we found φ est q = φ est s as a function of α. This shows as expected that the quark basis is a good basis for having a single mixing angle. This method can be used with future more precise lattice QCD calculations, to determine the optimal basis of interpolating operators, such that the amplitudes can be described by a single mixing angle.
In Fig. 10 we plot the mixing angle in the quark basis (φ fit ), obtained using the fit method, as a function of the square of the pion mass. We also include some experimental numbers for the mixing angle from the summary in [4]. We see that our results from staggered fermions are qualitatively in agreement with the results of the other two lattice groups and experiment.

VI. CONCLUSIONS
We have reported a lattice QCD calculation of the masses and mixing of the η-η ′ mesons. Our key results for the masses (Fig. 8) and mixing angles (Fig. 10), and the ratio plot (Fig. 4 ) do not give any grounds for concern over the validity of the staggered fermion formulation. Our results are qualitatively in agreement with experiment and with lattice analyses using other formulations.
Using a technique specific to staggered fermions [57] we were able to compute the disconnected diagrams at a We also include data from the UKQCD/RBC [32], Hadron Spectrum Collaboration [33], and UKQCD collaborations [61]. The UKQCD/RBC number is quoted in the SU(3) basis, so we have used equation 14 to convert it to the quark basis.
cost roughly 10 times the cost of the connected correlators. However, as always, higher statistics are required for lattice QCD calculations which require disconnected diagrams.
The main shortcoming of the calculation is that we have used only one quark mass at each of two different lattice spacings. It is clearly important to extrapolate the masses of the η and η ′ mesons to the physical quark masses and to take the continuum limit. To do this, it may be better to use a staggered fermion formalism with reduced flavor symmetry breaking. The MILC collaboration [68] and HOT collaboration [69] have started to simulate the HISQ [70] fermion action, which uses an additional level of fat links over the ASQTAD action so reducing the mass splitting between the masses of the (γ 5 × 1) and (γ 5 × γ 5 ) pions.
The stout link staggered action used by Aoki et al. [71,72] could also be used, since this formalism has reduced taste splitting compared to ASQTAD and the quark masses are within 10% of the physical light and strange masses.
For the phenomenology of η and η ′ it is important to study the decay constants and mixing angle of the η and η ′ mesons. Although staggered fermions are computationally cheap and have a Ward identity for the Goldstone pion operator, it is not clear whether we can obtain a strong enough signal for the (γ 4 γ 5 × 1) operator to be able to extract decay constants [64].
Even more challenging would be to study the mix-ing of the the η ′ and η with the pseudoscalar glueball. This requires an extension of the formalism for mixing described in section V. In principle, the mixing angle between the η and η ′ and the pseudoscalar glueball can be obtained by fitting various branching ratios. However, a clear picture from the phenomenological approach has yet to emerge [4,5,73,74], so a first principle calculation from lattice QCD would be valuable. The KLOE experiment [8] has analysed their branching fraction data for vector meson radiative decay to pseudoscalar mesons with a model and obtained an estimate of η ′ -glueball mixing which is three sigma from zero. The KLOE-2 experiment plans to reduce the experimental errors on the analysis [8]. Given that we have already published [30] a study of the pseudo-scalar glueball on these ensembles, we could at the very least include the gluonic and fermionic pseudoscalar interpolating operators in a single variational calculation (as was done for scalar operators in [75]). This calculation would be interesting from a nonperturbative QCD perspective, but also would be important for exploring CP violation in the decays of B and B s mesons to final states that include the η and η ′ .
Lattice QCD calculations of the η and η ′ mesons involve an interesting mixture of: toplogy, anamolies, mix-ing via quark loops, and have important applications to phenomenology. The relatively high cost of lattice QCD calculations of the η and η ′ mesons has meant that they are not as highly developed as the calculation of flavor nonsinglet hadrons [14]. Since the computational costs are becoming less exorbitant, we expect that lattice calculations of the η and η ′ mesons may soon become precision studies.

VII. ACKNOWLEDGMENTS
We thank Steven Miller, Zbyszek Sroczynski for participating in the early part of this project. We thank Chulwoo Jung and Mike Clark for providing optimized code used in this project. We also thank Alistair Hart, Andreas Kronfeld, Konstantin Ottnad and Carsten Urbach for their helpful comments. We used the North West grid, Scotgrid and Bluegene/P at Daresbury. The measurements used the Chroma software system [76]. The configurations were generated on the QCDOC [77]. EBG was supported at UCY by Cyprus Research Promotion Foundation grant ∆IEΘN HΣ/ΣT OXOΣ/0308/07.