Determination of the charm quark mass in lattice QCD with 2 + 1 flavours on fine lattices

We present a determination of the charm quark mass in lattice QCD with three active quark flavours. The calculation is based on PCAC masses extracted from Nf = 2 + 1 flavour gauge field ensembles at five different lattice spacings in a range from 0.087 fm down to 0.039 fm. The lattice action consists of the O(a) improved Wilson-clover action and a tree-level improved Symanzik gauge action. Quark masses are non-perturbatively O(a) improved employing the Symanzik-counterterms available for this discretisation of QCD. To relate the bare mass at a specified low-energy scale with the renormalisation group invariant mass in the continuum limit, we use the non-pertubatively known factors that account for the running of the quark masses as well as for their renormalisation at hadronic scales. We obtain the renormalisation group invariant charm quark mass at the physical point of the three-flavour theory to be Mc = 1486(21) MeV. Combining this result with five-loop perturbation theory and the corresponding decoupling relations in the MS¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{\mathrm{MS}} $$\end{document} scheme, one arrives at a result for the renormalisation group invariant charm quark mass in the four-flavour theory of Mc(Nf = 4) = 1548(23) MeV, where effects associated with the absence of a charmed, sea quark in the non-perturbative evaluation of the QCD path integral are not accounted for. In the MS¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{\mathrm{MS}} $$\end{document} scheme, and at finite energy scales conventional in phenomenology, we quote mcMS¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {m}_{\mathrm{c}}^{\overline{\mathrm{MS}}} $$\end{document}(mcMS¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {m}_{\mathrm{c}}^{\overline{\mathrm{MS}}} $$\end{document}; Nf = 4) = 1296(19) MeV and mcMS¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {m}_{\mathrm{c}}^{\overline{\mathrm{MS}}} $$\end{document}(3 GeV; Nf = 4) = 1007(16) MeV for the renormalised charm quark mass.


Introduction
Apart from the strong coupling constant, quark masses are the only fundamental parameters of quantum chromodynamics (QCD), the theory of the strong interaction. As such, their precise knowledge is not only of general interest, but also essential in the search for new physics beyond the Standard Model of particle physics. In particular, heavy quark masses serve as key parametric inputs for tests of the Standard Model, via hadronic contributions by QCD matrix elements to weak decays and ensuing constraints on CKM matrix elements [1], as well as in Higgs branching ratios to charm and bottom quarks; see for instance refs. [2,3]. The extraction of these and other Standard Model parameters that have their origin in the strong interaction requires to evaluate the QCD path integral. Since quarks are confined in bound states and cannot propagate freely, their properties have to be investigated in the hadronic regime at low energies where perturbative methods fail. We base our calculations on lattice QCD, where the theory is formulated on a discretised space-time grid to enable, in conjunction with Monte Carlo simulation techniques, an ab initio and thereby non-perturbative solution of the underlying highly non-linear equations with, in principle, full control of statistical uncertainties and systematic effects.

JHEP05(2021)288
Our determination of the charm quark mass is performed on a subset of the gauge configuration ensembles, generated by the Coordinated Lattice Simulations (CLS) cooperative effort of European lattice QCD teams, with N f = 2 + 1 flavours of non-perturbatively O(a) improved Wilson quarks and the tree-level Symanzik-improved gauge action [4,5]. The sea quark content in these discretised QCD configurations comprises a doublet of light degenerate quarks (with masses to be identified with the mean up and down quark mass) and a third, single flavour representing the strange quark. The masses employed in our numerical calculations differ from their physical counterparts in nature, which requires approaching the physical point through a chiral extrapolation to the physical quark mass values corresponding to QCD with 2 + 1 flavours.
In practice, the tuning of the parameters of the theory is guided by the computation of meson masses, where the physical point is defined by the pion and kaon (as well as some charmed meson) masses taking their physical values. As the sum of sea quark masses is held constant in our setup, the strange quark mass reaches its physical value from below. However, we would like to anticipate already now that no significant dependence on the light quark masses is observed in our analysis of observables in the charm sector under study.
The charm quark is included as partially quenched valence quark in our calculation, which means that it is not incorporated as a dynamical degree of freedom in the lattice action governing the path integral of the theory. The size of the systematic uncertainty induced by this approximation has been examined in ref. [6] and was found to be, owing to the decoupling of heavier quarks, of the order of a few per cent at most. On the other hand, the comparison of lattice QCD charm quark determinations relying on 2 + 1 [7][8][9][10][11] and 2 + 1 + 1 dynamical flavours [12][13][14][15][16][17], collected and reviewed in ref. [18] as far as available at that time, demonstrates an excellent agreement. As for the effects of QED and isospin breaking, the statistical error on our final result will turn out to be substantially larger than the level, at which one expects their influence to affect current lattice QCD computations; this justifies to neglect them here.
Lattice QCD offers a sound, first principles definition of quark masses that derives from a lattice version of the PCAC, viz., partially conserved axial current relation, holding in the chiral symmetry realm of QCD. Once bare quark masses based on the PCAC relation on the lattice are extracted (as we do it here in case of the charm quark mass), they need to be properly renormalised for their continuum limit to exist. A convenient object to consider is the renormalisation group invariant (RGI) quark mass, as it is a renormalisation scheme-and scale-independent quantity (if the renormalisation conditions are imposed at zero quark mass [19]), and the determination of the RGI charm quark mass is thus the main goal of this paper. For connecting the bare (PCAC) mass at a specified low-energy scale with the RGI mass in the continuum limit, we can draw on the study exposed in ref. [20], where the renormalisation group running of the quark mass over a wide range of scales for the present three-flavour discretisation of QCD was performed with full non-perturbative precision in the Schrödinger functional scheme, following the strategy of ref. [21]. The results of ref. [20] split into a factor accounting for the universal (i.e., regularisation and flavour independent) running of the quark masses and an associated renormalisation factor at a low-energy scale in an appropriate hadronic scheme, defined via the Schrödinger func-JHEP05(2021)288 tional. Moreover, our computations involve non-perturbatively determined improvement coefficients and renormalisation constants from refs. [20,[22][23][24][25], implementing Symanzik improvement for removal of discretisation effects from correlation functions and chiral Ward identities for finite normalisations of quark bilinears (such as the axial vector current) in the O(a) theory. Therefore, leading discretisation errors of O(a 2 ) in the bulk (except for O(g 2 0 a) uncertainties at the time boundaries) are encountered when taking the continuum limit of correlation functions and renormalised quark masses deduced from them.
The combination of non-perturbative O(a) improvement and the range of lattice spacings from ≈ 0.087 fm down to very fine resolutions of ≈ 0.039 fm covered in our work (where the calibration of the lattice spacing to its physical value was done through the intermediate hadronic gradient flow scale t 0 in ref. [26]) allow us to carry out a controlled continuum limit extrapolation, together with the chiral extrapolation in a joint global fit along a trajectory of constant physics towards the physical point of the theory. In this context, the use of different lattice definitions of the charm quark mass, which distinguish themselves considerably in their cutoff effects, proves to be very valuable, because we are able to confirm that no significant systematic uncertainties are introduced by the global extrapolation to the joint chiral and continuum limit. We find that linear cutoff effects are absent (as expected), but those of O(a 2 ) and higher give relevant contributions already at lattice spacings around 0.06 fm.
The structure of this paper is as follows. In section 2, we introduce the three different definitions of the bare quark mass on the lattice employed and explain their renormalisation pattern in the O(a) improved theory. Section 3 specifies the set of gauge field configuration ensembles and how they realise, in lattice parameter space, the approach to the point of physical quark masses (fixed, in practice, by physical values of pion, kaon and suitable charmed meson masses) along a chiral trajectory and in the continuum limit. There, we also detail the computation of two-point correlation functions, as well as the extraction of meson and bare quark masses from them, and provide the expressions from which the renormalised charm quark mass is obtained and eventually translated into the associated RGI mass. The joint global chiral-continuum extrapolation of the resulting estimates at still unphysical meson masses and finite lattice spacings are then thoroughly described in section 4, where we also elucidate the model averaging procedure that is adopted to arrive at our final result and to quantify its systematic uncertainties. Our final result on the RGI charm quark mass is presented in section 5, together with a careful discussion of the various contributions from the single steps of our calculation to its full error budget. In section 6, we conclude by summarising our findings and compare our final value for the charm quark mass with existing ones from other lattice QCD determinations, see also appendix A, where the conversion through perturbative five-loop running to other conventional energy scales in the MS renormalisation scheme is outlined, which represents the customary reference scheme in phenomenological applications. Lastly, results from several intermediate steps of the analysis, such as meson and (bare) quark masses in lattice units, are tabulated in appendix B.
A preliminary account of our work has been given in ref. [27]. The determination of the mass-degenerate light and strange quark masses on (2 + 1)-flavour CLS ensembles, along the same renormalisation strategy and partly similar analysis methods, was presented in ref. [28].

JHEP05(2021)288 2 Quark masses on the lattice with Wilson fermions
Our determination of the physical charm quark mass is based on various different definitions of the quark mass in the formulation of lattice QCD with Wilson quarks. We begin our discussion of these definitions with the bare subtracted quark mass of flavour i, which is defined as where κ i is the corresponding hopping parameter and κ cr is the critical hopping parameter defined from vanishing current quark masses in the SU(N f ) symmetric limit. The additive quark mass renormalisation via m cr is required due to the loss of chiral symmetry in the Wilson fermion action. The average of two subtracted quark masses of flavours i and j is Since composite operators built from Wilson fermion fields show considerable cutoff effects that are enhanced at the scale of the charm quark, O(a) improvement à la Symanzik drastically improves their approach towards the continuum limit. A detailed discussion of the underlying improvement pattern in the case of three flavours of quarks can be found in ref. [29]. We repeat the expression for the renormalised and improved subtracted quark mass, with the sea quark mass matrix M q = diag(m q,1 , m q,2 , . . . , m q,N f ) , (2.5) and the mass renormalisation constant Z m that not only depends on the bare coupling, but also explicitly on the renormalisation scheme and associated renormalisation scale aµ. 1 We note in passing that, for achieving O(a) improvement in a mass independent scheme, the bare coupling g 2 0 has to be improved with a mass dependent term,g 2 0 ≡ g 2 0 (1+b g 1 N f aTr [M q ]) as introduced in ref. [30]. Yet the chiral trajectory in the quark mass plane associated with our ensembles (see section 3 below) is constructed such as to imply a constantg 2 0 at fixed inverse coupling β = 6/g 2 0 [4], up to negligible violations [26], irrespective of the knowledge of the improvement coefficient b g . Therefore, it is legitimate to utilise renormalisation constants as functions of g 2 0 . As apparent from eq. (2.3), the subtracted quark mass receives a shift from non-zero sea quark masses which scales with (r m − 1), where r m is the ratio of flavour singlet and non-singlet scalar density renormalisation constants. The cutoff effects at order O(a) are mass dependent and get fully cleared if the coefficients of the corresponding terms are fixed non-perturbatively.

JHEP05(2021)288
As an alternative definition of the quark mass on the lattice we employ the bare current quark mass entering the PCAC (partially conserved axial current) relation. For two distinct flavours i = j it reads , (2.6) in terms of the improved axial vector current and the pseudoscalar density P ij , which in the presence of non-degenerate masses are offdiagonal bilinear fields. The standard choice for the discretised derivatives is given by the forward and backward differences while∂ denotes the average of these two. Following refs. [24,31,32], we also adopt another realisation of discretised derivatives based on the substitutions Since the discretisation effects introduced by these derivatives acting on smooth functions are of O(a 4 ), we will refer to this second choice as improved derivatives. Upon renormalisation, eq. (2.6) leads to the renormalised O(a) improved current quark mass m R,ij , which in the theory with mass non-degenerate quarks is equivalent to the average of two renormalised current quark masses of flavour i and j. Here, Z A (g 2 0 ) and Z P (g 2 0 , aµ) label the renormalisation constants of the axial current and the pseudoscalar density, respectively, the latter carrying the renormalisation scheme and scale dependences of the quark mass as does Z m in eq. (2.3) above.
We can now equate both expressions for the renormalised quark mass, eqs. (2.3) and (2.10), to derive where Z is the scale independent combination of renormalisation constants , (2.12) and employ this relation to eliminate the bare subtracted quark mass in eq. (2.10) in favour of PCAC masses to obtain Note that in eq. (2.13) all dependence on m cr has been removed. Yet another definition of the renormalised quark mass with different cutoff effects may be constructed along the so-called ratio-difference method [33], where we again combine the two definitions for renormalised quark masses, eqs. (2.3) and (2.10), but in a different way than before. Namely, we introduce a ratio of current quark masses, r, and a difference of subtracted quark masses, d, through where i and i label two distinct but mass degenerate quark flavours. This setup is chosen such that the multiplicative renormalisation of the current quark masses cancels in the ratio, while the additive renormalisation of the bare subtracted quark masses (via m cr ) and the leading-order dependence on Tr [M q ] drop out in the difference. As explained in detail in ref. [34], the O(a) improved versions of r ij and d ij can be written as From the latter, an estimator for the renormalised quark mass is formed by with the mass renormalisation factor Z m already introduced in eq. (2.3). For the improvement and renormalisation of the axial current we use the improvement coefficient c A calculated in ref. [22] and the values of Z A (g 2 0 ) determined within the chirally rotated Schrödinger functional [23]. The scale dependent renormalisation constant Z P was computed in ref. [20], together with the universal factor that relates quark masses at finite renormalisation scale to their RGI counterparts. The renormalisation constant Z defined in eq. (2.12) was calculated for our setup in ref. [24]. It can also be extracted exploiting the numerical results on Z S /Z P presented in ref. [34].
The non-perturbative determination of the improvement coefficients (b A − b P ) and b m has been performed in ref. [24], along with Z, for two different lines of constant physics (LCP) realised in the Schrödinger functional framework for the range of lattice spacings relevant here. More precisely, there is a set of coefficients computed on an LCP adapted for use in a fully massless setup (called 'LCP-0' in [24]), whereas, by contrast, a further set of estimates (labelled 'LCP-1') for Z, (b A − b P ) and b m refers to an LCP at fixed heavy quark mass in the valence sector. Therefore, particularly the latter appears to be predestined for applications with heavy (valence) quarks, because mass dependent cutoff JHEP05(2021)288 effects of all orders are expected to be accounted for in the improvement coefficients and renormalisation constant from this LCP. In fact, the investigations in refs. [35,36] have shown that the scaling of meson observables involving heavy quarks towards the continuum limit benefits from including these mass dependent cutoff effects in the LCP definition of the improvement coefficients and renormalisation factor in question. This conjecture will be tested later, when we discuss the chiral-continuum extrapolation of our results.
Non-perturbative results for r m for the full range of couplings recently became available from the determination in ref. [25] (earlier findings for a subset of the CLS ensembles are available in ref. [37]). We employ these for the corresponding improvement terms in eqs. (2.13) and (2.18) and observe only a marginal influence on heavy quark masses at finite lattice spacing.
Finally, for the improvement coefficient (b A −b P ) andb m multiplying cutoff effects depending on the sea quark masses, non-perturbative results with errors of O(100%) were obtained in refs. [38,39] using coordinate-space methods. Given the large uncertainties in these determinations and the fact that the coefficients are of O(g 4 0 ) in perturbation theory (as opposed to O(g 2 0 ) for the coefficients b X mentioned before), we ignore the corresponding terms of O(aTr [M q ]) in our improvement procedure. We note that the associated contributions from the sea quark masses are highly suppressed compared to the heavy valence quark mass effects that are removed non-perturbatively.

Computational details on the quark mass extraction
The numerical study described in the following has been carried out on the (2 + 1)-flavour CLS gauge field configuration ensembles described in refs. [4,5]. Gluons are implemented via the tree-level improved Lüscher-Weisz gauge action [40] and quarks are treated as clover improved Wilson fermions [41]. All ensembles employed in our computations feature open boundary conditions in time direction to prevent a freezing of the topological charge [42,43]. Twisted mass reweighting [44] has been applied in the simulations for the doublet of light sea quarks and the RHMC algorithm [45,46] has been used for the single strange quark in the sea.
Charm quarks are introduced in the valence sector only, rendering our theory a partially quenched approximation of QCD. Whereas the effect of a dynamical charm quark in lowenergy purely gluonic (and light flavour) observables is expected to be suppressed below the available accuracy [47], this may be different for observables at the scale of the charm quark, and therefore a systematic uncertainty is present owing to neglecting a dynamical charm in the sea in our final result for the physical charm quark mass. A conservative upper bound of 5% for the relative size of this uncertainty was estimated in refs. [6,48] from calculations with two dynamical charm quarks.
Ensembles with five different lattice spacings in a range from 0.087 fm down to 0.039 fm, together with non-perturbative O(a) improvement, allow for a controlled continuum limit extrapolation. Moreover, these ensembles cover pion masses in the range [200 MeV, 420 MeV] such that the dependence of our results for the charm quark mass on the light quark masses may be carefully investigated. In general, however, the effect of JHEP05(2021)288 unphysically large sea quark masses on observables at the scale of the charm quark is expected to be small. We perform our calculations on a set of ensembles, for which the sum of the bare sea quark masses Tr [M q ] is held fixed. In such a setup, the strange quark mass approaches its physical value from below when the light quark masses are lowered towards their physical values.
The strategy to fix an LCP trajectory in bare lattice parameter space within the Tr [M q ] = const. CLS ensembles and its exact specification are explained in great detail in refs. [26,49]. Since a fixed value of the bare trace of the quark mass matrix does not imply a constant renormalized trace of the quark mass matrix, the generated gauge field configuration ensembles do not lie strictly on an LCP. Therefore, and since the deviations of Tr [M R,q ] from a constant value have been observed to be larger than what would be expected from O(a) cutoff effects [26], the chiral trajectory has been redefined in terms of a constant value of where the gluonic quantity t 0 is defined from the Wilson flow [50] and m π and m K denote the masses of the lightest mesons composed of light and light-strange quarks, respectively.
To approach the physical point in a chiral extrapolation, the value of φ 4 on each ensemble has to be fixed to the physical value based on the physical values of pion and kaon masses in isospin symmetric QCD as detailed in ref. [51] and the physical value of t 0 determined in ref. [26]. In this appropriately adapted setup, the dependence of physical observables on the sea quark masses can now be parametrised via on each ensemble, we resort to small quark mass shifts onto the chiral trajectory by means of a Taylor expansion of the entering correlation functions as described in ref. [26]. Following this reference, the derivative of an observable O with respect to a change in the quark mass m q,i of flavour i is given by In order to carry out this shift at the level of the correlation functions employed in our work, their partial derivatives were numerically evaluated (simultaneously with the calculation of the correlators themselves) utilising the mesons program package [52], augmented by sets of measurements for the derivative terms of the action with respect to the quark mass in (3.5) that we had to extend compared to what was already available as a result of [26].

JHEP05(2021)288
The shift of observables depending on N f = 2 + 1 quark masses is then accomplished by the replacement and we choose the same value of ∆m q,i for all three quark flavours. In practice, the value of ∆m q,i on each ensemble was obtained in an iterative procedure that stops when the shifted φ 4 matches its physical value well below the statistical precision.

Correlation functions and fixing of the physical charm quark mass
As already indicated above, we have performed correlator measurements on 15 ensembles on the Tr [M q ] = const. trajectory from the (2 + 1)-flavour CLS gauge configuration data base, using the mesons code [52]. The specifications of these ensembles and their associated statistics in the number of molecular dynamic units (MDU) are summarised in table 1. To extract ground state meson masses as well as quark masses through the PCAC relation, we consider zero-momentum two-point correlation functions defined by where r and s are flavour indices. y 0 denotes the time coordinate of the source, i.e., the time slice where the source has a non-vanishing norm, and x 0 is the time coordinate of the sink. The operators O rs are quark bilinear covariant fields composed as in terms of Γ, representing the combination of Euclidean gamma matrices that read γ 0 γ 5 for the time component of the axial current, A 0 , and γ 5 for the pseudoscalar density, P . 2 Since translational invariance in time direction is broken by the open boundary conditions, we fix the temporal source positions y 0 to the two time slices at a and T − a as proposed in ref. [53]. To achieve good statistical accuracy, we average over 16 U(1) noise sources per source position and exploit time reversal symmetry.
The correlation functions f rs OO were evaluated for all combinations of light and strange as well as for two choices of heavy valence quark flavours. For the latter, the hopping parameters have been chosen such that they encompass the supposed hopping parameter of the physical charm quark κ c sufficiently close. An overview of these hopping parameters is given in table 2. Distance preconditioning [54], in its variant as sketched in ref. [55] (and implemented in the mesons measurement code [52]), has been employed for the calculation of the heavy quark propagators, in order to ensure the stability and accuracy of the associated numerical inversions of the Dirac operator across the complete temporal extent of the lattice.
We follow two different strategies to tune to the physical charm quark mass. Firstly, since Tr [M q ] is held constant in our framework, the flavour averaged meson mass We list the geometry (T × L 3 ) and the bare parameters of the action as well as the resulting approximate lattice spacings and meson masses. The statistics included is quantified by the number of molecular dynamics units N MDU . We performed measurements every 4 MDU except for J303 and J500, where the spacing between two measurements is 8 MDU. The exponential autocorrelation time is defined from t 0 as given by eq. (3.10).
built from D and D s mesons, is approximately constant as well. At the same time, the statistical precision of the corresponding effective masses is convincingly good. Secondly, we consider the effective mass of the connected part of the η c meson, called m ηc from now on. The signal of this observable is clean, and the effect of neglecting the disconnected part is expected to be subleading [56], particularly so for singlet mesons composed of heavy quarks. Hence we have performed all measurements including heavy quarks for both charmed hopping parameters and, to eventually obtain the physical charm quark mass, interpolate the meson mass results to their respective physical, i.e., experimental values jointly with a chiral-continuum extrapolation procedure, as will be explained in detail in the next section. Let us mention that we also had tested the spin-flavour averaged meson mass constructed from heavy-light and heavy-strange pseudoscalar and vector mesons [27] to fix the charm quark mass. 3 However, even though one is able to extract the ground states of the vector mesons (despite the correlators being considerably noisier than in the pseudoscalar channel) on all ensembles, we did not see any improvement in the approach of JHEP05(2021)288  Table 2. Hopping parameters of the partially quenched heavy valence quarks in the charm region serving to implicitly fix κ c , provided by [59] in the context of determining leptonic decay constants of D and D s mesons [55,60]. our results to the continuum limit compared to the other two methods. Moreover, due to the still relatively large statistical uncertainties of the vector meson masses, there is also no gain in overall statistical precision such that we did not include the spin-flavour averaged meson mass in our final analysis. We carry out our statistical error analysis using the Γ-method [61] deriving from autocorrelation functions, supplemented by automatic differentiation for the error propagation [62], and take into account the effect of the exponential tails in the autocorrelation functions owing to slow modes [63]. As estimate for the required exponential autocorrelation time on the CLS ensembles we recourse to the formula quoted in ref. [4] for our set of ensembles.

Extraction of meson and bare PCAC quark masses
Ground state masses of pseudoscalar mesons are extracted from the effective mass , (3.11) which forms a plateau for sufficiently large source-sink separations. The computation of pion and kaon masses on lattices with open boundaries has been discussed in refs. [4,26,53], and we employ the same methods to determine the plateau region, i.e., fits to the exponential corrections originating from boundary effects and short source-sink separations.
Including the dominant contributions amounts to model the effective masses with an ansatz where m PS gives the mass of the pseudoscalar meson in question, E 1 = m − m PS the difference between ground and first excited state energies, and E 2PS 2m PS . The energies and the parameters c i are chosen to be free fit parameters. The same procedure is applied to the effective mass of the η c meson. For heavy-light and heavy-strange mesons we restrict ourselves to a region where the statistical fluctuations are small such that effects of the

JHEP05(2021)288
boundary at maximal source-sink separation do not affect the quality of the plateau. We show representative effective masses for these mesons, entering into the interpolations to fix the physical charm quark mass, in the left part of figure 1; there the exponential corrections for small source-sink separations and close to the boundaries are clearly visible.
To specify the optimal fit window, we vary the minimal source-sink distance included in the fit and monitor its quality through the so-called 'corrected' χ 2 , cf. eq. (4.8), of an uncorrelated fit. χ 2 corrected is based on the correlations present in the data, as motivated and elaborated in ref. [64]. The beginning of the plateau range is set to the time slice where the contribution of the excited state corrections is four times smaller than the statistical error on the effective mass. The upper end of the fit interval is taken from a fit to the boundary corrections in case of mesons with constant signal to noise ratio, while it is a bound of 3% on the relative error of the correlation function in case of heavy-light mesons. The meson mass can then be defined either as the weighted plateau average within the plateau range thus determined or directly from the fit to the plateau and its exponential corrections according to eq. (3.12). Because of their compatibility with each other, in the final analysis we just opted for the latter. Note that thanks to the use of distance preconditioning the data on heavy quark propagators over the full temporal range can be incorporated in the plateau fits. This effect is actually enhanced on the J ensembles close to the continuum limit, which have a maximal source-sink separation of 189a.
As an identity on the operator level, the PCAC relation is valid on every time slice and, consequently, current quark masses exhibit a plateau, too. Contact terms lead to deviations from this plateau behaviour for small source-sink separations, as do cutoff effects close to the boundary. Following ref. [26], we estimate the time slice extent of these deviations with the help of the same methods as for the effective meson masses above, i.e., modeling the corrections by exponential fits. In identifying plateaux, we ensure to extract the mass values in the bulk of the lattice where O(g 2 0 a) cutoff effects from the open boundaries are absent. For heavy-light and heavy-strange quark masses, the signal is lost in noise for large source-sink separations. Typical current quark mass plateaux are displayed in the right part of figure 1.
In appendix B we tabulate the ensemble-wise results for the meson masses and the different bare current quark masses. In each table we list results before and after the mass shifting procedure along the prescription in eq. (3.6). Accounting for these mass shifts increases the statistical uncertainties of the meson and bare quark masses considerably. This is due to the fact that the target for the mass shift, φ phys 4 , is only known with about 1 % precision, dominated by the uncertainty of the physical lattice scale. It should also be noted that the target value of φ phys 4 is the same for all ensembles. For this reason the shifted masses from different gauge field ensembles cannot be considered uncorrelated.

Quark mass combinations
The PCAC definition (2.6) [41,94] hl, [40,54] 0.582 0.584 hh , [61,158] 0 20 40 [28,102] hl, [30,58] 0.1365 0.1366 hh , [45,166]  where c and c denote two distinct but mass degenerate quark flavours and the appropriate renormalisation and improvement pattern is given in eq. (2.13). This definition always provides clean signals and therefore small statistical errors of the plateau averages. The R,c on the light quark masses (encoded in its pion mass dependence) is small. Since the mass dependent cutoff effects are at the scale of the charm quark mass, we expect them to be sizeable, although effects of O(a) are cancelled non-perturbatively in our lattice formulation of the theory.
The size of these cutoff effects is reduced when heavy-light and heavy-strange correlation functions are employed from which estimators for the mass of the charm quark may then be obtained via As evident in our data, the definitions in eq. (3.14) have a relevant dependence on the sea quark masses which, however, turns out to be significantly reduced in the associated flavour-averaged charm quark mass because of the average sea quark mass being held constant along our trajectory in lattice parameter space towards the physical point. The signal of heavy-light and heavy-strange quark masses deteriorates for large source-sink separations. Accordingly, the statistical error on the charm quark mass definition m R,c is generally larger than on m (c) R,c . As a third option to arrive at a sensible prescription to compute the charm quark's mass, we adopt the ratio-difference method, eq. (2.19). Here, the first flavour is chosen to be the heavy valence quark and the second flavour is one of the sea quarks. Adhering to the same argument in favour of the linear combination in (3.15), we deduce the renormalised charm quark mass from the corresponding flavour average also in this case, i.e., with the improved ratio r I,ij and difference d I,ij already given in eqs. (2.17) and (2.18), respectively. The anticipated advantage of this estimator is that systematic uncertainties owing to an imprecise determination of the chiral point, which are naturally associated with subtracted quark masses, are cancelled in the quark mass from the ratio-difference method.

Renormalised quark masses
As a central result of this paper, we wish to quote the renormalisation group invariant (RGI) value of the charm quark mass. For a given quark flavour, the RGI quark mass, M , is defined through the formally exact expression where it is assumed that renormalisation conditions are imposed at zero quark mass, which suggests calling those schemes mass independent. This implies ratios of renormalised quark masses for different flavours to become scale and scheme independent constants. In the JHEP05(2021)288 same way as the QCD Λ-parameter associated with the running of the renormalised coupling g R (µ), M belongs to the RGI quantities whose total dependence on the renormalisation scale µ vanishes. In addition we recall that, even though perturbative expansions for the RG group functions β(g) and τ (g) exist, the Λ-parameter and the RGI quark masses are generically defined independent of perturbation theory and, particularly, the RGI quark masses M are renormalisation scheme independent. Therefore, we consider them ideal for comparisons with results from either experiments or other theoretical determinations.
Conventions regarding the RG β-and τ -functions as well as their respective perturbative coefficients b i and d i can be found, e.g., in refs. [21,65]. The renormalised coupling g R (µ) and the continuum RG functions of the running coupling and quark mass entering eq. (3.17) were accurately obtained for the three-flavour theory at hand by applying the non-perturbative step scaling approach within the Schrödinger functional scheme [20,66,67]. Hence, we have all the ingredients available to merge the above definition of M with one of the estimators for the renormalised charm quark mass introduced in the previous subsection and calculate the RGI charm quark mass from it. Explicitly, on the basis of the expressions in eqs. (3.13), (3.15) and (3.16), their relation to the RGI mass then reads: The renormalisation scale dependence of m R,c is inherited from the corresponding scale dependence of the respective Z-factor (see section 2), which turns the underlying bare quark mass into the renormalised one. Thanks to the determination of the universal, flavour independent ratio M m R (µ had ) = 0.9148(88) , (3.19) in the Schrödinger functional scheme for N f = 3 massless flavours at µ had = 233(8) MeV [20], the renormalisation factors and the ratio M/m R can be matched at the hadronic scale µ = µ had such that M c is readily returned. Further details concerning the non-perturbative computation of the running factor (3.19) can be found in ref. [20] and references therein.

Chiral-continuum extrapolations
After having extracted the renormalised quark masses as well as the relevant meson masses from every gauge field configuration ensemble, we proceed with a combined chiral and continuum limit extrapolation in order to obtain the value of the charm quark mass at the physical point, defined by φ 2 = φ phys 2 and zero lattice spacing. All gauge field ensembles listed in table 1 are included in the analysis, except for H200 with a spatial extent of ∼ 2 fm that we consider too small. 4 In the extrapolation, we include our different definitions (3.13), (3.15) and (3.16) of the renormalised charm quark mass, in conjunction with JHEP05(2021)288 eq. (3.19) and made dimensionless by attaching the factor √ 8t 0 . For the renormalisation constant Z, as well as the improvement coefficients b A − b P and b m entering these definitions, the numbers denoted as 'LCP-1' in ref. [24] are inserted, as their determination was specifically suited for the heavier quark regime and provides high statistical accuracy.
As consistency checks, we also tried using the 'LCP-0' values of these parameters from ref. [24], as well as those from the alternative Ward identity determination of Z in ref. [34], and found slight differences at finite lattice spacing but full agreement in the continuum limit. We utilise the results for r m from ref. [25] in which is guided by the following reasoning: having corrected the chiral trajectory by means of the mass derivatives in eq. (3.6), the light quark mass dependence is governed by the sea pion mass only, whereas the kaon mass is implicitly fixed by the Tr[M q ] = const. condition realised along the trajectory in lattice parameter space towards the physical point. Based on the Gell-Mann-Oakes-Renner relation [69], we thus assume the leading contribution to be proportional to φ 2 . Moreover, we model the interpolation to the physical value of the bare charm quark mass via where m H labels the mass of a charmed meson for which we either chose the flavour averaged D meson mass, mD, or the mass of the pseudoscalar charmonium meson, m ηc (neglecting all quark-disconnected contributions). Motivated by heavy quark effective theory, we expect the charm quark mass to be linearly dependent on the meson mass of choice [58].
For the continuum part we presume the leading cutoff effects to be of O(a 2 ) as all relevant renormalisation constants and improvement coefficients that propagate into our computation are known non-perturbatively. In addition, we also allow for terms describing the next two higher orders, O(a 3 ) and O(a 4 ). These cutoff effects may be quark mass independent or explicitly depend on the quark masses. Therefore, our general ansatz for the lattice spacing dependence of the charm quark mass is parametrised by all cutoff effects proportional to the light quark masses are already neglected here, because it turned out that these effects cannot be resolved in the fits. Furthermore, for now we also and particularly masses of heavy mesons [68] are expected to be unaffected by moderate violations of the conditions that typically qualify lattice QCD ensembles as resembling an infinite volume situation. 5 As argued at the end of section 2, the terms proportional tobm and (b A −b P ) are not expected to be significant and therefore ignored in our improvement procedure.

JHEP05(2021)288
neglect terms proportional to mixed powers of a √ 8t 0 and φ H , such as a 3 (8t 0 ) 3/2 φ 2 H or a 4 (8t 0 ) 2 φ 2 H . We will demonstrate below that the inclusion of these terms (which a priori cannot be excluded on theoretical grounds) has no effect on our final result. As explained at the end of section 2, we do not expect terms of O(g 4 0 aTr [M q ]) to contribute significantly to the cutoff effects. We have checked explicitly that the inclusion of such terms does not alter the continuum limit and found the corresponding fit parameters to vanish within errors.
To arrive at a combined model for the global fitting procedure, we can compose the chiral and continuum parts either linearly, by adding eqs. (4.1) and (4.3), viz. 4) or in a non-linear fashion by multiplication of the two: the c i in (4.1) and (4.3) are fit parameters. From now on we will refer to different model parametrisations with the following labelling scheme: • The first entry is either 'l' for the linear combination of the chiral and continuum parts or 'nl' for the non-linear variant.
• The second entry indicates the term which accounts for cutoff effects proportional to a 2 . Here, '1' stands for the term parametrised by c 1 , whereas '2' stands for the (heavy meson) mass dependent one, parametrised by c 2 . When both pieces, proportional to c 1 and c 2 , are included, the entry is labelled by '3'.
• The third entry represents the contribution accounting for cutoff effects of O(a 3 ).
Here, '0' specifies fits without any contribution of this order, whereas '1' does so for fits with the c 3 -term, '2' for fits with the c 4 -term and '3' for fits with both of them.
• The fourth entry labels the terms describing discretisation effects of O(a 4 ). In analogy to the third entry, '0' labels no term of this order, '1' the term proportional to c 5 , '2' the one proportional to c 6 and '3' both of them.
• Finally, the fifth entry indicates whether a pion mass dependence is included in the fit, i.e., '0' marks models without and '1' models with the term parametrised by the coefficient c L .
As an example, the model labelled by ('l', 3, 0, 1, 0) corresponds to the fit function (4.6) We always include at least one term describing cutoff effects of O(a 2 ) and restrict the maximum number of terms related to the cutoff effects to three. In total this amounts to K = 108 different models, with which we attempt to fit our data.

JHEP05(2021)288
In order to estimate the optimal parameters for a given model (generically designated p i and f in the formula below), we employ an extended χ 2 minimisation procedure, taking into account the errors of both, the independent and the dependent variables [70], viz.
The independent variables x a are promoted to fit parametersx a , which albeit stay constrained to their initial values x a and thereby only indirectly affect the number of degrees of freedom. In our case, the full correlation matrix is ill-conditioned and eludes safe inversion. 6 For this reason we ignore all correlations and set the off-diagonal elements of the error covariance matrix C explicitly to zero. While minimising the uncorrelated χ 2 yields reliable maximum likelihood estimators for our parameters (see ref. [71]), the minimal value of the uncorrelated χ 2 cannot serve as a meaningful indicator for the goodness of fit, as we do not know the true number of degrees of freedom. In practice, we would thus obtain values of χ 2 /(N − k) < 1 also for models that do not describe the data well. We circumvent this issue by evaluating χ 2 exp , the so-called 'expected' χ 2 proposed and advocated for use in such situations in ref. [64], 7 which equals the expectation value of χ 2 given normally distributed data with the full covariance matrix C. As explained in this reference, it qualifies to derive a corrected measure for the goodness of fit, in the sense that χ 2 corrected /(N − k) = 1 when the model describes the data.

Systematic effects
In order to address systematic effects in our computation of the charm quark mass quantitatively, we classify our data in different categories w.r.t. the various choices and technical ingredients that characterise the analysis: • For the bare current quark masses involved, either standard derivatives according to eq. (2.8) or its improved version in eq. (2.9) are used.
• To fix the physical charm quark mass, we either employ the flavour-averaged D meson mass, mD, or the mass of the pseudoscalar charmonium groundstate, m ηc (neglecting the quark-disconnected contributions). 6 Whereas Monte Carlo data from different gauge field ensembles is uncorrelated by construction, the full covariance matrix does not assume a block diagonal form, but receives additional contributions induced by common renormalisation and improvement parameters as well as the mass shifting procedure. 7 Note that χ 2 exp was already successfully applied before in a fit analysis of HQET correlators in the context of extracting semi-leptonic form factors [72].

JHEP05(2021)288
• To ensure that cutoff effects are properly described by the fits, we either include all gauge field ensembles, i.e., covering a range of lattice spacings 0.039 ≤ a fm ≤ 0.087, or a cut of 0.039 ≤ a fm ≤ 0.076 is imposed on the data, i.e., all ensembles at β = 3.4 are excluded.
This amounts to 24 categories in total. For each category we perform fits with all 108 fit forms specified in the previous subsection, resulting in a total of 2592 fits. In the following we discuss some of these fits in more detail. As measure for the best fit in each category we adopt the Akaike information criterion (AIC) [73] based on the above-introduced χ 2 corrected , where k is the number of parameters in the fit function.
In figure 2 we compare the best fits in different representative categories. The top part illustrates the effect of using improved derivatives in the PCAC masses on the definitions of the charm quark mass. Starting with m R,c , we see that improvement of the derivative leads to larger cutoff effects compared to the standard one. However, both variants agree very well in the continuum limit. This is also the case for m

JHEP05(2021)288
general tend to slightly larger values at finite lattice spacing. However, this ambiguity vanishes in the continuum limit as expected.
Since the fit distinguished as the best one by the AIC (or any other sensible measure) could in fact be also an outlier, we opt for performing a model averaging procedure within each category, similar to what has been proposed in ref. [75]. In analogy to this reference, we define the model average via (4.10) with the respective model weights w k and k indexing the alltogether K = 108 models (i.e., fit forms) within each of the 24 categories at disposal. For the systematic error arising from this model selection procedure we assume (4.11) For uniform weights, w k ≡ 1/K, this expression reduces to the variance of the different model estimates. Yet instead of setting the weights uniformly, we rather set the w k for each fit relying on the AIC, where the normalisation constant N is chosen to ensure K k=1 w AIC k = 1. In figure 3 the outcome of the model averaging procedure is summarised for two representative examples among the 24 categories. The individual circles correspond to the result of fitting the data to the respective model, labelled along the horizontal axis according to scheme introduced above. Their opacity indicates the associated weight in the average, the most opaque points having largest weight. The filled square and diamond mark the model averages with their statistical errors, while the shaded bands depict their systematic uncertainties emerging from the model averaging procedure. As a consequence of the definitions of the renormalised charm quark mass, eqs. (3.13), (3.15) and (3.16), which were suitably adapted to the position of the CLS ensembles in physical parameter space and the approach towards the chiral limit ensuing from it, we expect the (sea) pion mass dependence of M c to be small. In fact, the best fits in all 24 categories (and also the majority of fits with a relevant weight) do not incorporate a pion mass dependence. As a representative example we illustrate the pion mass dependence for the charm quark mass defined in eq. (3.15) in figure 4. Other quark mass definitions show the same flat behaviour along the chosen chiral trajectory for which φ 4 is kept constant.
As can be inferred from the top part of figure 3, almost all fits with non-negligible weight agree very well, which is reflected by the small systematic error. In this case, singling out only the best fit would give nearly an identical result compared to the model average. In the bottom part of figure 3, the best fit is ('nl', 2, 1, 1, 0)   large number of models that describe the data reasonably, outweigh the few outliers. The somewhat wider spread of the fits with more significant weight is revealed by the larger systematic model averaging error in this case. Figure 5 displays the mean values and errors from the model averages within each of the 24 categories specified above. In principle, the result in every single category represents a valid determination of the charm quark mass, which means that for infinite statistics we would expect all 24 determinations to coincide perfectly. Hence, to decrease the final error, one could even consider combined fits of different categories. However, as we will discuss in the next section, our error is not dominated by the statistical error arising from the underlying Monte Carlo data, but by the propagated error from external inputs. For this reason we have decided to conservatively take the spread of the individual categories' results as additional measure for the systematic error of our calculation.
For our final result of M c we chose the weighted average over all categories in figure 5, except (as argued before) for the ones with m (c) R,c and the standard derivative. The weights are given by the inverse of the respective errors (i.e., statistical and systematic ones added in quadrature). From the weighted standard deviation (4.11) of the model averages within these categories we estimate the size of the systematic error in our determination of the charm quark mass. However, we have convinced ourselves that an unweighted average leads to a number which is indistinguishable from the weighted average in relation to the total error for the final result quoted below. 9 In order to check the robustness of our analysis method, we have also induced different weights in the model averaging procedure, devised to set different penalties on the number of model parameters and thereby to favour over-or underfitting opposed to the standard procedure. Specifically, the Bayesian information criterion [76] as well as the corrected AIC [77] were tried, which give higher penalties for models with more fit parameters, while basing the weights on just χ 2 corrected introduces no penalty at all. The resulting variation on our final estimate of M c when using these different weights is by far smaller than our systematic error as assessed above. Furthermore, as already mentioned earlier, we investigated the effect of allowing for terms such as a √ 8t 0 (am H ) 2 or a 2 8t 0 (am H ) 2 in the generic fit ansatz for the joint chiral and continuum limit extrapolation, thus increasing the number of possible models in each category to K = 168. Again, also this extension of the number of models has a completely negligible effect on the final result. In summary, we therefore conclude that grounding our analysis on a large number of models and categories very consistently yields a reliable result that is fairly independent of human bias, because the manual input is only limited to a theoretically well justified pre-selection of considered categories and models. input data variation 8.2% Table 3. Error budget of our computation, quantified by the contributions to the squared error of our final estimate for M c . We group them into the three main categories 'Statistical error', 'Renormalisation group running' and 'Systematic uncertainty', which comprise the final (total) error, where the listed individual subcategories indicate how these uncertainties are distributed among the main categories.

Result
Following the data analysis presented in the foregoing sections, we quote for the RGI charm quark mass from our lattice QCD computation as final result which constitutes an unambiguous result in (the continuum limit of) the N f = 2 + 1 theory studied here. The first error is statistical, the second stems from the factor converting renormalised into RGI quark masses, available from ref. [20], while the third one represents our estimate of systematic effects. In table 3 we summarise the error budget of our computation. The most dominant contribution to the total error originates from the (universal) conversion factor M/m R (µ had ) that connects renormalised quark masses at finite scale via non-perturbative renormalisation group running to their RGI counterparts (cf. subsection 3.4) and was obtained for the three-flavour theory in ref. [20]. The next relevant contributions come from the scale setting of ref. [26], and the O(a) improvement b-coefficients multiplying the quark mass dependent pieces, non-perturbatively obtained in ref. [24]. For the systematic uncertainty, which is made up of two contributions, the one from the spread of outcomes of physical point extrapolations within different categories (i.e., characterised by variations of the input data as outlined in subsection 4.1) outweighs the other one from the model averaging procedure itself (i.e., involving the various fit forms). Note that the statistical error of the correlation functions entering the individual definitions of the charm quark mass is by far subleading in our calculation and that the uncertainties on the experimental values of the meson masses, which feed into the analysis to set the physical point, are completely irrelevant for the total error of our central result (5.1).

JHEP05(2021)288
As the final uncertainty of our charm quark mass determination is completely dominated by external input, there is hardly any room for a significant gain in precision, even if the underlying data base of gauge field configuration ensembles on the Tr[M q ] = const. trajectory in lattice parameter space were moderately larger. One key prerequisite for an improved overall precision on the RGI (charm) quark mass along the methodology presented in this work would be a more accurate knowledge of the non-perturbative quark mass running factor. A promising step in that direction is the idea of a renormalisation strategy proposed in ref. [78] that exploits the decoupling of heavy quarks from low-energy physics. Another substantial refinement in precision may be anticipated from an update of the physical scale setting parameter, 8t phys 0 , as it is currently being pursued within the CLS effort (see refs. [79,80]).
The effects of charm loops are absent in our calculation, as we performed it on gauge ensembles with 2 + 1 quark flavours in the sea. Still, for comparison to other results (from lattice and non-lattice QCD determinations) as well as later reference, it is customary to convert our result to the four-flavour theory through the well-established prescription relying on perturbative decoupling in the MS scheme. This procedure is sketched in some detail in appendix A, where we also provide values of the charm quark mass in the MS scheme at conventional energy scales. Particularly, for the four-flavour RGI mass we arrive at: Its error splits into the part derived from the uncertainty of the three-flavour RGI charm quark mass above, to which the subdominant part invoked by Λ

Summary and discussion
We have calculated the charm quark mass in QCD with N f = 2 + 1 dynamical quark flavours, on the basis of a large set of gauge field configuration ensembles generated by the Coordinated Lattice Simulations (CLS) initiative [4,5] and employing a partially quenched setup, consisting of non-perturbatively O(a) improved Wilson quarks in the sea supplemented by a quenched charm quark in the valence sector. Bare current quark masses were extracted from the PCAC relation. Notable features of our determination are: five (very) fine lattice spacings from 0.087 fm down to less than 0.04 fm, fully non-perturbative O(a) improvement of the quark mass and the use of the non-perturbative renormalisation group running, available from ref. [20], in order to arrive at the RGI (renormalisation group invariant, i.e., scheme and scale independent) estimate of the charm quark mass as central result. Moreover, our analysis exploits three different definitions of the quark mass on the lattice, which also allows systematic effects to be investigated in a controlled way. To ensure reliability and robustness of our result that, for given possible variations in the analysis, includes a large enough (albeit theoretically well-founded) number of model functions for the joint chiral and continuum extrapolation whose individual contributions are correctly JHEP05(2021)288 weighted by their fit to the actual data, we have adopted a model averaging procedure inspired by the prescription proposed in ref. [75]. In this regard, our final result can therefore be considered as model-independent and its error as accounting for all systematic effects of the present computation in the studied N f = 2 + 1 approximation of QCD.
The final numerical outcome of our work is to be taken from eqs. (5.1) and (5.2) of the previous section, where we quote the RGI charm quark mass in the three-flavour (i.e., N f = 2 + 1) and four-flavour theories, respectively, the latter being obtained by virtue of perturbative decoupling relations. It is worth emphasising that the precision of our result can only be substantially improved by reducing the uncertainties on external quantities (most notably, the renormalised-to-RGI quark mass conversion factor and the lattice scale in physical units) which enter our calculation. Furthermore, while the joint extrapolations of the charm quark mass to the physical point exhibit, in line with O(a) improvement, no corrections linearly, but quadratically and of higher order in the lattice spacing, they are not able to resolve any dependence on the sea pion mass; this stands in distinct contrast to the calculation of light quark masses on a subset of the gauge configuration ensembles at hand, reported in ref. [49], for which a decent decrease of errors can be expected from accommodating more chiral ensembles (especially with fine lattice spacings) and more statistics at light quark masses.
As soon as the relative error is smaller than about 1%, isospin splitting and QED effects may become noticeable. For instance, the impact of the inclusion of quenched QED was found to be around 0.2% in ref. [17], which is well below our final uncertainty. As apparent from figure 6 where we compare our result with other lattice QCD ones from the literature, no significant difference between 2 + 1 and 2 + 1 + 1 flavour results can be unravelled at the current level of precision of lattice QCD calculations.
After all, in appendix A, we explain how we translated eqs. (5.1) and (5.2) to associated charm quark mass values in the MS scheme, present our results at customary energy scales, and confront them with those of other lattice QCD collaborations. Here, we compare our four-flavour RGI mass estimate M c (N f = 4) = 1548 (23) MeV to the current FLAG average for calculations based on N f = 2 + 1 dynamical quarks [18], which is composed of the numbers from refs. [7][8][9]

A Perturbative quark mass running to conventional scales
In view of further applications or reference in phenomenological studies, the charm quark mass is often quoted in the MS renormalisation scheme and the four-flavour theory. In the case of heavy quarks, in particular, it is also common practice to do this in form of the scale invariant mass, which amounts to implicitly fixing the renormalisation (i.e., energy) scale of the running MS quark mass at its own value.
For the purpose of this conversion, starting from the RGI mass value, we first use the result Λ ; for the numerical evaluation, we proceed along the lines of appendix A of ref. [82]. Next, at the energy scale defined by the scale invariant charm quark mass, one has to employ (perturbative) decoupling relations, which account for the effect of the change in the number of active quark flavours when passing across the respective quark mass thresholds while running with energy, in order to consistently arrive at a result in the four-flavour theory; for this step, we utilise the RunDec software package described in refs. [83][84][85]. Then one can perturbatively run this number in the four-flavour theory, either to the scale invariant quark mass (which gets slightly shifted by the decoupling process) or upwards to the conventional energy scale µ = 3 GeV. In this way we provide results, for both 4-loop and 5-loop perturbative accuracy of the β-function and the quark mass anomalous dimension, 10  where the first error is the one from our determination discussed in the main text, whereas the second error stems from the uncertainty of Λ MS . After the switch to the four-flavour theory, we may employ perturbative running again in order to eventually obtain the four-flavour RGI charm quark mass, for which we quote    [86] (black triangle), the N f = 2 + 1 + 1 and N f = 2 + 1 FLAG averages [18] (grey vertical bands), as well as the lattice results that entered these averages: HPQCD 18 [16], FNAL/MILC/TUMQCD 18 [15], HPQCD 14A [14], ETM 14A [13], ETM 14 [12], JLQCD 16 [9], χQCD 14 [8], HPQCD 10 [7]. The upper part of the figure displays lattice results originating from simulations with four quark flavours in the sea (marked by diamonds), whereas the lower part refers to those obtained on gauge configurations with three dynamical quark flavours (marked by squares). and the 5-loop value already stated in eq. (5.2). Note that the error contribution from Λ (3) MS becomes negligible in the four-flavour RGI limit and that the 4-loop and 5-loop results agree up to less than 1 MeV. We thus conclude that any ambiguities, which arise from the perturbative running and decoupling in the MS scheme, are much smaller than the total error of our final result.
In figure 6 we compare our result for m MS c (m MS c ) in the four-flavour theory with those from other lattice QCD determinations that entered the averages in the FLAG review in ref. [18] based on gauge field configurations with 2+1 and 2+1+1 quark flavours in the sea. Except for the result from [13], labelled 'ETM 14A' in figure 6, ours is compatible with all of the shown numbers within 1σ uncertainties. Given that the underlying strategies to extract the physical quark mass, as well as the discretised actions, renormalisation procedures and further potential systematics, differ significantly among these calculations, this overall agreement appears very reassuring. For completeness, we mention that our result is in line with the findings of ref. [17] for 2 + 1 + 1 quark flavours and the result of ref. [11] (which is an extension of ref. [10]) for 2 + 1 flavours. Both references have been published after the closing date for the FLAG review and do not show tension with the results in figure 6.
Let us close this appendix with a general remark on the renormalisation strategy which we follow in this work. Indeed, one might argue that it is not ideally suited when one is sim-

JHEP05(2021)288
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.