Hadronic physics from a Wilson fermion mixed-action approach: Charm quark mass and $D_{(s)}$ meson decay constants

We present our first set of results for charm physics, using the mixed-action setup introduced in a companion paper. Maximally twisted Wilson valence fermions are used on a sea of non-perturbatively $O(a)$-improved Wilson fermions, made up by CLS $N_{\mathrm{\scriptstyle f}}=2+1$ ensembles. Our charm-sector observables are free from $O(am_c)$ discretisation effects, without need of tuning any improvement coefficient, and show continuum-limit scaling properties consistent with leading cutoff effects of $O(a^2)$. We consider a subset of CLS ensembles -- including four values of the lattice spacing and pion masses down to 200 MeV -- allowing to take the continuum limit and extrapolate to the physical pion mass. A number of techniques are incorporated in the analysis in order to estimate the systematic uncertainties of our results for the charm quark mass and the $D_{(s)}$-meson decay constants. This first study of observables in the charm sector, where the emphasis has been on the control of the methodology, demonstrates the potential of our setup to achieve high-precision results.


Introduction
Heavy flavour physics is a key frontline in the endeavour to test the limits of the Standard Model, and look for new fundamental physics.Ever-increasing precision for fundamental parameters such as quark masses and Cabibbo-Kobayashi-Maskawa matrix elements, as well as for weak matrix elements that control the low-energy hadronic contribution to weak decay amplitudes, is necessary to keep pace with experimental developments.
First-principles, systematically improvable computations performed in Lattice QCDpossibly, beyond a certain precision threshold, with QED corrections -are of course the basic source of input.When dealing with heavy quark physics, however, lattice computations face a non-trivial multiscale problem.Since computations involve both an ultraviolet cutoff -the inverse lattice spacing a −1 -and an infrared cutoff -the inverse size L −1 of the finite box computations are performed in -all physical scales should best fit comfortably between the cutoffs, lest control on their removal is compromised.A standard criterion for finite-volume effects to be sufficiently suppressed in typical hadronic quantities involves a constraint on the box size m π L ≳ 4; for the typical range of pion masses explored, which nowadays routinely reaches the physical point, this implies box sizes in the 3 to 7 fm ballpark.Having m c ≪ a −1 , and especially, m b ≪ a −1 , then requires values of the lattice size L/a that are close to or simply beyond current computational capabilities.This problem is much worsened by the extra difficulty to approach the very fine lattice spacings needed to accommodate heavy quark masses: the computational cost of typical simulations scales as ∼ a −7 [2], and for a ≲ 0.05 fm the algorithmic problem of topology freezing sets in, which in practice impedes simulations long enough to control statistical uncertainties reliably [3].
Facing these problems requires a specific toolset for heavy quark physics on the lattice, that, in particular, relies on input from effective theories to try and control the ultraviolet cutoff dependence: Symanzik effective theory [4][5][6][7] to understand and suppress the leading cutoff effects, heavy quark effective theory input to guide the construction of lattice actions or the extraction of physics, 1 etc.The resulting sophisticated frameworks often rely on assumptions about the systematic impact of the use of effective theory, and/or require the determination of ancillary quantities such as O(a) improvement coefficients.A full overview of lattice techniques and results for heavy quark physics can be obtained from the latest FLAG review [9].A main theme underpinning all studies in the charm and, especially, the B sector is that having results from a variety of approaches is essential to gain confidence on the systematic uncertainties affecting hadronic observables relevant for flavour physics.
The main motivation of the mixed-action setup used in this work, and fully discussed in [1], is to devise an optimal framework for heavy quark physics that bypasses many of the difficulties mentioned above.The first ingredient is the use of CLS N f = 2 + 1 ensembles obtained with non-perturbatively O(a) improved Wilson fermions [10] and open boundary conditions for the gauge field [11,12], which allows to enter the realm of very fine lattice spacings while keeping control on statistical uncertainties.The second ingredient is to compute heavy quark observables by means of a valence twisted-mass Wilson setup [13,14], which leads to automatic O(a) improvement [15].Working with a mixed action of course leads to new requirements, such a precise matching between the valence and sea sectors, and a careful analysis of the relative cutoff effects.This is discussed in the companion paper [1].Here we will focus on illustrating how the technique is able to obtain precise, reliable results for basic observables in the hadronic sector.Progress report of this long-term project have been given in [16][17][18][19][20][21][22][23].
In particular, we will focus on determining the value of the charm quark mass, and of the leptonic decay constants of the D and D s mesons.Our results are based on a subset of the available CLS ensembles which allow us to illustrate the properties of the setup.We also emphasise our development of a variant of the existing techniques to assess systematic uncertainties in lattice observables based on information criteria [24,25] applied to appropriate goodness-of-fit estimators [26].Still, despite the fact that our current results use a subset of the CLS ensembles, they are already at a point where they have competing precision in the context of the state-of-the-art determination of these quantities that enter current FLAG averages [27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45].Results with a larger set of CLS ensembles, including finer lattice spacings and physical pion mass ensembles, will be the object of upcoming publications.
Let us conclude this introduction by describing the organisation of the paper.Sec. 2 summarises the main aspects of our mixed-action approach, discussed at length in [1].Sec. 3 deals with our approach to matching the scale of our partially-quenched charm quark, and numerical aspects of computations in the charm sector.Secs. 4 and 5 discuss our determination of the charm mass and decay constants, respectively.Finally, Sec. 6 contains our conclusions and outlook.

Mixed-action setup
In this section we review the basic features of our setup, with an emphasis on their implications for heavy quark physics.We refer the reader to [1] for a fully detailed discussion of our approach.

Generalities
All our results stem from a mixed-action setup.In the sea sector we employ a tree-level improved gauge action [6,46], and a non-perturbatively O(a)-improved Wilson fermion action [47].This has indeed been used in the generation of the CLS N f = 2 + 1 ensembles [10,[48][49][50] that we employ.In the valence sector, on the other hand, we use a fully-twisted tmQCD [13] fermion action.Both actions include the same massless Wilson-Dirac operator [47,51] where ∇ µ and ∇ * µ are, respectively, the forward and backward covariant derivatives, and Fµν is the clover-leaf definition of the field strength tensor as spelled out in [7].
The mass term in the sea has the form while the tmQCD action is obtained by adding a mass term of the form [13,14] where m cr is the standard Wilson critical mass, and the signs have been set so that the values of the twisted masses µ f are implied to be positive.We will always work in the isospin limit, where the up and down quark masses take the same values both in the sea and in the valence (i.e., m d and µ u = µ d ≡ µ l ).The procedure to fully define the mixed action involves the matching between Wilson and tmQCD valence actions, and a specific prescription to define the critical mass used in our setup.To that purpose, for any given ensemble we first tune µ l , µ s and m cr such that the quantities ϕ 2 and ϕ 4 -depending on pion and kaon masses, as defined in Eq. (2.6)coincide for sea and valence actions, while imposing that the (u,d) standard PCAC quark mass -including all known O(a)-improvement counterterms -vanishes in the valence sector.This ensures equivalent physics and sets the twist angle to π/2, ensemble by ensemble.

Properties of the twisted valence sector
The most interesting property of this setup for the purpose of the results presented in this paper is that it results in automatic O(a) improvement of observables extracted from valence correlation functions [15], up to terms proportional to the trace of the subtracted sea quark mass matrix, atr{m (s) q } [1].Since the latter only involves up, down, and strange quarks, the value of the trace in lattice units is of O(10 −2 ) on our ensembles.Furthermore, these terms arise from loop effects, and their coefficient is thus formally of perturbative order α 2 s .Given our typical statistical uncertainties, the natural size of these atr{m (s) q } lattice artefacts therefore amounts to a subdominant contribution.This property can be furthermore verified a posteriori by inspecting the scaling of observables towards the continuum limit.This is very important for heavy quark observables, since we are then assured that the leading cutoff effects associated to a quark of mass µ h are of order (aµ h ) 2 , without need of fine-tuning improvement coefficients to ensure the elimination of linear effects, as would be the case in the standard O(a) improved setup.
Note that automatic O(a) improvement holds even in the absence of the clover term in the valence fermion action; we have however kept it for a number of reasons.First, it simplifies the matching between sea and valence, since the regularisations coincide in the chiral limit.Secondly, for the same reason, it allows to use non-perturbative renormalisation constants determined with standard methods -e.g., to obtain renormalised quark masses [53].Finally, it has been observed that keeping the clover term leads to a better control over the O(a 2 ) flavour-breaking effects induced by the twisted mass term, thus improving the overall scaling of the setup [54,55].
A second, more generic benefit is that the use of a twisted mass regularisation implies multiplicative renormalisation of (twisted) quark masses, and the possibility to determine decay constants without need of finite normalisation factors such as Z A .This is a result of the explicit chiral symmetry breaking pattern at full twist, which leaves exactly conserved axial currents.In the twisted quark field basis implicitly assumed when writing our valence mass terms, the relevant on-shell (x ̸ = 0) Ward-Takahashi identity reads where ∂ * µ is the backward lattice derivative; O is any gauge-invariant local operator; µ q,r are the Lagrangian twisted masses for the corresponding flavours q, r, that are here assumed to carry different signs in the twisted mass matrix µ of Eq. (2.3);3 P qr = ψq γ 5 ψ r is a non-singlet pseudoscalar density; and Ṽ qr µ is the point-split vector current4 Since the current is exactly conserved, there are two important consequences.First, current and Lagrangian quark masses coincide, and renormalise with Z µ = Z −1 P . 5Second, meson decay constants can be extracted from a two-point function of the pseudoscalar density, by setting O = P rq in Eq. (2.4) and using the fact that the l.h.s. of the Ward identity is exactly normalised.These will be the basis of our determinations of the charm quark mass in Sec. 3 and of f D (s) in Sec. 5.

Ensembles and line of constant physics
CLS ensembles have been generated along three different lines of constant physics.Our results are based on a subset of the ensembles generated at (approximately) constant value of tr{m (s) q }, which we list in Table 1.In order to define a precise line of constant physics, we use the quantities where t 0 is the gradient flow scale introduced in [57], and whose value in physical units has already been determined using CLS ensembles in [1,48,58,59].A renormalised line of constant physics can thus be fixed by setting ϕ 4 constant and equal to its physical value; extraction of the physics will then proceed by a combined continuum-chiral limit fit that hits the physical value of ϕ 2 .The condition that ϕ 4 is constant is well approximated by keeping tr{m (s) } fixed, since it is proportional to ϕ 4 at leading order in the effective chiral description of QCD dynamics.Small deviations from the correct value of ϕ 4 in each ensemble can be corrected by means of the mass shifting prescription introduced in [48], and incorporated into the fitting procedure -see [1] for technical details.Our renormalised chiral trajectory is ultimately set at ϕ phys 4 = ϕ isoQCD 4 = 1.101(7) (5), where the second error quoted is the systematic uncertainty coming from our Bayesian model averaging (see below), and the first error comprises the statistical uncertainty, the one associated to chiral-continuum extrapolations, and those related to input parameters -improvement coefficients, renormalisation constants, and the input pion and kaon masses.The values of the latter employed to fix ϕ 4 are those in the QCD isospin symmetric limit (isoQCD) given by [9] m isoQCD π = 134.9768(5)MeV, (2.7) In the remainder of this paper we will use the superscript "phys" for quantities defined in the isoQCD prescription for the continuum theory, as fixed above.
Id In this work we employ our determination of the physical scale from the gradient flow scale t 0 .To set the scale, we use the following combination of pion and kaon decay constants (2.9) At NLO in SU (3) χPT, this quantity remains constant up to logarithmic terms.From the chiral-continuum extrapolated value of √ 8t 0 f πK we eventually extract the flow scale t 0 in physical units by using as physical inputs the isoQCD values for f π and f K .Specifically, we use [9] f isoQCD π = 130.56 (13) MeV, (2.10) The full details of our scale setting procedure through a combination of the O(a)-improved Wilson results with the ones from the valence Wilson Twisted Mass regularisation can be found in [1].The resulting value of t 0 we will use to convert our results to physical units is where the uncertainty is split in the same way as described above for ϕ phys 4 .

Charm correlators and scale setting
In this section we discuss the technical details behind the computation of physical observables in the charm region from our mixed action setup.We introduce the GEVP setup used to extract meson masses and matrix elements throughout this work and explain our strategy to match the charm quark mass to its physical value.

Computation of correlation functions
To extract physical observables we have measured two-point correlation functions at zero momentum on CLS gauge configurations listed in Table 1.Fermionic two-point correlators have the form where y 0 and x 0 are, respectively, the source and sink time coordinates; q and r are flavour indices; and a trace over spin and colour is implicit.O q,r Γ are quark bilinear operators defined as where Γ is a spin matrix.The operator content will be denoted by subscripts in straightforward notation -we will refer to f PP when Γ = Γ ′ = γ 5 , f AP when Γ = γ 0 γ 5 and Γ ′ = γ 5 , and so on.In all computations in this work we have fixed the time position of the source at y 0 = T /2, to maximise the distance from the boundaries: when dealing with heavy-light and heavy-heavy flavour content in the operators O q,r Γ in Eq. (3.2), we observe that the region in which the signal for the considered two-point function is accessible lies entirely within the lattice bulk, and that the boundary effects are strongly suppressed.Ten time-diluted U (1) stochastic sources are employed in the computation of the quark propagators in each gauge field configuration.Moreover, the numerical inversion of the quark propagator in the charm region is performed using distance preconditioning techniques [60,61], in order to reduce signal deterioration and enhance accuracy at large Euclidean times.Error analysis and propagation are based on the Gamma method of [62] and automatic differentiation, as implemented in the ADerrors package [63].
Light and strange propagators are computed at the values of m cr , µ l and µ s determined to ensure maximal twist and pion and kaon masses matched to the sea (see Section 2).We note that this is a independent set of computations of the propagators with respect to those employed in the matching procedure [1], where a grid of values for the mass parameters is employed to accurately interpolate to the matching point.Moreover, this grid was also employed to compute the mass corrections to the renormalised chiral trajectory [1].Heavy propagators are computed at three different values of the twisted mass µ (i) c around the physical charm region (save for one ensemble where only two masses have been used), so that observables are interpolated at the physical value of the charm quark mass.In Table 2 we specify the twisted mass values and the critical hopping parameter κcr used to impose the maximal twist condition for each ensemble used in the analysis.

Extraction of meson masses
In our analysis meson masses are employed to fix the renormalised line of constant physics and match the quark masses to some target physical value.Light and strange quark masses are matched between the sea and valence sectors using ϕ 2 and ϕ 4 in Eq. ( 2 1.The critical value of the hopping parameter required to set the valence sector to maximal twist [1] is denoted by κcr .The values of aµ l and aµ s are the light and strange bare twisted quark masses, in lattice units, that match the corresponding sea quark masses [1].Finally, the last three columns contain the three values of heavy bare twisted quark masses in the charm region.In the case of the D200 ensemble two values that straddle the charm point were considered. partially quenched charm quark we use different combinations of mesons masses matched to their physical values, as explained in Section 3.3.The ground state meson masses are extracted from a generalised eigenvalue problem (GEVP) variational method defined as where C(t) is a matrix of Euclidean correlation functions of the form in Eq. (3.1), such that the indices i, j in C ij (t) correspond to different choices of Γ, Γ ′ and source/sink location, and t = x 0 − y 0 .This leads to the spectral expansion Here N denotes the matrix dimension, and we have assumed non-degenerate energy levels.The GEVP is solved in the regime t ref ≥ t/2, where a better control over excited state contributions is achieved [64].The matrix C(t) in our setup is built from pseudoscalar two-point functions f PP shifted in time as where τ is the value of the time shift.Several values of the time shift have been tested, and we observe a mild dependence on small values of τ for the extraction of eigenvalues and eigenvectors.We refer to Appendix A for a detailed discussion of our setup, together with sanity checks on the GEVP.In what follows we set τ = 3a.The ground state meson mass is extracted from the eigenvalues of the GEVP using Eq.(A.1).In order to assess the systematic effects and correctly identify the plateau region, we perform several uncorrelated χ 2 fits to a constant, by varying the time ranges of the fitting interval.Correlated fits are impractical due to the fact that sample covariance matrices display very small modes and thus have ill-behaved inverses.However, as the data is correlated, the uncorrelated χ 2 is not a suitable quantity to assess the goodness-of-fit; we therefore quantify the latter with the expectation value of χ 2 , denoted χ 2 exp , and the corresponding p-value, as introduced in [26].Through this procedure we assign a weight to each fit based on the χ 2 minimisation, and we eventually extract our ground state masses by means of the model averaging procedure described in Appendix B. An example of a GEVP plateau for the heavy-light pseudoscalar mass together with a summary of the model average procedure for an ensemble used in the analysis is shown in Figure 1.

Matching of the charm quark mass
In Section 2 we recalled the matching of the light sector worked out in [1], which ensures that physical observables involving only light and strange quarks computed in the valence and sea sectors coincide up to cutoff effects, so that unitarity is recovered in the continuum limit.A similar procedure is needed for the charm quark, designed to ensure that its physical value is obtained upon taking the continuum limit and performing chiral fits.Since the charm is partially quenched this matching procedure involves observables with only valence charm quark propagators.
In order to establish a connection with the physical point, we require that some charm-like observable O c matches its physical value.In this paper we studied three different charm scale settings based on three choices of O c , all in terms of meson masses; we will denote the latter as m (i) H , i = 1, 2, 3, and often express them in units of √ 8t 0 as ϕ H .The first possibility, corresponding to ϕ H , consists in using the flavour average meson mass combination m (1) built from heavy-light H and heavy-strange H s pseudoscalar meson masses with heavy-quark masses in neighbourhood of the charm.Since we require the considered CLS ensembles to hold a constant value of the flavour average combination of pion and kaon masses -denoted as ϕ 4 in Eq. (2.6) -we also expect the flavour average combination ϕ H to remain fairly constant along the chiral trajectory.The physical value of m   Our second strategy, corresponding to ϕ H , is to consider the mass-degenerate pseudoscalar meson mass m conn η h extracted from the quark-connected two-point correlation function made of heavy quark propagators with a mass in the neighbourhood of the charm mass, m The physical value for this mass, m (2),phys H , is set from the experimental value of the η c meson mass [65], m exp ηc = 2983.9(4)MeV, from which a correction of about 6 MeV, with 100% error, is subtracted to account for the absence of quark-disconnected diagrams and QED effects [43,[66][67][68][69]. Specifically, we employ, One potential advantage of this choice of matching observable is that the overall precision of the η conn c meson mass is substantially better than the one for heavy-light meson masses, as it does not suffer from the increase in noise-to-signal ratio with Euclidean time; this is illustrated in Figure 2, where we show the D, D s and η conn c pseudoscalar correlators for a one specific ensemble.Finally, as a third matching quantity we also tested the spin-flavour averaged mass combination m which involves a combination of heavy-light pseudoscalar m H (s) and vector m H * (s) meson masses in the charm region, and is motivated by heavy-quark symmetry.However, we observe that chiral-continuum fits coming from the spin flavour-averaged matching condition lead to worse χ 2 values, and as a result their weights are highly suppressed by our model average prescription.We interpret this finding as a reflection of relatively poor control of heavy-light vector states, whose masses are extracted with significantly larger errors than those of heavy-light pseudoscalar states.In the rest of the discussion we will therefore focus on the results coming from the other two matching conditions.
Any of these matching conditions can in principle be imposed ensemble by ensemble, even away from the physical point.However, by doing so we would as a result build in the charm quark mass a dependence on the value of the reference scale t phys 0 , as well as O(a 2 ) effects coming from the specific choice of O c .To avoid this, we have opted instead for setting the physical charm quark mass jointly with the chiral-continuum extrapolation, in a similar way as the one we employ to hit the physical point in the light and strange sector.What this means in practice is that the charm quark mass dependence of any given observable O is parameterised as O(a, ϕ 2 , ϕ ).This will be the procedure applied below in the determination of the physical value of the charm quark mass and of the decay constants f D and f Ds .
Note that, as a consequence of our matching procedure and of working on a line of constant physics where ϕ 4 is kept constant, it is non-trivial that by adopting any of our matching procedures the mass of any particular meson reaches its physical value in the chiral-continuum show the results of the highest weight fit contributing to the model average procedure and the corresponding plateau interval.We observe the expected increase of the statistical uncertainties at large time separations when increasing the mass-difference among the quarks propagators of the pseudoscalar two-point correlators.
limit; checking that it does is therefore a test of the robustness of our procedure.As an illustration, we show in Fig. 3 how the physical values of the D and D s meson masses arise when the charm scale is matched through either m D or m conn ηc .In either case we show results for the specific model of the lattice spacing, charm mass and pion mass dependence of the form where i = 1, 2 according to the notation introduced above and where c 1 and p j , j = 1, 2, 3, stand for the fit parameters.Note that the agreement is excellent, in spite of the different implications of the two setups for the specific case of m D (s) ; for instance, when m D is used for the matching cutoff effects are very small by construction, while the use of m conn ηc leads to sizeable cutoff effects which are however very well described by an O(a 2 ) term.

Renormalised charm quark masses
In Sec. 2 we have summed up the argument why renormalised quark masses can be easily retrieved from bare Lagrangian twisted masses.In our mixed-action setup, as discussed in detail in [1], the resulting O(a)-improved expression for the renormalised charm mass m c (µ) reads where Z P (g 2 0 , aµ) is a suitably defined renormalisation constant for the non-singlet pseudoscalar density at renormalisation scale µ.As we have already discussed, the improvement term H = m H , while those on the right use the mass-degenerate pseudoscalar meson mass, m H = m conn η h .The empty symbols correspond to the D (s) meson masses determined on a given ensemble, while the red square symbols show the extrapolated values at the physical point.Dashed lines show the fit forms projected to each individual lattice spacing, and the blue shaded bands are a projection to the continuum limit on the chiral plane.Data points are projected to the physical point ϕ (i),phys H . Finally, the green horizontal band shows the isoQCD input values for the corresponding masses in Eq. (3.7), in units of √ 8t 0 .
∝ tr{m (s) } can be neglected in practice, so O(a)-improved renormalised quark masses can be obtained by just applying the renormalisation constants to the exactly known Lagrangian masses.
In this work we will use the non-perturbative values of Z P computed in [53] in the Schrödinger Functional scheme, at a fixed renormalisation scale µ had = 233(8) MeV and for the range of values of g 2 0 covered by CLS.It will be used to obtain renormalised quark masses for each of our ensembles, that can then be used to determine the value of the charm quark mass in the continuum and at physical kinematics.Contact with other renormalisation schemes can then be made by computing the renormalisation group invariant (RGI) quark mass M RGI c , using the continuum (flavour-independent) ratio also computed in [53] M m(µ had ) = 0.9148 (88) .
Values of renormalised masses in, say, the MS scheme can then be obtained by using the perturbative value of m(µ) M at any convenient scale µ.

Charm quark mass chiral-continuum fits
Having determined the renormalised charm quark masses in the Schrödinger Functional scheme at the hadronic renormalisation scale µ had for all the ensembles listed in Table 1, we now describe our strategy to obtain results in the continuum limit and at the physical point, following the approach outlined in Sec. 3. The matching procedure of the light and strange sectors is already devised so that the physical value of the kaon mass is recovered at , where the physical value of ϕ 2 is computed with the isospin-symmetric values of the pion mass quoted in [9], and the physical scale t phys 0 is the one determined in [1].The charm scale is matched through the two different prescriptions described in Sec. 3.All quantities entering the fit are made dimensionless through the appropriate power of the factor √ 8t 0 , and physical units for the final result are restored by using our value for t phys 0 .We parameterise the continuum dependence of the renormalised charm quark mass on ϕ 2 and any of the ϕ Based on the heavy quark effective theory expansion [70] at lowest order, we expect a linear dependence of the charmed meson masses as a function of the the charm quark mass, hence the latter term in the ansatz.This assumption is supported by our data that show indeed a linear behaviour in the charmed meson masses, as illustrated in Figure 4.Note that this form is used only to describe the dependence within a short interval in mass values, and interpolate the charm scale from points close by.When considering the pion dependence of the charm quark mass, we assume that the leading order contributions exhibit a linear behaviour in ϕ 2 .With the current set of ensembles employed in this work we do not observe any deviations from the leading order term in the pion mass dependence.
Regarding the lattice spacing dependence of the charm quark mass, we assume the leading cutoff effects to be O(a 2 ), as discussed above.Corrections of odd order in a are generically expected to be highly suppressed at maximal twist, by way of the extension of the argument for automatic O(a) improvement; we thus include a 4 terms to account for deviations from linear behaviour in a 2 .Finally, we allow for terms proportional to m 2 π and to various powers of the charm mass.The generic ansatz to parameterise lattice spacing dependence thus take the following form In order to estimate the systematic effects arising from the model variation, we consider all the possible combinations where some of the c i coefficients vanish, save for c 1 which is always kept.Furthermore, following [45], we allow for cutoff effects to enter either linearly or non-linearly, viz., H in the chiral-continuum extrapolation have been projected to the physical point.The red square symbols indicate the continuum results at the physical value ϕ phys H .We observe a linear dependence of the charm quark mass on the different matching conditions used in this work.
We thus end up with a total of 64 functional forms for each of the two charm matching conditions, i.e., a total of 128 models.Fit parameters are estimated minimising an uncorrelated χ 2 where, however, the covariance between the independent variables and the data is taken into account.As previously discussed, the goodness-of-fit of fit can still be obtained in this case H the η conn h matching prescription in Eq. (3.9).The last column reports the combined result from these two matching procedures according to our model average prescription.The first error is statistical, while the second is the systematic uncertainty arising from the model variation.
from the measurement of the χ 2 exp and the associated p-value.The TIC result for each model is then fed into the model averaging procedure summarised in App.B, which finally allows to quote a systematic uncertainty that reflects the fluctuations engendered by the variety of fit ansaetze.
In Table 3 we report the results for µ R c in units of √ 8t 0 obtained with each of the two matching conditions independently, as well as for the combined model average.
In Figure 5 we summarise the model average procedure, showing some of the best fit results coming from the functional forms defined in Eq. (4.6) for the two matching conditions studied in this work.Each circle corresponds to a result coming from a particular model, and the opacity is associated to its weight determined from our Takeuchi's Information Criterion (TIC) as explained in App.B. We observe that for both matching conditions the majority of the models with relevant weights nicely agree, and as a result the systematic error is subleading with respect to the statistical uncertainty.Figure 6 shows a weighted histogram of our results coming from different fits.We observe that models cluster mainly around two values, which are adequately covered by our quoted systematic uncertainty.
Figure 7 illustrates typical fits for each of the matching conditions, chosen among those with higher weights according to the TIC prescription.The plot shows the continuum limit behaviour of the charm quark mass in units of √ 8t 0 .Results coming from the two matching strategies perfectly agree in the continuum, in spite of displaying a qualitatively different structure in cutoff effects.We observe a scaling of the charm quark mass in reasonable agreement with the O(a 2 ) leading order, confirming the automatic O(a)-improvement of our setup; nevertheless, we notice that given the current statistical accuracy, fits with O(a 4 ) terms are the preferred ones from the model average, since they allow to properly describe the curvature in our data.Note also the overall small size of scaling violations, which are at the few percent level.Finally, Figure 8 shows the pion mass dependence of the charm quark mass.As expected, we observe a mild dependence of the charm mass on the light quark masses.

Results for the charm quark mass
The renormalised charm quark mass µ R c can be obtained once we combine the results collected in Table 3 with our determination of t phys 0 in Eq. (2.12).As discussed at the beginning of this section, the knowledge of the renormalisation group running factors allows to quote results for the RGI and MS values of the charm quark mass.c in units of √ 8t 0 .We collect a subset of the best results according to the TIC procedure, coming from different models, for the flavour-averaged matching condition ϕ H .The opacity of each circle data point reflects the associated normalised weights W as given from the TIC.The yellow shaded band represents the systematic error computed with Eq. (B.4), while the left-most red square symbol corresponds to the result extracted from the model average procedure.The labels of the 32 models specified in the horizontal axis are related to the terms appearing in Eq. (4.5)characterising the lattice spacing dependence -in the following way: 'a2' corresponds to the term depending on the fit parameter c 1 .Similarly, 'a2l', 'a2h2', 'a4', 'a4h2', 'a4h4' refer to c 2 , . . ., c 6 , respectively.Given that the parameter c 1 is included in all the models, the associated label is not explicitly specified for all cases appearing in the horizontal axis.
After combining the results from our 128 fitting models through the model average procedure, and using the running factor in Eq. (4.2), we quote for the three-flavour theory the value for the RGI quark mass where the first error is statistical, including the uncertainty on t phys 0 , the second accounts for the systematic uncertainty, derived from the model average, the third is the error contribution from the RGI running factor in Eq. (4.2), and the last error in brackets is the total uncertainty.Figure 9 illustrates the relative contribution of various sources of error to the uncertainty of our determination of M RGI c .The dominant source of error comes from the renormalisation group running of Eq. (4.2), while the second most relevant contribution arises from the statistical error of the correlation functions computed in each ensembles.The error coming from the uncertainty on t phys 0 based on our scale setting procedure [1], as well as the systematic error from the model average are subleading contributions.We therefore expect that the inclusion in  this charm quark mass analysis of further ensembles -with finer lattice spacings and at physical pion masses -will only have a significant impact if combined with improved determinations of the RGI running factor and the scale setting procedure.
In order to quote results in the MS scheme, we use five-loop perturbation theory for the quark mass anomalous dimension [71][72][73] and the beta function [74][75][76].The matching between the N f = 3 and N f = 4 theories uses the four-loop decoupling effects [77] incorporated into the RunDec package [78][79][80].Renormalisation group equations are solved using as input the value Λ (3) MS = 341 (12) MeV from [81].The correlation arising from the fact that a common subset of gauge field configuration ensembles were employed in the computation of Λ (3) MS and the non-perturbative running factor in Eq. (4.2) is taken into account.We thus arrive to the following results for the RGI and MS-scheme charm quark masses in the 4-flavour theory, [13] GeV , (4.9) 11) Λ (5) trunc.[16] GeV , (4.10) where the first and second errors arise from the statistical and systematic errors, respectively, in the value of M RGI c (N f = 3) in Eq. (4.7), the third error is due to the non-perturbative running factor in Eq. (4.2), the fourth error is related to the uncertainty in Λ (3) MS , the fifth error is an estimate of the truncation uncertainty from the deviation between the 5-loop and 4-loop results, and the last error in brackets is the total error.We observe that at the lower renormalisation scale, µ = m c , the scale invariant MS charm mass, m c (µ = m c , N f = 4), Comparison of the continuum limit approach for the two charm matching prescriptions.Shown are two of the fits with highest weights from the TIC, projected onto the lattice spacing dimension.In yellow we show results for the η conn h matching condition, while the blue points illustrate the flavour-averaged matching.Each data-point in this plot is projected to the physical pion mass and the physical charm quark mass, in order to properly visualise the lattice spacing dependence.receives a large contribution to its error from the uncertainty of Λ In Figure 10 we compare our determinations of the charm quark mass in the MS scheme with the results from other lattice QCD calculations also based on N f = 2 + 1 dynamical simulations and with the corresponding FLAG average [9].We observe in particular a good agreement with the results from [45] which are also based on CLS ensembles but employ Wilson fermions in the valence sector.

Computation of decay constants
Along with the charm quark mass, in this paper we present a first computation of the D (s) meson decay constants within our setup.In the absence of electromagnetic interactions, the decay constant fully determines the leptonic decay amplitude of flavoured pseudoscalar mesons, and is given by the matrix element of the axial current as  .The dominant piece comes from the error in the non-perturbative determination of the renormalisation group running factor to the RGI mass quoted in Eq. (4.2).The label statistical plus χ-continuum limit stands for the error arising from the statistical accuracy of our data and the chiral-continuum extrapolation, while the scale setting piece comes from the physical value of the gradient flow scale t phys 0 .Finally, the model average piece illustrates the systematic error arising from the set of models considered in this work.
where the state |P qr ⟩ is the ground state for a pseudoscalar meson with flavour content qr, and m PS its mass.The factor 1/ √ 2m PS L 3 comes from the usual relativistic normalisation of one-particle states in finite volume.
With Wilson fermions, the computation of the above matrix elements requires the finite current normalisation factor Z A and, if O(a) effects are to be subtracted, a number of improvement coefficients.With our fully twisted valence sector this is completely bypassed: when qr belong in a twisted quark doublet -i.e., have different signs in the twisted mass matrix in Eq. ( 2.3) -the physical axial current, expressed in twisted quark variables, becomes a vector current, and the Ward identity in Eq. (2.4) allows to obtain it from the pseudoscalar two-point function.The resulting expression of the correctly normalised pseudoscalar decay constant reads We will extract the matrix element ⟨0|P qr |P qr (p = 0)⟩ from the normalised eigenvector v n (t, t 0 ) of the GEVP according to Eq. (A.3).In order to extract the large time plateau where excited state contributions are suppressed we perform several fits to constant behaviour by varying the fit ranges, and we assign a weight to each fit by means of the TIC prescription as described in App.B. The results for the ground state matrix element are then extracted through the model MS .Left: comparison for the m c (µ = 3 GeV, N f = 4).Right: comparison for m c (µ = m c , N f = 4).Starting from the bottom, results are taken from: PDG [65], HPQCD 08B [82], HPQCD 10 [27], χQCD [35], JLQCD 16 [37], Maezawa 16 [83], Petreczky 19 [42], ALPHA 21 [45].average given by Eq. (B.3).In Figure 11 we show a representative plateau for a heavy-light decay constant, together with a summary of the model average with different fit intervals.

Chiral-continuum fits and results for f D (s)
The chiral-continuum fits for the D (s) meson decay constants are performed similarly to the ones for the charm quark mass.By exploiting Chiral Perturbation Theory with heavy quarks [84,85] to construct appropriate fit functions, we extract the physical point observables trough a global fit of the f D and f Ds decays, and estimate the systematic effects by applying the model average procedure based on the TIC.
The quantities we fit to are combinations of meson masses and decay constants of the form for which a Heavy Quark Effective Theory (HQET) scaling law in powers of the inverse heavy quark mass exists.The general continuum heavy and light quark mass dependence can be expressed as the product of the individual contributions to arrive at the generic expression Here Φ χ governs the heavy-quark mass dependence while δΦ χPT controls the light quark behaviour as approaching the physical point.Finally the lattice spacing dependence describing cut-off effects is regulated by δΦ D (s) a .In the following, we analyse these terms independently to arrive at a final expression for the Φ D (s) approach to the physical point.
Models [t min /a, t max /a]  The continuum heavy-quark mass dependence, Φ χ , admits an expression in HQET of the form where ϕ H = √ 8t 0 m H monitors the heavy quark mass dependence with m H being the flavouraverage m H or the η conn h pseudoscalar meson masses.In general, this expression is not expected to have high accuracy in the charm mass region, due to it being at the limit of applicability of HQET.Furthermore, perturbative values for the matching factor C HQET (m h ) have notoriously poor convergence behaviour. 6However, we are not interested in modelling the heavy quark mass dependence in a wide region of masses -we rather want to interpolate to the charm point from the nearby values of the heavy masses we compute at.Therefore, we will simply take an expression with the same functional form for the m h power corrections, and a constant overall coefficient, as a convenient ansatz for the interpolation part of our fits.In HQET terms, this amounts to neglecting the small logarithmic dependence on m h in a short interval of values.
The light quark mass dependence term, following Heavy Meson χPT (HMχPT) considerations, reads [38,85] are fit parameters and g 2 is the H * Hπ coupling in the static and chiral limits, here treated as a free fit parameter alongside p (i) χ .In Eq. (5.7) we introduced the notation for the chiral logarithm corrections Here ϕ 2 and ϕ 4 are the usual hadronic combinations introduced in Eq. (2.6), which control the light and strange quark mass dependence.When working at NLO in the chiral expansion, the term ϕ f appearing in Eq. (5.7), which introduces the χPT scale, is here replaced by the continuum physical value of √ 8t 0 f πK , as determined from our setup [1] at full twist, with f πK given by 7 6 This is readily observed in the expression for the coefficient in the MS scheme [86,87], where, for QCD, γ0 = −4, γ1 = −254/9 − 56π 2 /27 + 20N f /9, while the perturbative coefficients of the β function have their usual values β0 = (11 − 2N f /3) and β1 = (102 − 28N f /3). 7We remind the reader that fπK is the quantity used to extract the physical scale t phys Finally, with similar arguments to the one discussed in the case of the charm quark mass, the lattice spacing dependence δΦ for the observables Φ D (s) can be parameterised as where p (0,1,2,... ) a are fit parameters.To summarise, for the continuum quark mass dependence of Φ D and Φ Ds we adopt the expressions obtained by combining the light and heavy quark dependencies δΦ χPT and Φ χ , respectively.Following Eq. (5.4), this then leads to the final ansatz for Φ D (s) of the form Since many fit parameters are shared between Φ D and Φ Ds , we opt for a global fit for determining the two quantities.Moreover, at the symmetric point, i.e., for those ensembles with degenerate light and strange quark masses µ l = µ s , the two decay constant coincide, and Φ D = Φ Ds .Therefore, a global fit also helps to constrain the parameters at the symmetric point.
Similarly to the case of the charm quark mass, we consider several specific forms of the fit ansatz, by setting some combination of fit parameters to zero.We furthermore again match the charm scale using the two different procedures described in Sec. 3. The result is a total of 57 different models for each matching condition, and we use our TIC criterion to extract a systematic uncertainty associated to the variation within the full set of fits.In this work, our current approach deliberately excludes fits involving cuts in β or pion masses, as with the current subset of ensembles they are significantly penalised by the TIC.As we look ahead to future updates with the complete set of ensembles we will incorporate cuts in the data within our analysis.
In Figure 12 we show the chiral extrapolations for f D and f Ds with larger weights in the model average.From our chiral-continuum extrapolations of Φ D and Φ Ds , we observe a mild dependence on the choice of the ϕ H used to match the charm scale.Therefore, in the Figures we illustrate the flavour-averaged matching condition only.We also notice that Φ D shows some curvature in ϕ 2 arising from the chiral logs, while Φ Ds presents a more linear behaviour while approaching the physical point.Figure 13 shows an illustration of the scaling towards the continuum limit of Φ D and Φ Ds .We observe that the continuum approach is very well H . Dashed lines refer to the mass dependence at finite values of the lattice spacing, while the blue band represents the projection to the continuum limit.Finally, the red square symbols indicate the physical point results.H .The last column reports the result of the combination of these two matching conditions.The first error is statistical while the second is the estimate of systematic uncertainty arising from the model averaging procedure.
described by leading cutoff effects of O(a 2 ), as expected for our valence action when it is tuned to maximal twist.
In Table 4 we show our determinations of Φ D and Φ Ds for each of the two procedures to match the charm scale, as well as the result from their combination.Using this combination we arrive at the following results for the the D (s) meson decay constants, f D = 211.3(1.9)(0.6)MeV, (5.15) f Ds = 247.0(1.9)(0.7)MeV, ( where the first error is statistical and the second the systematic uncertainty from the model average.The error for the D (s) decay constants is dominated by the statistical uncertainty of correlators and the error on chiral-continuum extrapolations.Therefore, we expect that a future addition of other ensembles with finer lattice spacing and physical pion masses will contribute to significantly reduce the uncertainty of our current determination.The different contributions to the variance of D (s) meson decay constants are shown in Figure 14.Finally, in Figure 15 we show a comparison between our results and other N f = 2 + 1 lattice QCD determinations.

Direct determination of f Ds /f D
In addition to the determination of f D and f Ds , we investigate the direct determination of the ratio f Ds /f D from a dedicated fit.This allows for a consistency check, since the ratio is dimensionless and thus does not require normalisation with a reference scale such as √ 8t 0 .One particular consequence is thus that this approach is only indirectly subject to the uncertainty of the lattice scale setting.Another advantage is that the ratio is exactly 1 by construction when m s = m l , i.e., the symmetric point of our ϕ 4 = constant trajectory, which is part of our line of constant physics.We can thus perform a fit that is highly constrained in the unphysical masses region, although at the price of reducing the total number of ensembles entering in the study of the approach to the physical point.
A first set of fit ansaetze is derived from the HMχPT expressions considered above for Φ D (s) .The generic form is  The label statistical plus χ-continuum limit represents the error arising from the statistical accuracy of our data and the chiral-continuum extrapolations.The scale setting label denotes the error coming from the physical value t phys 0 as determined within our setup [1], while the model average represents the systematic error arising from the model variation according to the TIC procedure.
Here δΦ χPT introduced in Eq. (5.7) labels the light quark mass dependence of the ratio, while δΦ D (s) a from Eq. (5.12) controls the continuum approach.It is worth noticing that at leading order the physical dependence on ϕ H , and also the lattice spacing dependence related to ϕ H , cancel out when expanding the ratio.Collecting all the terms entering in Eq. (5.17) from the previous section, we end up with In this expression we consider all the possible combinations of non-vanishing fit parameters, and perform our TIC-weighted model average among the different functional forms tested to quote a systematic uncertainty.Given that various terms cancel in the HMχPT expressions, we will further explore the systematic uncertainties by considering also functional forms based on a Taylor expansion of Φ D (s) .The generic expression then reads where Φ D (s) χ is the value in the chiral limit and at the physical value of the heavy-quark Only data points with filled symbols contribute to the FLAG averages.Starting from the bottom, results are taken from: HPQCD 10 [28], PACS-CS 11 [88], FNAL/MILC 11 [29], HPQCD 12A [30], χQCD 14 [35], RBC/UKQCD 17 [39], χQCD 20A [89].
mass.In this expansion, the heavy and light mass dependence terms read The lattice spacing dependence δΦ D (s) a can be parameterised in a similar fashion to that in Eq. (5.12).Combining these expressions into a functional form for the ratio of decay constants one then has Table 5: Results of the model average for f Ds /f D for the two charm-quark matching conditions.The last column reports the combined result.The first error is statistical while the second is the systematic uncertainty arising from the model variation procedure.simultaneously.In Table 5 we report our results for the ratio of decay constants from the model average separately for each charm matching condition, as well as their combination.Also for the ratio we observe good agreement for the two different ϕ H tested in this work.Finally, for the result combining the two matching conditions, we quote where the first error is statistical and the second is the systematic uncertainty based on the model average procedure.
In Figure 16 we show the HMχPT chiral-continuum fit of the Φ Ds /Φ D ratio with highest weight in the model averaging procedure.In particular the plot on the left shows the chiral approach to the physical point, while the plot on the right represents the lattice spacing dependence.The observed dependence on ϕ 2 shows only a mild curvature arising from the chiral logs, while cutoff effects appear to be highly suppressed at the current level of statistical precision of our data.
Figure 17 shows a summary of the model average procedure for the ratio Φ Ds /Φ D , displaying the fit results for the two matching conditions together with the associated weights, for the HMχPT and Taylor functional forms.
In Figure 18 we show the major error sources contributing to our final determination of the ratio, where we notice that the major contribution is given by the statistical and chiralcontinuum error.Finally, in Figure 19 we show a comparison between our result for f Ds /f D , the FLAG21 average and results from other collaborations.

Conclusions and outlook
In this work we have described our first computations of physical observables in the charm sector using the Wilson fermion mixed-action setup described in greater detail in [1].Emphasis is put in setting up our methodology and exhibiting the characteristics of the framework.Our results for the charm quark mass and the D (s) meson decay constants are based on a subset of CLS ensembles, yet they already sport a level of precision similar to that of several state-of-the-art results.We quote the values H and ϕ H .Each circular symbol represents the result of a specific functional form, and the opacity is associated to the normalised weight W of the model based on its TIC value.The yellow band represents the systematic uncertainty arising from the set of tested models, while the left-most red point is our final averaged result.The labels of the 20 models specified in the horizontal axis are related to the terms characterising the dependencies on the mass and lattice spacing in the following way: 'HMChPT' stands for the expression in Eq. (5.18) where only the leading terms depending on the fit parameters p the model selection, the third arises from the RGI running factor in Eq. (4.2), and the last one in brackets is the total error.For the decay constants f D , f Ds and their ratio f Ds /f D , the first error is statistical and the second is the systematic uncertainty from the model averaging, and the total error is given in brackets.
We foresee that these results could be improved in the future by means of a more extensive analysis including additional CLS ensembles with a finer value of the lattice spacing and physical pion mass simulations.This is expected to have a significant impact in reducing the statistical uncertainty of the decay constants.The error on the charm quark mass, on the other hand, is dominated by the uncertainty induced by the non-perturbative renormalisation group running and thus work on that front would be required to improve the precision significantly.In a related line of work, we are also applying our framework to the computation of semileptonic form factors for charmed meson decay, for which preliminary results have already been presented in [91,92].Together with the computations illustrated in this paper, they show how a comprehensive programme of precision heavy-flavour physics can be pursued in the framework of Wilson fermion regularisations, reaching an excellent compromise between the latter's advantages from the point of view of field-theoretical control and the aim of highprecision computations.[9].Only the results with filled symbols contribute to this average.Starting from the bottom, results are taken from: PACS-CS 11 [88], FNAL/MILC 11 [29], HPQCD 12A [30], RBC/UKQCD 17 [39], RBC/UKQCD 18A [90] Criteria (TIC) proposed in [25] to assign a weight to each model, which then allows to perform a weighted model average to arrive at a final result for the systematic uncertainty [24].Specifically, the value of the TIC assigned to each fitting model is To each model m in the complete set, consisting of M models, we assign a normalised weight W m defined as follows The result of the model average for an observable O that has been determined for each of the models is then given by ⟨O⟩ =

( 1 ) 7 )
,phys H is obtained by setting m H (s) to the following prescription for the isoQCD values of D (s) meson masses, The uncertainties in these isoQCD values are chosen to cover the deviation with respect to the experimental values [65] of the D ± and D ± s meson masses, m exp D ± = 1869.66(5)MeV and m exp D ± s = 1968.35(7)MeV, respectively.We observe that the larger uncertainty in the isoQCD inputs of the D and D s meson masses in Eq. (3.7) -as compared to the corresponding

Figure 1 :
Figure1: Illustration of the extraction of the ground-state mass after applying a GEVP analysis, illustrated for the ensemble J303.Top: Heavy-light pseudoscalar meson mass plateau showing the two fit intervals with higher weights W contributing to the model average.We also indicate the range of variations allowed for the interval in Euclidean time where the plateau is taken.Bottom: Summary of determinations of am H when considering variations over the fit intervals [t min /a, t max /a] together with the corresponding normalised weights W based on Takeuchi's Information Criterion (TIC), p-values and χ 2 /d.o.f..In the upper panel, the shaded blue band corresponds to the model average result.

H
), and we perform a global fit to obtain its physical value O(0, ϕ phys 2 , ϕ (i),phys H

Figure 2 :
Figure2: Illustration of the effective meson masses involved in the matching procedure to the physical charm scale for the ensemble J303.We show three cases where the effective mass of the pseudoscalar meson H is that of the D (left), D s (center ) and η conn c (right), normalised by the central value of the corresponding plateau averaged mass.The horizontal red bands show the results of the highest weight fit contributing to the model average procedure and the corresponding plateau interval.We observe the expected increase of the statistical uncertainties at large time separations when increasing the mass-difference among the quarks propagators of the pseudoscalar two-point correlators.

Figure 4 :
Figure 4: Heavy mass dependence of the renormalised charm quark mass µ R c in units of √ 8t 0 for the fits with larger weights according to the TIC criteria.Top: Results shown for the flavour-averaged matching condition ϕ (1) H = √ 8t 0 m H . Bottom: Results shown for the η conn h

Figure 5 :
Figure 5: Model average procedure for the renormalised charm quark mass µ Rc in units of √ 8t 0 .We collect a subset of the best results according to the TIC procedure, coming from different models, for the flavour-averaged matching condition ϕ

Figure 6 :
Figure 6: Weighted histogram illustrating the model average procedure for √ 8t 0 µ R c .The result from each of the 128 models -including both matching conditions ϕ (1) H and ϕ (2) H -parameterising the lattice spacing dependence is weighted by its normalised weight W based on the TIC.The vertical line represents the central value from the model average, while the vertical band shows the corresponding estimate of the systematic error.

( 3 )
MS and from the truncation error.These specific sources of uncertainty are less prominent in the RGI mass, M RGI c (N f = 4).

. 1 )φ phys 2 aFigure 8 :(N f = 3 )
Figure8: Pion mass dependence of the charm quark mass for one of the best fits according to the TIC criteria.Results are shown for the flavour-averaged matching condition.Each point corresponds to the value for a given ensemble, projected to the physical charm quark mass.The dashed lines represent the chiral trajectories at finite lattice spacing, while the blue shaded band is a projection to the continuum limit.The red point shows our final result extrapolated at the physical point in the continuum.

Figure 9 :
Figure 9: Relative contributions to the total variance of our final result for M RGI c

Figure 11 :
Figure 11: Illustration of the extraction of the heavy-light pseudoscalar decay constants, after applying a GEVP analysis, for ensemble J303.Top: plateau for the heavy-light pseudoscalar decay constant for the two fit intervals with higher weights in the model average.Bottom: summary of results from different fit ranges together with weights W , p-values and χ 2 /d.o.f..The shaded blue band represents the model average result.

Figure 12 :
Figure 12: Chiral behaviour of the best fits according to the TIC criteria applied to Φ D (top) and Φ Ds (bottom).Each point is projected to the physical charm quark mass, and results are shown for the flavour-averaged matching condition ϕ

Figure 13 :
Figure 13: Continuum limit extrapolation of the fits according to the TIC criteria applied to Φ D (top) and Φ Ds (bottom).Results are shown for the flavour-averaged matching condition ϕ (1) H .The blue band represents the projection to the physical ϕ 2 = ϕ phys 2 and ϕ H = ϕ phys H , while the red square denote the results in the continuum.

Figure 14 :
Figure14: Relative contributions to the total error of our determinations of f D (left) and f Ds (right).The label statistical plus χ-continuum limit represents the error arising from the statistical accuracy of our data and the chiral-continuum extrapolations.The scale setting label denotes the error coming from the physical value t phys

. 21 )f
Then, in order to arrive at a final determination of f Ds /f D we perform a model average among all the HMχPT and Taylor functional forms for the two different matching conditions ϕ Ds /f D 1.177(15)(6) 1.178(15)(6) 1.177(15)(5)

Figure 16 :
Figure 16: Illustration of the chiral-continuum extrapolation of the ratio Φ Ds /Φ D for the HMχPT model with highest TIC value.Results are shown for the flavour-averaged matching condition.Top: Chiral approach to the physical point.The dashed lines illustrate the chiral trajectories at finite lattice spacing, while the blue shaded band is a projection of the continuum fit.The red square symbol represents the physical result in the continuum.The black cross symbol corresponds to the symmetric point.Data points at finite lattice spacing are projected to the physical charm quark mass.Bottom: Lattice spacing dependence of Φ Ds /Φ D .The red square symbol indicates the continuum result, while the blue shaded band shows the fitted functional dependence on the lattice spacing.Points at finite lattice spacing are projected to the physical values of ϕ 2 and ϕ H .

14 WFigure 17 :
Figure 17: Summary of the model average procedure for the ratio Φ Ds /Φ D based on the combination of the two matching conditions, ϕ

aaχ
are considered .Similarly, 'taylor' refers to Eq.(5.18)where only the terms depending on the fit parameters p are kept.The labels 'p(2)' and 'p(4)' correspond to the addition of the higher order terms depending on the parameters p in Eq. (5.18), respectively, while 'pm(2)' denotes the addition of p (2) m from Eq. (5.21).Finally, 'p(3)' denotes the inclusion of the fit parameter p (3) a parameterising higher order lattice spacing dependence appearing in both the HMχPT and Taylor functional forms in Eq. (5.18) and Eq.(5.21).

f
Ds

Figure 18 :
Figure18: Left: Relative contributions to the total error on the determination of the ratio f Ds /f D .The label statistical plus χ-continuum limit represents the error arising from the statistical accuracy of our data and the chiral-continuum extrapolation.The scale setting label denotes the error coming from the physical value t phys 0 , while the model average represents the systematic error arising from the model variation according to the TIC procedure.Right: Details of the relative contributions to the statistical and chiral-continuum extrapolation error arising from specific gauge field configuration ensembles.

1 Figure 21 :
Figure 21: Illustration of the ground state and first excited state for a heavy-light pseudoscalar meson mass as extracted from the GEVP for the ensemble J303.We use t ref = t/2 such that the condition t ref ≥ t/2 is fulfilled.The shaded bands correspond to the plateaus of choice.Finally, to estimate the systematic uncertainty arising from the model variation we employ the weighted variance defined as follows

Table 1 :
List of CLS N f = 2 + 1 ensembles used in the present study.L/a and T /a refer to the spatial and temporal extent respectively of the lattice.The values κ l and κ s refer to the hopping parameters of the light and strange quark masses in the sea sector.Approximate values of the pion mass m π , the kaon mass m K , and of the product m π L are provided in the last three columns.

Table 2 :
List of run parameters for each ensemble in Table [9]parison of our charm quark mass determinations in the MS scheme with the FLAG average[9]and the results from other lattice QCD calculations based on N f = 2 + 1 dynamical simulations.In our results, shown in blue, we indicate both the total uncertainty and the error when excluding the uncertainty arising from Λ [89]CD 20A[89].Figure 20: Illustration of the ground-state effective masses determined from a GEVP analysis with three different ways of setting the value of t ref for the ensemble J303.The effective masses are normalised by the central value of the mass extracted from conservative plateau choices.The parameter t ref is either kept fixed, t ref /a = 1, 5, or varied by setting t ref = t/2 in such a way that the condition t ref ≥ t/2 is fulfilled.