The decay constants ${\mathbf{f_D}}$ and ${\mathbf{f_{D_{s}}}}$ in the continuum limit of ${\mathbf{N_f=2+1}}$ domain wall lattice QCD

We present results for the decay constants of the $D$ and $D_s$ mesons computed in lattice QCD with $N_f=2+1$ dynamical flavours. The simulations are based on RBC/UKQCD's domain wall ensembles with both physical and unphysical light-quark masses and lattice spacings in the range 0.11--0.07$\,$fm. We employ the domain wall discretisation for all valence quarks. The results in the continuum limit are $f_D=208.7(2.8)_\mathrm{stat}\left(^{+2.1}_{-1.8}\right)_\mathrm{sys}\,\mathrm{MeV}$ and $f_{D_{s}}=246.4(1.3)_\mathrm{stat}\left(^{+1.3}_{-1.9}\right)_\mathrm{sys}\,\mathrm{MeV}$ and $f_{D_s}/f_D=1.1667(77)_\mathrm{stat}\left(^{+57}_{-43}\right)_\mathrm{sys}$. Using these results in a Standard Model analysis we compute the predictions $|V_{cd}|=0.2185(50)_\mathrm{exp}\left(^{+35}_{-37}\right)_\mathrm{lat}$ and $|V_{cs}|=1.011(16)_\mathrm{exp}\left(^{+4}_{-9}\right)_\mathrm{lat}$ for the CKM matrix elements.


Introduction
The charmed D and D s mesons decay weakly into a lepton and a neutrino. The experimental measurement of the corresponding decay rates together with their prediction from within the Standard Model (SM) provides for a direct determination of the CKM matrix elements |V cd | and |V cs |. Leptonic D (s) decays have therefore been studied extensively by a number of experiments (CLEO-c [1][2][3][4][5][6][7], BES [8], Belle [9] and BaBar [10]). Together with the perturbative prediction of electroweak contributions in the SM this leads to the results [11] |V cd |f D + = 45.91(1.05)MeV , |V cs |f D + s = 250.9(4)MeV . (1.1) The reliable SM prediction for the decay constants where q = d, s , (1.2) hence allows for the determination of |V cd | and |V cs | respectively. Combined with calculations of other CKM matrix elements [11,12], a number of SM tests of, for instance, the unitarity of the CKM matrix can be devised (see e.g. refs. [13,14]). Surprisingly, relatively few state-of-the-art predictions for the decay constants with a reliable control of systematic uncertainties exist to date (see e.g. the discussion of the results of N f = 2+1 [15][16][17][18] and N f = 2+1+1 [19,20] in the 2016 FLAG report [21]). The computation of decay constants in lattice QCD is by now a well established exercise and many results with sub-percent precision within isospin symmetric QCD do exist for pions and kaons. This is less so however in the case of the D-and D s -meson decay constants. The major difficulty there lies in the fact that with the charm-quark mass slightly above 1 GeV, cut-off effects arising from the charm mass remain a serious concern in lattice simulations. Ensembles with dynamical quarks and sufficiently small lattice spacing to allow controlled continuum extrapolations have only become feasible in recent years.
Here we present the first calculation of the D-and D s -meson decay constants using the domain wall discretisation for the charm as well as the light and strange quarks on RBC/UKQCD's gauge ensembles with large volumes and physical values of the light-quark masses. Domain wall fermions (DWF) on the lattice provide chiral symmetry to a good approximation and as a result of this automatic O(a) improvement. In particular the latter is important when discretising charm quarks since no further work is required for tuning improvement coefficients in the action and for operators. Discretising both the light as well as the charm quark within the same framework will also allow to correctly reproduce GIM cancellation [28] in lattice computations of quantities such as K [29], ∆M K [30] and processes such as K → πl + l − [31,32] and K → πνν [33] which the RBC/UKQCD collaboration is pursuing.
Given the novel nature of domain wall fermions as heavy quark discretisation we have investigated their properties in detail in two preparatory publications [34,35]. We studied the continuum limit behaviour of heavy-strange decay constants over a wide range of lattice cut-offs (2-6 GeV) within quenched QCD. In this way we determined parameters of the domain wall discretisation which resulted in small discretisation effects for charmed meson decay constants. Refs. [34,35] show that DWF with suitably chosen parameters show mild cut-off dependence for charmed meson masses and decay constants. We expect this to hold also in the presence of sea-quarks on RBC/UKQCD's gauge ensembles, thus the study at hand.
This work summarises our computation within this setup of the the D-and D s -meson decay constants f D and f Ds , respectively, their ratio and the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements |V cd | and |V cs |. For convenience we anticipate the numerical results: where errors are statistical, systematic and experimental, respectively. Note that the quoted results for the decay constants are for isospin symmetric QCD. 1 Systematic errors are due to choices made when fitting and parameterising the lattice data, lattice scale setting, finite volume, isospin breaking and renormalisation, based on N f = 2 + 1 flavour simulations. Isospin breaking effects in the determination of CKM matrix elements are based on estimates. In the determination of the CKM matrix element the lattice statistical and systematic errors have been combined in quadrature. The paper is structured as follows: in section 2 we summarise our numerical set up and give details of all simulation parameters. Section 3 presents our complete analysis for the D (s) -meson decay constants and constitutes the main body of this paper. In section 4 we extract the corresponding CKM matrix elements before concluding in 5.

Numerical simulations
This report centres mainly around ensembles with physical light-quark masses in large volumes [37]. The data analysis will however also take advantage of information obtained on ensembles with unphysically heavy pions [38][39][40]. The gauge field ensembles we use (cf. table 1) represent isospin symmetric QCD with N f = 2 + 1 dynamical flavours at three different lattice spacings in the range 0.11 fm-0.07 fm (Coarse, Medium and Fine). All ensembles have been generated with the Iwasaki gauge action [41,42]. For the discretisation of the quark fields we adopt either the DWF action with the Möbius kernel [43][44][45] or the Shamir kernel [46,47]. The difference between both kernels in our implementation corresponds to a rescaling such that Möbius domain wall fermions (MDWF) are loosely equivalent to Shamir domain wall fermions (SDWF) at twice the extension in the fifth dimension [37]. Möbius domain wall fermions are hence cheaper to simulate while providing the same level of lattice chiral symmetry. Results from both formulations of domain wall fermions lie on the same scaling trajectory towards the continuum limit with cut-off effects starting at O(a 2 ). Even these O(a 2 ) cut-off effects themselves are expected to agree between our Möbius and Shamir formulations, with their relative difference at or below the level of 1% [37] for the finite values of L s used in our simulations.  Basic properties of all ensembles used in this work are summarised in tables 1 and 2. The lattice scale and physical light-quark masses have been determined for all ensembles bar F1 in [37] using m π , the m K and the m Ω as experimental input. The finest ensemble F1 was generated later as part of RBC-UKQCD's charm and bottom physics program and its properties are described in appendix A. We repeated the analysis of [37] after including also ensemble F1 and in this way determined the value of the lattice spacing and physical strange-quark mass on this ensemble.
In the valence sector we simulate light, strange and charm-like heavy quarks. The light-quark masses were chosen to be unitary. The strange-quark mass was slightly adjusted (partially quenched) to its physical value as determined in [37], in cases where this was known prior to running the measurements (C1, C2, M1 and M2). Otherwise the unitary value was chosen. The parameters of the light and strange sector are listed in table 2.
Besides the bare quark mass, DWF have two further input parameters that need to be specified in each simulation: the extent of the fifth dimension L s and the domain wall height parameter M 5 , respectively (for details see [35,37]). More specifically, M 5 is the negative mass parameter in the 4-dimensional Wilson Dirac operator which resides in the 5-dimensional DWF Dirac operator. The parameter M 5 can have some effect on both, the rate of exponential decay of the physical modes away from the boundary in the 5th dimension and the energy scale of unphysical modes that are not localised at the boundary. While changing L s mostly changes the magnitude of residual chiral symmetry breaking, the choice of M 5 in principle changes the ultraviolet properties of the discretisation. Hence, different choices of M 5 lead to in principle different scaling trajectories. The most important consequence is that calculations combining (quenched) charm quarks with a distinct M 5 from the light sea quarks must formally be treated as a mixed action when renormalising flavour off diagonal operators. Two observations in this context which we made in our quenched DWF studies [34,35] are crucial for understanding the choice of simulation parameters made here: studying the pseudoscalar heavy-heavy and strange-heavy decay constants we found cut-off effects to be minimal for M 5 ≈ 1.6. We find that the residual chiral symmetry breaking effects are well suppressed for L s = 12. At the same time we observed a rapid increase in discretisation effects as the bare input quark mass was increased above am h ≈ 0.4. Here we work under the assumption that these observations carry over to the case of dynamical simulations with N f = 2 + 1 flavours which motivates our choice of M 5 = 1.6 for charm DWF keeping am h ≤ 0.4. As we will see later this assumption is well justified. For the light quarks we use M 5 = 1.8 both in the valence-and sea-sector.
All simulated bare charm-quark masses are listed in table 3. Note, that we allow for one exception to the bound am h ≤ 0.4 by generating data for am h = 0.45 on ensemble C0. With this we tested whether the reach in the heavy quark mass for DWF with M 5 = 1.6 observed in the quenched theory [35] also persists in the dynamical case. The quantity that we monitor in this context is the residual quark mass am res , which provides an estimate of residual chiral symmetry breaking in the DWF formalism. It is defined in terms of the axial Ward identity (AWI) where ∆ − µ is the lattice backward derivative, am the bare quark mass in lattice units in the Lagrangian, P is the pseudoscalar density and J 5q is the pseudoscalar density in the centre of the 5th dimension. It motivates the definition . Figure 1 shows the behaviour of the residual mass on the C0 ensemble. As expected (see reference [35] for details) the residual mass does plateau and remains flat for am h 0.4, confirming the validity and upper bound of the mass point at am h = 0.4. The fact that this behaviour is not observed with am h = 0.45 indicates that cut-off effects change in nature when am h is pushed beyond 0.4, in agreement with our findings in ref. [35]. We therefore exclude any data with am h > 0.4 in the remainder of this paper.   Table 3. Möbius domain wall parameters for the heavy quarks on all ensembles. All quoted values for am h are bare quark masses in lattice units. As described in the text, the value indicated in red was only used to verify our assumptions about the applicability of the quenched pilot study to the dynamical case.

Data Analysis
In this section we describe how we make predictions for the decay constants starting from the evaluation of Euclidean two-point correlation functions on all ensembles and for all parameters discussed above. In particular, from fits to the simulated data we determine pseudoscalar masses m P and decay constants f P and ratios thereof, which depend on the simulation parameters (a, m l , m h , m s ). Of particular relevance for the analysis are P = π, K, D, D s and η c (the latter one corresponding to an unphysical cc state made of two distinct flavours of quarks having the same charm quark mass). We extrapolate the data for each observable O to physical light, strange and charm-quark masses as well as to vanishing lattice spacing and infinite volume, O(a = 0, m phys l , m phys h , m phys s ). Besides the decay constants and ratios thereof, for reasons that will become clear later, we will carry out the analysis of the decay constants in terms of the quantity f P √ m P and only remove the factor √ m P in the final step.

Correlation functions
In practice we determine the matrix element in equation (1.2) from the time dependence of Euclidean QCD zero-momentum two-point correlation functions, (3.1) We consider the cases i = π, K, D, D s or η c . The operators O i,M are interpolating operators with the quantum numbers of the desired mesons, e.g. O Ds,M =c Γ M s , where we consider Γ A = γ 0 γ 5 and Γ P = γ 5 , respectively. The sum on the r.h.s. of eq. (3.1) is over excited states k and in practice we will consider the ground and first excited state in the data analysis. The constants Z is the corresponding kth excited meson state.
When computing the quark propagators we use Z(2)×Z(2) stochastic wall sources [48][49][50] on a large number of time planes. Details of how many different source planes are used for the various ensembles are listed in the column "hits/conf" in table 1. The results on a given gauge configuration are averaged into one bin.
The calculation of the light and strange quark propagators were performed using the HDCG algorithm [51], reducing the numerical cost and hence making this computation feasible. For the heavy quark propagators a CG inverter was used and we monitored satisfactory convergence using the time-slice residual introduced in ref. [52].
Masses and decay constants have been determined by simultaneous multi-channel fits to the two-point correlation functions C i,AP and C i,P P for a given choice of i. We attempted the use of correlated fits, but found the correlation matrix to be too poorly estimated for a reliable inversion. We therefore carry out uncorrelated fits, i.e. assume the correlation matrix to be diagonal (compare ref. [37]).
The statistical precision of the ground state mass and matrix elements is improved by fitting the ground state as well as the first excited state, allowing for earlier time slices (with smaller statistical errors) to contribute to the fit. This was done for i = D, D s , where we are interested in the matrix elements and not merely in the masses as in the case for i = π, η c . During all these fits the χ 2 /d.o.f. were monitored. Tables with the bare results in lattice units on all ensembles can be found in appendix B. Figure 2 illustrates the correlator fit for the heaviest charm-quark mass on the C0 ensemble. The fit ranges were chosen by systematically varying t min and t max , fixing them in a region where no dependence is observed. Figure 3 shows the ground state results of varying t min and t max for the case of the heaviest mass point on the coarse ensemble C0.
A first impression of the range of ensembles (lattice spacing, light sea quark mass, charm-quark mass) for which we generated data is given in figure 4, representatively for the ratio of decay constant f Ds /f D plotted against the inverse of the measured η c mass.

Non perturbative renormalisation
To make contact between lattice regulated data and quantities in a continuum theory, the heavy-light current needs to be renormalised. Since we use a mixed action current, by  Since in the free field theory the modification of the action represents a modest change to irrelevant parameters, we might hope the impact on renormalisation constants is small. This is something we can verify. In fact, the estimate of the systematic error arising from this change in action is found to be below percent level, and this is discussed in more detail in section 3.2.2 along with a proposal for the determination of fully non-perturbatively   renormalised axial currents using appropriate ratios of off-shell mixed and unmixed action vertex functions. Since empirically this is indeed a small effect, in all of the following we will extract the renormalisation constants from a light-light unmixed action current and associate a systematic error devised from the non-perturbative renormalisation (NPR) study.
It is worth mentioning that we have recently developed a massive renormalisation scheme [53], RI/mSMOM, by extending the massless RI/SMOM scheme [54,55], with the renormalised composite fields defined away from the chiral limit. This includes finite masses in one, i.e. the mixed case, or both quark fields entering the bilinear operator. Using this scheme the renormalisation constant for the heavy-light axial current can be extracted non-perturbatively. For more details, refer to ref. [53] and appendix D.

Unmixed action axial current renormalisation constants
The light-light axial renormalisation constant can be found from fitting the time behaviour of the relation to a constant [37,56]. Here C AP (t) is the correlation function between the conserved point-split axial vector current defined on the links between the lattice sites [56] and the pseudoscalar density, whilst C AP (t) is the same correlation functions with the conserved axial vector current replaced by the local current one defined on the lattice sites. The folded time behaviour of the light-light current for all ensembles scaled to the interval [0, 1) is shown in figure 5. As expected we can see a lattice spacing dependence. The slight difference between C1 and C2 (M1 and M2) arises from the slightly different light-quark mass m l . We can also identify a plateau region to which we can fit a constant to obtain the values of the renormalisation constants. The results of these fits are listed in table 4 and are in good agreement with ref. [37].

Vertex functions of mixed action current
We use the non-exceptional Rome-Southampton renormalisation scheme, RI/SMOM, [54,55] to investigate the effects of a change in the action on quantities entering the axial current renormalisation constant. In particular we evaluate the variation of the projected, amputated vertex function for the axial current P[Λ A ], on each of the ensembles C2, M1 and F1. The details of the numerical computation of the amputated axial vertex function Λ A in the RI/SMOM scheme are discussed in appendix D. Tables 5, 6 Table 7. The ratios of projected amputated vertex function for the axial currents with different actions on the F1 ensemble. The quark mass for both fields is taken to be am l = 0.002144. ciently close to the chiral limit -is used. The data has been generated using ten gauge field configurations which leads to sufficiently precise results. We see that the ratio P[Λ A ](1.8,1.8) on each of the ensembles has the largest deviation from unity as compared to the other ratio combinations, which is expected since both the quark fields entering the bilinear have different actions between the numerator and the denominator. The main feature emerging from this study is that the deviation from unity is at most of order 0.4% across a range of momenta around 2 GeV. This is negligible on the scale of our other uncertainties and for the purposes of the present work we can simply include it as a sub-dominant systematic error.

Strange-quark mass correction
We determine the physical strange-quark mass on all ensembles considered here by repeating the global fit to light meson observables detailed in [37] but including our new fine ensemble F1. At the time of the data generation the values of am phys s were not yet known for ensembles with near-physical pion masses (C0 and M0) and the new finer ensemble F1. To correct for the resulting small mistuning we repeated the simulation of all charmed meson observables on C1 and M1 with both the unitary and the physical valence strangequark mass. From this we obtain information on the (small) corrections on C0, M0 and F1. We define the parameters α O (a, m i h ) for the different observables O using To obtain the values for F1, we need to extrapolate α O (a, m h ) measured on C1 and M1 to (a, m h ) appropriate for F1. Given the linear behaviour in the inverse heavy quark mass evident from the left panel of figure 6, we linearly extrapolate the values for each given lattice spacing to the corresponding m i h (F 1). This is then extrapolated to the lattice spacing of F1 by fitting the data to The results for this are shown by the green diamonds in the right-hand panel of figure 6 and summarised in

Fixing the physical charm quark
We have a number of possible choices for the meson H that fixes the charm-quark mass.
The ones we will consider are H = D, D s and η c . Each of these has slightly different advantages and disadvantages attached. The lattice data for the D meson is comparably noisy and has a strong light-quark-mass dependence, making it difficult to disentangle the extrapolation to physical light-quark masses from the interpolation to the physical charmquark mass. The D s is statistically cleaner and depends less on the light sea-quark mass than the D, but we need to correct for a mistuning in the valence strange-quark mass as discussed in section 3.3. Finally the η c is statistically the cleanest, but it differs by quarkdisconnected Wick contractions with respect to the corresponding physical particle listed by the Particle Data Group (PDG) [11]. However, this is assumed to be a small effect. 2 We will investigate all three choices and use the spread as an indication of potential systematic errors. The masses of these mesons, stated by the PDG [11] are (14) GeV, m ηc = 2.9836(6) GeV.

Global fit
Our fit ansatz corresponds to a Taylor expansion around the physical value of the relevant meson masses. It is given by and H = D, D s or η c . This means we simultaneously fit the continuum limit dependence (coefficients C CL ), the pion mass dependence (coefficients C χ ) and heavy quark dependence (coefficients C h ) as well as cross terms (coefficients linear in ∆1/m h , i.e. C 1 χ and C 1 CL ) in one global fit. The coefficients C 1 CL and C 1 χ capture mass dependent continuum limit and pion mass extrapolation terms. This arises by expanding C CL (m h ) and C χ (m h ) in powers of ∆m −1 H . When considering the individual decay constants (as opposed to their ratio) we use the quantity Φ D (s) = f D (s) √ m D (s) in the subsequent analysis. As can be seen in figure 7 this quantity is, within statistical resolution, linear in 1/m ηc . Also its dependence on the light (sea and valance) quark-mass, as shown exemplarily in figure 8, is linear irrespective of the heavy-quark mass. We estimate the systematic uncertainty by limiting the data entering the fit to the case of pion masses not larger than 450, 400 or 350 MeV in turn. Another variation we have already mentioned is the choice of the meson that fixed the charm-quark mass. Finally we can modify the fit form (3.6) by setting some of the parameters to zero by hand (e.g. C 1 CL and C 1 χ ), which we will do when the data is not sufficiently accurate to resolve them clearly.
One example for the global fit according to (3.6) for the case of the observable f Ds /f D . In the case presented here the charm-quark mass is fixed by the η c meson and a pion mass cut of m π < 400MeV is employed. The grey band shows the fit result at physical pion masses and vanishing lattice spacing. The coloured bands correspond to the fit projected to the given pion mass and lattice spacing for the corresponding ensembles. In this fit we ignore heavy mass dependent continuum and pion mass terms.

Global fit for ratio of decay constants
For the ratio of decay constants, fully correlated fits could be achieved. Figure 9 gives one example of such a fully correlated fit for the ratio of decay constants. The fit shown here has a pion mass cut of m π < 400MeV and uses the η c mass to fix the charm-quark mass. Furthermore, heavy mass dependent coefficients of the continuum limit and the extrapolation to physical pion masses are ignored (i.e. C 1 CL = 0 and C 1 χ = 0). Table 14 in appendix C summarises the results of all fit variations for f Ds /f D . The results of these are also shown in figure 10. The red (blue, green) data points correspond to  Finally the label at the x-axis describes which fit was used by stating the number of coefficients for the continuum limit (CL) and pion mass limit (χ) respectively. E.g. fits results labelled (2, 2) correspond to the fit form (3.6) whilst (2, 1) corresponds to keeping two coefficients for the continuum limit extrapolation, but only one coefficient for the pion mass extrapolation by setting C 1 χ to zero. Cases where one of the coefficients C 1 CL and C 1 χ is compatible with zero at the one sigma level are indicated by the corresponding data point being partially transparent.
From the results shown in figure 10 we can make a few observations. We find that the ratio of decay constants is insensitive to the way we fix the charm-quark mass. This is not surprising as the ratio of decay constants does not strongly depend on the heavy quark mass (compare figure 9). We find that a dependence is observed when including pions with m π > 400 MeV, for this reason we restrict ourselves to m π ≤ 400 MeV. We can also see that when allowing for heavy mass dependent pion mass and continuum extrapolation terms, these can not be resolved with the present data. They also do not significantly change the central value of the fit result but increase the statistical error. This is again not surprising, given the mild behaviour with heavy quark mass displayed by the data.
From this discussion we choose the highlighted fit (i.e. the one presented in figure 9) as our final fit result and as statistical error. We assign the systematic error arising from the fit form from the maximal spread in the central value of the fit results as we vary the parameters of the fit, maintaining m π ≤ 400 MeV. More precisely, we take the maximal difference between the central value of the preferred fit and the central values of all fits with m π ≤ 400 MeV displayed in figure 10. From this we quote where the first error is statistic and the second error captures the systematic error associated with the chiral-continuum limit as well as the way the charm-quark mass is fixed. Figure 11 shows the chosen fit results for Φ D (top) and Φ Ds (bottom) respectively. In both cases the heavy quark mass is fixed by the η c mass and a pion mass cut of m π ≤ 400 MeV is used. Contrary to the fit of the ratio of decay constants, correlated fits of Φ D and Φ Ds proved to be unstable. Uncorrelated fits were used instead which lead to slightly larger errors. Tables 15 and 16 in appendix D summarise the results of all fit variations for Φ D and Φ Ds respectively. Similar to the previous section we vary the fit parameters to determine the stability of the results. We find that we can consistently resolve the C 1 CL coefficient in the case of Φ Ds , whilst this is less clear in the case of Φ D . For this reason we choose (CL, χ) = (1, 1) for the case of Φ D and (CL, χ) = (2, 1) for the case of Φ Ds (see figure 12). Again, little dependence is observed in the case of m π ≤ 400 MeV so this pion mass cut is used. The dependence is larger in the case of Φ D than for Φ Ds in agreement with intuition. We see little dependence in the way the heavy quark mass is fixed, even though (contrary to the ratio of decay constants) the heavy mass dependence is now significant. Overall we see more variation in the results of the fit than we have for the ratio of decay constants. Following the same procedure to determine the systematic error associated with the fit as above we find Φ D = 0.2853(38)( +24 −18 ) fit GeV 3/2 , Φ Ds = 0.3457 (26)

Global fit for Φ D and Φ Ds
(3.8)

Systematic error analysis
So far we have discussed central values, statistical errors and the systematic errors due to the fit for the value of Φ D , Φ Ds and f Ds /f D and the non-perturbative renormalisation. We now discuss the remaining systematic error budget due to scale setting, mistuning of the strange-quark mass, finite volume and isospin breaking and continuum limit. The uncertainty in determining the lattice spacing has been propagated throughout the entire analysis by creating a bootstrap distribution with the width of the error quoted in table 1. The uncertainty in the physical strange-quark masses arising from ref. [37] has been treated in the same way. Both of these are therefore already included in the statistical error.
We have already discussed the systematic error arising from the correction of the mistuning of the strange-quark mass in section 3.3 and came to the conclusion that this yields an uncertainty of 9 × 10 −6 GeV 3/2 for Φ Ds and 3 × 10 −4 for f Ds /f D .
We estimate the finite size effects by comparing our values of m π L to a study the MILC collaboration has undertaken [20]. They studied the volume dependence of charmed meson For the decay constants the variations they found are < 0.3% and < 0.15%. Applied to our data, this leads to estimates of the systematic errors of δΦ D ≈ 0.001 GeV 3/2 and δΦ Ds ≈ 0.0006 GeV 3/2 . We expect f D and f Ds to be similarly affected by finite size effects, and therefore expect cancellations in the ratio. We conservatively take the larger relative error (0.3%) as an estimate for the ratio, yielding δ f Ds f D ≈ 0.0035. Given that the minimum value of m π L for our ensembles is 3.8, results derived from these numbers are a good conservative estimate.
In our simulations we treat the up and down quark masses as degenerate, which is   not the case in nature, and neglect electromagnetic effects. This affects in particular the masses of the mesons we consider. In principle these effects cannot be disentangled. In the determination of the decay constants we neglect electromagnetic effects since they are defined as pure QCD quantities. However, for the determination of the CKM matrix elements these effects will need to be taken into account [11,20].
We devise a systematic error associated to the way we fix the heavy quark mass by considering how much the fit result for Φ D changes when we replace the input mass m D ± = 1.86961(09) GeV by m D 0 = 1.86484(05) GeV [11]. We estimate the effect of this shift using the fit result of the coefficient C 0 h for the case of h = D and multiplying these by  Table 9. Summary of the systematic error budget for the quantities Φ D , Φ Ds and the ratio of decay constants. Details of the discussion leading to these results can be found in the text.
m π ± instead of m π 0 as input mass, i.e. calculating C 0 χ m 2 π ± − m 2 π 0 . From this we find δΦ D ∼ 0.00029 GeV 3/2 , δΦ Ds ∼ 0.00001 GeV 3/2 and δ f Ds f D ∼ −0.00080. Adding these two effects in quadrature we obtain the values listed in the column m u = m d in table 9.
Given that the continuum limit coefficient C 0 CL is compatible with zero for the fits chosen for Φ D (C 0 CL = −0.003(11) GeV 7/2 ) and f Ds /f D = 0.005(25) GeV 2 , we neglect higher order O(a 4 ) effects. For Φ Ds we find C CL 0 = −0.027(10) GeV 7/2 (the heavy mass dependent continuum limit term vanishes at the physical charm-quark mass). To estimate the impact of higher order discretisation effect terms we write Substituting the numbers for C 0 CL and the coarsest and finest lattice spacings we find C 0 CL a 2 /Φ Ds ∼ 0.026 and 0.010 respectively. Assuming D 0 CL /C 0 CL = (0.5 GeV) 2 (i.e. setting the scale such that discretisation effects grow as a/Λ with Λ = 500 MeV) we find D 0 CL /C 0 CL a 2 ∼ 0.008 and ∼ 0.003. So the residual discretisation effects are 8% (3%) of the leading discretisation effects, yielding at most 0.2% of the absolute value.
Combining these errors in quadrature we arrive at our final values for Φ D (s) . Using the masses of D ± and D ± s [11] (compare eq. (3.5)) we obtain values for the decay constants f D (s) , (3.10) We are now in a position to compare our results to the results found in the literature. Adding our results to those presented in the most recent FLAG report [21] we obtain the plots in figure 13. The smaller error bar presents the statistic error only, whilst the larger error bar shows the full error (statistic and systematic). In all cases the error budget is dominated by the statistical error. We find good agreement with the literature and have errors competitive with the other results displayed in figure 13 [15-20, 27, 57-66].  Figure 13. Superposition of our results (blue circles) to the data presented in the most recent FLAG report [21]. The small error bar shows the statistic error only, whilst the large error band includes both, the statistic and the systematic error.

CKM matrix elements
Having obtained the decay constants, we can make a prediction of the CKM matrix elements |V cd | and |V cs |. However, the values shown in (1.1) are obtained in nature and therefore we need to adjust these values to those of an isospin symmetric theory. In other words, the measured decay rate |V cq | f Dq does include electroweak, electromagnetic and isospin breaking effects, so before extracting |V cq | we need to correct the decay rate for these effects. Ref. [20] distinguishes between universal long-distance electromagnetic (EM) effects, universal short distance electroweak (EW) effects and structure dependent EM effects. All of these modify the decay rate to match the experimental value to the theory in which we simulate. The combined effect of the universal long-distance EM and shortdistance EW effects is to lower the decay rate by 0.7% [20,67,68]. We adjust the decay rates from (1.1) and then calculate the CKM matrix elements from this. We find (4.1) Again, we can superimpose our results to those obtained in the most recent FLAG report [21], shown in figure 14. This combines the results of refs. [15-20, 64, 69, 70]. Again we find good agreement between previous works and obtain a competitive error.  Figure 14. Superposition of our results (blue circles) to the data presented in the most recent FLAG report [21]. The smaller error bars of our results show the lattice error only, whilst the large error bands include both, the theoretical and the experimental errors, added in quadrature.

Conclusion and outlook
(f D ), 1.0% (f Ds ) and 0.7% (f Ds /f D ) the results are competitive and establish domain wall fermions as a powerful discretisation for heavy quarks. We hope that our results will provide useful input to a wide range of applications in (Beyond) Standard Model phenomenology. Looking ahead, we are exploring changes in the formulation of the domain wall action, such as gauge link smearing, which we found increases the reach in the heavy quark mass on a given ensemble before cut off effects become substantial [71]. This will allow us to do computations directly at the physical charm quark mass also on our coarsest ensemble.

A Properties of ensemble F1
Here we present details and properties of ensemble F1 which was generated in order to allow for a continuum limit with three lattice spacings. It has not appeared in any of RBC/UKQCD's previous analyses.
In table 10 we summarise basic simulation parameters and properties for ensemble F1. All integrated autocorrelation times were estimated using the technique described in ref. [73]. We have used the same implementation of the exact hybrid Monte Carlo algorithm for the ensemble generation as in ref. [37], with five intermediate Hasenbusch masses, (0.005,0.017, 0.07 , 0.18, 0.45), for the two-flavor part of the algorithm. A rational approximation was used for the strange quark determinant. See ref. [37] for more details. Figures 15 and 16 show the Monte Carlo evolution and histograms of the topological charge (measured with the GLU package [74]) and the Wilson flow scales w 0 and √ t 0 , respectively. The measured autocorrelation time motivates our choice to separate evaluations of observables on ensemble F1 in the main part of this paper by 40 molecular dynamics steps (for comparison: the separation is 40 on C1 and 20 on all remaining ensembles).

B Correlator fit results
Tables 11, 12 and 13 summarise the fit ranges and the fit results of the correlation function fits to the heavy-light, heavy-strange and heavy-heavy pseudoscalar mesons, respectively. Since the different channels may have a slightly different approach to the ground state (compare figure 2) we quote the smaller value of t min . In all cases t AP min ≤ t P P min . The t max quoted is the first time slice which is not included in the fit.           where Λ A,R is the amputated axial vertex function and the subscript R denotes a renormalised quantity. The momentum q out of the vertex satisfies the symmetric non-exceptional condition: The momenta are determined by for every lattice spacing L such that the magnitude of p is around 2 GeV for an integer n. Note that in order to reach the intermediate momenta we use twisting [75]: Eq. (D.1) can be written in terms of the amputated bare vertex function, the field renormalisation Z q and the axial operator renormalisation Z A as follows, Tr (q · Λ A ) γ 5/ q sym = 1, (D. 5) where the bare quantity is what is computed numerically on the lattice. Above we have denoted P[Λ A ] ≡ lim m→0 1 12q 2 Tr (q · Λ A ) γ 5/ q sym .
(D. 6) In principle, if the action is substantially changed for the heavy as compared to light quarks, it is possible in the massless renormalisation scheme, to apply the axial RI-SMOM condition (D.1) to the mixed action bilinear Again, the indices i = 1, 2 refer to the action of the first and second quark field entering the bilinear operator. In the chiral limit, it is possible to systematically eliminate the factors of the quark field renormalisation Z q from the SMOM condition [53]. However, since we can compute the corresponding unmixed vertex functions, and know how to determine Z 11 A and Z 22 A from the conserved domain wall current as in eq. (3.2), we are able to determine the axial current renormalisation, Z 12 A , from ratios of vertex functions as follows: (P[Λ A ](1.6, 1.6)) (P[Λ A ](1.8, 1.8)) (P[Λ A ](1.6, 1.8)) 2 = Z 12 This result is an interesting, and fully non-perturbative, analogue of the Fermilab [76] partially-perturbative approach to currents that are conserved in the continuum. The ratio of vertex functions in our case is very near unity since the difference in the actions of the quark legs is very small, as is seen in tables 5, 6 and 7. Furthermore, if the RI-SMOM conditions are modified we can form a consistent set of conditions in the massive case [53]. Since the axial current is partially conserved in the continuum, the difference between the schemes is necessarily only a lattice artefact, implying that either approach may be taken when we take the continuum limit, and yields a universal result. The scaling violations will of course differ between these approaches and for the present paper we have taken the simpler approach of using the near massless vertex function data to define the scaling trajectory for axial current matrix elements.