Charm physics with overlap fermions on 2+1-flavor domain wall fermion configurations

Charm physics provides a rich phenomenology and offers a valuable platform for precisely testing the Standard Model and understanding low-energy quantum chromodynamics (QCD). The decay constants of charmed pseudoscalar mesons combined with the experimental data of the relevant leptonic decays can be used to extract the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements $ V_{cd} $ and $ V_{cs} $ (see [1]), or serve as probes for new particles of new physics, such as charged Higgs bosons (see [2]). However, perturbative theory and heavy quark expansion are rather precarious in the charm sector, and our most reliable theoretical tool is lattice QCD. An overview of the decay constants of charmed pseudoscalar mesons can be found in Ref. [3] for lattice QCD calculations before 2022. Two recent $ 2+1 $-flavor calculations of these decay constants appeared in Refs. [4, 5].

For the pseudoscalar meson D or $ D_s $, its pure leptonic decay $ D_{(s)}\to \ell\nu_\ell $, of which the decay width is proportional to $ f^2_P|V_{cq}|^2 $ ($ q = d\text{ or }s $ with $ f_P $ the decay constant of $ D\text{ or }D_s $, respectively), provides a clean channel to determine the CKM matrix elements. Moreover, the semileptonic processes $ D_{(s)}\to h\ell\nu_\ell $, with h representing a pion or kaon, are induced by the vector current $ \bar{q}\gamma^\mu c $, while the pure leptonic processes are induced by the axial-vector current $ \bar{q}\gamma^\mu\gamma_5 c $. Comparing the CKM matrix elements extracted from these two kinds of processes allows us to test the $ V-A $ structure of weak interactions. The vector mesons $ D_{(s)}^* $ show approximately 100% decay in relation to the corresponding charmed pseudoscalar mesons, with relatively small leptonic decay branching ratios. The first measurement of $ f_{D_s^*} $ from the leptonic decay $ D_s^*\to e^+ \nu_e $ was recently reported by BESIII [6]. In the phenomenological analysis of the nonleptonic weak decays of charmed vector mesons, the decay constants denoted by $ f_V $ are necessary inputs [7, 8]. In each of the $ b\to c $ induced semileptonic or nonleptonic decays of bottom mesons, the decay constant $ f_V $ of the charmed meson in the final state, as well as the decay constant $ f_V^T $ defined by the tensor current, appears as nonperturbative inputs [9, 10]. Furthermore, the ratio of the decay constants of vector and pseudoscalar mesons $ f_V/f_P $ approaches 1 in the heavy quark limit, which can be used to study the breaking of heavy quark symmetry. Currently, unquenched lattice QCD calculations of decay constants of charmed vector mesons are relatively sparse [11−15]. The results from sum rules can be found in, for example, Refs. [16, 17].

For the charmonium system, the leptonic decay constant of $ J/\psi $ can be directly determined experimentally [18]. In this study, we calculated the ratios of decay constants $ f_{J/\psi}/f_{\eta_c} $ and $ f^T_{J/\psi}/f_{J/\psi} $, where the former can be directly used in the amplitude analysis of B meson decays into charmonium states, and the latter appears in $ b\to c $ decays induced by specific new physics operators. By combining the decay width $ \Gamma(J/\psi\to e^+e^-) $ measured experimentally with the ratio $ f_{J/\psi}/f_{\eta_c} $, we can determine the decay constant $ f_{\eta_c} $, which is used for calculating decay amplitudes $ \Gamma(\eta_c\to\gamma^*\gamma) $ and form factors appearing in $ \gamma^*\gamma^*\to\eta_c $ processes [19−21]. Theoretically, the decay constants $ f_{J/\psi} $ and $ f_{\eta_c} $ were calculated in quenched lattice QCD in [22]. Two-flavor lattice QCD calculations of $ f_{J/\psi} $ and $ f_{\eta_c} $ are given in [23, 24]. The HPQCD collaboration obtains $ f_{J/\psi} $ and/or $ f_{\eta_c} $ in $ 2+1 $-flavor [25, 26] and $ 2+1+1 $-flavor simulations [27]. The latter work also considers the quenched quantum electrodynamics (QED) effects of charm quarks.

Our previous $ 2+1 $-flavor calculation [28] of the decay constants of charmed mesons were performed at an inverse lattice spacing $ a^{-1} = 1.730(4)\, $GeV. The discretization error was estimated to be approximately 2%. We now present the results on a finer lattice with $a^{-1} = 2.383(9)\,$GeV to better understand the lattice cutoff effects.

In the following, we present the details of our calculation framework in Sec. II, which includes the lattice setup and the computation of two-point correlation functions. Sec. III.A and Sec. III.B show the analyses of the meson masses and decay constants, respectively. The discussions and summary can be found in Sec. IV.

The $ 2+1 $-flavor gauge configurations used in this study were generated by the RBC-UKQCD collaborations [29]. The dynamical quarks are domain-wall fermions with degenerate light (up and down) quark masses $ am_l^{\rm sea}=0.004 $, $ 0.006 $, $ 0.008 $ and strange quark mass $ am_s^{\rm sea}= $ 0.03 in lattice units. The inverse lattice spacing $ a^{-1} $ determined in Ref. [30] and the other parameters of the configurations are listed in Table 1. The spatial extension of the lattice is approximately $ La\sim2.7 $ fm. The light sea quark masses $ am_l^{\rm sea} $ mentioned above correspond to pion masses $ m_\pi^{\rm sea} $ at approximately $ 302 $, 360, and 412 MeV, respectively [29]. More information about the configurations can be found in Ref. [29].

$ a^{-1} $/GeV	Label	$ am_l^{\rm sea}/am_s^{\rm sea} $	Volume	$ N_{\rm conf}\times N_{\rm src} $	$ am_{\rm res} $
2.383(9)	$ \mathrm{f004} $	0.004/0.03	$ 32^3\times64 $	$ 628\times1 $	0.0006664(76)
	$ \mathrm{f006} $	0.006/0.03	$ 32^3\times64 $	$ 42\times16 $
	$ \mathrm{f008} $	0.008/0.03	$ 32^3\times64 $	$ 49\times16 $

Table 1.  Configurations used in this study. The residual mass of the dynamical fermion $ am_{\rm res} $ is in the two-flavor chiral limit from Ref. [1]. $ N_{\rm conf} $ is the number of configurations, and $ N_{\rm src} $ the number of point sources on each configuration.

The valence quarks used in this study are overlap fermions. The massless overlap Dirac operator [31] is defined as

where ε is the matrix sign function and $ D_{\rm w}(\rho) $ is the usual Wilson fermion operator, except for a negative mass parameter $ - \rho = 1/2\kappa -4 $ with $ \kappa_c < \kappa < 0.25 $ and $ \kappa_c $ corresponding to a massless Wilson operator. In practice, we use $ \kappa = 0.2 $, which corresponds to $ \rho = 1.5 $. The massive overlap Dirac operator is defined as

To accommodate the $ S U$(3) chiral symmetry, it is usually convenient to use the chirally regulated field $ \hat{\psi} = \left(1 - \dfrac{1}{2} D_{\rm ov}\right) \psi $ in lieu of ψ in the interpolation field and the currents. This is equivalent to leaving the currents unmodified and adopting the effective propagator instead

where $ D_c = \dfrac{\rho D_{\rm ov}}{1 - D_{\rm ov}/2} $ satisfies $ \{\gamma_5, D_c\}=0 $ [32].

The decay constant $ f_P $ of a pseudoscalar meson P is defined through the following matrix element of the axial-vector current $ A_\mu(x)=\bar{q}_1(x) \gamma_\mu \gamma_5 q_2(x) $,

where $ p_\mu $ is the four-momentum of the meson P, with $ q_{1,2} $ referring to the spinor fields of the constituent quarks present. For charmed mesons $ D_{(s)} $, the two quarks are charm and light (strange) quarks, and for charmonium $ \eta_c $, both of the constituent quarks are charm quarks. In lattice QCD, the matrix element in Eq. (4) can be extracted by calculating the two-point function involving $ A_\mu(x) $,

where $ t\equiv x_0-s_0 $ is the time displacement in lattice units between the source point s and sink point x. The interpolating operators are $ {\cal{O}}=\bar{q}_1\Gamma q_2 $, where $ \Gamma=\gamma_\mu \gamma_5 $ for $ {\cal{O}}=A_\mu $. Only the connected quark contractions are considered for two-point functions of charmonia.

In addition, one needs to determine a normalization constant $ Z_A $ for $ A_\mu(x) $, because with finite lattice spacings, $ A_\mu $ is no longer a conserved current in the chiral limit. The numerical determination of $ Z_A $ introduces additional uncertainties.

For chiral lattice fermions, such as domain-wall and overlap fermions, one can use the partially conserved axial vector current (PCAC) relation to avoid the computation of $ Z_A $. One can get $ f_P $ from the matrix element of the pseudoscalar density as

where $ m_{1,2} $ are the quark masses and $ m_P $ is the mass of the pseudoscalar meson. For overlap fermions the renormalization constants of the quark mass and pseudoscalar density $ P=\bar{q}_1\gamma_5 q_2 $ cancel each other ($ Z_m^{-1}=Z_P $), of which the numerical verification can be found in, for example, Ref. [33]. Then, we can use $ \Gamma=\gamma_5 $ to calculate the two-point function in Eq. (5) and obtain $ f_P $ from Eq. (6).

Vector mesons have two decay constants $ f_V $ and $ f_V^T $, which are defined by the matrix element of the vector and tensor currents, respectively, between the vacuum and a vector meson V as

where $ \epsilon_\mu(p,\lambda) $ is the polarization vector of meson V with helicity λ, and $ \sigma_{\mu\nu}=\frac{\text{i}}{2}[\gamma_\mu,\gamma_\nu] $. $ \Gamma=\{\gamma_i,\sigma_{0i}\} $ are used in Eq. (5) to extract the decay constants $ f_V $ and $ f_V^T $, respectively.

The local vector and tensor currents on the lattice should be renormalized by the constants $ Z_V $ and $ Z_T $, respectively. For overlap fermions, we have $ Z_V = Z_A $, which was verified numerically in Refs. [34−36] on various gauge ensembles. The renormalization constants $ Z_A $ and $ Z_T $ for this study have been calculated in Ref. [37] as

where we give values of $ Z_T $ in the commonly used $ {{\overline{{\rm{MS}}}}} $ scheme and at a scale of $ \mu=2 $ GeV.

To obtain the two-point functions for various mesons, we calculate quark propagators with a range of masses from the light to charm quark on three ensembles by using Z3-random point sources. The spatial locations of the point sources on ensembles $ \mathrm{f006} $ and $ \mathrm{f008} $ are randomly chosen to reduce correlation. On each of these two ensembles, 16 equally-distributed time slices per configuration are used to set the point sources. On ensemble $ \mathrm{f004} $, one source at the origin of the lattice on each configuration is used. In total, more than six hundred measurements are done on every ensemble (see Table 1). The valence quark masses $ am_{l,s,c} $ in lattice units are given in Table 2.

$ am_l $	0.00460, 0.00585, 0.00677, 0.00765, 0.00885
	0.01120, 0.01290, 0.01520, 0.01800, 0.02400
$ m_\pi $	$ \sim $ 220 $ - $ 500 MeV
$ am_s $	0.037, 0.040, 0.043, 0.046, 0.049, 0.052
$ am_c $	0.450, 0.492, 0.500, 0.550

Table 2.  Valence quark masses used in this work. The physical mass point of the valence charm quark is estimated to be approximately 0.492 in lattice units (see below).

The physical mass points of light, strange, and charm quarks are fixed by using the meson masses $ m_\pi^{\rm expt} $, $ m_K^{\rm expt} $, and $ m_{J/\psi}^{\rm expt} $, measured experimentally. However, the corrections from the difference $ m_d-m_u $ and the electromagnetic effects should be removed from the experimental values. The corresponding values are defined to be $ m_\pi^{\rm phys}\equiv m_{\pi^0}^{\rm expt} $ and ($ m_K^{\rm phys})^2\equiv \dfrac{1}{2}[(m_{K^+}^{\rm latt})^2+(m_{K^0}^{\rm latt})^2] $ in isospin-symmetric QCD with $ m_{K^+}^{\rm latt}= $ 491.405 MeV and $ m_{K^0}^{\rm latt}= $ 497.567 MeV [38−40]. We then use $ m_\pi^{\rm phys}= 134.98 $ MeV, $ m_K^{\text{phys}}=494.49 $ MeV, and $ m_{J/\psi}^{\rm phys}=3.0969 $ GeV to determine the physical mass points. Note that at the physical pion mass, $ D^* $ decays into the P-wave $ D\pi $ state through the strong interaction. In Table 3, we have listed some masses of $ D^{(*)} $, with different valence pion masses on our lattice. At each valence pion mass, the mass of the vector meson $ D^* $ is below the $ D\pi $ threshold. Therefore, in this study, $ D^* $ is stable on our lattice. Similarly, $ D_s^* $ is also a stable particle.

$ am_l $	0.00460	0.00585	0.00677	0.00765	0.00885	$ am_c $
$ M_{D^*} $	1.937	1.940	1.941	1.943	1.945	0.450
$ M_{D} $	1.767	1.769	1.770	1.772	1.774	0.450
$ m_{\pi} $	0.221	0.250	0.269	0.285	0.307

Table 3.  Masses of $ D^{(*)} $ in units of GeV, with the corresponding valence pion masses listed, at different valence quark masses.

As mentioned earlier, 16 sources separated by four time slices are used on each configuration of ensembles $ \mathrm{f006} $ and $ \mathrm{f008} $. We check the autocorrelations among the measurements by computing the normalized autocorrelation function of two-point functions $ C(t) $ at several chosen t. For example, we compute the autocorrelation function $ \rho(i;t) $ of the two-point function $ C(t) $ calculated on the $ \mathrm{f008} $ ensemble with $ {\cal{O}}=\bar{s}\gamma_5 c $ and $ am_s/am_c= $ 0.037/0.450

where $ i,j $ label the measurements. We choose $ t=0,16 $, and the variation of $ \rho(i;t) $ with i is shown in Fig. 1. It is observed that $ \rho(i;t) $ decays very rapidly with the measurement separation i, indicating little correlation among the measurements of the two-point function. When estimating the statistical errors in our data analyses using the Jackknife resampling method, removing 16 measurements from the same configuration on ensemble $ \mathrm{f006} $ and $ \mathrm{f008} $ for each resampling is sufficiently safe.

Figure 1.  (color online) The behavior of autocorrelation functions with measurement interval i. The gray lines are schematic exponential functions $\rm e^{-i} $.

We now turn to fitting the two-point function in Eq. (5), of which the numerical result can be obtained by using the quark propagators. After the intermediate state insertion, the spectral expression of the two-point function under the temporal boundary condition reads

where $ N_{\cal{O}} $ denotes the number of all hadron states with the same quantum numbers as the operator $ {\cal{O}} $. At large t, the contribution from the ground state dominates the two-point function. Thus, by fitting $ C(t) $ through the function form on the right hand side of Eq. (10), we can extract the mass $ m_0 $ of the ground state and the decay constant that is encoded in $ A_0 $ by the relation $ A_0\equiv \left|\langle0|{\cal{O}}|H_0\rangle\right|^2/(2m_0) $. In doing so, we perform a correlated minimal-$ \chi^2 $ fit, and the statistical errors are estimated through the Jackknife analysis. Practically, we fold the data along $ T/2=32 $ before the fitting and have taken the following steps to ensure the stability of the fitting result:

1. Vary $ t_{\text{min}} $ for the fitting range $ [t_{\text{min}}, T/2] $ to select a stable fitting result.

2. Ensure that $ \chi^2/\text{d.o.f}\lesssim 1.0 $.

3. Use a constrained multi-state fit [41] to obtain the ground state mass $ m_0 $ again, ensuring that it does not change with the increase of states in the fitting function.

4. Verify that the plateau of the effective mass $ m_{\text{eff}} $ obtained from $ C(t) $ is consistent with the fitted $ m_0 $ from both the one-state and multi-state fit.

Taking the fit for $ M_{D_s} $ on ensemble $ \mathrm{f004} $ as an example, the left panel of Fig. 2 shows the fitted mass $ aM_{D_s} $ in lattice units obtained by varying $ t_{\rm min} $ for the fitting range $ [t_{\text{min}}, T/2] $. Here, as $ t_{\rm min} $ decreases, the contamination from excited states becomes more and more significant. It can be seen that $ aM_{D_s} $ shifts upwards when $ t_{\text{min}} $ is less than 14. We have chosen the fitting range as $ 15\leq t \leq 32 $, with the corresponding $ \chi^2/\text{d.o.f}=0.9 $. The blue points on the middle panel of Fig. 2 represent the effective mass defined by

Figure 2.  (color online) Left: fitted $ aM_{D_s} $ as a function of $ t_{\text{min}} $ with fixed $ t_{\text{max}}=32 $. Middle: effective masses (blue points) from the two-point function, compared with $ aM_{D_s} $ from the fit (orange band) with fitting range $ 15\leq t \leq 32 $. Right: fitted $ aM_{D_s} $ as a function of the number of states from the multi-state fit.

while the orange band represents $ aM_{D_s} $ obtained from the single-state fit with the fitting range $ 15\leq t \leq 32 $. It can be observed that the two are consistent within the margin of error. The right panel of Fig. 2 shows $ aM_{D_s} $ from the multi-state fits as a function of the number of states n in the fitting function [41]. After n exceeds $ 4 $, the fitting results stabilize and agree with that from the one-state fit.

For the other charmed mesons and charmonia, we have adopted similar procedures to determine their fitting ranges. We found that different valence quark masses do not significantly affect the behavior of the effective mass with respect to the lattice time t. Therefore, for the same meson with different valence quark masses, we have chosen the same value of $ t_{\text{min}} $. The fitted results for $ aM_{D_s} $ from one-state fits on ensemble $ \mathrm{f004} $ and the statistical errors obtained by Jackknife analyses are listed in Table 4. These 24 different meson masses are not the final physical results, as the aim is to obtain the meson masses and decay constants at the physical mass point of valence quarks.

$ am_s $	0.037	0.040	0.043	0.046	0.049	0.052	$ am_c $
$ aM_{D_s} $	0.7733(9)	0.7763(9)	0.7793(8)	0.7824(8)	0.7854(8)	0.7884(8)	$ 0.450 $
	0.8179(10)	0.8208(9)	0.8238(9)	0.8267(9)	0.8297(8)	0.8327(8)	$ 0.492 $
	0.8263(10)	0.8293(9)	0.8322(9)	0.8352(9)	0.8381(9)	0.8411(8)	$ 0.500 $
	0.8791(11)	0.8820(10)	0.8849(10)	0.8878(10)	0.8907(9)	0.8935(9)	$ 0.550 $

Table 4.  Fitted masses of $ D_s $ at various valence quark masses, with statistical uncertainties estimated by Jackknife analyses.

We also obtained the masses of pion and kaon. Then, we used the aforementioned $ m_\pi^{\rm phys} $, $ (m_{ss}^{\rm phys})^2\equiv2(m_K^{\rm phys})^2- (m_\pi^{\rm phys})^2 $, and $ m_{J/\psi}^{\rm phys} $ as inputs to determine the physical meson masses and decay constants. From Table 5, $ (am_\pi)^2 $ shows a linear dependence on the light quark mass, while $ (am_{ss})^2 $ is independent of it but shows a linear dependence on the strange quark mass, which are consistent with the expectations of the chiral perturbation theory. Notably, another method exists for determining the physical mass point of the strange quark, which is conducted by calculating $ m_{\eta_s} $ directly, instead of $ m_{ss} $. Here $ \eta_s $ is a fictitious $ s\bar s $ pseudoscalar meson whose two-point functions are calculated by considering only the contributions from the QCD-connected contraction. By using $ m_{\eta_s}^{\rm phys}=689.89 $ MeV as an input [42, 43], we obtained $r_{\rm phys}\equiv (m_{\eta_s}^{\rm phys})^2/ (m_{ss}^{\rm phys})^2=1.004(16)$ and $ f_{\eta_s}= $ 179.6(4.1) MeV, which are consistent with the previous results [40, 44, 45], and found that the final results of the decay constants of charmed mesons change by less than 0.2%. Table 6 lists some fitted results of $ am_{\eta_s} $ and r on ensemble $ \mathrm{f004} $ as examples.

$ am_l $	0.0046	0.0585	0.0677	$ \cdots $	0.0152	0.018	0.024	$ am_s $
$ (am_{\pi})^2 $	0.0086(2)	0.0110(2)	0.0127(2)	$ \cdots $	0.0282(2)	0.0333(3)	0.0442(3)
$ (am_{ss})^2 $	0.0696(10)	0.0696(9)	0.0696(9)	$ \cdots $	0.0697(8)	0.0698(8)	0.0700(8)	$ 0.037 $
$ (am_{ss})^2 $	0.0751(10)	0.0751(9)	0.0751(9)	$ \cdots $	0.0753(8)	0.0754(8)	0.0757(8)	$ 0.040 $
$ (am_{ss})^2 $	0.0807(10)	0.0807(10)	0.0807(9)	$ \cdots $	0.0809(9)	0.0810(8)	0.0814(8)	$ 0.043 $
$ (am_{ss})^2 $	0.0862(10)	0.0863(9)	0.0863(9)	$ \cdots $	0.0866(9)	0.0867(9)	0.0871(8)	$0.046$
$ (am_{ss})^2 $	0.0918(11)	0.0919(10)	0.0919(10)	$ \cdots $	0.0922(9)	0.0924(9)	0.0928(9)	$ 0.049 $
$ (am_{ss})^2 $	0.0974(11)	0.0975(11)	0.0975(10)	$ \cdots $	0.0979(9)	0.0981(9)	0.0986(9)	$ 0.052 $

Table 5.  Fitted results of $ m_\pi^2 $ and $ m_{ss}^2 $ at various valence quark masses, with statistical uncertainties estimated by Jackknife analyses.

$ am_s $	0.037	0.040	0.043	0.046	0.049	0.052
r	0.996(8)	0.999(8)	1.002(8)	1.005(8)	1.008(7)	1.011(7)
$ am_{\eta_s} $	0.2632(7)	0.2739(7)	0.2844(7)	0.2945(7)	0.3043(7)	0.3140(7)

Table 6.  Fitted results of $ r\equiv m_{\eta_s}^2/(2m_K^2-m_{\pi}^2) $ and $ am_{\eta_s} $ at various strange quark masses on ensemble $ \mathrm{f004} $, with statistical uncertainties estimated by Jackknife analyses.

After obtaining the meson masses and decay constants at different valence quark masses on ensemble $ \mathrm{f004} $, we used linear interpolation (extrapolation) to obtain the results at the physical mass point. The Particle Data Group [1] and Fermilab Lattice and MILC collaborations [39] give

Combining the above with the inverse lattice spacing $ a^{-1}=2.383(9)\, $GeV, we have

Here, we regard the error from hadron masses as tiny and negligible, compared with that from our lattice spacing. The form of the linear interpolation function is as follows:

with $ \Delta am_{J/\psi}=am_{J/\psi}-(am_{J/\psi})_{\rm phys} $ and $ \Delta (am_{ss})^2=(am_{ss})^2- (am_{ss})^2_{\rm phys} $ as inputs, and $ b_1 $, $ b_2 $ and $ aM_{D_s}^{\rm fit} $ as fitting parameters. On ensemble $ \mathrm{f004} $, we can then obtain

where the statistical error also comes from Jackknife analyses, considering the correlation between different valence quark masses. An illustrative plot of the interpolation for the 24 data points of $ aM_{D_s} $ in Table 4 is shown in Fig. 3. The four rows of blue points in the figure correspond to the variation of $ aM_{D_s} $ with respect to $ (am_{ss})^2 $ for the four different charm quark masses. The red point represents the result of the interpolation at the physical mass point. Eq. (12) describes the data well, and the physical mass point of the valence charm quark is approximately $ am_c=0.492 $.

Figure 3.  (color online) Distribution of $ aM_{D_s} $ (blue dots) with respect to $ \Delta (am_{ss})^2 $ and/or $ \Delta am_{J/\psi} $ and linear interpolation of $ aM_{D_s} $ (gray lines) to the physical point (red dot).

The computed $ D^{(*)} $ masses on configuration $ \mathrm{f004} $ are listed in Table 7. Similarly, we used linear interpolation and extrapolation methods to calculate the $ D^{(*)} $ mass at the physical mass point. The fitting function for the interpolation and extrapolation is given by

$ am_l $	0.0046	0.0585	0.0677	$ \cdots $	0.0152	0.018	0.024	$ am_c $
$ aM_{D} $	0.7417(28)	0.7423(24)	0.7429(22)	$ \cdots $	0.7507(14)	0.7535(13)	0.7595(11)	$ 0.450 $
	0.7867(31)	0.7873(27)	0.7879(25)	$ \cdots $	0.7955(15)	0.7982(14)	0.8042(12)	$ 0.492 $
	0.7953(31)	0.7958(27)	0.7964(25)	$ \cdots $	0.8040(16)	0.8067(14)	0.8126(12)	$ 0.500 $
	0.8484(35)	0.8489(30)	0.8495(28)	$ \cdots $	0.8570(17)	0.8596(16)	0.8655(13)	$ 0.550 $
$ aM_{D^{*}} $	0.8129(56)	0.8140(49)	0.8147(46)	$ \cdots $	0.8211(28)	0.8236(25)	0.8292(21)	$ 0.450 $
	0.8541(58)	0.8553(51)	0.8560(47)	$ \cdots $	0.8625(29)	0.8650(26)	0.8706(22)	$ 0.492 $
	0.8620(59)	0.8631(52)	0.8639(48)	$ \cdots $	0.8704(29)	0.8729(26)	0.8785(22)	$ 0.500 $
	0.9112(62)	0.9125(54)	0.9133(50)	$ \cdots $	0.9199(31)	0.9223 (28)	0.9279(23)	$ 0.550 $

Table 7.  $ M_{D^{(*)}} $ with statistical uncertainties estimated by Jackknife analyses.

Here, the mass of $ D^{(*)} $ depends on the light quark mass, so we included $ b_3\Delta (am_{\pi})^2 $ in Eq. (14), where $\Delta (am_{\pi})^2= (am_\pi)^2-(am_\pi)^2_{\rm phys}$. The fitting results are

The relative error of $ aM_{D^{(*)}}^{\rm fit} $ is larger than that of $ aM_{D_s^*} $ because our valence pion masses ranging from 220 to 500 MeV are larger than the physical one, and thus an extrapolation in $ (am_\pi)^2 $ is needed.

We repeated the previous steps on ensemble $ \mathrm{f006} $ and $ \mathrm{f008} $, and combined the results at the physical mass points with those from $ \mathrm{f004} $, all of which are listed in Table 8. For the mass of $ \eta_c $, an interpolation only on $ am_{J/\psi} $ is performed.

Label	$ am_\pi^{\rm sea} $	$ aM_D^{\rm fit} $	$ aM_{D^*}^{\rm fit} $	$ aM_{D_s}^{\rm fit} $	$ aM_{D_s^*}^{\rm fit} $	$ aM_{\eta_c}^{\rm fit} $
$ \mathrm{f004} $	0.1269(4)	0.7837(30)	0.8521(57)	0.8258(10)	0.8914(13)	1.24967(59)
$ \mathrm{f006} $	0.1512(3)	0.7836(27)	0.8644(34)	0.8261(9)	0.8951(12)	1.24796(59)
$ \mathrm{f008} $	0.1727(4)	0.7843(25)	0.8684(32)	0.8264(8)	0.8964(18)	1.24884(52)

Table 8.  Fitted masses of charmed mesons and charmonia on ensembles $ \mathrm{f004} $, $ \mathrm{f006} $, and $ \mathrm{f008} $ at physical mass point of valence quarks, with $ m_\pi^{\rm sea} $ representing the pion masses corresponding to the light sea quark masses [29].

The sea quark masses $ m_l^{\rm sea} $ for none of the three ensembles are at the physical mass point. Therefore, we need to perform another round of linear extrapolations to obtain the final physical results (taking $ aM_{D_s} $ as an example),

where $ aM_{D_s}^{\rm fit} $ and $ (am_\pi^{\rm sea})^2 $ are taken from Table 8, $ m_\pi^{\rm sea} $ refers to the pion masses corresponding to the light sea quark masses [29], and $ aM_{D_s}^{\rm phys} $ is a fitting parameter. In Eq. (16), $ aM_{D_s} $ can be replaced by other observables.

Figure 4 shows the extrapolation of $ aM_{D_s} $ with respect to the sea quark masses by using Eq. (16). Meson masses do not get corrections from renormalization, so we can multiply the results in lattice units by the inverse lattice spacing $ a^{-1}=2.383(9)\, $GeV to obtain the physical results,

Figure 4.  (color online) Extrapolations of $ aM_{D_s}^{\rm fit} $ (blue) to the physical mass point of light sea quarks (red) with uncertainties (gray band). The horizontal coordinates are the pion masses squared $ (am_\pi^{\rm sea})^2 $ in lattice units corresponding to the sea quark masses on ensembles $ \mathrm{f004} $, $ \mathrm{f006} $, and $ \mathrm{f008} $.

Here, the first error is statistical and the second one is from the uncertainty of the lattice spacing. The statistical errors for $ D^{(*)} $ are much larger than those for $ D^{(*)}_s $ because the physical mass point of the light valence quark is implemented by extrapolation rather than interpolation.

As shown in Fig. 5, our meson masses are consistent with the values in the Particle Data Group within $ 1\sigma $ deviation:

Figure 5.  (color online) Comparisons of our meson masses with experimental values. The experimentally measured masses $ m_\pi^{\rm expt} $, $ m_K^{\rm expt} $, and $ m_{J/\psi}^{\rm expt} $ are used as inputs. The right panel shows the ratios of our meson masses to their experimental values. Our results are in good agreement with the experimental measurements.

The right panel of Fig. 5 shows the ratios of our meson masses to their experimental values, which are consistent with one. Compared to the results at a coarse lattice spacing $ a^{-1}=1.730(4)\, $GeV [28]:

Our current results are also consistent within $ 1\sigma $ deviation. The largest difference in the center values between the current and previous results is in the mass of $ D^* $, which is approximately 1%. Therefore, we estimate the discretization error of our masses of charmed mesons and charmonia to be approximately 1%. This confirms the estimation of the discretization error given in the previous work [28].

Before analyzing the decay constants, we first provide the renormalization constants for the vector and tensor currents. As mentioned before, we use chiral lattice fermions in this study. Thus, the decay constants $ f_P $ obtained from the two-point functions of pseudoscalar operators do not get corrections from renormalization. The tensor current is renormalized by $ Z_T $ in the $ \overline{\text{MS}} $ scheme at the scale of $ 2\, $GeV, and the renormalization constant for the vector current is the same as that of the axial-vector, viz. $ Z_V=Z_A $, of which the numerical results are copied from Ref. [37] in Eq. (8).

Following the fitting procedure introduced in the previous section, we use single-state fits, as in Eq. (10), for two-point correction functions $ C(t) $. The decay constant can be obtained from the amplitude $ A_0 $. For pseudoscalar mesons, the decay constant is given by

where $ q_1 $ and $ q_2 $ represent the quark components in the pseudoscalar operator. For vector mesons, the decay constant (before renormalization) is given by

from fits to the correlation functions with the vector current inserted.

In addition to calculating the decay constants themselves, we also performed joint fits of two-point functions and obtained the ratio of decay constants $ f_V/f_P $ and $f^{\rm{T}}_V/f_V$.

The subsequent steps for going to the physical mass point are similar to those in the previous section. First, we collected decay constants and their ratios on different ensembles at different valence quark masses. Some of the results are listed in Table 9. We then performed a linear interpolation (extrapolation) of the results on each ensemble. Taking $ f_{D_s} $ on $ \mathrm{f004} $ as an example, the interpolation function takes the form similar to Eq. (12)

$ am_s $	0.037	0.040	0.043	0.046	0.049	0.052	$ am_c $
$ af_{D_s} $	0.1030(9)	0.1040(9)	0.1050 (9)	0.1059 (9)	0.1067 (9)	0.1076 (9)	$ 0.45 $
	0.1041(9)	0.1050(9)	0.1060(9)	0.1069(9)	0.1078 (9)	0.1087 (9)	$ 0.492 $
	0.1042(9)	0.1052(9)	0.1062(9)	0.1071(9)	0.1080(9)	0.1089(9)	$ 0.50 $
	0.1054(10)	0.1064(10)	0.1073(10)	0.1083(10)	0.1092(9)	0.1102(9)	$ 0.55 $

Table 9.  Decay constants of $ D_s $ extracted from $ \chi^2 $-fit on ensemble $ \mathrm{f004} $, with statistical uncertainties estimated by Jackknife analyses.

The corresponding figure is shown in Fig. 6, where the abscissa values are of $ \Delta am_{J/\psi}+\Delta (am_{ss})^2 $. The figure shows that the linear fitting function describes our data well. Some decay constants and ratios fitted at the physical mass point are collected in Table 10.

Figure 6.  (color online) Distribution of $ af_{D_s} $ (blue dots) with respect to $ \Delta am_{J/\psi}+\Delta (am_{ss})^2 $, with linear interpolation to the physical mass point (red dot).

Label	$ am_\pi^{\text{sea}} $	$ (af_D)_{\text{fit}} $	$ (af_{D^*})_{\text{fit}} $	$ (af_{D^*}^T)_{\text{fit}} $	$ (f_{D^*}/f_D)_{\text{fit}} $	$ f^T_{D^*}/f_{D^*}(\text{fit}) $
$ \mathrm{f004} $	0.1269(4)	0.0907(18)	0.0944(38)	0.0797(24)	1.041(41)	0.845(28)
$ \mathrm{f006} $	0.1512(3)	0.0920(18)	0.1049(22)	0.0870(25)	1.139(26)	0.829(17)
$ \mathrm{f008} $	0.1727(4)	0.0915(18)	0.1079(19)	0.0918(15)	1.179(29)	0.850(9)
Label	$ am_\pi^{\text{sea}} $	$ (af_{D_s})_{\text{fit}} $	$ (af_{D_s^*})_{\text{fit}} $	$ (af_{D_s^*}^T)_{\text{fit}} $	$ (f_{D_s^*}/f_{D_s})_{\text{fit}} $	$ (f^T_{D_s^*}/f_{D_s^*})_{\text{fit}} $
$ \mathrm{f004} $	0.1269(4)	0.1063(10)	0.1103(11)	0.0937(10)	1.038(12)	0.8492(43)
$ \mathrm{f006} $	0.1512(3)	0.1061(9)	0.1126(12)	0.0963(11)	1.061(14)	0.8550(63)
$ \mathrm{f008} $	0.1727(4)	0.1054(8)	0.1131(16)	0.0962(14)	1.073(17)	0.8505(71)
Label	$ am_\pi^{\text{sea}} $	$ (af_{\eta_c})_{\text{fit}} $	$ (af_{J/\psi})_{\text{fit}} $	$ (af_{J/\psi}^T)_{\text{fit}} $	$ (f_{J/\psi}/f_{\eta_c})_{\text{fit}} $	$ (f^T_{J/\psi}/f_{J/\psi})_{\text{fit}} $
$ \mathrm{f004} $	0.1269(4)	0.1748(12)	0.1704(11)	0.1527(11)	0.9754(70)	0.8960(14)
$ \mathrm{f006} $	0.1512(3)	0.1753(9)	0.1760(10)	0.1573(10)	1.004(66)	0.8935(14)
$ \mathrm{f008} $	0.1727(4)	0.1731(8)	0.1699(11)	0.1521(10)	0.9818(57)	0.8950(16)

Table 10.  Decay constants before renormalization and their ratios on ensembles $ \mathrm{f004} $, $ \mathrm{f006} $, and $ \mathrm{f008} $ of charmed mesons and charmonia fitted at the physical mass point of valence quarks, with $ m_\pi^{\text{sea}} $ representing the pion masses corresponding to the sea quark masses [29].

Finally, we performed a linear extrapolation for the sea quark mass $ m_l^{\rm sea} $ to obtain the decay constants at the physical sea quark mass point, as shown in Fig. 7. We converted the decay constants in lattice units to physical units and multiplied them by the appropriate renormalization constants. The physical results are as follows (in units of MeV):

Figure 7.  (color online) Examples of extrapolations of decay constants or ratios (blue) to the physical mass point (red) of light sea quarks with uncertainties (gray band). The abscissa values are of pion masses squared $ (am_\pi^{\rm sea})^2 $ in lattice units corresponding to the light sea quark masses on ensembles $ \mathrm{f004} $, $ \mathrm{f006} $, and $ \mathrm{f008} $.

Here, for each of the decay constants, the first error includes statistical and extrapolation/interpolation uncertainties, the last error is from the uncertainty of the lattice spacing, and the error in the middle, if there is one, is due to the uncertainty in the renormalization constant. It is seen in Eq. (8) that the uncertainty of $ Z_V(=Z_A) $ can be regarded as negligible. The uncertainty of $ Z_T $ is approximately 1% and is much smaller than the first error.

For chiral fermions as used in this study, the discretization effects are proportional to the squared lattice spacing $ a^2 $. With three or more lattice spacings we would be able to do linear fittings in $ a^2 $ and obtain results in the continuum limit. The lattice spacing squared, corresponding to our current study ($ a^2\approx 0.007 $ fm$ ^2 $), happens to lie close to the middle point between the continuum limit and the previous work [28] ($ a^2\approx 0.013 $ fm$ ^2 $). A simple linear extrapolation in $ a^2 $ using the results at the two lattice spacings has no degree of freedom and will shift our current result by an amount almost equal to the difference between the two lattice spacings. Therefore, we think the discretization errors can be estimated by simply comparing the results of these two works. For $ f_{D_s^{(*)}} $, our results are consistent within $ 1\sigma $ deviation with the previous results, respectively, where $ f_{D_s}=249(5)\, $MeV and $ f_{D_s^*}= $ 274(5) MeV [28] and we have dropped the discretization error assigned to the previous work. For $ f_{D^{(*)}} $, our results are also consistent within $ 1\sigma $ deviation with the previous results, respectively, where $ f_{D}=213(2)\, $MeV and $ f_{D^{*}}=234(3)\, $MeV. The differences in the central values are approximately from 1%−4%. Therefore, we take an average value of 3% as the estimate of the discretization errors in our decay constants and get (in units of MeV),

where the first error is the square root of the quadratic sum of the errors in Eq. (21), and the second error is the 3% discretization error, for each of the decay constants.

The decay constant ratios of the charmed mesons are collected in Table 11 and Table 12, which are consistent with those obtained in our previous work on a coarser lattice. The first error in the table considers the contributions from statistics, interpolation/extrapolation, and renormalization. The second error is our estimate of the discretization error from the center value differences of this work and the previous work. The ratios involving $ D^{(*)} $ have larger errors than those from our previous work because of the extrapolation with respect to the light quark mass (both valence and sea). The result of $f_{D_s^*}^{\rm{T}}/f_{D_s^*}$ from this work has a smaller error because the uncertainty from the renormalization factor $ Z_T^{{{\overline{{\rm{MS}}}}}} $(2 GeV)$ /Z_A $ is now smaller.

This work	$ D^*/D $	$ D_s^*/D_s $	$ J/\psi/\eta_c $	[28]	$ D^*/D $	$ D_s^*/D_s $
$ f_V/ f_{P} $	1.045(83)(55)	1.097(30)(3)	1.060(15)	$ f_V/ f_{P} $	1.10(2)(2)	1.10(3)(2)
$ f_V^T/f_V $	0.872(47)(38)	0.909(14)(11)	0.961(10)	$ f_V^T/f_V $	0.91(3)(2)	0.92(3)(2)

Table 11.  Ratios of decay constants. The discretization errors of the results from [28] were estimated to be 2%.

	This work	[28]
$ f_{D_s}/ f_D $	1.185(45)(23)	1.163(14)(23)
$ f_{D_s^*} / f_{D^*} $	1.231(73)(61)	1.17(2)(2)

Table 12.  Ratios of decay constants for $S U$(3) flavor symmetry breaking effects.

In addition, we calculated the decay constants and their ratios for $ \eta_c $ and $ J/\psi $. We used the renormalization constants, as shown in Eq. (8). The decay constants themselves are 6%−9% higher than the values in the continuum limit obtained by other lattice groups, such as the HPQCD collaboration [27]. We think this is due to large discretization effects in the decay constants of charmonia, which comprise two heavy quarks. In the work of HPQCD [27], sizable discretization effects were also spotted in those decay constants at non-zero lattice spacings.

We observed smaller discretization effects in ratios of decay constants for charmonia. The physical results are

Here, the error includes statistical uncertainty and the one in renormalization constants. Combining these values with the experimental measured $ f_{J/\psi}^{\rm expt}=407(4)\, $MeV, which is obtained from its pure leptonic decay width $ \Gamma(J/\psi\to e^+e^-) $ [27], we can derive $ f_{\eta_c}=383.8(6.7)\, $MeV and $f^{\rm{T}}_{J/\psi}=391.1(5.4)\,$MeV. These results are consistent within $ 1.5\sigma $ deviation with the results $ f_{\eta_c}= $ 394.7(2.4) MeV [25], $f^{\rm{T}}_{J/\psi}=$ 392.7(2.7) MeV and $f^{\rm{T}}_{J/\psi}/ f_{J/\psi}=$ 0.9569(52) [46] from the HPQCD collaboration obtained in the continuum limit. The better agreement in the ratio demonstrates that the discretization effects shrunk by the cancellation between the numerator and the denominator. Because we do not have the results of decay constants of charmonia from the previous work to compare, we do not include discretization errors for those constants and ratios in this work.

In the study of heavy flavor physics, decay constants of mesons are fundamental and important quantities. They are essential input parameters of theoretical calculations of semileptonic decays, such as $ B_{(c)}\to D,\eta_c $, as well as pure or non-leptonic decays of D mesons, by using either the QCD-based factorization theory or phenomenological models, such as topological diagram methods. Such decay constants are also crucial for extraction of CKM matrix elements and precise tests of the Standard Model. Theoretical calculations of decay constants are essential and lattice QCD provides a systematic way to improve their precision. Depending on only fundamental parameters of QCD without any modelling assumptions, lattice QCD is considered the best method for computing decay constants.

In this study, we have computed decay constants and their ratios for D, $ D^* $, $ D_s $, $ D_s^* $, $ \eta_c $, and $ J/\psi $, using $ 2+1 $-flavor configurations. We have also calculated the decay constants with tensor currents of vector mesons and provided estimates of discretization errors. Our final results for the decay constants adding up all errors are given in Table 13 (in units of MeV), of which the ratios are given in Tables 11 and 12 for charmed mesons, and in Eq. (23) for charmonia. The ratios in Tables 11 and 12 reflect the magnitudes of heavy quark symmetry breaking and $S U$(3) flavor symmetry breaking, respectively. The precision of $f_{D_s^*}^{\rm{T}}/f_{D_s^*}$ is improved compared with that from our previous study [28]. We compare the $ f_{D^*_{(s)}} $ in this study and other lattice QCD calculations in Table 14.

	D	$ D^* $	$ D_s $	$ D_s^* $	$ \eta_c $	$ J/\psi $
$ f_{P/V} $	215(11)	224(18)	255.7(8.9)	277(11)	383.8(6.7)	−
$f_V^{\rm{T}}$	−	191(14)	−	252(10)	−	391.1(5.4)

Table 13.  Decay constants in units of MeV, among which the $ f_{\eta_c} $ and $f_{J/\psi}^{\rm{T}}$ require $ f_{J/\psi} $ as an additional input. The ones with tensor currents are renormalized in the $ {{\overline{{\rm{MS}}}}} $ scheme at the scale of 2 GeV.

	This work	[28]	[11]	[13]	[14]	[12]
$ f_{D^*} $	224(18)	234(6)	−	−	223.5(8.4)	278(17)
$ f_{D_s^*} $	277(11)	274(7)	274(6)	264(15)	268.8(6.6)	311(9)

Table 14.  Comparisons of $ f_{D_{(s)}^*} $ (in MeV).

By using our result of $ f_{D_s^*}=277(11) $ MeV, one can obtain the decay width of the pure leptonic decay of $ D_s^* $, which is $ \Gamma(D_s^*\to \ell\nu_\ell)|_{\ell=e,\mu}=2.5(2)\times10^{-6} $ keV, as HPQCD did in [11]. Combining this result with the total decay width of $ D_s^* $, which is $ \Gamma_{\rm tot}(D_s^*)= $0.0587(54) keV [47], we then find the branching ratio

This value can be confirmed with future experiments.

Discretization effects in the decay constants of charmed mesons are found to be larger than those in the meson masses. As for the decay constants of charmonia, the discretization effects are even larger. However, the lattice artefacts are much smaller for the decay constant ratios.

Currently, our estimation of discretization errors are based on the analyses of results from two sets of lattices with inverse lattice spacings of $ 1.730\, $ GeV and $ 2.383\, $GeV. To remove the discretization errors, we need to repeat our calculations on other sets of lattice with different spacings. Additionally, QED corrections and contributions from the QCD-disconnected contractions for charmonia are necessary for high-precision studies in the future.

We thank the RBC-UKQCD collaboration for sharing the domain wall fermion configurations. The GWU code [48, 49] is acknowledged. The computations were performed on the HPC clusters at the Institute of High Energy Physics (Beijing).