Flavor-singlet meson decay constants from $N_f=2+1+1$ twisted mass lattice QCD

We present an improved analysis of our lattice data for the $\eta$--$\eta'$ system, including a correction of the relevant correlation functions for residual topological finite size effects and employing consistent chiral and continuum fits. From this analysis we update our physical results for the masses $M_\eta=557(11)_\mathrm{stat}(03)_{\chi\mathrm{PT}}\,\mathrm{MeV}$ and $M_{\eta'}=911(64)_\mathrm{stat}(03)_{\chi\mathrm{PT}}\,\mathrm{MeV}$, as well as the mixing angle in the quark flavor basis $\phi=38.8(2.2)_\mathrm{stat}(2.4)_{\chi\mathrm{PT}}^\circ$ in excellent agreement with other results from phenomenology. Similarly, we include an analysis for the decay constant parameters, leading to $f_l=125(5)_\mathrm{stat}(6)_{\chi\mathrm{PT}}\,\mathrm{MeV}$ and $f_s=178(4)_\mathrm{stat}(1)_{\chi\mathrm{PT}}\,\mathrm{MeV}$. The second error reflects the uncertainty related to the chiral extrapolation. The data used for this study has been generated on gauge ensembles provided by the European Twisted Mass Collaboration with $N_f=2+1+1$ dynamical flavors of Wilson twisted mass fermions. These ensembles cover a range of pion masses from $220\,\mathrm{MeV}$ to $500\,\mathrm{MeV}$ and three values of the lattice spacing. Combining our data with a prediction from chiral perturbation theory, we give an estimate for the physical $\eta,\eta' \rightarrow \gamma\gamma$ decay widths and the singly-virtual $\eta,\eta'\rightarrow\gamma\gamma^*$ transition form factors in the limit of large momentum transfer.


Introduction
The axial anomaly and the topological nature of Quantum Chromodynamics (QCD) are able to explain the large experimentally observed mass of the η meson. This understanding was possible due to perturbative arguments leading to the Witten-Veneziano [1,2] formula. This was recently confirmed non-perturbatively using lattice simulations in Ref. [3].
Beyond masses, there is also the phenomenon of mixing: η and η mesons are not flavor eigenstates, but represent a mixing of an octet and a singlet state. In contrast to the ω-φ meson mixing in the vector channel, the mixing in the pseudoscalar channel is ideal at the SU(3) symmetric point with mixing angle φ = 54.7 • (in the quark flavor basis), but exhibits significantly smaller φ-values at physical quark masses [4]. The precise knowledge of the mixing angle is important for several phenomenological applications, most notably for improving the theoretical estimate of the hadronic contribution to the anomalous magnetic moment of the muon [5].
In this paper we extend our previous studies [4,16,17] of properties of η and η mesons in two ways: first by an improved analysis: on an enlarged number of Monte Carlo ensembles we perform consistent chiral and continuum extrapolations. This leads to slightly changed results when compared to Ref. [4], mostly within the quoted error bars, with the exception being the mixing angle φ. The larger change in φ comes from the fact that we are now able to resolve the lattice spacing dependence in φ, too, due to more ensembles at the finest lattice spacing value. Moreover, by using derivatives of correlation functions instead of the correlation function themselves, we are able to remove systematic effects from not optimally sampled topological sectors.
Second, for the first time we estimate flavor-singlet pseudoscalar decay constants from lattice QCD. We use these to determine also physical η, η → γγ decay widths and the singly-virtual η, η → γγ * transition form factors in the limit of large momentum transfer. A first account of this work can be found in Ref. [17].
For the determinations of the aforementioned decay constants we rely on chiral perturbation theory instead of determining them directly from flavor-singlet axial-vector matrix elements. The reason for this procedure is an unfavorable signal-to-noise ratio in some of those matrix elements, which prevents a meaningful analysis. The results are compared to phenomenological determinations mostly summarized in Ref. [21].
The results presented here are based on gauge configurations produced by the European Twisted Mass Collaboration (ETMC) with active up/down, strange and charm quarks. Based on three values of the lattice spacing and pion masses ranging from 220 to 500 MeV, the ETMC ensembles allow us to perform a controlled chiral and continuum extrapolation. Dedicated ensembles with varied strange quark masses let us control also the strange quark mass dependence.
In the twisted basis the fermionic action containing a mass-degenerate, light quark doublet χ l = (χ u , χ d ) T reads [27][28][29] S l [χ l , χ l , U ] = a 4 x while for a non-degenerate, heavy quark doublet χ h = (χ c , χ s ) T we have [29,30] where the Pauli matrices τ i , i = 1, 2, 3 act in flavor space. The massless Wilson Dirac operator depends implicitly on the gauge links U . The doublets χ l and χ h are related to doublet fields in the physical basis ψ l , ψ h via chiral rotations. The bare strange and charm quark masses m s , m c are given in terms of the bare input parameters µ δ and µ σ µ c,s = µ σ ± Zµ δ , where Z = Z P /Z S denotes the ratio of pseudoscalar and scalar flavor non-singlet renormalization factors. The renormalized quark masses require an additional factor of non- which is the same as for the light bare quark mass, i.e. µ r l = µ l /Z P .
We employ 17 gauge ensembles as detailed in Tab. 1 at three different values of β corresponding to three different values of the lattice spacing a, c.f. Tab. 2. Compared to previous studies of the η,η -system [3,4], we have added one more ensemble at the finest lattice spacing (D20.48) and significantly increased the statistic on the B55.32 ensemble. In general, all observables have been computed with the full statistic as given in Tab. 1 with exception of the kaon mass M K and the kaon decay constant f K , that have been computed only on a subset of configurations in many cases. However, this is sufficient for our purposes as we are neither interested in M K nor f K themselves in this study. The resulting errors for derived observables are always dominated by the statistical uncertainties in the flavor-singlet sector, anyways.
Tab. 2 also contains the results for Z at each value of β, which are needed for the computation of matrix elements. The labels M1 and M2 refer to the two different methods used in Ref. [31] for the determination of renormalization factors. In addition, we included the values for the Sommer scale r 0 at each value of the lattice spacing, r 0 /a, that were taken again from Ref. [31] together with the physical value which is required to set the scale in our study. Note that in an earlier publication in Ref. [4] slightly different values have been used for r 0 and r 0 /a. However, those are now superseded by the values from the final analysis in Ref. [  While the values of aµ σ and aµ δ defining the bare strange and charm masses are generally fixed for each choice of β, we include a few dedicated ensembles (A80.24s, A100.24.s and D45.32sc), which have been generated with different µ δ -and in case of D45.32sc -also different µ σ values. This allows us to explicitly resolve the dependence on the strange quark mass and obtain more stable results from chiral fits.
For the computation of quark-disconnected diagrams we employ stochastic volume sources. The corresponding number of stochastic samples N S is included in Tab. 1 and is chosen such that the final statistical errors are dominated by gauge noise. The statistical errors for all observables are computed using the blocked bootstrap with 10000 samples and blocklengths chosen such that the effective length of every block corresponds to at least 20 HMC trajectories. This has been found sufficient to deal with autocorrelations in an earlier study in Ref. [16].

Computation of masses and amplitudes
The extraction of masses and matrix elements for the η,η -system has already been discussed in detail in previous publications [3,4,16,17]. In the following we will first briefly summarize the relevant methods and then introduce a modification leading to systematic improvement of our existing analysis. This improvement concerns a possible contamination of the large-t behavior of the flavor-singlet correlations functions induced by imperfectly sampled topology. Finally, we detail the extraction of mixing parameters and further, derived observables such as decay widths.

Correlation function matrix
In the physical basis we consider the following three local pseudoscalar operators whereψ l (x), ψ l (x) andψ h (x), ψ h (x) denote degenerate light and non-degenerate heavy quark doublets, respectively. While the doublet structure of the fields is required for the rotation to the twisted mass basis, the flavor projector (1 ± τ 3 )/2 disentangles charm ("+") and strange ("−") contributions. Upon rotation to the twisted basis the operators read at maximal twist In Refs. [4,16] it has been shown that the charm quark operator essentially has no overlap with the η,η -states and can hence be neglected. Therefore, we drop it from the actual analysis and keep only light and strange quark operators, i.e. S tm l (x) ≡ S 3,tm l (x) and P tm s (x) ≡ P −,tm h (x). Considering renormalization, the operators can be written as Pulling out a factor Z 2 P from the resulting 2 × 2 correlation function matrix we have whereC contains operators that are projected to zero-momentum and renormalized up to a global factor of flavor non-singlet Z P , i.e.S tm The mixing of flavor non-singlet pseudoscalar and scalar currents in the heavy quark basis remains manifest in the corresponding ratio of renormalization constants Z. Note that the construction of these operators requires only the ratio Z instead of Z P and Z S separately. We will show later that renormalization of the correlation function matrix up to a global factor Z 2 P is sufficient for the calculation of the mixing parameters, as this factor will be absorbed by the renormalization of corresponding factors of quark masses.
As first proposed in Ref. [32] and subsequently used in Refs. [4,9], we replace the quark-connected pieces in the correlation functions by the ground state contribution. This allows to extract the η and η states from the resulting principal correlators from the earliest available timeslice t = 2a on, leading to a significant improvement in the signal-to-noise ratio. For further technical details on this procedure, including fit parameters and results, we refer to appendix 7 and tables therein.

Correlator improvement for topological effects
In Ref. [33] it has been pointed out that the large-t behavior of quark-disconnected contribution C disc 2pt (t) to pseudoscalar flavor-singlet correlation functions in finite volume and for fixed (or imperfectly sampled) topology differs from zero. This has been further investigated numerically in Ref. [18], where the effect on the correlation functions has been shown explicitly for different topological charge sectors. In fact, the leading term in the 1/V expansion contributing to at large-t at fixed topological charge Q behaves as where χ t denotes the topological susceptibility and c 4 the kurtosis of the topological charge distribution. While we do not find deviations from a zero topological charge in the gauge average on any of the ensembles used in this study, it is still to be expected that an imperfectly sampled topological charge distribution leads to deviations from the infinite volume result. Only on very few ensembles we find a shift at the level of single correlation functions, which is not compatible with zero within errors. This is most notably the case for the light quark correlation functionC ll (t) = S tm l (t)S tm l (0) on D45.32, which yields the dominant contribution to the η principal correlator that is shown in the left panel of Fig. 1. This ensemble exhibits the smallest physical volume of all ensembles in this study.
In order to remove any constant shift from our correlation functions we replace the correlators in Eq. (15) by a naive time-derivative, i.e. the difference of two adjacent timeslices before solving the generalized eigenvalue problem (GEVP) where n = η, η denotes the two states for the 2 × 2 problem and using fixed t 0 /a = 1. Taking into account periodic boundary conditions, this changes the functional form of the correlation functions in Eq. (15) from an even cosh-like to an odd sinh-like behavior where A n i (A n j ) denotes the physical amplitudes for the n-th state with respect to the i-th (j-th) element of the basis containing N operators. The asymptotic behavior for the principal correlators then takes the same form which is fitted to the lattice data to extract the energies. Including only correlation functions projected to zero momentum, we have E n = M n , which yields M η and M η for the two lowest states. The information on physical amplitudes can be retrieved from eigenvectors in the standard way [34] Note that amplitudes computed in this way from the correlation function matrix in Eq. (15) are renormalized only up to a factor of Z P , which turns out sufficient for our purposes.
In addition to removing a constant shift in the η principal correlator, the derivative correlator has much smaller point errors compared to the standard approach, which can be seen comparing the left and the right panel of Fig. 1. The data in the right panel has been generated using a time shift of ∆t/a = 1. In fact, the actual choice of ∆t has little impact on the final results from correlated fits within errors, hence we use ∆t/a = 1 throughout our analysis. The only effect of values ∆t/a > 1 is less reduction in correlation between adjacent timeslices inC (t) but in the final fits this is compensated by the behavior of individual point errors.

Extraction of mixing parameters and decay widths
The general definition of meson decay constants employs the axial-vector matrix element in the physical basis where P = π, K, η, η , ... denotes the desired meson state (with momentum p) and the index a is used to distinguish different flavor structures for the axial-vector current.
First we consider the charged meson sector for a light quark doublet imposing exact isospin symmetry. In this case the physical axial-vector current transforms into the vector current in the twisted basis at maximal twist. By virtue of the PCVC relation this leads to a convenient expression for pion decay constant which can be used to compute f π in tmLQCD without the need for any renormalization and to high statistical precision due to the pseudoscalar current [35][36][37]. A similar relation holds in the heavy-light meson sector for the kaon . While there is again no overall renormalization factor needed, the ratio Z is required due to mixing between scalar and pseudoscalar currents and also implicitly in µ s as defined in Eq. (5)).
In the flavor-singlet sector, there exist two popular choices for the basis of two local operators made from degenerate light quark fields and a strange quark field. The first one is the so-called octet-singlet basis which is the preferred basis in the formulation of (χPT). A second choice is the quark flavor basis, defined by In any case, the most general parametrization of the decay constant parameters f a P for two local operators made from degenerate light quark fields and a strange quark field, reads: where a = 8, b = 0 or a = l, b = s, for octet-singlet and quark flavor basis, respectively. For the octet-singlet basis it is found in χPT that the leading contribution to the difference of the two mixing angles |φ 0 − φ 8 | is a purely SU (3) F -breaking effect, while at the same order in the chiral power counting OZI-violating terms contribute only to the flavor-singlet decay constant parameter f 0 [38,39]. Therefore |φ 0 − φ 8 | cannot be expected to be small. On the other hand in the quark flavor basis the corresponding difference φ l − φ s is proportional only to an OZI-violating term. Since in the SU (3) F -symmetric theory, the mixing angles fulfill which has been confirmed numerically in a previous lattice study [4]. Therefore, it is reasonable to define a simplified scheme in the quark flavor basis (the so-called FKS scheme [39]), employing only a single mixing angle, i.e. rewriting Eq. (31) as For a more detailed discussion of the relation between the two schemes, we refer to the review given in [21].
In principle, the l.h.s of Eq. (33) could be computed directly from the lattice. However, we find that axial-vector interpolating operators do not give a sufficient signal in practice. Therefore, in the physical basis of QCD we resort to pseudoscalar matrix elements where again P = η, η and m a for a = l, s denotes the quark mass for the light and strange quark, respectively. The pseudoscalar flavor-singlet operators are the ones given in Eqs. (8,9), i.e. P l = P 0,phys l for light quarks and P s = P −,phys h (x) for the strange component. Applying χPT to the same order as for the splitting of the mixing angles, one obtains the desired relation between the mixing parameters from the axial-vector case in Eq. (33) and the matrix elements h P a [21] h η On the lattice we work in the twisted basis, thus the pseudoscalar currents in the physical basis of QCD need to be replaced by their twisted counterparts. Considering renormalization the matrix elements that are actually computed on the lattice are given by for P = η, η . They are obtained from Eq. (23) solving the GEVP for the correlation function matrix in Eq. (15) and multiplying by a factor of the twisted bare light or strange quark mass. Note that the factor Z P , which would otherwise be needed to renormalize the operatorsS tm l ,P tm s , is canceled by the factor 1/Z P required for the renormalization of the quark masses µ l , µ s . We point out that it is an intrinsic advantage of the twisted-mass formulation that the computation of the flavor-singlet mixing parameters does not require knowledge of the singlet pseudoscalar renormalization factor Z 0 P at all and that even the non-singlet Z P is not required explicitly. The latter is similar to f π , which can be computed in the twisted mass formulation without the need for renormalization, or f K , which involves only the ratio of non-singlet renormalization factors Z = Z P /Z S .
Finally, we note that the mixing angle φ is always invariant under renormalization as it is computed from the (double-) ratio of matrix elements unlike the decay constant parameters f l , f s , that depend on Z.
Using the values for the mixing parameters it is possible to derive estimates and constraints for further physical observables which are driven by the chiral anomaly. First of all, there are relations to the η and η decay widths derived from effective field theory, given by [39] Secondly, effective field theory yields relations for the pseudoscalar transition form factors F ηγγ (q 2 ) and F η γγ (q 2 ) in the limit of large Euclidean momentum transfer Q 2 and the mixing parameters in the quark lim where we have introduced the shorthand notationF P γγ * , P = η, η for later use.

Results
Since our simulations are performed at unphysical values of the quark masses and finite lattice spacing, a chiral extrapolation is required to obtain physical results. To this end we employ a fit ansatz for each observable O inspired by leading order in χPT where we defined  as leading order proxies for the light and strange quark mass, respectively. In the above fit model n is an integer such that r n 0 O is dimensionless and m denotes the power of the observable required for the chiral expansion, i.e. m = 2 for masses and m = 1 for decay constant parameters and the mixing angle φ. The same values of m are used for ratios of the respective quantities, e.g. m = 2 for fitting M η /M K . The first term on the r.h.s of Eq. (43) is included as a free parameter only for observables that do not vanish in the SU (3) F chiral limit and take a non-trivial value, such as M η , decay constants or ratios thereof. For observables with an analytically known value (e.g. M η /M K → 1, φ → arctan √ 2) the parameterO is replaced by the respective value.
The constants L l,s,β are always free parameters and determined from the fit. They are used to perform chiral and continuum extrapolations, as well as to correct our lattice data for unphysical values of the quark masses and possibly lattice artifacts in plots, e.g. in the right panel of Fig. (2). Note that the resulting point errors are highly correlated. The O(a 2 ) term ∼ L β has been included to parametrize the leading lattice artifact.
The actual observables O considered in our fits are listed in Tab. 4 together with the resulting fit parameters and values for χ 2 /dof. Note that for the decay constant parameters f l and f s we have to resort to fitting ratios f l /f π and f s /f k , which cancel most of the lattice artifacts and m s -dependence (in case of f s ) that otherwise prevents reasonably fitting the above model.
The lattice results for the masses M π , M K , M η and M η are listed together with the mixing angle φ in Tab. 3. Additional information on the fits of the asymptotic form in Eq. (22) to the data for the η,η principal correlators can be found in the appendix in Tab. 11. For M η , M η and the mixing angle φ we find that the data are well described by this fit ansatz. In particular the light quark mass dependence for M η ,  For observables with analytically known / trivial value in the chiral limit, the respective parameter (r n 0O ) m has been fixed to this value (values without error) and is not a free parameter in the fit. For the decay constant ratios we include results for both renormalization methods, while all other (invariant) results are for fits to data using Z from M2. The values for the fit parameters of φ are obtained assuming that φ andφ are given in radian measure. Additionally, we include the reduced χ 2 values and degrees of freedom. Errors are statistical only.
M η is mild over the full range of available pion masses, as can be seen in Fig. 2. However, the mass of the η depends strongly on the strange quark mass, which is expected from χPT. Therefore, we consider the ratio with the kaon mass M η /M K , which cancels most of the m s dependence and assign a systematic error from the difference of the two central values. For the remaining observables we assess the uncertainty related to our fitting procedure by performing a cut in the pion mass for the data entering the fits. Including only ensembles with M π < M cut π = 390 MeV reduces the number of data points in the fit from 17 to 11, which is still large enough to obtain reliable fits. We assign a systematic uncertainty to each observable from the difference to the central values from a fit with the aforementioned pion mass cut M cut π = 390 MeV, that should reflect the uncertainty related to the chiral extrapolation. Our final results for the η and η masses at the physical point read: where we have used the experimental values for M π and M K to set the light and strange quark mass to their physical values and the Sommer parameter r 0 in Eq. (7) to set the scale. Both results are in good agreement with experiment and compatible with the previous result in Ref. [4]. The statistical errors are slightly increased compared to the old results, which is due to the additional degrees of freedom in the now fully consistent, combined chiral and continuum fits.
As mentioned before, we have computed the ratio M η /M K to assess the uncertainty of the chiral and continuum extrapolation for the η mass. This ratio has indeed been found to cancel most of the strange quark mass dependence in M η [4,16]. The combined leading order chiral and continuum extrapolation to the physical point using the fit ansatz in Eq.
Plugging in the neutral kaon mass gives M η = 0.554(15) stat , which confirms our result from the direct fit and extrapolation of M η . We point out that we use neutral meson masses (M exp π0 , M exp K0 ) to set the quark masses and define the physical point, as we do not include electromagnetic effects in our simulations. Using charged meson masses leads to an ambiguity of a few MeV, which is still below the statistical uncertainty even for M η . π corrected to physical value of the strange quark mass using LO chiral extrapolation. Physical value from full LO fit (see text). The extrapolations resulting from chiral fits at fixed lattice spacing and for physical strange quark mass are shown as red, pink and gray error bands, corresponding to β-values of 1.90, 1.95 and 2.10, respectively. Blue error band and solid line correspond to continuum limit extrapolation. (b) Same, but data also corrected for continuum limit. The SU(2)-chiral extrapolation band is shown only for continuum limit and at physical strange quark mass.
In Fig. 3 we show results from the improved analysis for the mixing angle φ in the quark flavor basis as defined in Eq. (38). The blue band in both panels represents the chiral extrapolation in M 2 π in the continuum limit and at physical strange quark mass as obtained from the fit model in Eq. 43. While the data in the right panel has been corrected also for the mismatch in the strange quark mass and the continuum limit, the data in the left panel are shown at finite values of the lattice spacing together with an error band from the corresponding chiral extrapolation at fixed lattice spacing. The final result for the mixing angle at physical quark masses and in the continuum reads in excellent agreement with results from phenomenology [5,21,40]. We remark that the value quoted above is lower by about three σ than what we found in Ref. [4]. The reason for this discrepancy is that in Ref. [4] we were not sensitive to lattice artifacts. Due to the improved analysis and the additional ensembles, the a 2 dependence can be resolved now, which is responsible for a 6 • decrease in the central value.
In Fig. 4 we show the lattice data for the decay constant parameters f l and f s for both choices of Z, i.e. the plots in the left and right columns show data for Z from method M1 and M2, respectively. While f l is essentially unaffected by the choice of renormalization, the impact on f S turns out to be very significant. Although formally of O(a 2 ), the difference from the choice of renormalization dominates the results and would lead to substantial systematic uncertainties in any attempt of a chiral and continuum extrapolation. Moreover, f s is rather sensitive to the strange quark mass, as can be inferred from comparing results for ensembles A80.24, A100.24 with their counterparts A80.24s and A100.24s, which have a lighter strange quark mass. This is not surprising, because the matrix element in Eq. (37) that determines f s is directly proportional to µ s . However, µ s as defined in Eq. (5) itself depends explicitly on Z. This seems to enhance the effect on f s of different  choices for Z. In fact, we cannot exclude that even terms of higher order (e.g. a term ∼ a 2 m s ) are numerically large for a chiral and continuum extrapolation of f s , hence it is not reasonable to attempt a leading order fit.
We find that taking ratios instead of fitting f l , f s individually allows to circumvent most of these issues and greatly improve the quality of the fits. In particular, we find that forming the ratio of f s with the kaon decay constant f K leads to a milder dependence on the choice of Z and it prevents extreme lattice artifacts as observed for f s using Z from method M1. This is immediately evident from comparing the plots in Fig. 5, which shows results for the ratio f s /f K for Z from both methods, with the corresponding ones for f s in Fig. 4.
Similarly, taking f l /f π leads to better fits than considering f l itself. Still, we find that values for Z computed from method M2 result in generally smaller cut-off effects compared to method M1, as can be seen from the fitted values fir L β in Tab. 4. Besides, method M2 gives a somewhat better fit for f s /f K . Therefore, we prefer to choose to use Z from method M2 to compute the final results in our analysis.
In Tab. 5 we collect the lattice results for f l , f s , f l /f π and f s /f K from method M2 that enter the final fits. Fig. 6 shows the chiral and continuum fit to f l /f π and f s /f K together with the extrapolated lattice data, which appear to be rather well described by the fit ansatz. The final physical results for the decay constant parameters read where we have used the experimental values f exp π = 130.50 MeV, f exp K = 155.72 MeV to extract f l and f s , respectively [41].
Plugging the physical values for M η , M η , f l , f s and φ into Eqs. (39,40) we can finally compute the η, η → γγ decay widths, leading to Γ phys η→γγ =0.71(9) stat (7) sys keV , The large statistical error for Γ phys η →γγ is dominated by the error on the η mass, which enters to third power in Eq. (40). Fig. 7 shows Γ phys η→γγ and Γ phys η →γγ together with the corresponding results at unphysical quark masses and finite lattice spacing, computed on the individual ensembles. Since Eqs. (39,40) become rigorous only in the chiral limit, it is expected hat their should be corrections at finite quark masses. Indeed, this is clearly observed for Γ phys η→γγ , which scales strongly with m s and exhibits also a residual light quark mass dependence. This might explain why the result differs by more than 2σ from the PDG value Γ exp η→γγ = 0.52(2) keV [41]. The situation is different for Γ phys η →γγ , which is essentially a constant in m l and m s , albeit with larger point errors. In fact, the data is compatible with a constant fit over the entire range in M π , confirming the applicability of the formula at least for the η . The result of such a fit is Γ phys η →γγ = 5.5(1.2) stat keV with χ 2 /dof = 20.5/16. In any case, Γ phys η →γγ is in agreement with the value Γ exp η→γγ = 4.4(2) keV within its large error.
Similarly, for the transition form factors at large momentum transfer in Eqs. (41,42) we obtain The much smaller relative statistical error forF η γγ * is caused by anticorrelation, which leads to cancellation of statistical fluctuations in the sum of the two terms on the r.h.s. of Eq. (42). In a similar way statistical fluctuations are enhanced in the difference of terms forF ηγγ * , while the absolute value of the result is smaller.
Nevertheless, even forF phys ηγγ * the relative statistical precision is better than 10%, whereas the systematic uncertainty due to neglecting higher orders in the chiral fits is clearly dominating. The situation seems slightly better forF η γγ * , but also in this case it is impossible to fully assess the systematics arising due to the use of Eq. 42 in the current setup. Therefore, any further improvement beyond the current precision must be subject to a future, dedicated study of η, η → γγ transition form factors, which ultimately should allow to map out the momentum dependence of the transition form factors as well.

Summary
In this paper we have for the first time presented results for flavor singlet, pseudoscalar decay constants using lattice QCD. Thanks to the gauge ensembles with N f = 2 + 1 + 1 dynamical quarks provided by the ETMC collaboration, we could study the continuum and chiral extrapolations in a controlled way. In particular, dedicated ensembles with varied bare strange quark mass allow to resolve the strange quark mass dependence.
For determining the decay constants f and f s , we had to rely on χPT in order to be able to extract them from pseudoscalar matrix elements. This was necessary, because the axial-vector matrix elements turned out to be too noisy. While masses and mixing angle(s) are independent of any renormalization constants, Z S and Z P are needed for f and f s , which have been determined in two different ways in Ref. [31]. Depending on which way (M1 or M2) we follow we observe large lattice artifacts in particular in f s . This we understand, because in f s the strange quark mass dependence is largest and the strange quark mass at separate β-values is strongest influenced by the way the renormalization constants are implemented. Therefore, we decided to rely on only method M2 with significantly smaller lattice artifacts.
Moreover, it turned out that the chiral extrapolation is most conveniently performed by using ratios f s /f K and f /f π . In those ratios most of the quark mass dependencies cancel out. Our final estimates for f and f s are in very good agreement to existing phenomenological determinations [21]. Having f and f s at hand we can also -again relying on χPT -estimate η, η decay widths. Within the large statistical uncertainties we observe reasonable agreement to the PDG values.
Compared to Ref. [4], we have also updated the results for M η , M η and the mixing angle φ in the quark flavor basis. This update was necessary because we have additional ensembles, more statistics and an improved analysis available. The results are, however, compatible. The somewhat larger difference in φ isas mentioned already before -due to the fact that we can now resolve the lattice spacing dependence.
Currently, we are working on improving the signal obtained from the axial-vector matrix elements. Moreover, we are working on extracting η and η directly at the physical point. These steps should allow us to cross-check the χPT formulae used in this paper and reduce the corresponding systematic uncertainties.   Figure 6: (a) Lattice data for f l /f π which has been corrected for the mismatch of the strange quark mass and lattice artifacts from leading order chiral fits. The physical value from the LO fit is shown together with a solid line and gray error band for the chiral extrapolation in M 2 π at m s = m s,phys and in the continuum limit. (b) Same but for f s /f K . The data in both panels has been generated using Z from method M2. Errors are statistical only and highly correlated due to the correction for quark masses and the continuum limit.

Appendix -Supplementary tables
In the following we collect tables containing additional results and information on the fits performed for this study. Tables 6,7,8,9,10 contain the value of the lower bound of the fit ranges [t ij 1 /a, .., t ij 2 /a], t ij 2 = T /2 − 2a to the quark-connected derivative correlators, which are required for replacing them by the respective ground state correlator. Each table corresponds to one matrix elements as indicated in the captions by the indices i, j = 0, ..., 2, which are used as a short-hand for labeling the elements of the original 3 × 3 correlation function matrix made from the operators in Eqs. (10,11). For technical reasons we perform the replacement of connected contributions at this level, before applying the ratio of renormalization factor Z and rotating to the final operator basis given in Eqs. (12,13) up to a factor Z P . The charm quark operator is only dropped after this for the actual calculation (i.e. before solving the GEVP), reducing the problem to a 2 × 2 matrix. Note that only five out of nine matrix elements involve a quark-connected contribution to the full correlator. These are the ones entirely made up from either only light or only heavy quarks fields. The contractions for matrix elements, which contain both types of quarks yield only quark-disconnected diagrams. In addition, we include the correlated, reduced χ 2 -values (χ 2 /dof) for each fit, together with the corresponding p-value and the value of the mass aM ij conn and its error, which where used for the identification of the final plateau from scanning different values of t 1 /a. Table 11 contains the upper bound for the final fit ranges [t η,η 1 /a, ..., t η,η 2 /a] to the η and η principal correlators from solving the GEVP in Eq. (20). Besides, the resulting correlated χ 2 /dof and p-values are given.
Finally, in Tab. 12 we have included the numerical data for the decay widths and transition form factors for each ensemble, which have been used for Fig. 7 and Fig. 8 (39,40) and the transition form factors in the large momentum limit as given in Eqs. (41,42). Result are obtained using Z factors from method M2, cf. Tab. 2.
Errors are statistical only.