Fermilab Lattice and MILC Collaborations

#  $D$-meson semileptonic decays to pseudoscalars from
four-flavor lattice QCD

Alexei Bazavov  Department of Computational Mathematics, Science and
Engineering, and Department of Physics and Astronomy, Michigan State
University, East Lansing, Michigan 48824, USA  Carleton DeTar  Department of
Physics and Astronomy, University of Utah,
Salt Lake City, Utah 84112, USA  Aida X. El-Khadra  Department of Physics,
University of Illinois, Urbana, Illinois, 61801, USA  Illinois Center for
Advanced Studies of the Universe, University of Illinois,
Urbana, Illinois, 61801, USA  Elvira Gámiz  CAFPE and Departamento de Física
Teórica y del Cosmos,
Universidad de Granada, E-18071 Granada, Spain  Zechariah Gelzer  Department
of Physics, University of Illinois, Urbana, Illinois, 61801, USA  Steven
Gottlieb  Department of Physics, Indiana University, Bloomington, Indiana
47405, USA  William I. Jay<EMAIL_ADDRESS>Center for Theoretical Physics,
Massachusetts Institute of Technology,
Cambridge, MA 02139, USA  Hwancheol Jeong  Department of Physics, Indiana
University, Bloomington, Indiana 47405, USA  Andreas S. Kronfeld  Theory
Division, Fermi National Accelerator Laboratory, Batavia, Illinois, 60510, USA
Ruizi Li  Department of Physics, Indiana University, Bloomington, Indiana
47405, USA  Andrew T. Lytle  Department of Physics, University of Illinois,
Urbana, Illinois, 61801, USA  Illinois Center for Advanced Studies of the
Universe, University of Illinois,
Urbana, Illinois, 61801, USA  Paul B. Mackenzie  Theory Division, Fermi
National Accelerator Laboratory, Batavia, Illinois, 60510, USA  Ethan T. Neil
Department of Physics, University of Colorado, Boulder, Colorado 80309, USA
Thomas Primer  Department of Physics, University of Arizona, Tucson, Arizona
85721, USA  James N. Simone  Theory Division, Fermi National Accelerator
Laboratory, Batavia, Illinois, 60510, USA  Robert L. Sugar  Department of
Physics, University of California, Santa Barbara, California 93106, USA  Doug
Toussaint  Department of Physics, University of Arizona, Tucson, Arizona
85721, USA  Ruth S. Van de Water  Theory Division, Fermi National Accelerator
Laboratory, Batavia, Illinois, 60510, USA  Alejandro Vaquero  Department of
Physics and Astronomy, University of Utah,
Salt Lake City, Utah 84112, USA  Departmento de Física Teórica, Universidad de
Zaragoza,
50009 Zaragoza, Spain

(August 28, 2024)

###### Abstract

We present lattice-QCD calculations of the hadronic form factors for the
semileptonic decays $D\to\pi\ell\nu$, $D\to K\ell\nu$, and $D_{s}\to
K\ell\nu$. Our calculation uses the highly improved staggered quark (HISQ)
action for all valence and sea quarks and includes $N_{f}=2+1+1$ MILC
ensembles with lattice spacings ranging from $a\approx
0.12\leavevmode\nobreak\ \text{fm}$ down to $0.042\leavevmode\nobreak\
\text{fm}$. At most lattice spacings, an ensemble with physical-mass light
quarks is included. The HISQ action allows all the quarks to be treated with
the same relativistic light-quark action, allowing for nonperturbative
renormalization using partial conservation of the vector current. We combine
our results with experimental measurements of the differential decay rates to
determine $|V_{cd}|^{D\to\pi}=0.2238(11)^{\rm Expt}(15)^{\rm QCD}(04)^{\rm
EW}(02)^{\rm SIB}[22]^{\rm QED}$ and $|V_{cs}|^{D\to K}=0.9589(23)^{\rm
Expt}(40)^{\rm QCD}(15)^{\rm EW}(05)^{\rm SIB}[95]^{\rm QED}$. This result for
$|V_{cd}|$ is the most precise to date, with a lattice-QCD error that is, for
the first time for the semileptonic extraction, at the same level as the
experimental error. Using recent measurements from BES III, we also give the
first-ever determination of $|V_{cd}|^{D_{s}\to K}=0.258(15)^{\rm
Expt}(01)^{\rm QCD}[03]^{\rm QED}$ from $D_{s}\to Kl\nu$. Our results also
furnish new Standard Model calculations of the lepton flavor universality
ratios $R_{\mu/e}^{D\to\pi}=0.98671(17)^{\rm QCD}[500]^{\rm QED}$,
$R_{\mu/e}^{D\to K}=0.97606(16)^{\rm QCD}[500]^{\rm QED}$, and
$R_{\mu/e}^{D_{s}\to K}=0.98099(10)^{\rm QCD}[500]^{\rm QED}$, which are
consistent within $2\sigma$ with experimental measurements. Our extractions of
$|V_{cd}|$ and $|V_{cs}|$, when combined with a value for $|V_{cb}|$, provide
the most precise test of second-row CKM unitarity, finding agreement with
unitarity at the level of one standard deviation.

††preprint: MIT-CTP/5513, FERMILAB-PUB-22-943-T

###### Contents

  1. I Introduction
  2. II Definitions
  3. III Simulation details
  4. IV Correlator Analysis
    1. IV.1 Definitions
    2. IV.2 Statistical analysis
    3. IV.3 Nonperturbative renormalization
  5. V Chiral-continuum analysis
    1. V.1 Description of the chiral-continuum fit function
    2. V.2 Chiral-continuum fits
    3. V.3 Alternative constructions of $f_{+}$ and $f_{0}$
    4. V.4 Model-independent $z$ expansion
    5. V.5 Spectator dependence
    6. V.6 Comparison with existing results in the literature
  6. VI Systematic Error Analysis
    1. VI.1 Chiral-continuum fits: stability analysis
    2. VI.2 Chiral-continuum fits: error breakdown
    3. VI.3 Finite-volume corrections
    4. VI.4 Nonequilibrated topological charge
  7. VII Phenomenology
    1. VII.1 Experimental measurements
    2. VII.2 Systematic uncertainty from strong isospin effects and QED
    3. VII.3 CKM matrix elements
    4. VII.4 Tests of CKM unitarity
    5. VII.5 Lepton flavor universality
  8. VIII Conclusions
  9. A Analysis of staggered correlation functions
  10. B Discretization errors for HISQ
  11. C Shrinkage of covariance and correlation matrices
    1. C.1 Linear shrinkage
    2. C.2 Nonlinear shrinkage
    3. C.3 Numerical examples of shrinkage
  12. D Fits: additional details and figures
    1. D.1 Correlator fits
    2. D.2 Chiral-continuum fits: Results for $D\to K$ and $D_{s}\to K$
    3. D.3 Chiral-continuum fits: Stability plots for $D\to K$ and $D_{s}\to K$
    4. D.4 Chiral continuum fits: Error breakdowns for $D\to K$ and $D_{s}\to K$
    5. D.5 $z$-expansion fits: Joint fits to lattice-QCD form factors and experimental data

##  I Introduction

Historically, measurements in quark-flavor physics have a strong precedent of
anticipating the direct discovery of new particles. To name one instance,
consider the charm quark, decays of which are the subject of this paper. Its
existence was conjectured on the basis of symmetry [1, 2], and its mass was
predicted to explain the rates of strangeness-changing neutral-current
processes [2, 3]. The discovery of the $J/\psi$ [4, 5] was then immediately
interpreted as charmonium [6, 7, 8, 9]. Another example is the measurement in
1987 of large mixing in neutral $B$ mesons by the ARGUS Collaboration [10],
which suggested the unusually large mass for the top quark (see, e.g., Ref.
[11]), eight years before its direct observation at the Tevatron in 1995 [12,
13, 14]. In light of several anomalies in measurements of $B$-meson decays and
tension in several tests of the Standard Model (SM) flavor structure [15, 16],
one can speculate that this area of particle physics is again pointing toward
something new. To illuminate the situation, it is timely to improve the
theoretical ingredients in confronting experiment with the Standard Model for
other quark-flavor processes. In this paper, we report on lattice-QCD
calculations relevant to the second row of the Cabibbo-Kobayashi-Maskawa (CKM)
matrix, enabling stringent tests of second-row CKM unitarity.

Within the Standard Model (SM), charged-current flavor-changing processes are
described by the CKM matrix

$V_{\rm CKM}=\begin{pmatrix}V_{ud}&V_{us}&V_{ub}\\\ V_{cd}&V_{cs}&V_{cb}\\\
V_{td}&V_{ts}&V_{tb}\end{pmatrix},$ (1)

which describes the mismatch between the propagating mass eigenstates and the
flavor eigenstates which participate in the weak interaction. By construction,
the CKM matrix is unitary, so each row and column should have unit norm.
Deviations from this expectation can arise if $V_{\rm CKM}$ is a $3\times 3$
submatrix in an extended flavor sector or if non-SM processes contribute to
measured decay and mixing rates. It is important to test the CKM paradigm
using independent determinations from multiple processes, for example,
comparing leptonic and semileptonic decays with the same flavor charge.
Improved precision for the individual matrix elements leads directly to more
stringent tests of the Standard Model. Any statistically significant deviation
from the predictions of CKM-unitarity would constitute evidence for new
physics beyond the Standard Model.

The strongest test of unitarity comes from the first row, where the matrix
elements are determined most precisely, with the exception of $|V_{ub}|$,
which plays a negligible role in the first row unitarity relation at the
current level of precision. Either taking the most precise value of $|V_{ud}|$
that comes from superallowed $\beta$ decays [17] 111Recent calculations of the
universal electroweak radiative corrections relevant for superallowed $\beta$
decays in Refs. [18, 19, 20, 21, 22] found larger values than those estimated
before, shifting the central value of $|V_{ud}|$ and increasing the tension
with unitarity. In addition, further, previously unaccounted, nuclear-
structure uncertainties in the inner radiative correction have considerably
increased the error for earlier determinations [19, 23]. and $|V_{us}|$ as
extracted from semileptonic $K_{\ell 3}\equiv K\to\pi\ell\nu$ decays, or using
only inputs from kaon and pion decays (i.e., $|V_{us}|$ from semileptonic
decays and $|V_{us}|/|V_{ud}|$ from the ratio of leptonic decays, $K_{\ell
2}\equiv K\to\ell\nu$ over $\pi_{\ell 2}\equiv\pi\to\ell\nu$ [24]), the
combination $|V_{ud}|^{2}+|V_{us}|^{2}+|V_{ub}|^{2}$ is in tension with
unitarity at the $3\sigma$ level [25]. There is also a $\sim 3\sigma$ tension
between the semileptonic and the leptonic determinations of $|V_{us}|$ [25],
where the leptonic determination uses $|V_{ud}|$ from superallowed decays as
an external input. In those tests, the relevant QCD nonperturbative inputs for
semileptonic and leptonic decays, the form factor $f_{+}^{K\pi}(0)$ [26, 27,
28, 29] and the ratio of decay constants $f_{K}/f_{\pi}$ [30, 31, 32, 33, 33,
34, 35], respectively, are calculated using lattice QCD with uncertainties
that have reached the $\sim 0.18$% level [36]. Experimental data for the decay
widths of $K_{\ell 3}$ and $K_{\ell 2}/\pi_{\ell 2}$ decays are similarly
precise [37, 38], leaving electromagnetic corrections as an important source
of uncertainty in the extraction of the corresponding CKM matrix elements.
Pioneering work addressing the calculation of structure-dependent QED
corrections both for pion and kaon leptonic decays using lattice techniques
[39, 40] and kaon semileptonic decays [41, 42, 43] including lattice
calculations of the $\gamma W$-box contribution, have been recently performed,
opening the door to an important reduction of the electromagnetic uncertainty.

Similarly precise tests for the CKM matrix elements in the second row have
been limited both by theory and experimental uncertainties. On the theory
side, the error for the decay constants $f_{D}$ and $f_{D_{s}}$ (roughly
$0.35$–$0.2\%$ [36]) are now subleading in the extraction of $|V_{cd}|$ and
$|V_{cs}|$, respectively, from leptonic decays thanks to the progress made by
lattice calculations in the last years [32, 31]. However, the situation is
very different for semileptonic extractions of those CKM matrix elements.
Since the decay rates are not suppressed by the lepton mass, experimental
measurements are more precise. For leptonic decays, the HFLAV world averages
for $f_{D_{s}}|V_{cs}|$ and $f_{D}|V_{cd}|$ have fractional errors of roughly
$1\%$ and $2\%$, respectively [44]. The corresponding semileptonic decay-rate
measurements are roughly a factor of two more precise in each case, with the
fractional errors in $f_{+}^{D\to K}(0)|V_{cs}|$ and
$f_{+}^{D\to\pi}(0)|V_{cd}|$ around $0.5\%$ and $1\%$, respectively [44].
Lattice-QCD calculations of semileptonic form factors (including both
normalization and shape), while more complex than for decay constants for
leptonic decays, have a long history in lattice QCD [45, 46, 47, 48, 49, 50,
51, 52, 53, 54, 55, 56, 57, 58]. Now, however, the current experimental errors
and the forthcoming improvements by BES III motivate further reducing the
lattice-QCD errors to the level of experimental precision.

In this work, we leverage the same theoretical tools that were successfully
employed in the calculation of decay constants and the $K_{\ell 3}$ form
factor [59, 32, 27]: the same highly improved relativistic lattice actions and
gauge-field ensembles with physical quark masses and small lattice spacings.
In particular, we compute the hadronic form factors for the semileptonic
decays $D\to\pi\ell\nu$, $D\to K\ell\nu$, and $D_{s}\to K\ell\nu$ in lattice
QCD, with the goal of obtaining percent-level determinations of $|V_{cd}|$ and
$|V_{cs}|$ when combined with experimental data. Our values for $|V_{cd}|$ and
$|V_{cs}|$ provide a stringent test of unitarity and their precision allows a
commensurate comparison with leptonic determinations. As a key aspect of our
analysis, we report the correlations between the hadronic form factors in the
different decay channels as well as between the final values for $|V_{cd}|$
and $|V_{cs}|$ (see Sec. VII). Preliminary results for the present calculation
of the form factors have been presented in Refs. [60, 61]. We note that the
HPQCD collaboration has recently presented a precise lattice-QCD calculation
of the form factors for $D\to K$ decay [62, 63], with a quoted lattice-QCD
uncertainty close to the experimental one in the extraction of $|V_{cs}|$
[62]. On the other hand, this paper yields the first percent-level
determination of $|V_{cd}|$ and enables the first stringent test of second-row
CKM unitarity from semileptonic $D$-meson decays.

With the hadronic form factors for a given decay in hand, it is
straightforward to construct the lepton flavor universality (LFU) ratios
$R_{\mu/e}$, which are defined as the ratio of the branching fractions to muon
versus electron final states; see Sec. VII.5. These ratios are expected to be
close but not identically equal to unity in the SM, with differences coming
from lepton-mass, isospin-breaking, and QED effects. Lattice QCD calculations
offer a theoretically clean method for determining the SM prediction to high
precision (up to QED corrections), contributing to stringent LFU tests in
those channels.

The rest of this article is organized as follows. Section II reviews the
definitions and formalism for relating experimentally measured decay rates to
the hadronic form factors we calculate. Section III gives details related to
the lattice-QCD simulation. Section IV reports the statistical analysis of
Euclidean correlation functions which yields renormalized form factors.
Section V describes the final chiral-continuum fit, which interpolates the
form factors to the physical hadron masses and extrapolates to the continuum
limit. Section VI analyzes the uncertainties in our calculation and summarizes
the complete statistical and systematic error budget for the form factors.
Section VII discusses applications to phenomenology, including determinations
of the CKM matrix elements and the LFU ratios in each channel. Finally, Sec.
VIII gives some concluding remarks. Four appendices provide additional
technical information. Appendix A contains useful formulae appearing in the
statistical analysis of staggered correlation functions. Appendix B presents
useful information about staggered fermions and heavy quark effective theory
when the bare lattice quark mass is large. Appendix C describes linear and
nonlinear shrinkage techniques for correlation and covariance matrix, the
latter of which is a novel aspect of the correlator analysis presented in this
work. Appendix D provides supporting details and figures regarding various
fits, which exceed the scope of the main text but illustrate the robustness of
our analysis.

##  II Definitions

The differential decay rate for the semileptonic decay $H\to L\ell\nu$ of a
heavy pseudoscalar meson $H\in\\{D,D_{s}\\}$ to a light pseudoscalar meson
$L\in\\{K,\pi\\}$ is given by

$\displaystyle\begin{split}\frac{d\Gamma}{dq^{2}}=&\frac{G_{F}^{2}}{24\pi^{3}}\eta_{\rm
EW}^{2}|V_{cx}|^{2}(1-\epsilon)^{2}(1+\delta_{\rm EM})\times\\\
&\left[\left|\bm{p}\right|^{3}\left(1+\frac{\epsilon}{2}\right)\left|f_{+}(q^{2})\right|^{2}+\left|\bm{p}\right|M_{H}^{2}\left(1-\frac{M_{L}^{2}}{M_{H}^{2}}\right)^{2}\frac{3\epsilon}{8}\left|f_{0}(q^{2})\right|^{2}\right],\end{split}$
(2)

where $\epsilon=m_{\ell}^{2}/q^{2}$ (with $m_{\ell}$ the lepton mass),222 In
our notation, $m_{\ell}$ with a cursive subscript always refers to the lepton
mass in the decay $H\to L\ell\nu$. The light-quark mass is denoted $m_{l}$.
$q$ is the momentum transfer, $M_{H}$ and $M_{L}$ are the masses of the heavy
initial and light final mesons, and $\bm{p}$ is the three-momentum of the
final-state meson in the rest frame of the initial hadron. Short-distance
electroweak corrections to $G_{F}$ are contained in $\eta_{\rm
EW}=1+(\alpha_{\text{QED}}/\pi)\ln(M_{Z}/\mu)|_{\mu=M_{D}}=1.009(2)$ [64],
where the error is an estimate of the scale uncertainty from a factor-of-two
variation around $\mu=M_{D}$.333Physically, the scale dependence of $\eta_{\rm
EW}$ should cancel against that of the structure-dependent electromagnetic
corrections which, though calculable in principle, have never been computed
for these decays. Computing these corrections exceeds the scope of the present
work.  Long-distance and structure-dependent electromagnetic corrections are
described by $\delta_{\rm EM}$.444Systematic uncertainties from neglected
electromagnetic corrections and strong isospin breaking are discussed in Sec.
VII.2.

The form factors $f_{+}(q^{2})$ and $f_{0}(q^{2})$ encapsulate the
nonperturbative hadronic structure of the decay. They arise in the usual way
from the Lorentz-covariant decomposition of the relevant transition matrix
elements,

$\displaystyle\langle L|\mathcal{V}^{\mu}|H\rangle$
$\displaystyle\equiv\sqrt{2M_{H}}\left[v^{\mu}f_{\parallel}(q^{2})+p_{\perp}^{\mu}f_{\perp}(q^{2})\right],$
(3) $\displaystyle\equiv
f_{+}(q^{2})\left(k^{\mu}+p^{\mu}-\frac{M_{H}^{2}-M_{L}^{2}}{q^{2}}q^{\mu}\right)+f_{0}(q^{2})\frac{M_{H}^{2}-M_{L}^{2}}{q^{2}}q^{\mu},$
(4) $\displaystyle\langle L|\mathcal{S}|H\rangle$
$\displaystyle=\frac{M_{L}^{2}-M_{H}^{2}}{m_{h}-m_{x}}f_{0}(q^{2}).$ (5)

In these expressions, $k^{\mu}$, and $p^{\mu}$ refer to the four-momentum of
the heavy initial and light final mesons; $m_{h}$ and $m_{x}$ refer to the
masses of the heavy and light quarks in the current; $v^{\mu}=k^{\mu}/M_{H}$
is the four-velocity of the heavy meson; $p_{\perp}^{\mu}=p^{\mu}-(p\cdot
v)v^{\mu}$ is the component of the final-state hadron’s momentum orthogonal to
$v$; and $q^{\mu}=k^{\mu}-p^{\mu}$ is the momentum transfer. The same form
factor $f_{0}$ appears in Eqs. (4) owing to partial conservation of the vector
current (PCVC), namely
$\partial_{\mu}\mathcal{V}^{\mu}=(m_{h}-m_{x})\mathcal{S}$ as an operator
identity.

In lattice gauge theory, we introduce bilinears of lattice fermion
fields—$J=V^{0}$, $V^{i}$, and $S$—and associated matching factors $Z_{J}$,
such that $Z_{J}J$ and the corresponding $\mathcal{J}$ have the same matrix
elements (up to controlled uncertainties). In this notation, and in the rest
frame of the decaying meson, the relations between form factors and matrix
elements take the following forms:

$\displaystyle f_{\parallel}(q^{2})$
$\displaystyle=Z_{V^{0}}\frac{\left<L\vphantom{V^{0}H}\right|V^{0}\left|H\vphantom{LV^{0}}\right>}{\sqrt{2M_{H}}},$
(6) $\displaystyle f_{\perp}(q^{2})$
$\displaystyle=Z_{V^{i}}\frac{1}{p^{i}}\frac{\left<L\vphantom{V^{i}H}\right|V^{i}\left|H\vphantom{LV^{i}}\right>}{\sqrt{2M_{H}}},$
(7) $\displaystyle f_{0}(q^{2})$
$\displaystyle=Z_{m}Z_{S}\frac{m_{h}-m_{l}}{M_{H}^{2}-M_{L}^{2}}\left<L\vphantom{SH}\right|S\left|H\vphantom{LS}\right>.$
(8)

[No sum is implied in Eq. 7.] Using the preceding equations, the vector form
factor is given by a linear combination of $f_{\perp}(q^{2})$ and
$f_{0}(q^{2})$,

$\displaystyle f_{+}(q^{2})$
$\displaystyle=\left(\frac{M_{H}-E_{L}}{\sqrt{2M_{H}}}\right)\left(1-\frac{E_{L}^{2}-M_{L}^{2}}{(M_{H}-E_{L})^{2}}\right)f_{\perp}(q^{2})+\left(\frac{M_{H}^{2}-M_{L}^{2}}{M_{H}-E_{L}}\right)\frac{f_{0}(q^{2})}{2M_{H}},$
(9)

which will be useful below.

In momentum space, PCVC implies the following condition for the lattice
currents:555In Minkowski space, the basic momentum-space operator relation
reads
$iq_{\mu}\left\langle\mathcal{V}^{\mu}(q)\right\rangle=(m_{h}-m_{l})\left\langle\mathcal{S}\right\rangle$.
Equation 10, in which all terms come with a positive sign and without factors
of $i$, amounts to a definition of the sign convention for Wick rotation and
the phase convention for the lattice currents.

$\displaystyle
Z_{V^{0}}(M_{H}-E_{L})\left<L\vphantom{V^{0}H}\right|V^{0}\left|H\vphantom{LV^{0}}\right>+Z_{V^{i}}\bm{q}\cdot\left<L\vphantom{\bm{V}H}\right|\bm{V}\left|H\vphantom{L\bm{V}}\right>=Z_{m}Z_{S}(m_{h}-m_{x})\left<L\vphantom{SH}\right|S\left|H\vphantom{LS}\right>.$
(10)

which can be used to extract the renormalization factors for the temporal and
spatial components of the vector current, $Z_{V^{0}}$ and $Z_{V^{i}}$ [48], as
explained in detail in Sec. IV.3. With the present treatment of all valence
quarks with the highly improved staggered quark (HISQ) action [65], the local
scalar density enjoys absolute normalization, $Z_{m}Z_{S}=1$ [66, 67].
Furthermore, PCVC allows one to express any single matrix element in terms of
the other two involved in the relation in Eq. 10, for example,

$\displaystyle f_{+}^{\rm alt}(q^{2})$
$\displaystyle=\frac{1}{\sqrt{2M_{H}}}\left[f_{\parallel}(q^{2})+(M_{H}-E_{L})f_{\perp}(q^{2})\right],$
(11) $\displaystyle f_{0}^{\rm alt}(q^{2})$
$\displaystyle=\frac{\sqrt{2M_{H}}}{M_{H}^{2}-M_{L}^{2}}\left[(M_{H}-E_{L})f_{\parallel}(q^{2})+(E_{L}^{2}-M_{L}^{2})f_{\perp}(q^{2})\right],$
(12)

with $f_{\parallel}$ and $f_{\perp}$ computed using Eqs. 6 and 7.666 Another
expression for $f_{+}$ in terms of $f_{0}$ and $f_{\parallel}$ exists but
involves a delicate numerical cancellation near $q^{2}_{\rm max}$. For this
reason it is excluded from the subsequent discussion. These alternative
constructions will be used to check for systematic errors in our analysis; see
Sec. V.3.

##  III Simulation details

Table 1:  A summary of the lattice spacings, lattice spatial ($N_{s}$) and temporal ($N_{t}$) sizes, valence quark masses, intermediate scale-setting parameters, number of source times and configurations $N_{\rm src}\times N_{\rm configs}$, source-sink separations $T/a$, and approximate Goldstone (pseudoscalar taste) pion masses used in our calculation. The text describes the ensembles’ sea and valence masses in more detail. The gauge ensembles were generated by the MILC collaboration [68, 69, 32]. The values for the gradient-flow scale $w_{0}/a$ have been calculated previously [70, 71]. The simulation program is described in detail in Ref. [69] and was later extended to smaller lattice spacings ($a\approx 0.042\leavevmode\nobreak\ \text{fm}$), as described in Ref. [32]. The number of source times $N_{\rm src}$ refers to the number of loose-solve source times employed in the truncated solver method; on each configuration one corresponding fine solve is used. Values for the sea- and valence-quark masses are given in Table 2.  $\approx a$ | $N_{s}^{3}\times N_{t}$ | $m_{l}$ | $m_{h}/m_{c}$ | $w_{0}/a$ | $N_{\rm src}\times N_{\rm configs}$ | $T/a$ | $\approx M_{\pi,P}$
---|---|---|---|---|---|---|---
[fm] |  |  |  |  |  |  | [MeV]
0.12 | $48^{3}\times 64$ | physical | $0.9,1.0,1.4$ | 1.4168(10) | $32\times 1352$ | {12, 13, 14, 16, 17} | 135
0.088 | $64^{3}\times 96$ | physical | $0.9,1.0,1.5,2.0$ | 1.9470(13) | $24\times 980$ | {16, 17, 19, 22, 25} | 130
0.088 | $48^{3}\times 96$ | $m_{s}/10$ | $0.9,1.0,1.5,2.0$ | 1.9299(12) | $24\times 697$ | {16, 19, 22, 25} | 224
0.057 | $96^{3}\times 192$ | physical | $0.9,1.0,1.1,2.2$ | 3.0119(19) | $32\times 877$ | {25, 28, 30, 34, 37} | 134
0.057 | $64^{3}\times 144$ | $m_{s}/10$ | $0.9,1.0,2.0$ | 2.9478(31) | $36\times 916$ | {23, 30, 34, 37} | 231
0.057 | $48^{3}\times 144$ | $m_{s}/5$ | $0.9,1.0,2.0$ | 2.8956(33) | $36\times 823$ | {23, 30, 34, 37} | 325
0.042 | $64^{3}\times 192$ | $m_{s}/5$ | $0.9,1.0,2.0$ | 3.9222(29) | $24\times 1008$ | {34, 39, 45, 52} | 308
Table 2: Sea- and valence-quark masses in lattice units for the ensembles used in this calculation. The first two columns specify the ensemble by the approximate lattice spacing and the ratio of light- and strange-quark masses. The next three columns give the sea-quark masses. The final three columns contain the valence-quark masses. $\approx a$ | $m_{l}/m_{s}$ | $(am_{l})^{\rm sea}$ | $(am_{s})^{\rm sea}$ | $(am_{c})^{\rm sea}$ | $(am_{l})^{\rm valence}$ | $(am_{s})^{\rm valence}$ | $(am_{h})^{\rm valence}$
---|---|---|---|---|---|---|---
0.12 | 1/27 | 0.001907 | 0.05252 | 0.6382 | $(am_{l})^{\rm sea}$ | $(am_{s})^{\rm sea}$ | {0.5744, 0.6382, 0.8935}
0.088 | 1/27 | 0.0012 | 0.0363 | 0.432 | $(am_{l})^{\rm sea}$ | $(am_{s})^{\rm sea}$ | {0.389, 0.432, 0.648, 0.864}
0.088 | 1/10 | 0.00363 | 0.0363 | 0.43 | $(am_{l})^{\rm sea}$ | $(am_{s})^{\rm sea}$ | {0.389, 0.432, 0.648, 0.864}
0.057 | 1/27 | 0.0008 | 0.022 | 0.26 | $(am_{l})^{\rm sea}$ | $(am_{s})^{\rm sea}$ | {0.257, 0.286, 0.572}
0.057 | 1/10 | 0.0024 | 0.024 | 0.286 | $(am_{l})^{\rm sea}$ | $(am_{s})^{\rm sea}$ | {0.257, 0.286, 0.572}
0.057 | 1/5 | 0.0048 | 0.024 | 0.286 | $(am_{l})^{\rm sea}$ | $(am_{s})^{\rm sea}$ | {0.257, 0.286, 0.572}
0.042 | 1/5 | 0.00316 | 0.0158 | 0.188 | 0.00311 | 0.01555 | {0.164, 0.1827, 0.365}
Figure 1:  Summary of the lattice spacings and light-quark valence masses used
in the present calculation. The sizes of the colored circles are proportional
to the number of configurations in each ensemble. Quantitative details are
given in Table 1.

Our calculation uses ensembles generated by the MILC Collaboration using a
one-loop Symanzik improved gauge action and $N_{f}=2+1+1$ flavors of dynamical
sea quarks with the HISQ action [68, 69, 32]. Table 1 and Fig. 1 summarize the
ensembles used in this work. Lattice spacings range from $a\approx 0.12$ fm
down to $a\approx 0.042$ fm. An ensemble with physical-mass light quarks
appears for all lattice spacings but $a\approx 0.042\leavevmode\nobreak\
\text{fm}$. For the finer lattice spacings, we also include ensembles with
heavier-than-physical light quarks with $m_{l}\approx m_{s}/10$ and
$m_{l}\approx m_{s}/5$.

The masses of the valence light and strange quarks generally match those in
the sea. In all ensembles the charm and strange quarks in the sea have (close
to) physical masses. The heavy valence quarks used in this study range from
around nine-tenths to around twice the physical charm mass. The precise values
for the sea- and valence-quark masses are given in Table 2.

Although the primary targets of this work are the dimensionless form factors
$f_{0}$ and $f_{+}$, many intermediate quantities (e.g., $f_{\parallel}$ and
$f_{\perp}$ and masses) are dimensionful. Throughout this work, the scale is
set on each ensemble using previously calculated values for the gradient-flow
scale $w_{0}/a$ [70, 71], also listed in Table 1. Details of the intermediate
scale-setting scheme in the chiral-continuum analysis are discussed below in
Sec. V.2.

##  IV Correlator Analysis

###  IV.1 Definitions

Figure 2: Schematic picture of the three-point functions defined in Eqs. 16,
18 and 17. The final-state hadron $L\in\\{\pi,K\\}$ is created with momentum
$\bm{p}$ at the time $t_{\rm src}$. An external current $J$ is inserted at
time $t+t_{\rm src}$. The initial-state hadron $H\in\\{D,D_{s}\\}$ is
destroyed at rest at time $T+t_{\rm src}$.

To access the matrix elements
$\left<L\vphantom{SH}\right|S\left|H\vphantom{LS}\right>$,
$\left<L\vphantom{V^{0}H}\right|V^{0}\left|H\vphantom{LV^{0}}\right>$, and
$\left<L\vphantom{V^{i}H}\right|V^{i}\left|H\vphantom{LV^{i}}\right>$, we
compute the following two- and three-point correlation functions:

$\displaystyle C_{H}^{P}(t)$ $\displaystyle=\sum_{\bm{x},\bm{y}}\left\langle
P_{H}(t_{\rm src},\bm{x})P_{H}(t+t_{\rm src},\bm{y})\right\rangle,$ (13)
$\displaystyle C_{H}^{A^{0}}(t)$
$\displaystyle=\sum_{\bm{x},\bm{y}}\left\langle A^{0}_{H}(t_{\rm
src},\bm{x})A^{0}_{H}(t+t_{\rm src},\bm{y})\right\rangle,$ (14) $\displaystyle
C_{L}^{P}(t,\bm{p})$
$\displaystyle=\sum_{\bm{x},\bm{y}}e^{i\bm{p}\cdot(\bm{x}-\bm{y})}\left\langle
P_{L}(t_{\rm src},\bm{x})P_{L}(t+t_{\rm src},\bm{y})\right\rangle,$ (15)
$\displaystyle C_{H\to L}^{S}(t,T,\bm{p})$
$\displaystyle=\sum_{\bm{x},\bm{y},\bm{z}}e^{i\bm{p}\cdot(\bm{x}-\bm{y})}\left\langle
P_{L}(t_{\rm src},\bm{x})S(t+t_{\rm src},\bm{y})P_{H}(T+t_{\rm
src},\bm{z})\right\rangle,$ (16) $\displaystyle C_{H\to
L}^{\bm{V}}(t,T,\bm{p})$
$\displaystyle=\sum_{\bm{x},\bm{y},\bm{z}}e^{i\bm{p}\cdot(\bm{x}-\bm{y})}\left\langle
P_{L}(t_{\rm src},\bm{x})\bm{V}(t+t_{\rm src},\bm{y})P_{H}(T+t_{\rm
src},\bm{z})\right\rangle,$ (17) $\displaystyle C_{H\to
L}^{V^{0}}(t,T,\bm{p})$
$\displaystyle=\sum_{\bm{x},\bm{y},\bm{z}}e^{i\bm{p}\cdot(\bm{x}-\bm{y})}\left\langle
P_{L}(t_{\rm src},\bm{x})V^{0}(t+t_{\rm src},\bm{y})A^{0}_{H}(T+t_{\rm
src},\bm{z})\right\rangle,$ (18)

where the labels denote the initial heavy hadron $H\in\\{D,D_{s}\\}$ and the
final light hadron $L\in\\{\pi,K\\}$. The schematic structure of the three-
point correlators in Eqs. 16, 17 and 18 is depicted in Fig. 2 and the spin-
taste structure of the operators in our simulations specified in Table 3. The
operators used for the scalar current $S$ and temporal vector current $V^{0}$
are local, but the spatial vector current $\bm{V}$ is the taste-singlet one-
link operator. The tastes of the meson creation and annihilation operators are
chosen so that the correlation functions are overall taste singlets. For
three-point functions involving $S$ and $\bm{V}$, Eqs. 16 and 17, it therefore
suffices to use local pseudoscalar operators $P$, corresponding to Goldstone
pseudoscalar mesons, at the source and sink. For three-point functions
involving $V^{0}$, Eq. 18, we use the local axial vector operator $A^{0}$,
corresponding to a non-Goldstone pseudoscalar meson, at either the source or
the sink. Our choice in Eq. 18 is to use $A^{0}$ for the initial-state hadrons
($D$ and $D_{s}$) and $P$ for the final-state hadrons ($\pi$ and $K$). To
reduce statistical noise, APE smearing [72] is applied to the gauge field
appearing in the one-link vector current, with 20 iterations and staple weight
$0.05$.

Table 3: The spin-taste structure of the staggered operators used in this work. Pseudoscalar mesons are created and annihilated using $P$ ($\pi$, $K$, $D$, and $D_{s}$) and $A^{0}$ ($D$ and $D_{s}$). Transitions between these states are induced by the currents $S$, $V^{0}$, and $V^{i}$. The operator $A^{0}$ is necessary to conserve taste in Eq. 18. Operator | Spin $\otimes$ Taste | Locality
---|---|---
$P$ | $\gamma^{5}\otimes\xi_{5}$ | Local
$A^{0}$ | $\gamma^{0}\gamma^{5}\otimes\xi_{0}\xi_{5}$ | Local
$S$ | $1\otimes 1$ | Local
$V^{0}$ | $\gamma^{0}\otimes\xi_{0}$ | Local
$V^{i}$ | $\gamma^{i}\otimes 1$ | One-link

In Eqs. 13, 14, 15, 16, 18 and 17, we work in the rest frame of the heavy
initial hadron $H$ and compute the recoiling light hadron $L$ with eight
different lattice momenta $\bm{p}=2\pi\bm{n}/N_{s}a$, where $N_{s}$ is the
spatial extent of the lattice, and $\bm{n}$ is $(0,0,0)$, $(1,0,0)$,
$(1,1,0)$, $(2,0,0)$, $(2,1,0)$, $(3,0,0)$, $(2,2,2)$, or $(4,0,0)$. For each
choice of heavy-quark mass in Table 1 and momentum above, we compute the
three-point function for several different source-sink separations $T$, given
in Table 1. The final-state light-quark and spectator-quark propagators are
computed using random corner-wall sources [73]. The heavy-quark propagators
are computed sequentially from the end of the spectator-quark propagator at
time $T+t_{\rm src}$ as shown in Fig. 2. For the light- and strange-quark
propagators, we employ the truncated solver method [74, 75], using a single
fine solve together with 24 to 36 loose solves on each configuration (see
Table 1 for details). To reduce autocorrelation in Monte Carlo time, the
source locations for the fine and loose solves are precessed in Euclidean time
from one configuration to the next.

As usual, states with both positive and negative parities contribute to the
staggered correlation functions. For the operators considered here, the
negative-parity states decay smoothly with Euclidean time, while the positive-
parity states oscillate while decaying in Euclidean time. The spectral
decompositions of Eqs. 13, 14, 15, 16, 18 and 17 take the following forms:

$\displaystyle C_{L}^{\mathcal{O}}(t,\bm{p})$
$\displaystyle=\sum_{n=0}(-1)^{n(t+1)}\frac{\left|\left<\emptyset\vphantom{\mathcal{O}_{L}n}\right|\mathcal{O}_{L}\left|n\vphantom{\emptyset\mathcal{O}_{L}}\right>\right|^{2}}{2E_{L}^{(n)}(\bm{p})}\left(e^{-E_{L}^{(n)}(\bm{p})t}+e^{-E_{L}^{(n)}(\bm{p})(N_{t}-t)}\right),$
(19) $\displaystyle C_{H}^{\mathcal{O}}(t)$
$\displaystyle=\sum_{m=0}(-1)^{m(t+1)}\frac{\left|\left<\emptyset\vphantom{\mathcal{O}_{H}m}\right|\mathcal{O}_{H}\left|m\vphantom{\emptyset\mathcal{O}_{H}}\right>\right|^{2}}{2M_{H}^{(m)}}\left(e^{-M_{H}^{(m)}t}+e^{-M_{H}^{(m)}(N_{t}-t)}\right),$
(20) $\displaystyle\begin{split}C_{H\to
L}^{J}(t,T,\bm{p})&=\sum_{m,n}(-1)^{n(t+1)}(-1)^{m(T-t-1)}\frac{\left<\emptyset\vphantom{\mathcal{O}_{L}n}\right|\mathcal{O}_{L}\left|n\vphantom{\emptyset\mathcal{O}_{L}}\right>\left<n\vphantom{Jm}\right|J\left|m\vphantom{nJ}\right>\left<m\vphantom{\mathcal{O}_{H}\emptyset}\right|\mathcal{O}_{H}\left|\emptyset\vphantom{m\mathcal{O}_{H}}\right>}{4E_{L}^{(n)}(\bm{p})M_{H}^{(m)}}\\\
&\phantom{000}\times\left(e^{-E_{L}^{(n)}(\bm{p})t}+e^{-E_{L}^{(n)}(\bm{p})(N_{t}-t)}\right)\left(e^{-M_{H}^{(m)}(T-t)}+e^{-M_{H}^{(m)}(N_{t}-T+t)}\right),\end{split}$
(21)

where $\mathcal{O}\in\\{P,A^{0}\\}$ is the appropriate interpolating operator
and $|\emptyset\rangle$ denotes the QCD vacuum state. In the final line, the
ground-state term contains the transition matrix elements,
$\left.\left<n\vphantom{Jm}\right|J\left|m\vphantom{nJ}\right>\right|_{n=m=0}\equiv\left<L\vphantom{JH}\right|J\left|H\vphantom{LJ}\right>$,
from which one can extract the desired form factors via Eqs. 6, 7 and 8.

For the sake of visualization, certain ratios of correlation functions prove
useful:

$\displaystyle R_{\parallel}(t,T,\bm{p})$ $\displaystyle=\frac{\bar{C}_{H\to
L}^{V^{0}}(t,T,\bm{p})\sqrt{2E_{L}}}{\sqrt{\bar{C}_{L}^{A^{0}}(t,\bm{p})\bar{C}_{H}^{P}(T-t)e^{-E_{L}t}e^{-M_{H}(T-t)}}},$
(22) $\displaystyle R_{\perp}(t,T,\bm{p})$
$\displaystyle=\frac{\sqrt{2E_{L}}}{p^{i}}\frac{\bar{C}_{H\to
L}^{V^{i}}(t,T,\bm{p})}{\sqrt{\bar{C}_{L}^{P}(t,\bm{p})\bar{C}_{H}^{P}(T-t)e^{-E_{L}t}e^{-M_{H}(T-t)}}},$
(23) $\displaystyle R_{0}(t,T,\bm{p})$
$\displaystyle=2\sqrt{M_{H}E_{L}}\left(\frac{m_{h}-m_{l}}{M_{H}^{2}-M_{L}^{2}}\right)\frac{\bar{C}_{H\to
L}^{S}(t,T,\bm{p})}{\sqrt{\bar{C}_{L}^{P}(t,\bm{p})\bar{C}_{H}^{P}(T-t)e^{-E_{L}t}e^{-M_{H}(T-t)}}},$
(24)

where the bars (e.g., $\bar{C}^{P}_{L}$) denote the time-slice-averaged
correlators defined in Eqs. 83 and 85. Up to discretization effects (and
renormalization), these ratios asymptotically approach the form factors at
large Euclidean times:

$\displaystyle R_{\parallel}(t,T,\bm{p})$
$\displaystyle\stackrel{{\scriptstyle 0\ll t\ll
T}}{{\longrightarrow}}Z_{V^{0}}^{-1}f_{\parallel}(\bm{p}),$ (25)
$\displaystyle R_{\perp}(t,T,\bm{p})$ $\displaystyle\stackrel{{\scriptstyle
0\ll t\ll T}}{{\longrightarrow}}Z_{V^{i}}^{-1}f_{\perp}(\bm{p}),$ (26)
$\displaystyle R_{0}(t,T,\bm{p})$ $\displaystyle\stackrel{{\scriptstyle 0\ll
t\ll T}}{{\longrightarrow}}f_{0}(\bm{p}).$ (27)

The subsequent analysis of statistical and systematic uncertainties was
conducted in a blinded fashion. More precisely, all of our three-point
correlation functions were multiplied by a blinding factor $X\in(0.95,1.05)$,
which was chosen randomly and held fixed across all ensembles, momenta,
currents, and heavy-quark masses in the three-point functions. The blinded
results were carried all the way through the phenomenological applications
described in Sec. VII. The blinding factor was removed only after the estimate
of systematic errors was complete and the analysis was frozen.

###  IV.2 Statistical analysis

Figure 3: Pion two-point correlation functions
$C_{\pi}^{P}(t,\bm{p}=2\pi\bm{n}/N_{s}a)$ and effective masses on the
physical-mass $a\approx 0.088\leavevmode\nobreak\ \text{fm}$ ensemble. To
reduce the visual impact of opposite parity states, the effective mass is
computed separately for even and odd times using Eq. 30 and plotted using
circles and triangles, respectively. After folding the data around the
midpoint of the lattice, the correlator is defined for
$t/a\in[0,N_{t}/2]=[0,48]$. Because of the form of Eq. 30 involves offsets by
2, the effective mass is defined on times $t/a\in[2,N_{t}/2-2]=[2,46]$.
Figure 4: Kaon two-point correlation functions
$C_{K}^{P}(t,\bm{p}=2\pi\bm{n}/N_{s}a)$ and effective masses on the physical-
mass $a\approx 0.088\leavevmode\nobreak\ \text{fm}$ ensemble. To reduce the
visual impact of opposite parity states, the effective mass is computed
separately for even and odd times using Eq. 30 and plotted using circles and
triangles, respectively. After the data is folded around the midpoint of the
lattice, the correlator is defined for $t/a\in[0,N_{t}/2]=[0,48]$. Because of
the form of Eq. 30 involves offsets by 2, the effective mass is defined on
times $t/a\in[2,N_{t}/2-2]=[2,46]$. Figure 5: $D$-meson two-point correlation
functions $C_{D}^{P/A^{0}}(t,\bm{p})$ at $\bm{p=0}$ and effective masses with
the heavy-quark mass near its physical value on the physical-mass $a\approx
0.088\leavevmode\nobreak\ \text{fm}$ ensemble. The corresponding results for
the $D_{s}$ meson are very similar. To reduce the visual impact of opposite
parity states, the effective mass is computed separately for even and odd
times using Eq. 30 and plotted using circles and triangles, respectively. Open
symbols indicate values of the correlation function that are negative and, for
ease of visualization, have been multiplied by $-1$. These negative points are
responsible for the behavior of the odd-site effective mass for
$C_{D}^{A^{0}}$ at early times. After folding the data around the midpoint of
the lattice, the correlator is defined for $t/a\in[0,N_{t}/2]=[0,48]$. Because
the form of Eq. 30 involves offsets by 2, the effective mass is defined on
times $t/a\in[2,N_{t}/2-2]=[2,46]$.

The statistical analysis consists of two stages. First, two-point functions
are analyzed in isolation. Second, the two- and three-point functions are
analyzed together to extract the form factors. Several features are common to
the fits in both stages. To avoid possible contamination from autocorrelation
in Monte Carlo time, the data are binned by 10 configurations prior to
fitting. The amount of binning was chosen by looking for stability and
saturation of errors in the fit results for the masses and form factors.

Our analysis employs standard Bayesian fits, which can be described generally
as least-squares regression to a model function $f(\bm{a})$ with parameters
$\bm{a}$ for some data set $D$. The likelihood function
$\text{pr}(D|\mathbf{a})\propto\exp\left[-\frac{1}{2}\chi^{2}\right]$ is
written in terms of the augmented chi-squared function $\chi^{2}=\chi^{2}_{\rm
data}+\chi^{2}_{\rm prior}$, with

$\displaystyle\chi^{2}_{\rm data}$
$\displaystyle=(\bar{y}-f(\bm{a}))^{T}\Sigma^{-1}(\bar{y}-f(\bm{a})),$ (28)
$\displaystyle\chi_{\rm prior}^{2}$
$\displaystyle=(\bm{a}-\tilde{\bm{a}})^{T}\tilde{\Sigma}^{-1}(\bm{a}-\tilde{\bm{a}}),$
(29)

where $\bar{y}$ is a vector with the data means, $\Sigma$ is the covariance
matrix, $\tilde{\bm{a}}$ is a vector with the prior values, and
$\tilde{\Sigma}$ is the prior covariance matrix. These expressions are
standard [76, 77, 78]. In the present analysis, $\bar{y}$ and $\Sigma$
correspond to the measured means and covariance matrices of the correlation
functions. The function $f(\bm{a})$ corresponds to the spectral decomposition
of Eqs. 19, 20 and 21, with the energies and matrix elements serving as the
parameters $\bm{a}$. The choices for the priors $\tilde{\bm{a}}$ are discussed
below. Instead of using the usual binned-sample covariance matrix in Eq. 28,
we used an improved estimator $\hat{S}_{n}$ employing nonlinear shrinkage,
which corrects for finite-sample-size effects [79]; for technical details, see
Appendix C.777To avoid possible confusion, we emphasize our correlator fits
use _nonlinear_ shrinkage. The chiral-continuum fits described below use
linear shrinkage, since it combines data from different ensembles, each with a
different statistical size; see the discussion in Sec. V.2. The general
procedure is as follows. First, we compute the binned variances $\sigma$.
Second, we compute the correlation matrix $C_{n}$ using the full (unbinned)
data. Third, we compute the shrinkage estimator $\hat{C}_{n}$, taking the
effective sample size to be the ratio of the total configurations to the bin
size. Finally, we construct
$\hat{S}_{n}=\operatorname{diag}(\sigma)\hat{C}_{n}\operatorname{diag}(\sigma)$.
Apart from the usage of shrinkage, a similar procedure has been employed in
the past, e.g., in Ref. [32]. In all cases, statistical uncertainties in the
fit parameters are determined via bootstrap resampling with $500$ draws. For
fits on bootstrap-resampled pseudoensembles, the covariance matrix is held
fixed to the binned-sample covariance matrix with shrinkage ($\hat{S}_{n}$)
for the full ensemble [80].

As mentioned above, the statistical analysis begins with two-point functions.
Figures 3, 4 and 5 display representative two-point functions and effective
masses for the pion, the kaon and the $D$ meson, respectively, on the
physical-mass $a\approx 0.12\leavevmode\nobreak\ \text{fm}$ ensemble with a
heavy-quark mass near its physical value for the $D$ meson. For the
correlation functions themselves, dramatic oscillations from opposite-parity
states are present only for the heavy mesons (see Fig. 5). When plotted in the
usual way, oscillations are visible in all the effective masses aside from the
zero-momentum pion. To reduce the distraction of opposite-parity states and
bring out the approach to the ground state, the effective mass is constructed
separately on even and odd time slices using

$am_{\rm
eff}(t)\equiv\frac{1}{2}\operatorname{arcCosh}\left[(C(t+2)+C(t-2))/2C(t)\right].$
(30)

In the effective mass plots, the triangle and circle markers correspond to the
even and odd time slices, respectively. As expected, the statistical noise
grows exponentially for correlators with nonzero momentum. High-momentum
correlators therefore become noisy at large times, especially those considered
here with $\bm{n}=(2,2,2)$ or $(4,0,0)$. Even so, clear plateaus spanning
several time slices are typically present in the effective mass at each
momentum. For $D$ mesons, contributions from excited states are visibly larger
when the interpolating operator $A_{0}$ is used. This observation informs
certain analysis choices below. The behavior of two-point functions is similar
to the ones shown in Figs. 3, 4 and 5 for other masses and lattice spacings.

Table 4: Preferred analysis settings for fits of two-point functions to Eq. 19 and Eq. 20. The same settings are applied uniformly across all ensembles. The larger $t_{\rm min}$ cut for $C_{H}^{A_{0}}(t)$ is taken to avoid the excited-state contributions visible in Fig. 5.  Correlator | $N_{\rm decay}+N_{\rm osc}$ |  $t_{\rm min}$ [fm] |  $t_{\rm max}$ cut
---|---|---|---
$C_{\pi}^{P}(t,\bm{p}=\bm{0})$ | 3 + 0 | $\approx 0.5$ | Noise $\leq 30\%$
$C_{K}^{P}(t,\bm{p}=\bm{0})$ | 3 + 1 | $\approx 0.5$ | Noise $\leq 30\%$
$C_{L}^{P}(t,\bm{p}\neq\bm{0})$ | 3 + 1 | $\approx 0.5$ | Noise $\leq 30\%$
$C_{H}^{P}(t)$ | 3 + 2 | $\approx 0.5$ | Noise $\leq 30\%$
$C_{H}^{A_{0}}(t)$ | 3 + 2 |  $\approx 0.75$–$0.85$ | Noise $\leq 30\%$
Table 5:  Priors used for the energies in fitting two-point functions to the spectral decomposition. All energy values are in MeV. At each lattice spacing, the values are converted to lattice units. Internally, the actual fit parameters are (the logarithm of) energy differences [76]. For instance, the splitting between the first and second excited states for the pion is $1700(400)-1300(400)=400(566)\leavevmode\nobreak\ \text{MeV}$. Priors for the amplitudes are discussed in the main text. For ensembles with heavier-than-physical pions, the central values for the priors for $E_{\pi}$ and $E_{K}$ are increased using the tree-level expectation from chiral perturbation theory for the quark-mass dependence of the hadron mass (see Sec. V.1 below).  $n$ | $E_{\pi}$ | $E_{\pi}^{\rm osc}$ | $E_{K}$ | $E_{K}^{\rm osc}$ | $E_{D}$ | $E_{D}^{\rm osc}$ | $E_{D_{s}}$ | $E_{D_{s}}^{\rm osc}$
---|---|---|---|---|---|---|---|---
0 | $135(50)$ | $500(300)$ | $498(100)$ | $800(300)$ | $1865(200)$ | $2300(700)$ | $1968(200)$ | $2317(200)$
1 | $1300(400)$ | — | $1460(400)$ | — | $2565(700)$ | $3000(700)$ | $2300(400)$ | $2713(400)$
2 | $1700(400)$ | — | $1860(400)$ | — | $3200(700)$ | — | $2700(400)$ | —

Each two-point correlator is fit to the corresponding spectral decomposition,
Eq. 19 or Eq. 20, using the choices in Table 4. We have verified that our
results are stable under reasonable variations of these choices, such as
including more states or changing the value of $t_{\rm min}$. The preferred
number of states is roughly the minimal number required to achieve
statistically significant fits (with, say, $\chi^{2}/{\rm DOF}\lesssim 1$ or
$p\gtrsim 0.1$ for goodness of fit888In this work, we use the augmented
$\chi^{2}$ when quoting reduced $\chi^{2}/{\rm DOF}$. Throughout the analysis,
judgements about fit quality are insensitive to the precise definition used,
and indistinguishable results are obtained for other reasonable definitions,
e.g., the alternative quality-of-fit metrics defined in Appendix B of Ref.
[81].), while the cuts on $t_{\rm min}$ and $t_{\rm max}$ are designed to
retain as much of the data as possible without undue contamination from
excited states at early times or statistical noise at late times. The choices
for the number of states and $t_{\rm min}$ are broadly similar to Fermilab-
MILC work on decay constants [32]. The main difference from Ref. [32] is that
the present analysis includes more states for the pion and kaon, e.g., $3+1$
versus $1+1$, in order to include data from shorter Euclidean times, which is
advantageous for the subsequent analysis with three-point functions. Table 5
summarizes the priors used for the energies in fitting the two-point
functions. For the amplitudes, we choose broad priors in lattice units:
$0.50(20)$ for decaying states and $0.1(1.0)$ for the oscillating states.

As an example, Fig. 6 shows the stability of the ground-state mass extracted
from $C_{\pi}^{P}(t,\bm{p}=\bm{0})$ on the physical-mass $a\approx
0.12\leavevmode\nobreak\ \text{fm}$ ensemble using fits with different choices
for the number of states and $t_{\rm min}/a$. Consistency with expectations
from the effective mass is demonstrated in Fig. 7. Similar studies inform the
other choices in Table 4.

Figure 6: Stability of the ground-state mass in multi-exponential fits to
$C_{\pi}^{P}(t,\bm{p}=\bm{0})$ on the physical-mass $a\approx
0.12\leavevmode\nobreak\ \text{fm}$ ensemble. The different colors show the
posterior values for the ground-state mass using different numbers of states.
The blue points indicating 1-state fits are aligned with the value of $t_{\rm
min}/a$ on the horizontal axis. The corresponding results for 2- and 3-state
fits at the same $t_{\rm min}/a$ are offset slightly to the right. The size of
the markers is proportional to the $p$-value of the fit. As described in Table
4, the preferred fit uses 3 states and $t_{\rm min}/a=4$ and is indicated by
the green star.  Figure 7:  Comparison of the ground-state mass from the
preferred fit and the effective mass for $C_{\pi}^{P}(t,\bm{p}=\bm{0})$ on the
physical-mass $a\approx 0.12\leavevmode\nobreak\ \text{fm}$ ensemble. To
reduce the visual impact of opposite parity states, the effective mass is
computed separately for even and odd times using Eq. 30 and plotted using
circles and triangles, respectively.

Figure 8:  Left: The pion dispersion relation. Within statistical
uncertainties, $E^{2}=m^{2}+\bm{p}^{2}$. Right: The behavior of the overlap
factor
$\left<\emptyset\vphantom{P^{(\pi)}\pi}\right|P^{(\pi)}\left|\pi\vphantom{\emptyset
P^{(\pi)}}\right>$, normalized by the value at zero momentum. Within
statistical uncertainties, the overlap factor is constant.

Figure 9:  Left: The kaon dispersion relation. Within statistical
uncertainties, $E^{2}=m^{2}+\bm{p}^{2}$. Right: The behavior of the overlap
factor
$\left<\emptyset\vphantom{P^{(\pi)}\pi}\right|P^{(\pi)}\left|\pi\vphantom{\emptyset
P^{(\pi)}}\right>$, normalized by the value at zero momentum. Within
statistical uncertainties, the overlap factor is constant.

The second stage of the analysis combines data from two-point and three-point
functions to extract $f_{0}$, $f_{\parallel}$, and $f_{\perp}$. The basic
procedure consists of simultaneous correlated fits to the spectral
decompositions, Eqs. 19, 20 and 21, for a particular value of the heavy-quark
mass and the current $J$ using the choices for the numbers of states and the
fit ranges in Table 4. For instance, a simultaneous correlated fit to
$C_{D}^{P}(t)$, $C_{\pi}^{P}(t,\bm{p})$, and $C_{D\to\pi}^{S}(t,T,\bm{p})$
furnishes $\left<\pi\vphantom{SD}\right|S\left|D\vphantom{\pi S}\right>$. For
consistency between the two-point and three-point functions, the fit window
for the three-point functions is taken to be $t\in[t_{\rm min}^{\rm
src},T-t_{\rm min}^{\rm snk}]$, where $t_{\rm min}^{\rm src}$ and $t_{\rm
min}^{\rm snk}$ are the values of $t_{\rm min}$ associated with the source and
sink operators, which in general differ. For example, when the $A_{0}$ sink
operator is used $t_{\rm min}^{\rm src}<t_{\rm min}^{\rm snk}$; see Table 4.
The Bayesian priors used in these fits incorporate knowledge about the ground-
state energies and overlap factors coming from the two-point fits. Let $M_{\rm
2pt}\pm\delta M_{\rm 2pt}$ denote the posterior value of the ground-state
energy emerging from a fit to Eq. 19 or Eq. 20, and let $E_{\rm
2pt}(\bm{p}^{2})\equiv\sqrt{M_{\rm 2pt}^{2}+\bm{p}^{2}}$ denote the value of
the energy obtained by boosting the central value. Similarly, let $A_{\rm
2pt}\pm\delta A_{\rm 2pt}$ denote the posterior value of the ground-state
amplitude from the same fit.999 At large times, a generic two-point function
is $C(t)=A_{\rm 2pt}^{2}\left(e^{M_{\rm 2pt}t}+e^{M_{\rm
2pt}(T-t)}\right)+\cdots$, so the amplitude $A_{\rm 2pt}$ contains the
momentum-dependent relativistic normalization of states in the denominator.
For the joint fits to the two- and three-point functions at zero momentum, the
central values for the amplitude and energy priors are taken to match the two-
point posterior central values ($M_{\rm 2pt}$ and $A_{\rm 2pt}$), while the
prior widths are taken to be ten times the posterior widths ($10\times\delta
M_{\rm 2pt}$ and $10\times\delta A_{\rm 2pt}$). For nonzero momentum, the
prior central values are obtained by boosting the corresponding ground-state
results assuming the continuum relativistic dispersion relation; the
fractional prior widths are taken to match the expected size of discretization
effects, e.g., $\delta E_{\rm 2pt}/E_{\rm
2pt}=\mathcal{O}(\alpha_{s}a^{2}\bm{p}^{2})$. Table 6 summarizes the choices
of these priors. For the excited states, the priors for the energy differences
are the same as those in Table 5, and the priors for the amplitudes are as
above.

For generic transition matrix elements, a prior of $V_{nm}=0.1(10)$ in lattice
units is used, where

$V_{nm}\equiv\frac{\left<\emptyset\vphantom{\mathcal{O}_{L}n}\right|\mathcal{O}_{L}\left|n\vphantom{\emptyset\mathcal{O}_{L}}\right>\left<n\vphantom{Jm}\right|J\left|m\vphantom{nJ}\right>\left<m\vphantom{\mathcal{O}_{H}\emptyset}\right|\mathcal{O}_{H}\left|\emptyset\vphantom{m\mathcal{O}_{H}}\right>}{4E_{L}^{(n)}(\bm{p})M_{H}^{(m)}}.$
(31)

For the special case of the ground state, the central value of $V_{00}$ is
estimated from the plateau in ratios following Eq. 108 below, and the width is
taken to be 50%.

Table 6:  Summary of how priors for simultaneous fits to two- and three-point functions incorporate information from the two-point fits. Values for $\alpha_{s}$ are given in Table 8.  Momentum | Quantity | Prior value
---|---|---
$\bm{p}^{2}=0$ | Energy | $M_{\rm 2pt}\pm 10\times\delta M_{\rm 2pt}$
| Amplitude | $A_{\rm 2pt}\pm 10\times\delta A_{\rm 2pt}$
$\bm{p}^{2}>0$ | Energy | $E_{\rm 2pt}\times\left(1\pm\alpha_{s}a^{2}\bm{p}^{2}\right)$
| Amplitude | $A_{\rm 2pt}\sqrt{M_{\rm 2pt}/E_{\rm 2pt}}\times\left(1\pm\alpha_{s}a^{2}\bm{p}^{2}\right)$

Figure 10:  Visual tests involving the ratio $R_{0}^{D\to\pi}$ in Eq. 24 to
form factors coming from the spectral decomposition. Data and results are
taken from the physical-mass $0.12\leavevmode\nobreak\ \text{fm}$ ensemble
with the charm-quark mass near its physical value. (Top) The approach to the
asymptotic plateau region for the ratio $R_{0}^{D\to\pi}$ at fixed momentum as
a function of the source-sink separation. The right panel shows the maximum
point from each color set of points on the left,
$\max_{t}R_{0}^{D\to\pi}(t,T)$. The horizontal black line in the top panels
shows the form factor’s posterior value, taken from the joint fit to the
spectral decomposition. (Bottom) The form factor’s momentum dependence. The
left panel shows the ratio $R_{0}^{D\to\pi}(t,T_{\rm max},\bm{p})$ at fixed
source-sink separation, with each color corresponding to a different momentum.
Horizontal lines show the central value for the form factors coming from the
fits (including all source-sink separations $T$). The right panel shows the
smooth momentum dependence of $f_{0}^{D\to\pi}(q^{2})$.

Once statistically acceptable fits (e.g., $\chi^{2}/{\rm DOF}\lesssim 1$ or
$p\gtrsim 0.1$) are obtained, a variety of visualizations give confidence that
the bare form factors have been extracted reliably. For instance, the fits
must reproduce the data visually with reasonable uncertainties and give
results for the ground-state masses and overlap factors that agree with the
initial analysis of two-point functions in isolation. As the priors in Table 6
suggest, energies are expected to satisfy the continuum dispersion relation,
$E^{2}=(M^{2}+\bm{p}^{2})(1+\mathcal{O}(\alpha_{s}a^{2}\bm{p}^{2}))$, and
overlap factors are expected to be constant, since only point-like
interpolators were used for the source and sink operators. Figures 8 and 9
demonstrate that both conditions are well satisfied. The blue points
correspond to the posterior (“best-fit”) results, while the dashed lines show
the size of the priors for $\bm{p}^{2}>0$, as defined in Table 6. As the
figure shows, the posteriors typically are much narrower than the priors. We
have verified that statistically consistent results, with similar statistical
precision, are obtained if the priors are relaxed by inflating the width by a
factor of ten. Figure 10 shows representative results for joint fits for
$D\to\pi$. The top rows show the approach to the asymptotic ($T/a\to\infty$)
plateau region. In the top left panel, data are plotted at fixed momentum
$\bm{p}=2\pi(1,0,0)/N_{s}a$, with each color corresponding to a different
source-sink separation $T$. The top right panel shows the approach to the
asymptotic plateau versus $T/a$, with each point corresponding to the maximum
point in the curves on the top left:
$\max_{t}R_{0}^{D\to\pi}(t,T,\bm{p}=2\pi(1,0,0)/N_{s}a)$. The horizontal black
line in the top panels shows the form factor’s posterior value, taken from the
joint fit to the spectral decomposition. The bottom panels shows the form
factor’s momentum dependence. In the bottom left panel, the data correspond to
the ratio $R_{0}^{D\to\pi}(t,T_{\rm max},\bm{p})$, with each color
corresponding to a different momentum. In each case, only the largest source-
sink separation $T_{\rm max}$ is displayed. Horizontal lines denote the
posterior central values for the form factors, coming from fits including all
source-sink separations $T$. The bottom right panel shows the smooth momentum
dependence of $f_{0}^{D\to\pi}(q^{2})$. Additional details, along with similar
figures for the decays $D\to K$ and $D_{s}\to K$ are given in Sec. D.1.

###  IV.3 Nonperturbative renormalization

Bare matrix elements are renormalized nonperturbatively by imposing the PCVC
relation, Eq. 10. Figure 11 shows the matrix elements entering this
expression, before and after renormalization, for the physical-mass $a\approx
0.12\leavevmode\nobreak\ \text{fm}$ ensemble with the charm-quark mass near is
physical value. The black points show the quantity
$(m_{c}-m_{q})\left<L\vphantom{SH}\right|S\left|H\vphantom{LS}\right>-q^{\mu}\left<L\vphantom{V^{\mu}H}\right|V^{\mu}\left|H\vphantom{LV^{\mu}}\right>$
with $L\in\\{\pi,K\\}$, $H\in\\{D,D_{s}\\}$, and $m_{q}\in\\{m_{l},m_{s}\\}$.
The fact that the open black circles differ slightly from zero gives a visual
indication that the renormalization factors $Z_{V^{0}}$ and $Z_{V^{i}}$ are
necessary to satisfy PCVC. The closed black squares, statistically consistent
with zero, show the precision with which the PCVC relation is satisfied after
renormalization.

Figure 11: Matrix elements (in arbitrary units) entering the PCVC relation,
Eq. 10, before and after renormalization for the physical-mass $a\approx
0.12\leavevmode\nobreak\ \text{fm}$ ensemble with the charm-quark mass near
its physical value, for the decays $D\to\pi$ (top left), $D\to K$ (top right),
and $D_{s}\to K$ (bottom) The open circles denote bare quantities, while the
filled squares are renormalized.

In principle, much freedom exists for extracting the vector-current
renormalization factors. The present analysis fits the bare matrix elements as
a function of momentum to Eq. 10 for each ensemble and choice of $m_{h}$,
treating $Z_{V^{0}}$ and $Z_{V^{i}}$ as free parameters. Recall $Z_{m}Z_{S}=1$
for the local staggered scalar current. When constructing the renormalized
matrix elements, correlations between the bare matrix elements and $Z_{V^{0}}$
and $Z_{V^{i}}$ are incorporated via the bootstrap resampling discussed above.

Figure 12: Vector-current renormalization factors $Z_{V^{0}}$, in all cases
with the charm-quark mass near its physical value. (Left) Results for the
quark-level transition $c\to l$ appearing in the decays $D\to\pi$ and
$D_{s}\to K$. As discussed in the text, the data are taken from fits to the
$D_{s}\to K$ data. (Right) Results for the quark-level transition $c\to s$
appearing in the decay $D\to K$. The light gray points show the
renormalization factors computed by HPQCD on the same ensembles, using the
same local $V^{0}$ current but slightly different valence masses [62]. The two
sets of values agree at $1$–$2\sigma$. Figure 13: Vector-current
renormalization factors $Z_{V^{i}}$, in all cases with the charm-quark mass
near its physical value. (Left) Results for the quark-level transition $c\to
l$ appearing in the decays $D\to\pi$ and $D_{s}\to K$. As discussed in the
text, the data are taken from fits to the $D_{s}\to K$ data. (Right) Results
for the quark-level transition $c\to s$ appearing in the decay $D\to K$.

Figures 12 and 13 show the results for renormalization factors for the
temporal and spatial components of the vector current, respectively, in all
cases for the data with the charm-quark mass near its physical value. The
transition $c\to l$ appears in both $D\to\pi$ and $D_{s}\to K$ decays,
differing only by the spectator quark. The data for the latter decay are
statistically more precise, which in turns yields more precise values for the
renormalization factors of the $\bar{l}c$ currents. We thus use the
renormalization factors extracted from the $D_{s}\to K$ data to renormalize
both $D_{s}\to K$ and $D\to\pi$ matrix elements.

At a given lattice spacing, uncertainties both in the bare matrix elements
(coming from the correlator fits) and in the renormalization factors
contribute to the total error budget for the form factors. The relative
importance of the renormalization error depends both on the form factor
($f_{\parallel}$ or $f_{\perp}$) and the momentum. For instance, for the
$D\to\pi$ decay on the physical-mass $a\approx 0.12\leavevmode\nobreak\
\text{fm}$ ensemble, the renormalization error in $f_{\parallel}$ from
$Z_{V^{0}}$ is $\lesssim 0.1\%$. For the same decay and lattice spacing, the
renormalization error in $f_{\perp}$ from $Z_{V^{i}}$ is around $1\%$. For
comparison, the individual statistical errors in both $f_{\parallel}$ and
$f_{\perp}$ (neglecting the renormalization error) range from around $1\%$ at
low momentum to around $8\%$ at large momentum. These observations are
consistent with the expectation that renormalization with PCVC should enable
sub-percent determinations of form factors.

##  V Chiral-continuum analysis

This section describes our chiral-continuum analysis, yielding results for
$f_{+}(q^{2})$ and $f_{0}(q^{2})$ at physical quark mass and in the continuum
limit. Section V.1 describes the fit function used in the analysis and its
connection to effective field theory (EFT). Section V.2 presents the results
of the fits and describes our definition of the physical point in isospin-
symmetric QCD. Section V.3 presents a cross check on our results by
constructing $f_{+}$ and $f_{0}$ in different ways. Section V.4 re-expresses
our results in a compact form using the model-independent $z$ expansion.
Section V.5 considers the spectator dependence of the form factors by
comparing our results for $D\to\pi$ and $D_{s}\to K$. Finally, Sec. V.6
compares our form factors with published results in the literature.

###  V.1 Description of the chiral-continuum fit function

Together, the bare matrix elements and renormalization factors calculated in
Sec. IV furnish the form factors $f_{\parallel}$, $f_{\perp}$, and $f_{0}$ at
four different lattice spacings, three different pion masses, and several
values of the heavy-quark mass. These results are extrapolated to the
continuum limit and interpolated to the physical point using guidance from
effective field theory.

We treat the light-quark mass dependence of the form factors $f_{\parallel}$
and $f_{\perp}$ with SU(2) heavy-meson rooted staggered chiral-perturbation
theory [82, 83]. Following earlier work [84, 85], we use the version for a
hard final-state hadron [86, 87, 88], hereafter referred to as “hard SU(2)
$\chi$PT.” We include the complete set of chiral logarithms and analytic
corrections through next-to-leading order (NLO) in the chiral expansion. To
account for truncation errors, we also include all analytic terms consistent
with the power-counting scheme of Ref. [82, 83] through next-to-next-leading
order (NNLO). These choices amount to considering the following functional
form for $P\in\\{\parallel,\perp,0,+\\}$:

$\begin{split}w_{0}^{d_{P}}f_{P}(E)=\frac{c_{0}}{w_{0}\left(E+\Delta_{xy,P}\right)}\times\Big{[}1&+\delta
f_{\rm logs}+c_{l}\chi_{l}+c_{s}\chi_{s}+c_{H}\chi_{H}+c_{E}\chi_{E}\\\
&+c_{l^{2}}\chi_{l}^{2}+c_{ls}\chi_{l}\chi_{s}+c_{s^{2}}\chi_{s}^{2}\\\
&+c_{lH}\chi_{l}\chi_{H}+c_{lE}\chi_{l}\chi_{E}+c_{sH}\chi_{s}\chi_{H}+c_{sE}\chi_{s}\chi_{E}\\\
&+c_{H^{2}}\chi_{H}^{2}+c_{HE}\chi_{H}\chi_{E}+c_{E^{2}}\chi_{E}^{2}\\\
&+\delta f_{\rm artifacts}^{(a^{2}+h^{2})}\Big{]},\end{split}$ (32)

where the exponent $d_{P}\in\\{1/2$, $-1/2$, $0$, $0\\}$ for
$P\in\\{\parallel,\perp,0,+\\}$ (respectively) and $c_{0}$ is a dimensionless
constant. Although, in principle, the function describing the chiral
logarithms, $\delta f_{\rm logs}$, depends on the form factor, in hard SU(2)
$\chi$PT it is the same for all $P\in\\{\parallel,\perp,0,+\\}$, as discussed
below.

The leading pole factor (in terms of the final-state hadron’s energy $E$)
arises from the exchange of a virtual $W$ boson, which couples to an excited
meson $D_{x}^{*}$ composed of $c$ and the final-state quark $x$, contributing
a factor proportional to

$\frac{1}{E+\Delta_{xy,P}}=\frac{2M_{D_{y}}}{M_{D^{*}_{x}}^{2}-q^{2}}.$ (33)

The intrinsic angular momentum and parity of the $D_{x}^{*}$ are those of the
virtual $W$ boson, which is $J^{P}=1^{-}$ for $f_{+}$ and $J^{P}=0^{+}$ for
$f_{0}$. According to the leading-order expectations of the heavy-quark
expansion [89], the same pole arises pairwise for $f_{\perp}$ as $f_{+}$, and
similarly for the pair $f_{\parallel}$ and $f_{0}$ (cf. Table 7). Equation 33
implies that location of the pole in the energy can be written as

$\Delta_{xy,P}=\frac{M_{D^{*}_{x}(J^{P})}^{2}-M_{D_{y}}^{2}-M_{L}^{2}}{2M_{D_{y}}},$
(34)

where $y$ is the spectator quark and $L\in\\{\pi,K\\}$ is the final-state
hadron. Values for the $\Delta_{xy,P}$ are collected in Table 7 for the decays
of interest.

Table 7: Approximate pole locations $\Delta_{xy,P}$ appearing in the decays $D\to\pi$, $D_{s}\to K$, and $D\to K$. Decay | $c\to x$ | $J^{P}$ | $D_{x}^{*}(J^{P})$ | $D_{y}$ | $L$ | $\Delta_{xy,P}$ | $(\Delta_{xy,P})^{\rm PDG}$
---|---|---|---|---|---|---|---
$D\to\pi$ | $c\to l$ | $1^{-}$ | $D^{*}$ | $D$ | $\pi$ | $\Delta_{ll,+/\perp}$ | $140\leavevmode\nobreak\ \text{MeV}$
|  | $0^{+}$ | $D^{*}_{0}(2300)$ | $D$ | $\pi$ | $\Delta_{ll,0/\parallel}$ | $480\leavevmode\nobreak\ \text{MeV}$
$D_{s}\to K$ | $c\to l$ | $1^{-}$ | $D^{*}$ | $D_{s}$ | $\pi$ | $\Delta_{ls,+/\perp}$ | $-25\leavevmode\nobreak\ \text{MeV}$
|  | $0^{+}$ | $D^{*}_{0}(2300)$ | $D_{s}$ | $\pi$ | $\Delta_{ls,0/\parallel}$ | $300\leavevmode\nobreak\ \text{MeV}$
$D\to K$ | $c\to s$ | $1^{-}$ | $D_{s}^{*}$ | $D$ | $K$ | $\Delta_{sl,+/\perp}$ | $200\leavevmode\nobreak\ \text{MeV}$
|  | $0^{+}$ | $D_{s0}^{*}(2317)$ | $D$ | $K$ | $\Delta_{sl,0/\parallel}$ | $440\leavevmode\nobreak\ \text{MeV}$

The $\chi_{n}$ are dimensionless expansion parameters defined according to

$\displaystyle\chi_{l}$ $\displaystyle=\frac{(M_{\pi}^{\rm
sim})^{2}}{8\pi^{2}f_{\pi}^{2}},$ (35) $\displaystyle\chi_{s}$
$\displaystyle=\frac{(M_{K}^{\rm sim})^{2}-(M_{K}^{\rm
PDG})^{2}}{8\pi^{2}f_{\pi}^{2}},$ (36) $\displaystyle\chi_{E}$
$\displaystyle=\frac{\sqrt{2}E}{4\pi f_{\pi}},$ (37) $\displaystyle\chi_{H}$
$\displaystyle=\frac{\Lambda_{\rm HQET}}{M_{H_{(s)}}^{\rm
sim}}-\frac{\Lambda_{\rm HQET}}{M_{D_{(s)}}^{\rm PDG}},$ (38)

where $f_{\pi}$ is the physical pion decay constant and $\Lambda_{\rm HQET}$
is the scale of heavy-quark effective theory. As in Ref. [32], we take
$\Lambda_{\rm HQET}=800\leavevmode\nobreak\ \text{MeV}$. The parameters
$\chi_{l}$ and $\chi_{E}$ describe the analytic dependence on the light-quark
mass $m_{l}$ (via the leading-order expression $M_{\pi}^{2}=2\mu m_{l}$) and
the final-state hadron energy $E$, respectively. Their normalization is such
that, according to the typical $\chi$PT power counting, the corresponding
coefficients in the fit function $c_{l}$ and $c_{E}$ are expected to be of
order 1. The parameter $\chi_{H}$ describes the heavy quark mass mistuning
through the difference between the simulated heavy meson mass
$M_{H_{(s)}}^{\rm sim}$ and the physical $M_{D^{0}}$ or $M_{D_{s}}$ from Ref.
[25]. This term allows for a simultaneous description of results across
several different heavy quark masses. Finally, $\chi_{s}$ describes the
strange-quark mass mistuning.

In hard SU(2) $\chi$PT, the chiral logarithms for $f_{\parallel}$ and
$f_{\perp}$ in Eq. 32 have the common form [84, 85],

$\begin{split}\delta f_{\rm logs}^{\text{SU(2)}}&=\frac{1}{(4\pi
f_{\pi})^{2}}\left(-\frac{1}{16}\sum_{\xi}\mathcal{I}_{1}(M_{\pi,\xi})+\frac{1}{4}\mathcal{I}_{1}(M_{\pi,I})+\mathcal{I}_{1}(M_{\pi,V})-\mathcal{I}_{1}(M_{\eta,V})+[V\to
A]\right)\\\ &\phantom{=}\times\begin{cases}(1+3g^{2}),\quad D\to\pi\\\
3g^{2},\phantom{(1+{})}\quad D\to K\\\ 1,\phantom{({}+3g^{2})}\quad D_{s}\to
K\end{cases},\end{split}$ (39)

where $\mathcal{I}_{1}(M)\equiv M^{2}\ln(M^{2}/\Lambda^{2})+4I_{1}^{\rm
FV}(M,ML)$, with $I_{1}^{\rm FV}(M,ML)$ being a calculable finite-volume
correction to the chiral logarithm which vanishes exponentially for large
volumes; see Sec. VI.3 below. Hard SU(2) $\chi$PT enjoys the further
simplification that nonanalytic self-energy corrections vanish for all three
decays considered here. These expressions [Eq. 39 and the self energies] were
originally derived for a non-staggered heavy quark, but because the heavy-
quark taste is conserved in all-staggered $\chi$PT, they hold in the present
case too [90]. In heavy-meson $\chi$PT, compact expressions are available for
$f_{\perp}$ and $f_{\parallel}$ [91, 92], while the corresponding results for
$f_{+}$ and $f_{0}$ follow as linear combinations. Because of their simple
connection to heavy-meson $\chi$PT, previous lattice calculations have
historically worked primarily in terms of $f_{\parallel}$ and $f_{\perp}$.
However, since the chiral logarithms have the same functional form for
$f_{\parallel}$ and $f_{\perp}$ in hard SU(2) $\chi$PT [see Eq. 39], the same
functional form also describes the chiral logarithms for $f_{0}$ and $f_{+}$.
In other words, Eq. 32 may be used directly for all four form factors, with a
$1^{-}$ pole for $f_{+,\perp}$ or a $0^{+}$ pole for $f_{0,\parallel}$.

Following Ref. [83], the arguments of the chiral logarithms involve the masses
of mesons with different tastes $\xi\in\\{I,P,V,A,T\\}$, that can be expressed
as

$\displaystyle M_{\pi,\xi}^{2}$ $\displaystyle=M_{uu,\xi}^{2}=M_{dd,\xi}^{2},$
(40) $\displaystyle M_{ij,\xi}^{2}$
$\displaystyle=\mu(m_{i}+m_{j})+\Delta_{\xi},$ (41) $\displaystyle
M_{\eta,V(A)}^{2}$
$\displaystyle=M_{uu,V(A)}^{2}+\frac{1}{2}\delta^{\prime}_{V(A)},$ (42)
$\displaystyle\bar{\Delta}$
$\displaystyle=\frac{1}{16}\sum_{\xi}\Delta_{\xi}.$ (43)

The low-energy constant $\mu$ and the taste splittings $\Delta_{\xi}$ have
been tabulated for these ensembles in Ref. [27]. At NLO in the chiral
expansion, the taste splittings $\Delta_{\xi}$ and the hairpin parameters
$\delta^{\prime}_{V,A}$ both scale like $\alpha_{s}^{2}a^{2}$, so their ratio
remains approximately constant as the lattice spacing changes. We follow Ref.
[32] and take $\delta^{\prime}_{A}/\bar{\Delta}=-0.88(09)$ and
$\delta^{\prime}_{V}/\bar{\Delta}=+0.46(23)$.

Chiral logs described above include the dominant discretization effects coming
from the taste-symmetry breaking of staggered fermions at NLO in the chiral
expansion. We also remove the leading-order (tree-level) heavy-quark
discretization effects in the form factors prior to fitting by applying a
multiplicative normalization factor $Z^{\rm HQET,LO}_{hx}$, described in
Appendix B. Because of the tree-level improvement of the HISQ action, the
remaining discretization effects are expected to arise at order
$\alpha_{s}(a\Lambda)^{2}$ or $\alpha_{s}(am_{h})^{2}$, where $\Lambda$ is the
scale of generic discretization effects. They are thus expected to be well
described by an expansion in terms of the parameters $x_{a^{2}}$ and $x_{h}$:

$\displaystyle x_{a^{2}}$
$\displaystyle=\frac{a^{2}\bar{\Delta}}{8\pi^{2}w_{0}^{2}f_{\pi}^{2}},$ (44)
$\displaystyle x_{h}$ $\displaystyle=\frac{2}{\pi}am_{h}.$ (45)

The quantity $x_{a^{2}}$ gives a dimensionless measure of order
$\alpha_{s}(a\Lambda)^{2}$ discretization corrections, while $x_{h}$ is the
natural expansion parameter for heavy-quark discretization effects. The HISQ
action was designed specifically to control lattice artifacts for charm
physics, and the leading heavy-quark corrections are suppressed both by
$\alpha_{s}$ and the velocity $v\approx\sqrt{1/10}$ of the charm quark within
the heavy hadron. Our preferred model thus takes the following simple _Ansatz_
for the discretization effects

$\displaystyle\delta f_{\rm artifacts}^{(a^{2}+h^{2})}$
$\displaystyle=c_{a^{2}}x_{a^{2}}+\alpha_{s}vc_{h^{2}}x_{h}^{2}.$ (46)

Values for $\alpha_{s}$ are given in Table 8.

To check for truncation effects with high-order discretization effects, we
also consider variations:

$\displaystyle\delta f_{\rm artifacts}^{(a^{2}+h^{2}+a^{4})}$
$\displaystyle=c_{a^{2}}x_{a^{2}}+\alpha_{s}vc_{h^{2}}x_{h}^{2}+c_{a^{4}}x_{a^{2}}^{2},$
(47) $\displaystyle\delta f_{\rm artifacts}^{(a^{2}+h^{2}+h^{4})}$
$\displaystyle=c_{a^{2}}x_{a^{2}}+\alpha_{s}vc_{h^{2}}x_{h}^{2}+vc_{h^{4}}x_{h}^{4}.$
(48)

With the HISQ action, discretization effects of order $x_{h}^{4}$ are
suppressed by $v$ (at the tree level) or by $\alpha_{s}$. These suppression
factors are numerically similar enough that the last term in Eq. 48 tests
both.

Table 8: Values for the strong coupling constant, which are based on the continuum value of $\alpha_{s}(5\leavevmode\nobreak\ \text{GeV},N_{f}=4)$ from Ref. [93]. Continuum perturbation theory is used to convert to the $\alpha_{V}$ scheme and to run to the scale $2/a$ [94].  $\approx a$ [fm] | $\alpha_{s}(2/a)$
---|---
0.15 | 0.3509
0.12 | 0.3091
0.088 | 0.2646
0.06 | 0.2236
0.042 | 0.2036

###  V.2 Chiral-continuum fits

We perform correlated fits, using the methodology described around Eqs. 28 and
29, to Eq. 32 with $\delta f_{\rm artifacts}^{(a^{2}+h^{2})}$ in Eq. 46 for
each of the above form factors for $D\to\pi$, $D\to K$, and $D_{s}\to K$
including all of the ensembles and heavy-quark masses described in Table 1.

In our preferred fits, the input data for the form factors
$f_{\perp/\parallel/0}$ are defined using a single matrix element each via
Eqs. 6, 7 and 8. The chiral-continuum fit results using the alternative
constructions $f^{\rm alt}_{\rm P}$ of Eqs. 11 and 12 are considered below in
the analysis of systematic effects.

Table 9: Summary of the priors used in the chiral-continuum fits to Eq. 32. Values for $(\Delta_{xy,P})^{\rm PDG}$ in the different decays are given in Table 7.  Parameter | Value
---|---
$c_{0}$ | $1\pm 10$
$c_{n}$ | $0\pm 1$
$g$ | $0.5\pm 0.2$
$\Delta_{xy,P}$ | $(\Delta_{xy,P})^{\rm PDG}\pm 200\leavevmode\nobreak\ \text{MeV}$

The free parameters varied in the fits are the coefficients $c_{n}$, the
coupling $g$, and the mass splittings $\Delta_{xy,\rm P}$. The Bayesian priors
for these parameters are given in Table 9. The leading coefficient $c_{0}$ is
well determined by the data, so the preferred analysis uses a broad prior (the
results are insensitive to the central value, and any reasonable variation
gives indistinguishable results). The chiral-continuum fit function is based
on power-counting arguments from effective field theory, according to which
the coefficients $c_{n}$ are expected to be of order unity. The preferred
analysis therefore uses priors of $0\pm 1$ for the parameters $c_{n}$. The
dimensionless (reduced) “$DD^{*}\pi$” coupling appearing as a coefficient of
the chiral logarithms is expected to be $g\approx 0.5$, both from experimental
measurement [95, 96, 97] and previous lattice-QCD calculations [98, 99, 100,
101, 102, 103, 104]. For compatibility with these results, our fits take a
prior of $0.5\pm 0.2$. Because a broad width is used for $\Delta_{xy,\rm
P}^{*}$, and since the fits are insensitive to the precise value, the priors
do not distinguish between the $J^{P}=0^{+}$ and $1^{-}$ states. The other
inputs to the fits are the measured initial- and final-state hadron masses on
each ensemble, the staggered parameters described in the previous section, and
the pion decay constant (which is held fixed to its physical value in Table
10).

A few words are in order regarding our choice of intermediate scale setting
using $w_{0}/a$ and its role in the chiral-continuum fit function. The
dimensionless expansion parameters $\chi_{l}$ and $\chi_{E}$ contain factors
of the pion decay constant in the denominator. Using the continuum values for
$f_{\pi}$ and $w_{0}$ in Table 10, we express the denominator as a
dimensionless number. For the numerator, we use the measured values of
$aM_{\pi}$ and $w_{0}/a$ to construct the dimensionless product. In other
words, on each ensemble $\chi_{l}$ is constructed as

$\displaystyle\chi_{l}\rightarrow\left(\frac{(w_{0}/a)^{\rm
sim}}{w_{0}}\right)^{2}\frac{(aM_{\pi}^{\rm sim})^{2}}{8\pi^{2}(f_{\pi}^{\rm
PDG})^{2}},$ (49)

and similarly for $\chi_{E}$ and $\chi_{H}$. In the continuum limit and at the
physical point, the $w_{0}$ dependence cancels in all the analytic terms.
Since $f_{+}$ and $f_{0}$ are dimensionless, the only residual scale-setting
dependence enters through the pole term, where the energy $E$ (in physical
units) must be converted to $w_{0}$ units.

On each ensemble, the input data for the form factors and meson masses and
energies are correlated using the results of the bootstrap fits from Sec. IV.
The resulting correlation matrices tend to be near-singular, with small,
poorly determined eigenvalues posing a difficult challenge for the fits. For a
given form factor, the dominant source of these eigenvalues is highly
correlated data at nearby heavy valence masses (e.g., $0.9\,m_{c}$,
$1.0\,m_{c}$, and $1.5\,m_{c}$). A common solution to this problem is SVD
cuts, which have recently been used in another lattice-QCD analysis of $D\to
K$ form factors [62] and which are summarized lucidly in Ref. [105]. Another
solution is shrinkage of the eigenvalue spectrum, as described in Appendix C.

In the analysis of correlation functions in Sec. IV, we could use nonlinear
shrinkage, which has the desirable feature of not involving any tunable
parameters. As described in Appendix C, however, the amount of shrinkage
applied to the eigenvalue spectrum is controlled by the concentration ratio
(the ratio of the number of random variables to the number of independent
statistical samples). Since the chiral-continuum extrapolation combines data
from different ensembles, there is no clear-cut concentration ratio. For the
chiral-continuum fits, we therefore employ linear shrinkage, which entails a
parameter $\lambda$. We find that $\lambda=0.1$ is large enough to regulate
the small eigenvalues (thus giving good fits) without discarding correlations
unnecessarily. As with SVD cuts [105], linear shrinkage improves the quality
of fits and tends to increase the uncertainty in the posterior values. The
quantitative effect of linear shrinkage is discussed alongside other
systematic effects in Sec. VI.1. A qualitative comparison of nonlinear
shrinkage, linear shrinkage, and SVD cuts is given in Appendix C.

These fits deliver the form factors in the continuum limit and at the physical
point. The continuum limit of Eq. 32 is defined by setting $\delta f_{\rm
artifacts}^{(a^{2}+h^{2})}$ equal to zero and setting the taste splittings to
zero in $\delta f_{\rm P,logs}$. The physical point is defined by setting the
input meson masses equal to their physical values, given in Table 10. By
construction, all quantities involving $\chi_{H}$ also vanish identically at
the physical mass of the decaying heavy meson. Our simulations and chiral
analysis are both done in the isospin limit (i.e., with a pair of degenerate
quarks with mass $m_{l}=(m_{u}+m_{d})/2$), so the final results for the form
factors correspond to QCD in the isospin limit. The physical meson masses in
Table 10 are chosen accordingly, following the prescription in Ref. [36]. The
systematic uncertainty with the isospin-symmetric approximation is discussed
below in Sec. VII.2.

Figure 14:  The result of the chiral-continuum fit for the $D\to\pi$ form
factors constructed using Eqs. 7, 6 and 8 in units of the gradient-flow scale
$w_{0}$. For visual clarity, only the physical-mass ensembles with heavy
valence masses $m_{h}/m_{c}\in\\{0.9,1.0,1.1\\}$ are shown, although all
ensembles and heavy valence masses in Table 1 are included in the fit.
Different colors label different lattice spacings and different shapes
correspond to the different values of the heavy-quark mass.

The results for the $D\to\pi$ form factors are shown in Fig. 14. To avoid
plotting many overlapping data and curves, the figures restrict to the three
ensembles with physical-mass pions and heavy valence masses with
$m_{h}/m_{c}\in\\{0.9,1.0,1.1\\}$. In all cases, the nearly coincident data
around the physical charm mass ($m_{h}\approx m_{c}$) suggest a mild
dependence on the lattice spacing. The black band denotes the result in the
continuum limit and at the physical point. The results for $D\to K$ and
$D_{s}\to K$ are quite similar and given in Sec. D.2. Table 11 summarizes the
fit quality for the preferred fits.

Table 10:  External inputs used to define the physical point in isospin-symmetric QCD using Eq. 32. As described in the text, the experimentally measured values of the heavy mesons are also used implicitly as inputs in Eq. 32. Quantity | Value | Reference
---|---|---
$f_{\pi}$ | $130.2(8)\leavevmode\nobreak\ \text{MeV}$ | Ref. [36]
$M_{\pi^{0}}$ | $134.9768(5)\leavevmode\nobreak\ \text{MeV}$ | Ref. [106]
$M_{K^{0}}$ | $497.611(13)\leavevmode\nobreak\ \text{MeV}$ | Ref. [106]
$w_{0}$ | $0.17177(67)\leavevmode\nobreak\ \text{fm}$ | Ref. [36]
$M_{D^{0}}$ | $1864.83(05)\leavevmode\nobreak\ \text{MeV}$ | Ref. [106]
$M_{D_{s}}$ | $1968.34(07)\leavevmode\nobreak\ \text{MeV}$ | Ref. [106]
Table 11:  Summary of the reduced $\chi^{2}$ values and associated degrees of freedom (in brackets) for the preferred fits to Eq. 32 for all decays and form factors. | $D\to\pi$ | $D\to K$ | $D_{s}\to K$
---|---|---|---
$f_{0}$ | 0.91 [126] | 0.48 [128] | 1.31 [134]
$f_{\parallel}$ | 0.59 [112] | 0.41 [123] | 0.88 [128]
$f_{\perp}$ | 0.64 [110] | 0.32 [111] | 0.66 [113]
$f_{+}$ | 0.59 [106] | 0.29 [109] | 0.60 [111]

###  V.3 Alternative constructions of $f_{+}$ and $f_{0}$

Our default construction for $f_{0}$ is given by Eq. 8 and obtained in the
preceding section. In an analogous way we construct $f_{+}$ from the
continuum-limit results for $f_{0}$ and $f_{\perp}$ in the preceding section
following Eq. 9. As discussed in Sec. II, the PCVC relation in Eq. 10 provides
the alternative constructions given in Eqs. 11 and 12. Additional freedom
exists in whether the linear combinations in Eqs. 11, 9 and 12 are taken
before or after the chiral-continuum limit. A comparison of the different
constructions is given in Fig. 15 for $D\to\pi$ and $D_{s}\to K$ ($D\to K$ is
similar), where excellent stability is observed throughout the kinematic
range. In the legend, the notation $CL$ specifies whether the continuum limit
is taken before or after computing the linear combination
[$CL(f_{\perp})+CL(f_{0})$ versus $CL(f_{\perp}+f_{0})$, respectively]. In all
cases, $f_{\parallel}$ and $f_{\perp}$ are directly related to vector matrix
elements via Eqs. 6 and 7. Because the results are statistically consistent,
our preferred analysis takes the results with the best statistical precision
(our default analysis). We take $f_{0}$ from Eq. 8. We construct $f_{+}$ via
Eq. 9, using results for $f_{\perp}$ and $f_{0}$ given by Eqs. 7 and 8, each
separately extrapolated to the continuum limit.

Figure 15:  Comparison of the form factors coming from different continuum-
limit constructions for the decays $D\to\pi$, and $D_{s}\to K$. Points are
offset horizontally for readability; the same value of $q^{2}$ is used in each
grouping. The notation used in the legend is explained in the main text. The
black points denote the preferred results with the best statistical precision.
Similar agreement was also found for $D\to K$.

###  V.4 Model-independent $z$ expansion

The previous section gave results for $f_{+}(q^{2})$ and $f_{0}(q^{2})$ in the
continuum limit and at the physical point. To facilitate comparison with
experimental measurements and other theoretical calculations, it is convenient
to re-express our results using the $z$ expansion. To start, consider the
decays $D\to\pi$ and $D\to K$. The $z$ expansion leverages the known analytic
structure of the form factors in the complex $q^{2}$-plane to write the form
factors as a rapidly convergent expansion in the variable $z(q^{2},t_{0})$,

$\displaystyle
z(q^{2},t_{0})=\frac{\sqrt{t_{+}-q^{2}}-\sqrt{t_{+}-t_{0}}}{\sqrt{t_{+}-q^{2}}+\sqrt{t_{+}-t_{0}}},$
(50)

where $t_{+}=(M_{D}+M_{L})^{2}$ denotes the start of the multiparticle cut,
$L\in\\{\pi,K\\}$, and $t_{0}\in[0,t_{+}]$ can be chosen for convenience. This
map sends the branch cut onto the unit circle, $|z(q^{2},t_{0})|=1$, and the
rest of the first Riemann sheet onto the open unit disk, $|z(q^{2},t_{0})|<1$.
Note that

$\displaystyle z(t_{+},t_{0})$ $\displaystyle=-1,$ (51) $\displaystyle
z(t_{0},t_{0})$ $\displaystyle=0,$ (52) $\displaystyle z(-\infty,t_{0})$
$\displaystyle=+1.$ (53)

Further, Eq. 50 maps the physical region for semileptonic decay onto an
interval on the real axis. Similar considerations apply for the decay
$D_{s}\to K$, except that the multiparticle cut begins at
$t_{+}=(M_{D}+M_{\pi})^{2}$ [and _not_ at $(M_{D_{s}}+M_{K})^{2}$]. Below, we
take $t_{0}=0$, so $q^{2}\in[0,q^{2}_{\rm max}]$ is mapped to $z\in[0,-z_{\rm
max}]$. Because $-z_{\rm max}\approx 0.332$, $0.190$, and $0.192$, for
$D\to\pi$, $D\to K$, and $D_{s}\to K$, respectively, one expects a series
expansion in $z$ to converge within our precision in roughly four or fewer
terms.

The form factors can be expressed in $z$ in various ways [107, 108]. We follow
Bourrely, Caprini, and Lellouch [108] as,

$\displaystyle f_{0}(z)$
$\displaystyle=\frac{1}{1-q^{2}(z)/M^{2}_{0^{+}}}\sum_{n=0}^{M-1}b_{m}z^{m},$
(54) $\displaystyle f_{+}(z)$
$\displaystyle=\frac{1}{1-q^{2}(z)/M^{2}_{1^{-}}}\sum_{n=0}^{N-1}a_{n}\left(z^{n}-\frac{n}{N}(-1)^{n-N}z^{N}\right).$
(55)

In these expressions, $M_{J^{P}}$ refers to a possible sub-threshold
($M_{J^{P}}^{2}<t_{+}$) pole, which requires explicit removal. For the scalar
or vector form factors, the pole corresponds to any sub-threshold particle
with quantum numbers $J^{P}=0^{+}$ or $1^{-}$, respectively, corresponding to
the helicity of the virtual $W$ boson. Such poles are present for the decays
$D\to K$ and $D_{s}\to K$ with $J^{P}=1^{-}$. No sub-threshold poles are
present for $D\to\pi$, but the fits are more stable if the nearby poles are
nevertheless included, as shown previously for $D\to\pi$ [57].

Table 12:  Pole masses and cut positions used in Eqs. 54 and 55. The closest pole and the start of the cut are the same for both $D\to\pi$ and $D_{s}\to K$, since they both involve the same $c\to d$ quark-level transition.  Decay | $\sqrt{t_{+}}$ | pole | $J^{P}=1^{-}$ | pole | $J^{P}=0^{+}$
---|---|---|---|---|---
$D\to K$ | $M_{D}+M_{K}$ | $D_{s}^{*}$ | $2112.2(4)\leavevmode\nobreak\ \text{MeV}$ | $D_{s0}^{*}$ | $2317.8(5)\leavevmode\nobreak\ \text{MeV}$
$D\to\pi$, $D_{s}\to K$ | $M_{D}+M_{\pi}$ | $D^{*}$ | $2006.85(05)\leavevmode\nobreak\ \text{MeV}$ | $D_{0}^{*}$ | $2300(15)\leavevmode\nobreak\ \text{MeV}$

Since the input data from the continuum results in Sec. V.3 spans the full
kinematic range of the decay, the $z$ expansion amounts to a convenient change
of variables. To carry out this procedure, we evaluate each form factor at
four evenly spaced points spread throughout the physical $q^{2}$-region:
$[0.1,0.37,0.63,0.9]\times q^{2}_{\rm max}$. We then perform a joint
correlated fit of these synthetic data to Eqs. 54 and 55, imposing the
kinematic constraint $f_{+}(0)=f_{0}(0)$ by taking a common coefficient for
$n=0$: $a_{0}\equiv b_{0}$. The pole masses entering Eqs. 54 and 55 are given
in Table 12. Table 13 reports the correlated posterior values for $a_{n}$ and
$b_{m}$ emerging from the preferred fits for the three decays analyzed. The
preferred fits have $N=M=4$ terms for all three decays. As shown in Fig. 16,
the posteriors for the coefficients stabilize with these choices. In all
cases, statistical uncertainties in the fit parameters are determined via
bootstrap resampling with $500$ draws. These bootstrap fits also furnish
estimates of the $21\times 21$ correlation matrix associated with the full set
of form factors ($f_{+}$ and $f_{0}$ for all three decays). The block-diagonal
correlations for each decay are also given in Table 13, while the full
correlation matrix is given in the supplementary material. The results for
$f_{+}(q^{2})$ and $f_{0}(q^{2})$ coming directly from the chiral-continuum
fits (before applying the $z$ expansion) are compared with those from the $z$
expansion in Fig. 17.

An alternative, common form of the $z$ expansion uses [107]

$\displaystyle f_{+}(q^{2})$
$\displaystyle=\frac{1}{P(q^{2})\phi(q^{2})}\sum_{n=0}^{N-1}a_{n}z^{n},$ (56)

with $P(q^{2})=1$ for $D\to\pi$ or $z(q^{2},M_{D_{s}^{*}}^{2})$ for $D\to K$
and outer function $\phi(q^{2})$ given by

$\displaystyle\phi(q^{2})$
$\displaystyle=\sqrt{\frac{\pi}{3}}m_{c}\left(\frac{z(q^{2},0)}{-q^{2}}\right)^{5/2}\left(\frac{z(q^{2},t_{0})}{t_{0}-q^{2}}\right)^{-1/2}\left(\frac{z(q^{2},t_{-})}{t_{-}-q^{2}}\right)^{-3/4}\left(\frac{t_{+}-q^{2}}{(t_{+}--t_{0})^{1/4}}\right),$
(57)

where $t_{0}=t_{+}(1-\sqrt{1-t_{-}/t_{+}})$ and
$m_{c}=1.25\leavevmode\nobreak\ \text{GeV}$. For comparison with the
experimental determination of the shapes in Sec. VII.3 below, we use Eq. 56
together with the refitting procedure described in Ref. [62].

Figure 16:  Stability analysis for the fit parameters appearing in the $z$
expansion for the decays $D\to\pi$ and $D\to K$. The preferred fit uses
$N=M=4$ terms (i.e., up to and including $a_{3}$ and $b_{3}$), at which point
the coefficients’ central values and errors have stabilized and higher-order
terms are expected to contribute negligibly at the current level of precision.
Results for $D_{s}\to K$ are qualitatively similar.  Table 13:  Correlated
posterior values for $a_{n}$ and $b_{m}$ for the coefficients of the $z$
expansion for the decays $D\to\pi$, $D\to K$, and $D_{s}\to K$. The
simultaneous fit to Eqs. 55 and 54 constrains $a_{0}\equiv b_{0}$. The pole
masses used in the fits are given in Table 12. The full correlation matrix is
given in the supplementary material.

$D\to\pi$ | $a_{0}\equiv b_{0}$ | $a_{1}$ | $a_{2}$ | $a_{3}$ | $b_{1}$ | $b_{2}$ | $b_{3}$
---|---|---|---|---|---|---|---
| 0.6300(51) | -0.610(99) | -0.20(30) | 0.30(19) | 0.330(51) | -0.31(25) | -1.90(39)
| 1.0000 | 0.5670 | 0.5189 | -0.2018 | 0.7547 | 0.3473 | 0.0861
|  | 1.0000 | 0.8912 | -0.2826 | 0.5148 | 0.2529 | 0.0747
|  |  | 1.0000 | -0.1482 | 0.5082 | 0.2782 | 0.1162
|  |  |  | 1.0000 | -0.1728 | -0.0496 | 0.0354
|  |  |  |  | 1.0000 | 0.8277 | 0.6066
|  |  |  |  |  | 1.0000 | 0.9442
|  |  |  |  |  |  | 1.0000

$D\to K$ | $a_{0}\equiv b_{0}$ | $a_{1}$ | $a_{2}$ | $a_{3}$ | $b_{1}$ | $b_{2}$ | $b_{3}$
---|---|---|---|---|---|---|---
| 0.7452(31) | -0.948(97) | 0.14(40) | 0.07(12) | 0.776(62) | 0.14(34) | 0.03(13)
| 1.0000 | -0.0332 | 0.0747 | -0.0201 | 0.7753 | 0.4920 | -0.0189
|  | 1.0000 | 0.3272 | -0.1586 | -0.0909 | -0.1090 | 0.0420
|  |  | 1.0000 | -0.7543 | 0.2071 | 0.2565 | 0.1457
|  |  |  | 1.0000 | -0.0594 | -0.1119 | -0.2259
|  |  |  |  | 1.0000 | 0.9087 | 0.1012
|  |  |  |  |  | 1.0000 | 0.2126
|  |  |  |  |  |  | 1.0000

$D_{s}\to K$ | $a_{0}\equiv b_{0}$ | $a_{1}$ | $a_{2}$ | $a_{3}$ | $b_{1}$ | $b_{2}$ | $b_{3}$
---|---|---|---|---|---|---|---
| 0.6307(20) | -0.562(65) | -0.19(20) | 0.33(29) | 0.347(27) | 0.44(18) | -0.21(43)
| 1.0000 | 0.1825 | 0.2612 | -0.0266 | 0.8467 | 0.5197 | 0.0973
|  | 1.0000 | 0.9274 | -0.2432 | 0.1899 | 0.1915 | 0.1180
|  |  | 1.0000 | -0.0514 | 0.3243 | 0.3065 | 0.1764
|  |  |  | 1.0000 | -0.0551 | -0.1580 | -0.2098
|  |  |  |  | 1.0000 | 0.8344 | 0.4260
|  |  |  |  |  | 1.0000 | 0.8442
|  |  |  |  |  |  | 1.0000

Figure 17:  Final results for $f_{+}$ and $f_{0}$ for the decays $D\to\pi$, $D\to K$, and $D_{s}\to K$ in the continuum limit and at the physical point before and after fitting to $z$ expansion. The solid curves show the results after the chiral-continuum fit, while the hatched curves show the result of the $z$ expansion. The left column shows the product $(1-q^{2}/M_{J^{P}}^{2})f_{0,+}$ as a function of $z$, while the right column shows the form factors versus $q^{2}$. Table 14:  Final results for $f_{+}(0)=f_{0}(0)$, $f_{+}(q^{2}_{\rm max})$, and $f_{0}(q^{2}_{\rm max})$ for the decays $D\to\pi$, $D\to K$, and $D_{s}\to K$, together with comparisons with existing $N_{f}=2+1+1$ results in the literature from HPQCD [62, 63] and ETMC [57]. The results of the present work, denoted “FNAL-MILC”, are all given at the physical point and in the continuum limit in isospin-symmetric QCD. Included in these results are all systematic errors discussed in Sec. VI and summarized in Table 15. Not included are additional systematic uncertainties associated with QED, isospin, and electroweak corrections (these effects are estimated in Sec. VII.2). The different groups use slightly different conventions to define the isospin-symmetric point. Shifts from these differences are expected to be small. Figure 24 suggests that the largest differences, perhaps amounting to a few percent, will be present near $q^{2}_{\rm max}$. process | collaboration | $f_{0}(0)$ | $f_{0}(q^{2}_{\rm max})$ | $f_{+}(q^{2}_{\rm max})$
---|---|---|---|---
$D\to\pi$ | FNAL/MILC | 0.6300(51) | 1.2783(61) | 3.119(57)
$D\to\pi$ | ETMC 17 | 0.612(35) | 1.134(49) | 2.130(96)
$D\to K$ | FNAL/MILC | 0.7452(31) | 1.0240(21) | 1.451(17)
$D\to K$ | HPQCD 22 | 0.7441(40) | 1.0136(36) | 1.462(16)
$D\to K$ | HPQCD 21 | 0.7380(40) | 1.0158(41) | 1.465(20)
$D\to K$ | ETMC 17 | 0.765(31) | 0.979(19) | 1.336(54)
$D_{s}\to K$ | FNAL/MILC | 0.6307(20) | 0.9843(18) | 1.576(13)

###  V.5 Spectator dependence

From the hadronic perspective, the decay channels $D\to\pi$ and $D_{s}\to K$
are quite similar, differing only by the mass of the valence spectator quark.
As illustrated in Fig. 18, we find that the vector and scalar form factors for
these two transitions agree with with each other at the level of $\lesssim
2\%$ throughout the full kinematic range of the $D_{s}\to K$ decay. The first
experimental measurement of the decay $D_{s}\to K$ by BES III [109] confirms
this picture within experimental uncertainties while old, unpublished results
by the HPQCD collaboration [52, 53] are also consistent with our findings.

Figure 18:  Comparison of the vector and scalar form factors between the
decays $D\to\pi$ and $D_{s}\to K$, which differ only by the mass of the
valence spectator quark. The form factors agree at the level of $\lesssim 2\%$
throughout the full kinematic range of the $D_{s}\to K$ decay. The long dashed
lines extending to $q^{2}\approx 3\leavevmode\nobreak\ \text{GeV}^{2}$
correspond to $D\to\pi$, while the shorter solid lines correspond to $D_{s}\to
K$.

###  V.6 Comparison with existing results in the literature

The form factors under consideration have been computed previously using
lattice QCD with $N_{f}=2+1+1$ flavors of dynamical fermions by ETMC [57, 58]
(for both $D\to\pi$ and $D\to K$) and by HPQCD (for $D\to K$) [62, 63]. The
more recent HPQCD calculation [63] includes the same set of $D\to K$
correlators as the earlier one [62], but they are analyzed together with
tensor-current three-point functions, data for heavier-than-charm quark
masses, and $D_{s}\to\eta_{s}\ell\nu$ form factor data [110]. Both the
correlator fits and the description of the heavy-quark-mass dependence and
discretization effects are thus different. Our $D\to\pi$ results for the form
factors and the semimuonic differential decay rate are compared with those of
ETMC in Fig. 19. At large $q^{2}$, our results for the form factors are
significantly larger than those in Ref. [57]. Due to phase-space suppression,
the difference is less visibly pronounced in the differential decay rate
$d\Gamma/dq^{2}$. In the low $q^{2}$ region, which is most relevant for
extractions of $|V_{cd}|$, good agreement is observed at the level of $\approx
1\sigma$. Similarly, our $D\to K$ results are compared with those of ETMC and
HPQCD in Fig. 20. Mild tension, at the level of $\approx 2\sigma$, is observed
between our results and ETMC. Good agreement with HPQCD is observed throughout
the kinematic range. Our results for $f_{+}(0)=f_{0}(0)$, $f_{0}(q^{2}_{\rm
max})$, and $f_{+}(q^{2}_{\rm max})$ are summarized in Table 14 alongside the
published results of Refs. [57, 58, 62, 63].

Figure 19:  Comparison of our results for the $D\to\pi$ form factors and
semimuonic differential decay rate
$(d\Gamma/dq^{2})(24\pi^{3}/G_{F}^{2})/|V_{cd}|^{2}$ with published results
from ETMC [57]. No QED or electroweak corrections [cf. $\eta_{\rm EW}$ in Eq.
2] or errors have been included. To account for differences in defining the
physical isospin-symmetric point, the errors in our curves have been inflated
with an estimate of SIB effects; see Sec. VII.2 below.

Figure 20:  Comparison of our results for the $D\to K$ form factors and
semimunoic differential decay rate
$(d\Gamma/dq^{2})(24\pi^{3}/G_{F}^{2})/|V_{cs}|^{2}$ with published results
from ETMC [57] and HPQCD [62]. No QED or electroweak corrections (cf.
$\eta_{\rm EW}$ in Eq. 2) or errors have been included. To account for
differences in defining the physical isospin-symmetric point, the errors in
our curves have been inflated with an estimate of SIB effects; see Sec. VII.2
below.

##  VI Systematic Error Analysis

The fits to the $z$ expansion described in Sec. V.4 and given in Table 13
provide our final results for the pure-QCD form factors at the physical point
in isospin-symmetric QCD. In this section, we examine and quantify the various
statistical and systematic uncertainties contributing to the calculations. The
complete final error budget is summarized in Table 15 for $f_{+}(0)=f_{0}(0)$,
$f_{+}(q^{2}_{\rm max})$, and $f_{0}(q^{2}_{\rm max})$ for all decay modes. As
discussed in Secs. VI.3 and VI.4, the very small corrections for the leading
finite-volume shifts ($\lesssim 0.01\%$) and the effect of nonequilibrated
topological charge (relevant for $a\approx 0.042\leavevmode\nobreak\
\text{fm}$ only) have been applied to the form factors prior to fitting and
thus are not included as separate errors. Systematic errors associated with
isospin breaking effects and QED corrections, which are external to our
calculation in isospin-symmetric QCD but necessary for comparison with
experimental results, are discussed in Sec. VII.2.

Table 15:  Complete statistical and systematic error budget for the vector and scalar form factors at $q^{2}=0$ and $q^{2}_{\rm max}$ for the decays $D\to\pi$, $D\to K$, and $D_{s}\to K$. All values are given in percent. The breakdown of the chiral-continuum fit errors is discussed in Sec. VI.2. Corrections for finite-volume and topological-charge effects, discussed in Secs. VI.3 and VI.4, are applied prior to the chiral-continuum fit and are negligibly small ($<0.01\%$). Experimental uncertainties on the meson masses are also negligible at our current level of precision. Decay | $D\to\pi$ |  |  | $D\to K$ |  |  | $D_{s}\to K$ |  |
---|---|---|---|---|---|---|---|---|---
Source | $f_{+}(0)$ | $f_{+}(q^{2}_{\rm max})$ | $f_{0}(q^{2}_{\rm max})$ | $f_{+}(0)$ | $f_{+}(q^{2}_{\rm max})$ | $f_{0}(q^{2}_{\rm max})$ | $f_{+}(0)$ | $f_{+}(q^{2}_{\rm max})$ | $f_{0}(q^{2}_{\rm max})$
Statistics $f_{\perp}$ | 0.21 | 1.46 | 0.01 | 0.07 | 0.95 | 0.01 | 0.07 | 0.73 | 0.02
Statistics $f_{0}$ | 0.70 | 0.39 | 0.40 | 0.39 | 0.36 | 0.22 | 0.29 | 0.18 | 0.12
Ctm. $w_{0}$ | 0.31 | 0.09 | 0.15 | 0.24 | 0.14 | 0.17 | 0.32 | 0.27 | 0.27
$\chi$EFT $f_{\perp}$ | 0.12 | 0.46 | 0.01 | 0.03 | 0.50 | 0.00 | 0.03 | 0.39 | 0.01
$\chi$EFT $f_{0}$ | 0.24 | 0.13 | 0.13 | 0.11 | 0.14 | 0.05 | 0.05 | 0.03 | 0.02
Discr. $f_{\perp}$ | 0.08 | 0.73 | 0.00 | 0.06 | 0.41 | 0.01 | 0.03 | 0.41 | 0.01
Discr. $f_{0}$ | 0.05 | 0.02 | 0.07 | 0.02 | 0.02 | 0.02 | 0.00 | 0.00 | 0.01
$f_{\pi}^{\rm PDG}$ | 0.16 | 0.10 | 0.13 | 0.06 | 0.04 | 0.03 | 0.12 | 0.04 | 0.09
Total error | 0.87 | 1.84 | 0.48 | 0.49 | 1.29 | 0.28 | 0.46 | 0.99 | 0.30

###  VI.1 Chiral-continuum fits: stability analysis

The results in Sec. V are the product of several choices. In this section, we
examine the stability of the results under reasonable variations to these
choices for the fiducial point $q^{2}=0$. First, the model for the EFT is
varied. The staggered chiral logarithms are replaced with their continuum
counterparts, setting the known taste splittings to zero by hand. Another
alternative is simply dropping the chiral logarithms $\delta f_{\rm P,logs}$
in Eq. 32. This variation is reasonable, since the ensembles with physical-
mass pions reduce the approach to the physical point from an extrapolation to
an interpolation. The final EFT variation consists of augmenting the analytic
terms in Eq. 32 to include all the N3LO terms (i.e., terms cubic in the
$\chi_{\ell}$, $\chi_{H}$, and $\chi_{E}$). Second, we consider variations to
the model for discretization effects as given in Eqs. 47 and 48. Third, the
widths of our Bayesian priors are increased, and the fits are rerun. In one
variation, the widths of the priors for the coefficients of the leading-order
analytic terms ($c_{l}$, $c_{E}$, and $c_{H}$) are increased by a factor of
ten. In another variation, the widths of all the priors are increased by a
factor of two. Fourth, the choice of the linear shrinkage parameter is tested
by fits varying it by a factor of 2 from its fiducial value ($\lambda=0.1$).
Finally, the choice of data used in the fits is varied, rerunning after
dropping the coarsest ensemble ($a\approx 0.12\leavevmode\nobreak\ \text{fm}$)
and after dropping the finest ensemble ($a\approx 0.042\leavevmode\nobreak\
\text{fm}$).

As Fig. 21 shows, for $D\to\pi$, that all variations are statistically
consistent with the preferred fit at the level of one standard deviation.
Stability plots for $D\to K$ and $D_{s}\to K$ are similar and given in Figs.
47 and 48 in Sec. D.2.

Figure 21: Stability of the $D\to\pi$ form factors $f_{\perp/\parallel/0}$ at
$q^{2}=0$ under variations to the EFT model, the model for discretization
effects, to the choice of data included in the fit, and other analysis choices
as described in the main body. The central values have been normalized by the
central value of preferred fit in green. All variations are statistically
consistent with the preferred fit, highlighted by the green band in each
panel. The statistical significance of the fits is indicated by the marker
size, with larger points denoting better fits.

The discussion in Sec. V.3 demonstrates good agreement for the physical form
factors constructed in different ways, while the discussion above shows that
alternative discretization models, as well as continuum-$\chi$PT fit functions
(without taste splittings in the chiral logarithms), give consistent results.

###  VI.2 Chiral-continuum fits: error breakdown

The form-factor results coming out of the chiral-continuum fits contain
several sources of uncertainty: statistical errors in the form factor on each
ensemble (the correlated uncertainty from the bare form factors and
renormalization constants), scale-setting errors coming from the continuum
value of $w_{0}$, choices in the fit function and chiral interpolation,
discretization effects, and errors in the input parameters (physical meson
masses and $f_{\pi}$ in Table 10). The different sources of error are
entangled in the total fit uncertainty; in particular, the fit function,
chiral interpolation and discretization errors are rather difficult to
separate unambiguously. Nevertheless, an estimate of each error can be
obtained using the package gvar [111] following the methodology described in
Ref. [112]. The discretization error is defined to be the error coming from
the parametric uncertainty in $\delta f_{\rm artifacts}^{(a^{2}+h^{2})}$ from
$c_{a^{2}}$ and $c_{h^{2}}$. The combined uncertainty from all other fit
parameters in Eq. 32 is defined to be the error in the fit function and chiral
interpolation. This error includes the uncertainty from the $DD^{*}\pi$
coupling, $g$, which turns out to have a small influence on the final results.
The experimentally measured values of the meson masses also contribute
negligibly to the total error.

Numerical results for the error breakdown are shown in Table 15 for $q^{2}=0$
and $q^{2}_{\rm max}$, and Figs. 22, 49 and 50 show the error budgets through
the full kinematic range for $D\to\pi$, $D\to K$, and $D_{s}\to K$,
respectively, after fits to the $z$ expansion. The colored curves sum in
quadrature to give the total error in black. Not shown are contributions from
uncertainties less than $0.01\%$; this includes the experimental values for
the input meson masses. Since the lattice data span the full kinematic range
in $q^{2}$, errors from the $z$ expansion are also negligible.

Several important qualitative features are evident in the error budgets. For
all three decays, $f_{+}$ has the largest errors near $q^{2}_{\rm max}$, since
this kinematic region involves an extrapolation $(\bm{p}\to\bm{0})$. Second,
because the $z$-expansion analysis uses a correlated joint fit to $f_{0}$ and
$f_{+}$, the final errors in each case include contributions from statistical
uncertainties in both $f_{0}$ and $f_{+}$. Third, because the $f_{\perp}$ term
vanishes in Eq. 9 at $q^{2}=0$, the contributions from statistical errors in
$f_{\perp}$ decrease for small $q^{2}$. Fourth, although the form factors are
dimensionless, the scale-setting uncertainty is significant and tends to
decrease for large $q^{2}$. At the physical point, the scale-setting
uncertainties vanish identically for the chiral logarithms and analytic terms.
The full uncertainty comes from the leading-order term in Eq. 32: since the
posterior values for $c_{0}$ and $\Delta_{xy,\rm P}$ are both implicitly in
intermediate units of $w_{0}$, so must the energy be. The associated scale-
setting uncertainty thus decreases when the energy is small.

The error budgets for $D\to K$ and $D_{s}\to K$ are qualitatively similar as
shown in Figs. 49 and 50. Over the whole kinematic range, statistics is the
dominant source of error for all three channels, except for $f_{0}^{D_{s}\to
K}$ near $q^{2}_{\rm max}$, where the scale-setting uncertainty dominates.

Figure 22:  Final error budget for the form factors $f^{D\to\pi}_{+}$ and
$f^{D\to\pi}_{0}$ after the fit to the $z$ expansion. Contributions less than
$0.01\%$ are not shown.

###  VI.3 Finite-volume corrections

In principle, the finite volume of our simulations is a systematic effect
influencing the results for the form factors. Within chiral perturbation
theory, the leading corrections amount to replacing loop integrals by discrete
sums [113, 114]. The basic infinite-volume loop integral appearing in the
present analysis is

$i\mu^{\epsilon}\int\frac{d^{4-\epsilon}q}{(2\pi)^{4-\epsilon}}\frac{1}{q^{2}-M^{2}}=\frac{1}{16\pi^{2}}I_{1}(M),$
(58)

with $I_{1}(M)=M^{2}\ln(M^{2}/\Lambda^{2})$ as in Eq. 39. In a finite volume,
this integral becomes the discrete sum

$\mathcal{I}_{1}(M)=\frac{1}{L^{3}}\sum_{\bm{q}}\int\frac{dq^{0}}{2\pi}\frac{1}{q^{2}-m^{2}+i\epsilon}\equiv
I_{1}(m)+I_{1}^{\rm FV}(m),$ (59)

where $I_{1}^{\rm FV}(M)$ is the finite-volume correction that vanishes
exponentially for large volumes. The correction has the explicit form

$I_{1}^{\rm FV}(M)=\frac{1}{4\pi^{2}}M^{2}\sum_{|\bm{n}|\neq
0}\frac{K_{1}(nML)}{nML},$ (60)

with the sum running over all nonzero lattice vectors
$\bm{n}\in\mathbb{Z}^{3}$ in the finite volume, and where $K_{1}$ is a
modified Bessel function of the second kind. As described in Sec. V, the
effect of this correction has already been included explicitly in our fits to
Eq. 32. To quantify the overall size of the finite-volume effect, it is useful
to compute the dimensionless ratio:

$\frac{I_{1}^{\rm
FV}(m)}{I_{1}(m)}=\frac{4}{\ln\left(M^{2}/\Lambda^{2}\right)}\sum_{|\bm{n}|\neq
0}\frac{K_{1}(nML)}{nML}.$ (61)

As shown in Table 16, the finite-volume corrections amount to $\lesssim 2\%$
shifts in $I_{1}(m)$. In the chiral-continuum fits to Eq. 32, the overall
contribution from the chiral logarithms enter at the level of a few percent.
The total size of finite-volume corrections to the form factors may be
estimated to be at the few permyriad level, $\mathcal{O}(0.01)\%$. Since the
leading correction to the chiral logarithm has already been included in our
fits to Eq. 32, and since the effect is so small, we do not include any
additional error for residual finite-volume effects in our final systematic
error budget.

Table 16: Finite-volume corrections to the chiral logarithm $I_{1}(M_{\pi})$ for the ensembles given in Table 1. $\approx a$ [fm] | $m_{l}/m_{s}$ | $L/a$ | $M_{\pi}L$ |  $I_{\rm FV}(M_{\pi})/I_{1}(M_{\pi})$ [%]
---|---|---|---|---
0.120 | 1/27 | 48 | 3.9 | 1.32
0.088 | 1/10 | 48 | 4.7 | 0.65
0.088 | 1/27 | 64 | 3.7 | 2.06
0.057 | 1/5 | 48 | 4.5 | 1.31
0.057 | 1/10 | 64 | 4.3 | 1.25
0.057 | 1/27 | 96 | 3.7 | 1.91
0.042 | 1/5 | 64 | 4.3 | 1.70

###  VI.4 Nonequilibrated topological charge

Efficiently sampling regions with different topological charges $Q$ in
lattice-QCD simulations becomes slow in standard algorithms, which use a
continuous updating procedure for the gauge fields. Brower et al. [115]
realized that chiral perturbation theory can be used to study the
$Q$-dependence of observables, and they showed how to extract physical results
from numerical data at fixed topology. Their calculations confirmed the
theoretical expectation that, due to locality and cluster decomposition, the
effects from fixed topology should be suppressed for large volumes. Subsequent
calculations by Bernard and Toussaint [116] extended these ideas to heavy-
light decay constants and meson masses in the context of heavy-meson chiral-
perturbation theory. The analysis was extended to light form factors in Ref.
[27].

Following those works, we account for the effect of the difference between the
correct $\left\langle Q^{2}\right\rangle$ and the simulation $\left\langle
Q^{2}\right\rangle_{\rm sample}$ in the extraction of heavy-light form factors
by applying a correction factor $\Delta_{Q}f_{\rm P}$, independent on $q^{2}$,
valid for all form factors considered in this work, and given by

$\displaystyle\Delta_{Q}f_{\rm P}$ $\displaystyle\equiv f_{\rm
P,corrected}-f_{\rm P,sample},$ (62)
$\displaystyle=-\frac{1}{2\chi_{T}V}\left.\frac{\partial^{2}f_{\rm
P}}{\partial\theta^{2}}\right|_{\theta=0}\left(1-\frac{\left\langle
Q^{2}\right\rangle_{\rm sample}}{\chi_{T}V}\right),$ (63)

with $f_{\rm P,sample}$ the simulation value of a given form factor at any
value of $q^{2}$, $\theta$ the vacuum angle, $\chi_{T}$ the topological
susceptibility, and $V$ the four-dimensional lattice volume. The second
derivative of the form factors with respect to the vacuum angle is obtained
using LO heavy-light $\chi$PT with $\theta\neq 0$

$\displaystyle\left.\frac{\partial^{2}f_{\rm
P}}{\partial\theta^{2}}\right|_{\theta=0}$
$\displaystyle=-\frac{1}{4}\left(\frac{m_{l}m_{s}}{m_{x}(m_{l}+2m_{s})}\right)^{2}\left.f_{\rm
P}(\theta)\right|_{\theta=0}.$ (64)

where $m_{l,s}$ are the light and strange quark masses respectively, and
$m_{x}$ is the mass of the spectator quark in the transition, i.e., $m_{l}$
for $D\to\pi(K)\ell\nu$, $m_{s}$ for $D_{s}\to K\ell\nu$. The value of
$\left\langle Q^{2}\right\rangle_{\rm sample}$ is understood to be the
measured value from the simulation. For the chiral susceptibility, we take the
prediction from leading-order staggered $\chi$PT [117],

$\chi_{T}=\frac{1}{4}f_{\pi}^{2}M^{2}$ (65)

where $1/M^{2}=2/M_{ll,I}^{2}+1/M_{ss,I}^{2}$ involves the taste-singlet non-
Goldstone states. At leading order, the only change from the familiar result
[118] is the replacement of $M_{\pi}^{2}$ by $M^{2}$. The masses of the taste-
singlet mesons are calculated using Eq. 41. Of the ensembles considered in
this work (cf. Table 1), the effects of nonequilibrated topological charge are
relevant only for the finest ensemble ($a\approx 0.042\leavevmode\nobreak\
\text{fm}$), for which $\left\langle Q^{2}\right\rangle_{\rm sample}=27.59$
[116]. The resulting corrections, $(\Delta_{Q}f_{\rm P})/f_{\rm P}\lesssim
0.0003$, are applied to the form factor data on the $a\approx
0.042\leavevmode\nobreak\ \text{fm}$ ensemble prior to the chiral-continuum
fit in Sec. V. Having accounted for the effect explicitly, and given the
smallness of the correction, no further systematic error is assigned for
nonequilibrated topological charge.

##  VII Phenomenology

The analysis of the preceding sections yields the semileptonic form factors
for $D\to\pi$, $D\to K$, and $D_{s}\to K$ in the idealized case of isospin-
symmetric QCD. For phenomenological applications, we have to consider the
effects of strong isospin and QED, and then combine corrected results with
experimental data. In this section, we first (Sec. VII.1) explore several
options for combining the experimental results with the lattice-QCD form
factors and next (Sec. VII.2) estimate QED and strong isopin-breaking effects.
We are then in a position to determine via Eq. 2 the CKM matrix elements
$|V_{cd}|$ and $|V_{cs}|$ with a full error budget (Sec. VII.3) and to carry
out tests of CKM unitarity (Sec. VII.4). We also compute the Standard Model
predictions for the LFU ratios $R_{\mu/e}$.

###  VII.1 Experimental measurements

The differential decay rates $d\Gamma/dq^{2}$ for semileptonic decays of
$D_{(s)}$ mesons to pseudoscalar light mesons have been measured by FOCUS
(shape only) [119], Belle [120], BaBar [121, 122], CLEO [123], and BES III
[124, 125, 126, 109, 127]. Table 17 summarizes the published measurements
according to decay channel.

Table 17: Summary of published measurements of semileptonic decays of $D$ mesons to pseudoscalar light mesons. FOCUS 2005 [119] obtained shape information only and is omitted.  Decay | Measurements | Notes
---|---|---
$D^{0}\to\pi^{-}e^{+}\nu$ | BaBar 2015 [122] |
| Belle 2006 [120] |  $e^{+}$ and $\mu^{+}$ averaged
| BES III 2015 [124] |
| CLEO 2009 [123] |
$D^{+}\to\pi^{0}e^{+}\nu$ | BES III 2017 [126] |
| CLEO 2009 [123] |
$D^{0}\to\pi^{-}\mu^{+}\nu$ | Belle 2006 [120] |  $e^{+}$ and $\mu^{+}$ averaged
| BES III 2018 [127] |
$D^{+}\to\pi^{0}\mu^{+}\nu$ | BES III 2018 [127] |
$D^{0}\to K^{-}e^{+}\nu$ | BaBar 2007 [121] |
| Belle 2006 [120] |  $e^{+}$ and $\mu^{+}$ averaged
| BES III 2015 [124] |
| CLEO 2009 [123] |
$D^{+}\to\bar{K}^{0}e^{+}\nu$ | BES III 2017 [126] |
| CLEO 2009 [123] |
$D^{0}\to K^{-}\mu^{+}\nu$ | Belle 2006 [120] |  $e^{+}$ and $\mu^{+}$ averaged
| BES III 2019 [128] |
$D^{+}\to\bar{K}^{0}\mu^{+}\nu$ | BES III 2016 [125] | total rate only
$D_{s}^{+}\to K^{0}e^{+}\nu$ | BES III 2019 [109] |

Due to the experimental challenge of reconstructing muons in the final state,
more measurements exist for the electron channels. The only published data
available for the semimuonic final states are from BES III, which measured the
rates for $D\to\pi\mu\nu$ [127] and $D\to K\mu\nu$ [125, 128]. Although Belle
measured both the semielectronic and semimuonic final states [120], numerical
values for the rate were not reported; instead, only values for the product
$|V_{cx}|f_{+}^{D\to\pi/K}(q^{2})$ averaged over the lepton final state are
available [129, 130], without any correlation information. Since both
experimental data and lattice-QCD form factors have now reached a level of
precision where the effects of the scalar form factor (which are proportional
to $m_{\ell}^{2}$) are no longer negligible, as discussed below, we exclude
Belle data from our subsequent analysis.

Besides the experimental difficulties associated with semimuonic final states,
the extraction of the CKM matrix elements $|V_{cd}|$ and $|V_{cs}|$ poses an
additional complication. Contributions to the differential decay rate from the
scalar form factor enter Eq. 2 with a factor of $m_{\ell}^{2}$. As Fig. 23
shows, the scalar form factor is negligible for semielectronic final states
everywhere except the lowest $q^{2}$ bin, where its contribution is roughly
$1\%$. The situation for semimuonic final states is entirely different, where
contributions from $f_{0}$ are roughly $(m_{\mu}/m_{e})^{2}\approx 10^{5}$
times larger and, thus, contribute at the few-percent level throughout the
full kinematic range. Many extractions of $|V_{cd}|$ and $|V_{cs}|$ have
neglected the contributions of the scalar form factor. But, with errors of
$\lesssim 1\%$ both from experiment and from the results of this paper,
determinations of $|V_{cd}|$ and $|V_{cs}|$ require the inclusion of both
terms in Eq. 2.

Figure 23:  Fractional contributions from the scalar and vector form factors
to the differential decay rate $d\Gamma/dq^{2}$ for $D^{0}\to\pi^{-}\ell\nu$
for electronic and muonic final states. Scalar contributions enter Eq. 2 with
a factor of $m_{\ell}^{2}$. In muonic decays, in general, scalar contributions
are thus a factor of $(m_{\mu}/m_{e})^{2}\approx 10^{5}$ larger than in the
corresponding electronic decays.

###  VII.2 Systematic uncertainty from strong isospin effects and QED

The form-factor results reported in Table 13 are computed in isospin-symmetric
QCD, i.e., in simulations with degenerate light quarks of mass
$m_{l}=(m_{u}+m_{d})/2$ in the sea and valence sectors. This theory is
slightly different from nature, which includes corrections from
electromagnetic effects and strong isospin breaking (SIB). An estimate of
these neglected effects is necessary before combination with experimental
data.

Consider first SIB. Isospin violation in the sea may be ignored at the current
level of precision. Because the matrix elements yielding the form factors are
symmetric under exchange of the up and down sea quarks ($m_{u}\leftrightarrow
m_{d}$), the leading contributions from SIB in the sea are of order
$(m_{d}-m_{u})^{2}$. This behavior appears in the $\chi$PT [83] expressions,
showing that sea SIB is smaller than the NNLO terms in the chiral expansion
[32]. To estimate the valence correction, we evaluate the form factors with a
different definition of the physical point, replacing the masses of the
neutral initial and final hadrons that define the physical point, see Table
10, with their charged counterparts, and then computing the fractional shift
$1-(f_{+,0}^{\rm neutral}/f_{+,0}^{\rm charged})$ as a function of $q^{2}$. To
account for this systematic effect, we increase our errors on the form factors
by $\pm(1-f_{+,0}^{\rm neutral}/f_{+,0}^{\rm charged})$, leaving the central
value unchanged. The systematic error profiles are shown as functions of
$q^{2}$ in Fig. 24. Although this treatment of SIB does not distinguish
between SIB in the sea and valence sectors, it is conservative insofar as both
sea and valence effects contribute the variation with the hadron masses.
Guidance from EFT calculations or dedicated simulations with $m_{u}\neq m_{d}$
would be useful to help quantify this effect more precisely. Due to phase-
space suppression at large $q^{2}$, isospin effects will turn out to be a
small (and sometimes neglible) contribution to the systematic error budgets
for quantities of phenomenological interest.

Figure 24:  Systematic error envelopes $1-f_{+,0}^{\rm neutral}/f_{+,0}^{\rm
charged}$ estimating the effects of isospin breaking from the variation in
form factors using an alternative definition of the physical point with the
charged initial and final hadron masses (cf. Table 10). The total errors on
the theoretical prediction for the form factor are increased, leaving the
central value unchanged. Although the systematic uncertainty from SIB
increases with $q^{2}$, its effect on $|V_{cx}|$ and $R_{\mu/e}$ ends up being
small due to phase-space suppression.

Some effects of QED are taken into account in the experimental measurements.
For instance, final-state radiation tends to degrade the momentum resolution,
which can lead to mis-measurement of the positron momentum if background
radiative events (e.g., $D^{0}\to\pi^{-}e^{+}\nu\gamma$) are not handled
correctly. Experimental groups correct for this effect using the Monte Carlo
tool PHOTOS [131, 132]. See Refs. [121, 123] for a discussion.

The long-distance electromagnetic corrections to the semileptonic decays
themselves ($\delta_{\rm EM}$ in Eq. 2) have not been calculated for the
decays $D_{(s)}\to K/\pi\ell\nu$. However, the analogous corrections to the
decay amplitudes for $K\to\pi\ell\nu$ have been computed in the framework of
$\chi$PT [133, 134] and more recently in a hybrid framework combining $\chi$PT
and Sirlin’s representation of SM radiative corrections [41, 42, 43]. The more
recent calculations confirm the older results but with smaller final
uncertainties. The overall picture, substantiated by Table 18, is that final
states with a charged hadron (e.g., $\pi^{-}e^{+}$) tend to have shifts of
$\delta_{\rm EM}\approx 1$–$1.5\%$, while the shifts for final states with a
neutral hadron (e.g., $\pi^{0}e^{+}$) are roughly a factor of 3–4 smaller.
Differences between the decays with an $e^{+}$ or a $\mu^{+}$ in the final
state are around an order of magnitude smaller. Since, as mentioned above, no
similar calculations exist for the decays at hand, we are unable to apply a
concrete correction $\delta_{\rm EM}$. Instead, using the results for
$K\to\pi\ell\nu$ as a rough guide, we include an additional systematic
uncertainty. For extractions of the CKM matrix elements, we add a conservative
error of $\pm 1\%$ to the final value $|V_{cd}|$ or $|V_{cs}|$. In all cases,
the uncertainty is inflated without shifting the central values.

Table 18:  Long-distance electromagnetic corrections for the $K_{\ell 3}$ decay amplitude, taken from Ref. [133, 134, 41, 42, 43]. Since the shifts are computed for the amplitude, the factor of two is necessary for use with the decay rate. Entries correspond to $\frac{1}{2}\delta_{\rm EM}$ in %.  Decay | Cirigliano et al. [133, 134] | Seng et al. [41, 42, 43]
---|---|---
$K^{0}\to\pi^{-}e^{+}\nu$ | $0.50\pm 0.11$ | $0.580\pm 0.016$
$K^{0}\to\pi^{-}\mu^{+}\nu$ | $0.70\pm 0.11$ | $0.77\pm 0.04$
$K^{+}\to\pi^{0}e^{+}\nu$ | $0.05\pm 0.13$ | $0.105\pm 0.024$
$K^{+}\to\pi^{0}\mu^{+}\nu$ | $0.08\pm 0.13$ | $0.25\pm 0.05$

In the analysis below, we also report values for the correlated ratio
$|V_{cd}|/|V_{cs}|$ as well as LFU ratios. A few additional remarks are
necessary concerning the QED uncertainty for these quantities.

Consider first the ratio $|V_{cd}|/|V_{cs}|$. As the results in Table 18 show,
QED corrections for $K^{+}$ decays are $\lesssim 0.25\%$, which suggests
similarly small corrections for $D^{+}$ decays. For the decays of $D^{0}$, the
QED corrections will be dominated by the Coulomb interaction between the
charged final-state particles. The Coulomb shift in the rate is approximately
given by $1+\pi\alpha/\beta$, where
$\beta=\sqrt{1-M_{L}^{2}m_{\ell}^{2}/(p_{L}\cdot p_{\ell})^{2}}$ is relative
velocity between the charged final-state particles [135, 136, 133, 137, 138].
For the decays considered here, the kinematics are such that $\beta\approx 1$.
Therefore, within the uncertainties of our calculation, the Coulomb
corrections for $D^{0}$ decays are essentially constant over the kinematic
range of the decays and would cancel in the ratio. Overall, we take a
conservative $0.5\%$ QED systematic uncertainty for the ratio
$|V_{cd}|/|V_{cs}|$.

Similar considerations apply for the LFU ratios. Again using $K\to\pi\ell\nu$
and Table 18 for guidance, the correlated [43] differences $\delta_{\rm
EM}(\mu)-\delta_{\rm EM}(e)$ are about 0.3–0.4%. As for the ratio of CKM
matrix elements, the Coulomb corrections are expected to introduce a factor,
$(1+\alpha\pi)$, which cancels, within the precision of our calculation, in
the ratio. For the same reasons as above, we thus take a conservative $0.5\%$
QED systematic uncertainty for the LFU ratios.

###  VII.3 CKM matrix elements

Our analysis extracts the CKM matrix elements using two different methods: the
joint $z$-expansion method and the binned method, discussed below.

First, the _joint $z$-expansion method_ fits experimental data for
$d\Gamma/dq^{2}$ together with synthetic data for our lattice-QCD form factors
$f_{+}(q^{2})$ and $f_{0}(q^{2})$. More precisely, the expected model for the
decay rate is given by Eq. 2 using the four-parameter $z$ expansions for both
$f_{+}(q^{2})$ and $f_{0}(q^{2})$ via Eq. 55 and Eq. 54. The CKM matrix
element $|V_{cx}|_{\rm joint}$ is treated as a free parameter in the fit which
serves as a floating relative normalization factor between the experimental
data for the rate and synthetic data for the form factors. The synthetic data
are computed using our results for $f_{+}(q^{2})$ and $f_{0}(q^{2})$ given in
Table 13, evaluated for $q^{2}$ at $[0.1,0.3,0.63,0.9]\times q^{2}_{\rm max}$.
The locations of these points are the same as the synthetic points used in
Sec. V.4. Since this method works directly with the full expression for the
differential decay rate, Eq. 2, it makes no assumptions about the relative
size of the vector and scalar contributions.

Figure 25:  The differential decay rates for $D\to\pi$ (top row) and $D\to K$
(bottom row) in the semielectronic (left) and semimuonic (right) channels. The
blue curves shows the result of evaluating Eq. 2 using our lattice-QCD form
factors, normalized by $|V_{cx}|^{2}_{\rm joint}$. The orange curves show the
result of the joint fit of experimental data and synthetic lattice-QCD data to
the $z$ expansion. The black data points indicate charged-hadron
($K^{-}/\pi^{-}\ell^{+}$) final states, while the green points indicate
experimental measurements for neutral-hadron ($K^{0}/\pi^{0}\ell^{+}$) final
states. In the top row, $D\to\pi$ results come from BaBar [122], CLEO [123],
and BES III [124, 126, 127]. In the bottom row, results come from BaBar [121],
CLEO [123], and BES III [124, 126, 128]. Results from different experiments
have are distinguished by different markers. The points have been slightly
offset horizontally for readability.

Figure 26:  (Left) The differential decay rate for $D_{s}\to K^{0}e^{+}\nu$.
The blue curve shows the result of evaluating Eq. 2 using our lattice-QCD form
factors, normalized by $|V_{cd}|^{2}_{\rm joint}$. The orange curve show the
result of the joint fit of experimental data and synthetic lattice-QCD data to
the $z$ expansion. The data points indicate experimental measurements from BES
III [109]. (Right) The binwise estimate of the CKM matrix element
$\left[|V_{cd}(q_{i}^{2})|\right]_{\rm Binned}$ from the decay $D_{s}\to K$.

Joint $z$-expansion fits have been carried out including all experimental data
for each decay process. The corresponding differential decay rates (orange
curves) are shown, together with the experimental data, in Fig. 25 ($D\to\pi$
and $D\to K$) and Fig. 26 ($D_{s}\to K$). For completeness, the fit posteriors
for the $z$-expansion coefficients are given in Sec. D.5. Measurements from
CLEO were reported including correlations between different decay channels
[123]; these correlations are included in our analysis. The published results
from BES III do not include correlations between different decays (e.g.,
$D^{0}\to\pi^{-}e^{+}\nu$, $D^{+}\to\pi^{0}e^{+}\nu$,