Charmed pseudoscalar decay constants on three-flavour CLS ensembles with open boundaries

TL;DR

This paper presents lattice QCD calculations of D and D_s meson masses and decay constants using CLS ensembles with open boundaries, exploring different lattice spacings and pion masses to improve continuum extrapolation accuracy.

Contribution

It introduces a detailed analysis of charmed meson decay constants on open boundary CLS ensembles, including implementation of distance preconditioning for heavy quark propagators.

Findings

Decay constants and masses for D and D_s mesons obtained at two lattice spacings.

Implementation of distance preconditioning improved the accuracy of heavy quark propagator calculations.

Ongoing efforts towards a stable continuum limit extrapolation.

Abstract

We determine the masses and pseudoscalar decay constants of D and D_s mesons employing lattice QCD with non-perturbatively O(a) improved Wilson quarks and a tree-level Symanzik-improved gauge action. Our analysis is based on the large-volume N_f=2+1 ensembles using open boundary conditions, generated within the CLS effort. The status of results presented here covers two lattice spacings, a ~ 0.0854 fm and a ~ 0.0644 fm, and pion masses varied from 420 to 200 MeV. We also report on our implementation of distance preconditioning for the calculation of heavy quark propagators and discuss the impact of the resulting accuracy improvements on the extraction of charmed meson masses and decay constants. This is part of a continuing analysis by the RQCD and ALPHA Collaborations, aiming at a stable continuum extrapolation using several lattice spacings. To extrapolate to the physical masses, we follow both, the (2*m_l+m_s)=const. and the m_s=const. line in parameter space.

Figures

Tables

Table 1. Table 1: Details of the ensembles analysed so far for the two trajectories to the physical point, keeping m ¯ ¯ 𝑚 \overline{m} fixed to the value at the symmetric point ( m ¯ = m sym ¯ 𝑚 subscript 𝑚 sym \overline{m}=m_{\rm sym} ) and keeping the renormalised strange quark mass ( m ^ s subscript ^ 𝑚 s \hat{m}_{\rm s} ) approximately equal to the physical value ( m ^ s phys superscript subscript ^ 𝑚 s phys \hat{m}_{\rm s}^{\rm phys} ). The light and strange quark hopping parameters are denoted κ l subscript 𝜅 l \kappa_{\rm l} and κ s subscript 𝜅 s \kappa_{\rm s} , respectively. The lattice volumes N t × N s 3 subscript 𝑁 𝑡 superscript subscript 𝑁 𝑠 3 N_{t}\times N_{s}^{3} , the pion ( m π subscript 𝑚 𝜋 m_{\pi} ) and kaon ( m K subscript 𝑚 𝐾 m_{K} ) masses and the statistics given by the number of molecular dynamics units ( N MD subscript 𝑁 MD N_{\rm MD} ) are also indicated.

trajectory	ensemble	${\kappa }_l$	${\kappa }_s$	$L{m}_{\pi }$	$N_t\times N_s^3$	$\frac{m_{\pi }}{\text{MeV}}$	$\frac{m_K}{\text{MeV}}$	$N_{\mathrm{MD}}$
$\beta =3.4$ [$a=0.0854\left(15\right)\text{fm}$],    $\sqrt{8t_0}/a=4.852\left(7\right)$
$\overline{m}=m_{\mathrm{sym}}$	H101	0.13675962	0.13675962	5.8	$96\times {32}^3$	422	422	8000
	H102	0.136865	0.136549339	4.9	$96\times {32}^3$	356	442	7988
	H105	0.13697	0.13634079	3.9	$96\times {32}^3$	282	467	11332
	C101	0.13703	0.136222041	4.6	$96\times {48}^3$	223	476	6208
${\hat{m}}_s={\hat{m}}_s^{\mathrm{phys}}$	H107	0.136945665908	0.136203165143	5.1	$96\times {32}^3$	368	549	6256
	H106	0.137015570024	0.136148704478	3.8	$96\times {32}^3$	272	519	6212
	C102	0.1370508458	0.136129062556	4.6	$96\times {48}^3$	223	504	6000
$\beta =3.55$ [$a=0.0644\left(11\right)\text{fm}$],    $\sqrt{8t_0}/a=6.433\left(6\right)$
$\overline{m}=m_{\mathrm{sym}}$	H200	0.137	0.137	4.4	$96\times {32}^3$	418	418	8000
	N203	0.13708	0.136840284	5.4	$128\times {48}^3$	345	441	6172
	N200	0.13714	0.13672086	4.4	$128\times {48}^3$	283	461	6800
	D200	0.1372	0.136601748	4.2	$128\times {64}^3$	199	479	4000
${\hat{m}}_s={\hat{m}}_s^{\mathrm{phys}}$	N204	0.137112	0.136575049	5.6	$128\times {48}^3$	351	544	2000
	N201	0.13715968	0.136561319	4.5	$128\times {48}^3$	284	522	6000
	D201	0.137207	0.136546436	4.1	$128\times {64}^3$	198	499	4312

Equations

$${\varphi }_4=8t_0\left(m_K^2+\frac{1}{2}{m}_{\pi }^2\right)\mathrm{and}{\varphi }_2=8t_0{m}_{\pi }^2,$$

$$⟨0\left|A_{\mu }^{\mathrm{lc}}\right|\mathrm{D}(p)⟩=\mathrm{i}{f}_{\mathrm{D}}{p}_{\mu },⟨0\left|A_{\mu }^{\mathrm{sc}}\right|{\mathrm{D}}_{\mathrm{s}}(p)⟩=\mathrm{i}{f}_{{\mathrm{D}}_{\mathrm{s}}}{p}_{\mu },$$

$$A_{\mu }^{\mathrm{qc},\mathrm{I}}=A_{\mu }^{\mathrm{qc}}+a{c}_{\mathrm{A}}\frac{1}{2}\left({\partial }_{\mu }+{\partial }_{\mu }^{\ast }\right)P^{\mathrm{qc}}.$$

$${\left(A_{\mu }^{qc}\right)}^{\mathrm{R}}=Z_{\mathrm{A}}\left[1+a\left(b_{\mathrm{A}}{m}_{\mathrm{qc}}+3{\tilde{b}}_{\mathrm{A}}\overline{m}\right)\right]A_{\mu }^{\mathrm{qc},\mathrm{I}}+\mathcal{O}(a^2)$$

$$m_{\mathrm{qc}}=\frac{1}{2}\left(m_{\mathrm{q}}+m_{\mathrm{c}}\right),3\overline{m}=m_{\mathrm{s}}+2m_{\mathrm{l}}.$$

$$C_{\mathrm{A}}(x_0,y_0)=-\frac{a^6}{L^3}\sum _{\vec{x},\vec{y}}⟨A_4^{\mathrm{qc},\mathrm{I}}(x){\left(P^{\mathrm{qc}}(y)\right)}^{†}⟩,C_{\mathrm{P}}(x_0,y_0)=-\frac{a^6}{L^3}\sum _{\vec{x},\vec{y}}⟨P^{\mathrm{qc}}(x){\left(P^{\mathrm{qc}}(y)\right)}^{†}⟩,$$

$$C_{\mathrm{A}}(x_0,y_0)\approx \frac{f_{\mathrm{qc}}^{\mathrm{bare}}}{2}A(y_0){\mathrm{e}}^{-m_{D_{\mathrm{q}}}(x_0-y_0)},C_{\mathrm{P}}(x_0,y_0)\approx \frac{{\left|A(y_0)\right|}^2}{2m_{D_{\mathrm{q}}}}{\mathrm{e}}^{-m_{D_{\mathrm{q}}}\left(x_0-y_0\right)},$$

$$f_{\mathrm{qq}}^{\mathrm{bare}}(x_0,y_0)\propto \sqrt{\frac{C_{\mathrm{A}}(x_0,y_0)C_{\mathrm{A}}(x_0,T-y_0)}{C_{\mathrm{P}}(T-y_0,y_0)}},$$

$$\Big{|}\sum_{y}(D[U]+m_{0})_{x,y}\,S^{n}(y)-\eta_{t}(x)\Big{|}<r_{\text{gl}}{}$$

$$P=\left(\begin{array}{ccccc}{p}_1& 0& \dots & \dots & 0\\ 0& p_2& 0& \dots & 0\\ 0& 0& \ddots & 0& 0\\ 0& \dots & 0& \ddots & 0\\ 0& \dots & \dots & 0& p_T\end{array}\right),$$

$$AS=\eta \mathrm{with}A=\sum _y(D[U]+m_0{)}_{x,y},$$

$$C_{\mathrm{A},\mathrm{P}}(x_0,y_0)=A_1(y_0){\mathrm{e}}^{-m_{{\mathrm{D}}_{\mathrm{q}}}(x_0-y_0)}+A_2(y_0){\mathrm{e}}^{-m^{\prime }(x_0-y_0)}+B_1(y_0){\mathrm{e}}^{-(E_{2{\mathrm{D}}_{\mathrm{q}}}-m_{{\mathrm{D}}_{\mathrm{q}}})(T-x_0+y_0)}+\dots ,$$

$$⟨O⟩=\frac{⟨O{W}_0{W}_1{⟩}_W}{⟨W_0{W}_1{⟩}_W}.$$

$$\frac{|A^{\prime }{|}^2{\mathrm{e}}^{-M^{\prime }(x_0^{\min }-y_0)}}{2M^{\prime }}<\frac{1}{4}\Delta C_{\mathrm{P}}(x_0^{\min },y_0),$$

Full text

1 Introduction

Recent experiments at the LHC, BEPC II and KEK supply us with an abundance of spectroscopic information about charmed systems with ever improving measurements of the leptonic decay rates of D and $\mathrm{D}_{\mathrm{s}}$ mesons. In order to fully exploit this progress on the experimental side and to convert it into an amended determination of the CKM matrix entries \big{|}V_{\rm cq}\big{|}, ${\rm q}={\rm d},{\rm s}$, it must be complemented on the lattice QCD side with a more precise computation of the associated low-energy hadronic matrix elements, characterised by the pseudoscalar decay constants $f_{\rm D}$ and $f_{\rm D_{s}}$, respectively. Charmed systems are challenging because of systematics, in particular, cutoff effects are generally significant. Therefore, all the more important is a fully controlled continuum extrapolation, $a\rightarrow 0$. In this proceedings we describe our ongoing efforts in calculating the charmed pseudoscalar decay constants.

2 CLS ensembles and general computational setup

Our analysis is based on the $N_{\rm f}=2+1$ Coordinated Lattice Simulations (CLS) ensembles generated with non-perturbatively ${\cal O}(a)$ improved Wilson-Sheikholeslami-Wohlert (clover) fermions and the tree-level improved Lüscher-Weisz gauge action, employing the open-source package openQCD [1].

Ensembles have been realised with lattice spacings in the range of $a\approx 0.0854-0.039\,{\rm fm}$ ($\beta=3.4-3.85$), aiming at a controlled extrapolation to the continuum limit [2, 3, 4]. Simulating with $a<0.05$ fm is possible through the use of open boundary conditions in the temporal direction, which avoids the problem of topological freezing as the continuum limit is approached [5, 6]. Readers, who are interested in the details regarding the algorithmic setup, are referred to Ref. [2]. Here, we only mention the use of twisted-mass reweighting for the two mass-degenerate light fermions [7], which prevents instabilities in the HMC algorithm due to near-zero modes of the Dirac operator, and the use of the rational approximation for simulating the strange quark. These algorithmic choices mean that the computation of observables, such as correlation functions, require reweighting, see Section 5.

So far we have analysed lattices with $\beta=3.4$ and $3.55$ corresponding to $a=0.0854(15)$ and $0.0644(11)\,{\rm fm}$, respectively [3]. For each $\beta$-value, the physical point is approached in the light and strange quark mass plane following 2 trajectories — (i) keeping the average lattice quark mass ($\overline{m}=\left(2m_{\rm l}+m_{\rm s}\right)/3$ 111Here we refer to the vector Ward identity masses, $m_{{\rm q}={\rm l(ight)},{\rm s}}=(1/\kappa_{\rm q}-1/\kappa_{\rm crit})/(2a)$, where $\kappa_{\rm crit}$ is the critical hopping parameter value at which the axial Ward identity mass in the symmetric limit, $m_{\rm l}=m_{\rm s}$, vanishes.) fixed such that the sum of the renormalised quark masses is constant up to ${\cal O}(a)$ effects, and (ii) keeping the renormalised strange quark mass approximately constant. Figure 1 displays the two trajectories for $\beta=3.4$, where the pion masses are varied from 422 down to 223 MeV. Simulations along the $\overline{m}$ constant line, first proposed by the QCDSF Collaboration [8], start from the $m_{\rm l}=m_{\rm s}=m_{\rm sym}$ symmetric point with the strange (light) quark mass becoming heavier (lighter) towards the physical point. The quark mass dependence of hadron observables is given by Gell-Mann-Okubo expansions in powers of $m_{\rm s}-m_{\rm l}$ or, alternatively, can be described by SU(3) chiral perturbation theory (ChPT). For the second trajectory, SU(2) ChPT applies. See Ref. [3] for details on how an almost constant renormalised strange quark mass was achieved in this case. As will be demonstrated in Section 6, having two trajectories available enables the extrapolation to the physical point to be tightly constrained.

Table 1 provides further details of the ensembles entering our analysis. They exhibit large spatial volumes with high statistics; in particular, $m_{\pi}L\gtrsim 4$ ($L=N_{s}a$) is satisfied by all ensembles, and more than 4000 molecular dynamic units $N_{\mathrm{MD}}$ have been generated in each case (with the exception of ensemble N204, which is still in production). For open temporal boundaries, one needs large physical time extents, since discretisation effects and the propagation of scalar particles into the bulk mean that, in general, timeslices close to the boundaries need to be discarded. The lattice spacing is set by the gradient (resp. Wilson) flow scale $t_{0}$, where we use $\sqrt{8t_{0}^{\rm phys}}=0.4144(59)(37)\mathrm{fm}$, as determined by BMW Collaboration [9], together with $\sqrt{8t_{0}}/a$-values from Refs. [2, 3]. Using the pion and kaon decay constants to set the scale, is discussed in Ref. [4]. The physical point is defined employing the combinations

[TABLE]

together with the pion and kaon masses in the isospin limit, $m_{\pi}=134.8(3)\,{\rm MeV}$ and $m_{K}=494.2(4)$ MeV, taken from the FLAG report [10]. We note that the starting value for the $\overline{m}=\mathrm{const.}$ trajectory at the symmetric point, $m_{\rm sym}$, was fixed by tuning $\phi_{4}=1.15$, slightly above the physical value of $\phi_{4}^{\rm phys}=1.117(38)$. This choice was motivated by expected discretisation and {\cal O}\big{(}(m_{\rm s}-m_{\rm l})^{2}\big{)} effects in $\phi_{4}$ along the trajectory. In addition, the computation of $t_{0}$ involved in this tuning procedure was performed at unphysical quark mass. In the final analysis, the trajectory may be found to slightly miss the physical point. Mistunings of this kind, however, can still be accounted for later by slightly adjusting the results through reweighting or via a Taylor expansion (following Ref. [4]). We leave these considerations for future study.

3 Definition of observables

Let us explain our notation and the observables analysed. The pseudoscalar decay constants $f_{\rm D}$ and $f_{\rm D_{s}}$ are defined via the matrix elements of the axial vector current between $\mathrm{D}$ and $\mathrm{D}_{\mathrm{s}}$ meson states $|\mathrm{D}(p)\rangle$ and $|\mathrm{D}_{\mathrm{s}}(p)\rangle$ at momentum $p$, respectively, and the vacuum, viz.

[TABLE]

where $A^{\rm qc}_{\mu}=\overline{q}\gamma_{\mu}\gamma_{5}c$ for quark flavours ${\rm q}={\rm l(ight)},{\rm s}$. In order to reduce leading discretisation effects to ${\cal O}(a^{2})$ in the lattice simulation, we employ the improved axial operator

[TABLE]

The pseudoscalar interpolator is given by $P^{\rm qc}=\overline{q}\gamma_{5}c$. The improvement coefficient $c_{\rm A}$ has been determined non-perturbatively in Ref. [11]. Due to the breaking of chiral symmetry by Wilson-type fermions, the axial current must be renormalised in order to ensure that the continuum axial Ward identity

[TABLE]

is satisfied, where the (bare) vector Ward identity quark mass combinations read:

[TABLE]

The renormalisation factor $Z_{\rm A}$ is taken from the non-perturbative evaluation of Ref. [12]. The coefficient $b_{\rm A}$ has been determined non-perturbatively in Ref. [13]. The same authors find $\tilde{b}_{\rm A}$ to be consistent with zero in a preliminary analysis at $\beta=3.4$. In order to perform a consistent analysis for both $\beta=3.4$ and 3.55, at present, we neglect the $\tilde{b}_{\rm A}\overline{m}$-term. Since the mass dependent corrections in Eq. (4) are anyway dominated by the term involving $m_{\rm qc}$, this choice does not introduce a significant systematic uncertainty.

In the following, we restrict the analysis to zero-momentum $\mathrm{D}$ and $\mathrm{D}_{s}$ mesons, for which only the temporal component of the axial current needs to be considered. The matrix elements in Eq. (2) are extracted from the following two-point correlation functions:

[TABLE]

with the source and sink interpolating operators inserted at timeslices $y_{0}$ and $x_{0}$, respectively. For sufficiently large time differences $x_{0}-y_{0}$ and $T-x_{0}$, for a lattice of physical time extent 222For lattices with $N_{t}$ points in temporal direction, $T=(N_{t}-1)a$ denotes the physical time extent. $T$, these correlators behave as [4]

[TABLE]

where $m_{D_{\rm q}}$ is the mass of a pseudoscalar meson containing a charm quark and an anti-quark of flavour ${\rm q}={\rm l},{\rm s}$. The (unrenormalised) lattice decay constant, $f^{\rm bare}_{\rm qc}$, corresponds to $\langle 0|A_{4}^{\rm qc,I}|\mathrm{D}_{\rm q}\rangle/m_{D_{\rm q}}$, while the source-time dependent amplitude $A(y_{0})$ encapsulates the matrix element of the pseudoscalar operator $\langle 0\left|P^{\rm qc}\right|\mathrm{D}_{\rm q}\rangle$ and effects arising from the boundary close to the source. For $y_{0}$ sufficiently far from the boundary, $A(y_{0})=\langle 0\left|P^{\rm qc}\right|\mathrm{D}_{\rm q}\rangle$ holds.

In this preliminary analysis we place the source at large enough times such that boundary effects can be neglected. In order to gain statistics, three sources are chosen corresponding to $y_{0}/a=30$, 47 and 65 for ensembles with $T/a=95$ and $\beta=3.4$ (see Table 1), while for $\beta=3.55$, the sources are placed at $y_{0}/a=27$, 47 and 68 for ensembles with $T/a=95$, and at $y_{0}/a=46$, 63 and 81 for those with $T/a=127$. Ground-state dominance at small source-sink separations is realised by employing Wuppertal smearing [14, 15] with APE-smoothed links [16] to the pseudoscalar source and sink interpolators. Figure 2 displays typical effective masses for the D meson derived from the correlation function $C_{\rm P}(x_{0},y_{0})$ for the different source positions as a function of the sink-source time separation. For small $x_{0}-y_{0}$, the results from the source in the center of the lattice ($y_{0}/a=47$) are compatible with those from the sources closer to the boundary (i.e., from the forward-propagating meson for $y_{0}=30a$ and the backward-propagating meson for $y_{0}=65a$), consistent with the absence of boundary effects. However, the latter are clearly evident for $y_{0}=47a$, as the sink time nears the boundary, while for $y_{0}=30a$ and 65$a$ the signal deteriorates before significant boundary effects are seen.

The pseudoscalar decay constant is extracted by fitting $C_{\rm A}(x_{0},y_{0})$ and $C_{\rm P}(x_{0},y_{0})$ to the functional forms given in Eq. (7), with $A(y_{0})$ independent of $y_{0}$. The fitting procedure is discussed in Section 5. Alternative approaches, for instance, employing the axial Ward identity (i.e., PCAC) quark mass [17], will also be considered in the future. We remark that for pseudoscalar mesons containing a quark and an anti-quark of the same flavour, for which the statistical noise of the correlators is constant with varying source-sink separation, one can also place the source close to the boundary and consider ratios such as

[TABLE]

see Ref. [18] for details. When extracting the decay constant from the correlation functions, we wish to maximise the time range for the fit. However, there is a well known problem with numerical accuracy towards large time separations, when computing correlation functions involving propagators of heavy quarks. This issue is discussed in the next section along with the distance preconditioning procedure of Ref. [19], which can be adopted to alleviate such difficulties.

4 Distance preconditioning of heavy quark propagators

4.1 Motivation for distance preconditioning

When numerically inverting the Dirac operator, one has to pay special attention to the quality of the obtained solutions. Internally, the solver routine checks if the condition

[TABLE]

is satisfied, where $D[U]$ is the discretized lattice Dirac operator, $m_{0}$ the bare quark mass in lattice units and $S^{n}(y)$ the approximate solution at the $n$-th iteration of the solver. $\eta_{t}(x)$ denotes a source located on a single timeslice $t$, while $r_{\text{gl}}$ is the global residuum, i.e., indicating the numerical accuracy one likes to achieve. Contributions to the norm above are negligible for heavy quarks, if one considers large time separations between the source $y_{0}$ and the sink $x_{0}$, because of the exponential decay of correlation functions $\propto\exp(-m_{0}|y_{0}-x_{0}|)$. Solutions for large time extents thus become increasingly inaccurate in this case.

In order to improve the overall precision in our charm quark propagator computations, we have implemented the so-called “Distance Preconditioning” (DP) technique, first proposed in [19]. Rather than modifying the Dirac operator directly, as it was done in the original paper by de Divitiis et al., we included this preconditioning in the Dirac operator inversion via multiplication with a diagonal matrix $P$,

[TABLE]

with $p_{i}=\exp{}\left(\alpha_{0}\cdot|{}y_{0}-x_{0}^{(i)}|{}\right)$, corresponding to timeslice $i$, and a (control) parameter $\alpha_{0}$. $P$ acts as unity matrix in spin, color and spatial coordinates. By this procedure, timeslices far away from the source are enhanced by a large exponential factor, which balances the rapid decay of the correlation function and ensures an accurate computation of the propagators even in the case of large quark masses (in lattice units) and very large time extents.

For our implementation, instead of numerically solving the un-preconditioned system

[TABLE]

we multiply it with $P$ from the left and rewrite $(PAP^{-1})(PS)=(P\eta)$ as $A^{\prime}S^{\prime}=\eta^{\prime}$. The inversion to calculate the solution $S^{\prime}=PS$ of the preconditioned system then becomes numerically much more stable, and to recover the original solution $S$, it suffices to scale $S^{\prime}$ back with $P^{-1}$.

4.2 Numerical tests

We performed several exploratory tests on three different subsets of the CLS ensembles H105, H200 and U101. These low-statistic analyses (usually 25 to 200 configurations were included in the correlator measurements) served to explore the validity of the distance preconditioning technique for two different values of $\beta$ and different quark masses. The U101 subset was chosen because of its larger temporal lattice extent of $T/a=127$.

Since there are contaminations from excited states near the source and excitations close to the boundary, the two-point functions defined in (6) have the general form

[TABLE]

where the first term mimics the ground state, the second one represents the first excited state with mass $m^{\prime}$ and the third term is the first excited boundary state with $E_{2{\rm D}_{\rm q}}\approx 2m_{{\rm D}_{\rm q}}$ (see Ref. [4]). As one moves farther away from the boundary, the excited states become strongly suppressed, and the behaviour of the two-point functions is governed by the first term.

First tests of the distance preconditioning method were performed employing the normal-equation variant of the conjugate gradient solver (CGNE), because of its simplicity and numerical robustness. Figure 3 illustrates the numerical difficulties, as the correlator $C_{\rm P}(x_{0},y_{0})$ deviates from the expected behaviour (i.e., linear on a logarithmic scale) for large time separations $|x_{0}-y_{0}|$. These results were obtained on ensemble H105 at $\beta=3.4$ with the correlation function averaged over two source positions $y_{0}/a=1$ and $T/a-1$ (DP-results were always generated utilising these source positions). The heavy-heavy correlator is evaluated for $\kappa=0.124503\approx\kappa_{\text{charm}}$. With a loose global residual of $r_{\text{gl}}=10^{-6}$ the results for the unmodified CGNE solver (in blue) show an irregular behaviour already around timeslice $\approx 26$. While this can still be partly improved via choosing a far smaller residual of order $10^{-12}$ (in green), it is obvious that brute-force methods become unfeasible on lattices with larger time extent. The outcome can, however, be greatly improved by distance preconditioning: While the results for the DP-modified CGNE solver with $r_{\text{gl}}=10^{-6}$ and control parameter $\alpha_{0}=0.4$ (red) already yield a modest improvement, those for $r_{\text{gl}}=10^{-12}$ (black) show a smooth and nearly perfect exponential decay over the whole range such that the effective pseudoscalar meson mass $am_{\text{eff}}(x_{0},y_{0})=\ln\left(\frac{C_{\rm P}(x_{0},y_{0})}{C_{\rm P}(x_{0}+a,y_{0})}\right)$ can be extracted along a wide plateau (the slight increase in $am_{\text{eff}}(x_{0},y_{0})$ at $x_{0}/a>85$ is likely due to boundary excited states back-propagating into the bulk).

In order to gain in computational speed, we also implemented the DP-method into another solver routine, namely, the domain-decomposition (SAP) preconditioned generalized conjugate residual solver (GCR). The left panel of Figure 4 shows the results for the local residual, defined as $r_{\text{loc}}(x_{0},y_{0})=\frac{|A(x_{0},y_{0})S(x_{0},y_{0})-\eta_{t}(x_{0},y_{0})|}{|S(x_{0},y_{0})|}$, for three different solver routines. While the CGNE solver (blue) shows substantial deviations between source $\eta_{t}(x_{0},y_{0})$ and solution $S(x_{0},y_{0})$, this is neither the case for the DP-modified CGNE solver residual (red) that stays almost constant over the whole time extent, nor the DP-modified SAP GCR solver one (black), which stays well below $r_{\text{loc}}=10^{-8}$. The right panel shows a representative result for H200 with $\beta=3.55$ and $\kappa_{l}=\kappa_{s}=m_{\text{sym}}$. If $\kappa_{\text{charm}}$ is chosen close to a value corresponding to $\frac{1}{2}am_{\text{eff}}^{\text{heavy-heavy}}$, one can afford to set the global residual to even larger values, which helps to decrease the computational cost of the inversion further.

Finally, Figure 5 highlights the influence of DP on $\mathrm{D}_{\rm s}$ meson (left panel) and D meson (right panel) effective masses. While the light and strange quark systems were solved in all cases by a SAP-preconditioned GCR solver with local deflation (DFL), the charm quark inversion was performed either in the same way or via the DP-modified SAP GCR solver (red) for comparison. In the $\mathrm{D}_{\rm s}$ meson case, an increase in accuracy is clearly visible, as the plateau region extends farther into the bulk. For the D meson, the benefits of the DP-technique are less pronounced. Further studies based on higher statistics and devoted to finding optimal choices for $r_{\text{gl}}$ and $\alpha_{0}$ are certainly needed, given the promising results of our first numerical tests.

5 Analysis details

Since the tests of the distance preconditioning technique and its final implementation within our correlator measurement setup are still on-going, the results we present below were generated by solving the un-preconditioned linear system (Eq. 11) with a decreased residual (Eq. 9) of $r_{\rm gl}=10^{-15}$ in the charm and $r_{\rm gl}=10^{-11}$ in the strange and light quark propagator computations. With this choice of (global) residual for the charm quark inversion, the numerical problems discussed in the previous section can be seen in charmonium two-point functions at distances of around 60 timeslices from the source.

As mentioned in Section 2, we use twisted-mass reweighting for the light quarks and the rational approximation for the strange quark in the HMC. To correct to the proper distribution, we must reweight the physical observables (in this case, the two-point functions) with the associated reweighting factors:

[TABLE]

The twisted-mass ($W_{0}$) and rational approximation reweighting ($W_{1}$) factors are defined in Ref. [2] (Eqs. (3.2) and (3.5), respectively).

The analysis of the two-point functions proceeds by first making use of time-reversal symmetry and averaging $C_{\rm X}(x_{0},y_{0})$ and $C_{\rm X}(T-x_{0},T-y_{0})$ (${\rm X}={\rm A},{\rm P}$, see Eq. 6) for the three source positions, e.g. $y_{0}=30a,47a$ and $65a$ for $T/a=95$ at $\beta=3.4$. In other words, we average the forward-propagating part of $y_{0}=30a$ with the backward-propagating branch of $y_{0}=65a$, while for correlators with $y_{0}=47a$, we average the forward and backward propagating mesons. To extract the lattice pseudoscalar decay constant, we fit the averaged correlators for $C_{\rm A}$ and $C_{\rm P}$ simultaneously (i.e., four correlators in total) to the single-exponential forms in Eq. 7. Of course, this functional form is only valid in a region, where contamination from excited states has fallen below the noise and boundary effects as well as the numerical problems associated with the charm quark inversion are absent. As discussed in Section 3, the values of $y_{0}$ have been chosen to be sufficiently far away from the boundary. Effects from the opposite boundary and the numerical inversion problems are not significant when fixing the end point of the fit range to $x^{\rm max}_{0}-y_{0}\lesssim 25a$ for the D meson and $x^{\rm max}_{0}-y_{0}\lesssim 30a$ for the $\mathrm{D}_{\mathrm{s}}$ meson. In order to select the starting point of the fit interval, $x_{0}^{\rm min}$, we first estimated the contribution of the excited states by fitting the correlators to a double-exponential form. $x_{0}^{\rm min}$ is then fixed to the point at which excited state contributions amount to less than one quarter of the statistical error of the correlation function. For example, for $C_{\rm P}$,

[TABLE]

where $A^{\prime}$ and $M^{\prime}$ denote the amplitude and the mass of the first excited state, respectively, as determined by the fit, and $\Delta C_{\rm P}$ denotes the statistical error of the correlator. Figure 6 illustrates this procedure. The lattice pseudoscalar decay constant is then extracted from a combination of the amplitudes of $C_{\rm A}$ and $C_{\rm P}$ (resulting from single-exponential fits between $x_{0}^{\rm min}$ and $x_{0}^{\rm max}$), in accordance with Eq. 7.

As the charm quark is partially quenched in our setup, we need to fix the corresponding hopping parameter ($\kappa_{\rm charm}$) in our analysis. In a first step, on a subset of the available statistics, $\kappa_{\rm charm}$ was estimated from the charmonium 1S spin-averaged mass $\overline{M}(1S)=\left(m_{\eta_{c}}+3m_{J/\Psi}\right)/4$. In order to maintain flexibility for possible adjustments later, the correlation functions were calculated with full statistics at two values of the hopping parameter, straddling the target $\kappa_{\rm charm}$; these values were chosen by allowing the lattice spacing to vary by $\pm 2\%$. $\kappa_{\rm charm}$ was then determined on this full set by interpolating $\overline{M}(1S)$ linearly in the inverse of the hopping parameter to the physical point, as given by the central value for $a$. A linear interpolation is valid, since both $\kappa$-values for the “heavier” and “lighter” charm quark are close enough to the physical value. Employing the charmonium 1S mass introduces an unknown (although expected to be small) uncertainty, since we neglect the disconnected diagrams and possible flavour mixing. As an alternative, free from these systematics, we also determined $\kappa_{\rm charm}$ utilising the spin-flavour-averaged 1S meson mass combination $M_{X}=\left(6m_{{\rm D}^{*}}+2m_{\rm D}+3m_{{\rm D}_{\rm s}^{*}}+m_{{\rm D}_{\rm s}}\right)/12$ along the $\overline{m}=\mathrm{const.}$ line, and the spin-averaged ${\rm D}_{\rm s}$ mass along the $\hat{m}_{\rm s}=\mathrm{const.}$ line. So far, the results for the interpolated decay constants, using $\kappa_{\rm charm}$ fixed via either $\overline{M}(1S)$ or $M_{X}$ and $m_{{\rm D}_{\rm s}}$, were consistent within one standard deviation. In the future, we plan to base our final determination of the hopping parameter of the (valence) charm quark on $M_{X}$ and $m_{{\rm D}_{\rm s}}$.

6 Preliminary results

The results presented in this section are preliminary and systematic errors arising from the uncertainty in the lattice spacing, fixing the charm quark mass, the renormalisation procedure, etc. are neglected. The statistical errors are estimated using a bootstrap analysis employing 500 samples. Let us first inspect the ratio of decay constants, $f_{\rm D_{s}}/f_{\rm D}$, for which the renormalisation factor ($Z_{\rm A}$) drops out and also the contributions from the ${\cal O}(a)$ improvement terms partially cancel (see Eqs. (3) and (4)). Figure 7 shows our results for the ratio as a function of $m_{\pi}^{2}$ for $\beta=3.4$ and 3.55 (corresponding to $a=0.0854\,{\rm fm}$ and $0.0644\,{\rm fm}$, respectively), consult Table 1 for the values of $m_{\pi}$ for each ensemble. We performed chiral extrapolations linear in $m_{\pi}^{2}$ along the trajectories of constant $\overline{m}$ and $\hat{m}_{\rm s}=\hat{m}_{\rm s}^{\rm phys}$ simultaneously, enforcing the value at the physical point to be the same in both cases. As the dominant discretisation artefacts are of ${\cal O}(a^{2})$, we show the results for the two $\beta$-values plotted against $a^{2}$ in Figure 8. The comparison with the recent FLAG averages for $N_{\rm f}=2+1$ and $N_{\rm f}=2+1+1$ [10] is promising.

The individual pseudoscalar D and ${\rm D}_{\rm s}$ meson decay constants are displayed in Figures 9 and 10. These depend on the renormalisation, and ${\cal O}(a)$ improvement has a big impact. If we focus on the $\mathrm{D}_{\mathrm{s}}$ meson, e.g., including the $b_{\rm A}m_{\rm qc}$-improvement term in Eq. (4) causes an upward shift of the physical-point value of $f_{\rm D_{s}}$ of about 35% for $\beta=3.4$ and 19% on the finer $\beta=3.55$ lattice. Similar shifts occur for $f_{\rm D}$. Most of the shift is due to the presence of the charm quark mass in $m_{\rm qc}$. In contrast, for the ratio $f_{\rm D_{s}}/f_{\rm D}$ there is only a 1% shift when adding this improvement term, in line with the expectation that the impact of the charm quark largely cancels here. In the future, the chiral extrapolations will be investigated in more detail and, once results at additional lattice spacings have been analysed, a careful continuum limit extrapolation will also be performed.

7 Conclusions and outlook

Our ongoing simulations in the charmed meson sector of $N_{\rm f}=2+1$ lattice QCD address the computation of the leptonic decay constants $f_{\rm D}$ and $f_{\rm D_{s}}$. To reach a precision high enough to be relevant for theory inputs to global analyses in flavour physics phenomenology, among the next necessary steps are: (i) increase of statistics, inclusion of further ensembles and correcting the results for possible small mismatches of the constant physics condition, enabling a stable joint chiral and continuum extrapolation; (ii) a full budget of (statistical and systematic) errors, also accounting for their correlations. This will be carried out in the future.

Acknowledgments. We thank Gunnar Bali, Tomasz Korzec, Stefan Schaefer, Rainer Sommer and Nazario Tantalo for useful discussions. This work is supported by the Deutsche Forschungsgemeinschaft (DFG) through the grants GRK 2149 (Research Training Group “Strong and Weak Interactions – from Hadrons to Dark Matter”, K. E. and J. H.) and SFB/TRR 55 (S. C., S. H. and W. S.). We are indebted to our colleagues in CLS for the joint production of the $N_{\rm f}=2+1$ gauge configurations. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. for granting computer time on SuperMUC at the Leibniz Supercomputing Centre. Additional simulations were performed on the Regensburg iDataCool cluster and on the SFB/TRR 55 QPACE 2 and QPACE B computers [20, 21]. Some of the two-point functions were computed using the Chroma [22] software package, along with the locally deflated domain decomposition solver implementation of openQCD [1], the LibHadronAnalysis library and the multigrid solver implementation of Ref. [23]; our implementation of distance preconditioning (Section 4) is based on [24].