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

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\approx0.12$ fm down to $0.042$ 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 K \ell\nu$. Our results also furnish new Standard Model calculations of the lepton flavor universality ratios $R^{D\to\pi}=0.98671(17)^{\rm QCD}[500]^{\rm QED}$, $R^{D\to K}=0.97606(16)^{\rm QCD}[500]^{\rm QED}$, and $R^{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.


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/ψ [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 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 CKM is a 3 × 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 β decays [17] 1 and |V us | as extracted from semileptonic K 3 ≡ K → π ν 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 2 ≡ K → ν over π 2 ≡ π → ν [24]), the combination |V ud | 2 + |V us | 2 + |V ub | 2 is in tension with unitarity at the 3σ level [25]. There is also a ∼ 3σ 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π + (0) [26][27][28][29] and the ratio of decay constants f K /f π [30][31][32][33][34][35], respectively, are calculated using lattice QCD with uncertainties that have reached the ∼ 0.18% level [36]. Experimental data for the decay widths of K 3 and K 2 /π 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 γ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 Ds (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 [31,32]. 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 1 Recent calculations of the universal electroweak radiative corrections relevant for superallowed β 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].
In this work, we leverage the same theoretical tools that were successfully employed in the calculation of decay constants and the K 3 form factor [27,32,59]: 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 → π ν, D → K ν, and D s → K ν 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 → 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 µ/e , which are defined as the ratio of the branching fractions to muon versus electron final states; see Sec. VII E. 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 → L ν of a heavy pseudoscalar meson H ∈ {D, D s } to a light pseudoscalar meson L ∈ {K, π} is given by where = m 2 /q 2 (with m the lepton mass), 2 q is the momentum transfer, M H and M L are the masses of the heavy initial and light final mesons, and 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 η EW = 1 + (α QED /π) ln(M Z /µ)| µ=M D = 1.009(2) [64], where the error is an estimate of the scale uncertainty from a factor-of-two variation around µ = M D . 3 Long-distance and structure-dependent electromagnetic corrections are described by δ EM . 4 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, 2) . (2.4) In these expressions, k µ , and p µ 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 µ = k µ /M H is the four-velocity of the heavy meson; p µ ⊥ = p µ − (p · v)v µ is the component of the final-state hadron's momentum orthogonal to v; and q µ = k µ − p µ is the momentum transfer. The same form factor f 0 appears in Eqs. (2.3) owing to partial conservation of the vector current (PCVC), namely ∂ µ V µ = (m h − m x )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 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: [No sum is implied in Eq. (2.6).] Using the preceding equations, the vector form factor is given by a linear combination of f ⊥ (q 2 ) and f 0 (q 2 ), which will be useful below.
In momentum space, PCVC implies the following condition for the lattice currents: 5 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 C. 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. (2.9), for example, with f and f ⊥ computed using Eqs. (2.5) and (2.6). 6 These alternative constructions will be used to check for systematic errors in our analysis; see Sec. V C.

III. SIMULATION DETAILS
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 [32,68,69]. Table I and Fig. 1 summarize the ensembles used in this work. Lattice spacings range from a ≈ 0.12 fm down to a ≈ 0.042 fm. An ensemble with physicalmass light quarks appears for all lattice spacings but a ≈ 0.042 fm. For the finer lattice spacings, we also include ensembles with heavier-than-physical light quarks with m l ≈ m s /10 and m l ≈ m s /5. 5 In Minkowski space, the basic momentum-space operator relation reads iq µ V µ (q) = (m h − m l ) S . Equation (2.9), 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. 6 Another expression for f + in terms of f 0 and f exists but involves a delicate numerical cancellation near q 2 max . For this reason it is excluded from the subsequent discussion. Table I. 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 src × N 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 [32,68,69]. 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 ≈ 0.042 fm), as described in Ref. [32]. The number of source times N 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 II.
[MeV] 0. 12  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 II.
Although the primary targets of this work are the dimensionless form factors f 0 and f + , many intermediate quantities (e.g., f and f ⊥ and masses) are dimensionful. Throughout this work, the scale is set on each ensemble using previously calculated values for the gradientflow scale w 0 /a [70,71], also listed in Table I. Details of the intermediate scale-setting scheme in the chiral-continuum analysis are discussed below in Sec. V B.

A. Definitions
To access the matrix elements L| S |H , L| V 0 |H , and L| V i |H , we compute the following two-and three-point correlation functions:   Table I.  Table III. The spin-taste structure of the staggered operators used in this work. Pseudoscalar mesons are created and annihilated using P (π, 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. (4.6).
Operator Spin ⊗ Taste Locality Fig. 2 and the spin-taste structure of the operators in our simulations specified in Table III. The operators used for the scalar current S and temporal vector current V 0 are local, but the spatial vector current 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 V , Eqs. (4.4) and (4.5), 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. (4.6), 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. (4.6) is to use A 0 for the initial-state hadrons (D and D s ) and P for the final-state hadrons (π 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. In Eqs. (4.1)-(4.6), we work in the rest frame of the heavy initial hadron H and compute the recoiling light hadron L with eight different lattice momenta p = 2πn/N s a, where N s is the spatial extent of the lattice, and 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 I and momentum above, we compute the three-point function for several different source-sink separations T , given in Table I. 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 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 I for < l a t e x i t s h a 1 _ b a s e 6 4 = " 3 o A 1 f y n 2 y C j s c v m I H u R K l 0 a x y U s j n a i / J z K m r B 2 q M O 9 U D P t 2 0 R u L / 3 n t F K O z I B M 6 S R E 0 n y 6 K U k k x p u M Y a F c Y 4 C i H O W H c i P x W y v v M M I 5 5 W C V f w y O P l W K 6 m / k s t K O 2 F 2 S + h A i f a M 2 j v h G 9 P j 7 N / Z a F a h y e t x j V X 3 J 7 X P d O 6 o 2 b R u 3 8 Y h Z j k e y R f X J I P H J K z s k V u S Y t w k l K n s k L e X X e n H f n 0 / m a t h a c 2 c w u m Y P z / Q N G a 6 g L < / l a t e x i t > T + t src < l a t e x i t s h a 1 _ b a s e 6 4 = " t c J H l Z p f 9 6 e g f 0 k w I z U y w 0 W 3 6 q 2 E v Z T n C j R y y a z t B H 6 G U c E M C i 5 h X A 5 z C x n j t 6 w P H U c 1 U 2 C j Y n r r m O 4 5 p U e T 1 L j S S K f q 9 4 m C K W u H K n a d i u H A / v Y m 4 n 9 e J 8 f k J C q E z n I E z u z z u 1 u t / w J 6 D T J C h J n Z S 4 6 G x 4 S 2 E 3 5 b k C j V w y a 9 u B n 2 F U M I O C S x h V w 9 x C x v i A 9 a D t q G Y K b F R M L h 3 R X a d 0 a Z I a V x r p R P 0 5 U T B l 7 V D F r l M x 7 N u / 3 l j 8 z 2 v n m B x H h d B Z j q D 5 9 6 I k l x R T O n 6 b d o U B j n L o C O N G u F s p 7 z P D O L p w q q G G e 5 4 q x X S 3 C F l s R + 0 g K k I J C T 7 Q e k B D I 3 p 9 f P j 1 W x G r k Q s v + B v V N L n e b w S H j Y P L g 3 r z p I x x k W y T H b J H A n J E m u S M X J A W 4 Q T I I 3 k i z 9 6 L 9 + q 9 e e / f r R W v n N k i v + B 9 f g E l k 6 J r < / l a t e x i t > J < l a t e x i t s h a 1 _ b a s e 6 4 = " s R u d I 5 O T d z P 1 g E Y 9 2 T 7 a l l p O x b I = " > A A A C I H i c b V B N S 8 N A E N 3 U z 9 Z v P X p Z L I K n k o i o x 6 I X j w q 2 C k 0 o m + 2 k X b q 7 C b s T p c T + A q / 6 B / w 1 3 s S j / h q 3 N Q e 1 P h h 4 v D f D z L w 4 k 8 K i 7 3 9 4 l b n 5 h c W l 5 W p t Z X V t f W N z a 7 t t 0 9 x w a P F U p u Y 2 Z h a k 0 N B C g R J u M w N M x R J u 4 u H 5 x L + 5 A 2 N F q q 9 x l E G k W F + L R H C G T r q 6 6 G 7 W / Y Y / B Z 0 l Q U n q p M R l d 8 u r h r 2 U 5 w o 0 c s m s 7 Q R + h l H B D A o u Y V w L c w s Z 4 0 P W h 4 6 j m i m w U T G 9 d E z 3 n d K j S W p c a a R T 9 e d E w Z S 1 I x W 7 T s V w Y P 9 6 E / E / r 5 N j c h o V Q m c 5 g u b f i 5 J c U k z p 5 G 3 a E w Y 4 y p E j j B v h b q V 8 w A z j 6 M K p h R r u e a o U 0 7 0 i Z L E d d 4 K o C C U k + E D r A Q 2 N 6 A / w 4 d d v R a z G L r z g b 1 S z p H 3 Y C I 4 b R 1 d H 9 e Z Z G e M y 2 S V 7 5 I A E 5 I Q 0 y Q W 5 J C 3 C C Z B H 8 k S e v R f v 1 X v z 3 r 9 b K 1 4 5 s 0 N + w f v 8 A i I b o m k = < / l a t e x i t > H Figure 2. Schematic picture of the three-point functions defined in Eqs. (4.4)-(4.6). The final-state hadron L ∈ {π, K} is created with momentum p at the time t src . An external current J is inserted at time t + t src . The initial-state hadron H ∈ {D, D s } is destroyed at rest at time T + t src . 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. (4.1)-(4.6) take the following forms: where O ∈ {P, A 0 } is the appropriate interpolating operator and |∅ denotes the QCD vacuum state. In the final line, the ground-state term contains the transition matrix elements, n| J |m | n=m=0 ≡ L| J |H , from which one can extract the desired form factors via Eqs. (2.5)-(2.7). For the sake of visualization, certain ratios of correlation functions prove useful: , (4.10) , (4.12) where the bars (e.g.,C P L ) denote the time-slice-averaged correlators defined in Eqs. (A1) and (A3). Up to discretization effects (and renormalization), these ratios asymptotically approach the form factors at large Euclidean times: 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 ∈ (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. C P π (t, p)  Figure 3. Pion two-point correlation functions C P π (t, p = 2πn/N s a) and effective masses on the physical-mass a ≈ 0.088 fm ensemble. To reduce the visual impact of opposite parity states, the effective mass is computed separately for even and odd times using Eq. (4.18) and plotted using circles and triangles, respectively. After folding the data around the midpoint of the lattice, the correlator is defined for t/a ∈ [0, N t /2] = [0, 48]. Because of the form of Eq. (4.18) involves offsets by 2, the effective mass is defined on times t/a ∈ [2, N t /2 − 2] = [2,46].

B. Statistical analysis
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 leastsquares regression to a model function f (a) with parameters a for some data set D. The likelihood function pr(D|a) ∝ exp − 1 2 χ 2 is written in terms of the augmented chi-squared function χ 2 = χ 2 data + χ 2 prior , with whereȳ is a vector with the data means, Σ is the covariance matrix,ã is a vector with the prior values, andΣ is the prior covariance matrix. These expressions are standard [76][77][78].
In the present analysis,ȳ and Σ correspond to the measured means and covariance matrices of the correlation functions. The function f (a) corresponds to the spectral decomposition of Eqs. (4.7)-(4.9), with the energies and matrix elements serving as the parameters a. The choices for the priorsã are discussed below. Instead of using the usual binned-sample covariance matrix in Eq. (4.16), we used an improved estimatorŜ n employing nonlinear shrinkage,   Figure 4. Kaon two-point correlation functions C P K (t, p = 2πn/N s a) and effective masses on the physical-mass a ≈ 0.088 fm ensemble. To reduce the visual impact of opposite parity states, the effective mass is computed separately for even and odd times using Eq. (4.18) 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 ∈ [0, N t /2] = [0, 48]. Because of the form of Eq. (4.18) involves offsets by 2, the effective mass is defined on times t/a ∈ [2, N t /2 − 2] = [2,46]. which corrects for finite-sample-size effects [79]; for technical details, see Appendix C. 7 The general procedure is as follows. First, we compute the binned variances σ. Second, we compute the correlation matrix C n using the full (unbinned) data. Third, we compute the shrinkage estimatorĈ n , taking the effective sample size to be the ratio of the total configurations to the bin size. Finally, we constructŜ n = diag(σ)Ĉ n diag(σ). 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 (Ŝ n ) for the full ensemble [80].
As mentioned above, the statistical analysis begins with two-point functions. Figures 3-5 display representative two-point functions and effective masses for the pion, the kaon and the D meson, respectively, on the physical-mass a ≈ 0.12 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 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 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-5 for other masses and lattice spacings. Each two-point correlator is fit to the corresponding spectral decomposition, Eq. (4.7) or Eq. (4.8), using the choices in Table IV. 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 min . The preferred number of states is roughly the minimal number required to achieve Table IV. Preferred analysis settings for fits of two-point functions to Eq. (4.7) and Eq. (4.8). The same settings are applied uniformly across all ensembles. The larger t min cut for C A 0 H (t) is taken to avoid the excited-state contributions visible in Fig. 5. 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) 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 π 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 A below). statistically significant fits (with, say, χ 2 /DOF 1 or p 0.1 for goodness of fit 8 ), while the cuts on t min and t 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 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 V 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 P π (t, p = 0) on the physical-mass a ≈ 0.12 fm ensemble using fits with different choices for the number of states and t min /a. Consistency with expectations from the effective mass is demonstrated in Fig. 7. Similar studies inform the other choices in Table IV.
The second stage of the analysis combines data from two-point and three-point functions to extract f 0 , f , and f ⊥ . The basic procedure consists of simultaneous correlated fits to the spectral decompositions, Eqs. (4.7)-(4.9), for a particular value of the heavy-quark (1 + 0) States (2 + 0) States (3 + 0) States Figure 6. Stability of the ground-state mass in multi-exponential fits to C P π (t, p = 0) on the physical-mass a ≈ 0.12 fm ensemble. The different colors show the posterior values for the groundstate mass using different numbers of states. The blue points indicating 1-state fits are aligned with the value of t min /a on the horizontal axis. The corresponding results for 2-and 3-state fits at the same t 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 IV, the preferred fit uses 3 states and t min /a = 4 and is indicated by the green star. mass and the current J using the choices for the numbers of states and the fit ranges in Table IV. For instance, a simultaneous correlated fit to C P D (t), C P π (t, p), and C S D→π (t, T, p) furnishes π| S |D . For consistency between the two-point and three-point functions, the fit window for the three-point functions is taken to be t ∈ [t src min , T − t snk min ], where t src min and t snk min are the values of t min associated with the source and sink operators, which in general differ. For example, when the A 0 sink operator is used t src min < t snk min ; see Table IV. 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 2pt ± δM 2pt denote the posterior value of the ground-state energy emerging from a fit to Eq. (4.7) or Eq. (4.8), and let E 2pt (p 2 ) ≡ M 2 2pt + p 2 denote the value of the energy obtained by boosting the central value. Similarly, let A 2pt ± δA 2pt denote the posterior value of the ground-state amplitude from the same fit. 9 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 2pt and A 2pt ), while the prior widths are taken to be ten times the posterior widths (10×δM 2pt and 10×δA 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., δE 2pt /E 2pt = O(α s a 2 p 2 ). Table VI summarizes the choices 9 At large times, a generic two-point function is C(t) = A 2 2pt e M2ptt + e M2pt(T −t) + · · · , so the amplitude A 2pt contains the momentum-dependent relativistic normalization of states in the denominator. . Comparison of the ground-state mass from the preferred fit and the effective mass for C P π (t, p = 0) on the physical-mass a ≈ 0.12 fm ensemble. To reduce the visual impact of opposite parity states, the effective mass is computed separately for even and odd times using Eq. (4.18) and plotted using circles and triangles, respectively.  Right: The behavior of the overlap factor ∅ P (π) π , normalized by the value at zero momentum. Within statistical uncertainties, the overlap factor is constant.
Table VI. Summary of how priors for simultaneous fits to two-and three-point functions incorporate information from the two-point fits. Values for α s are given in Table VIII.

Momentum Quantity
Prior value of these priors. For the excited states, the priors for the energy differences are the same as those in Table V, 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 For the special case of the ground state, the central value of V 00 is estimated from the plateau in ratios following Eq. (D1) below, and the width is taken to be 50%. Once statistically acceptable fits (e.g., χ 2 /DOF 1 or p 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 VI suggest, energies are expected to satisfy the continuum dispersion relation, E 2 = (M 2 +p 2 )(1+O(α s a 2 p 2 )), and overlap factors are expected to be constant, since only point-like interpolators were used for the source and sink operators.  points correspond to the posterior ("best-fit") results, while the dashed lines show the size of the priors for p 2 > 0, as defined in Table VI. 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 → π. The top rows show the approach to the asymptotic (T /a → ∞) plateau region. In the top left panel, data are plotted at fixed momentum p = 2π(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 D→π 0 (t, T, p = 2π(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 D→π 0 (t, T max , p), with each color corresponding to a different momentum. In each case, only the largest source-sink separation T 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 D→π 0 (q 2 ). Additional details, along with similar figures for the decays D → K and D s → K are given in Appendix D 1.

C. Nonperturbative renormalization
Bare matrix elements are renormalized nonperturbatively by imposing the PCVC relation, Eq. (2.9). Figure 11 shows the matrix elements entering this expression, before and after renormalization, for the physical-mass a ≈ 0.12 fm ensemble with the charm-quark mass near is physical value. The black points show the quantity 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.
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. (2.9) 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. 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 → l appears in both D → π and D s → 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 thē lc currents. We thus use the renormalization factors extracted from the D s → K data to renormalize both D s → K and D → π 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 or f ⊥ ) and the momentum. For instance, for the D → π decay on the physical-mass a ≈ 0.12 fm ensemble, the renormalization error in f from Z V 0 is 0.1%. For the same decay and lattice spacing, the renormalization error in f ⊥ from Z V i is around 1%. For comparison, the individual statistical errors in both f and f ⊥ (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. Matrix elements

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 A describes the fit function used in the analysis and its connection to effective field theory (EFT). Section V B presents the results of the fits and describes our definition of the physical point in isospin-symmetric QCD. Section V C presents a cross check on our results by constructing f + and f 0 in different ways. Section V D re-expresses our results in a compact form using the model-independent z expansion. Section V E considers the spectator dependence of the form factors by comparing our results for D → π and D s → K. Finally, Sec. V F compares our form factors with published results in the literature.

A. Description of the chiral-continuum fit function
Together, the bare matrix elements and renormalization factors calculated in Sec. IV furnish the form factors f , f ⊥ , 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 and f ⊥ with SU(2) heavymeson 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) χ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 ∈ { , ⊥, 0, +}: artifacts , where the exponent d P ∈ {1/2, −1/2, 0, 0} for P ∈ { , ⊥, 0, +} (respectively) and c 0 is a dimensionless constant. Although, in principle, the function describing the chiral logarithms, δf logs , depends on the form factor, in hard SU(2) χPT it is the same for all P ∈ { , ⊥, 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 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 ⊥ as f + , and similarly for the pair f and f 0 (cf . Table VII). Equation (5.2) implies that location of the pole in the energy can be written as 3) Table VII. Approximate pole locations ∆ xy,P appearing in the decays D → π, D s → K, and D → K.
where y is the spectator quark and L ∈ {π, K} is the final-state hadron. Values for the ∆ xy,P are collected in Table VII for the decays of interest. The χ n are dimensionless expansion parameters defined according to where f π is the physical pion decay constant and Λ HQET is the scale of heavy-quark effective theory. As in Ref. [32], we take Λ HQET = 800 MeV. The parameters χ l and χ E describe the analytic dependence on the light-quark mass m l (via the leading-order expression M 2 π = 2µm l ) and the final-state hadron energy E, respectively. Their normalization is such that, according to the typical χPT power counting, the corresponding coefficients in the fit function c l and c E are expected to be of order 1. The parameter χ H describes the heavy quark mass mistuning through the difference between the simulated heavy meson mass M sim H (s) and the physical M D 0 or M Ds from Ref. [25]. This term allows for a simultaneous description of results across several different heavy quark masses. Finally, χ s describes the strange-quark mass mistuning.
In hard SU(2) χPT, the chiral logarithms for f and f ⊥ in Eq. (5.1) have the common form [84,85], see Sec. VI C below. Hard SU(2) χPT enjoys the further simplification that nonanalytic selfenergy corrections vanish for all three decays considered here. These expressions [Eq. (5.8) and the self energies] were originally derived for a non-staggered heavy quark, but because the heavy-quark taste is conserved in all-staggered χPT, they hold in the present case too [90]. In heavy-meson χPT, compact expressions are available for f ⊥ and f [91,92], while the corresponding results for f + and f 0 follow as linear combinations. Because of their simple connection to heavy-meson χPT, previous lattice calculations have historically worked primarily in terms of f and f ⊥ . However, since the chiral logarithms have the same functional form for f and f ⊥ in hard SU (2) χPT [see Eq. (5.8)], the same functional form also describes the chiral logarithms for f 0 and f + . In other words, Eq. (5.1) may be used directly for all four form factors, with a 1 − pole for f +,⊥ or a 0 + pole for f 0, .
Following Ref. [83], the arguments of the chiral logarithms involve the masses of mesons with different tastes ξ ∈ {I, P, V, A, T }, that can be expressed as The low-energy constant µ and the taste splittings ∆ ξ have been tabulated for these ensembles in Ref. [27]. At NLO in the chiral expansion, the taste splittings ∆ ξ and the hairpin parameters δ V,A both scale like α 2 s a 2 , so their ratio remains approximately constant as the lattice spacing changes. We follow Ref. [32] and take δ A /∆ = −0.88(09) and δ V /∆ = +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 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 α s (aΛ) 2 or α s (am h ) 2 , where Λ 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 : The quantity x a 2 gives a dimensionless measure of order α s (aΛ) 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 α s and the velocity v ≈ 1/10 of the charm quark within the heavy hadron. Our preferred model thus takes the following simple Ansatz for the discretization effects  Table VII.
Values for α s are given in Table VIII.
To check for truncation effects with high-order discretization effects, we also consider variations: With the HISQ action, discretization effects of order x 4 h are suppressed by v (at the tree level) or by α s . These suppression factors are numerically similar enough that the last term in Eq. (5.17) tests both.

B. Chiral-continuum fits
We perform correlated fits, using the methodology described around Eqs. (4.16) and (4.17), to Eq. (5.1) with δf (a 2 +h 2 ) artifacts in Eq. (5.15) for each of the above form factors for D → π, D → K, and D s → K including all of the ensembles and heavy-quark masses described in Table I. In our preferred fits, the input data for the form factors f ⊥/ /0 are defined using a single matrix element each via Eqs. (2.5)-(2.7). The chiral-continuum fit results using the alternative constructions f alt P of Eqs. (2.10) and (2.11) are considered below in the analysis of systematic effects.
The free parameters varied in the fits are the coefficients c n , the coupling g, and the mass splittings ∆ xy,P . The Bayesian priors for these parameters are given in Table IX. 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 ± 1 for the parameters c n . The dimensionless (reduced) "DD * π" coupling appearing as a coefficient of the chiral logarithms is expected to be g ≈ 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 ± 0.2. Because a broad width is used for ∆ * xy,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 X).
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 χ l and χ E contain factors of the pion decay constant in the denominator. Using the continuum values for f π and w 0 in Table X, we express the denominator as a dimensionless number. For the numerator, we use the measured values of aM π and w 0 /a to construct the dimensionless product. In other words, on each ensemble χ l is constructed as and similarly for χ E and χ 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 → 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 λ. We find that λ = 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 A. 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. (5.1) is defined by setting δf (a 2 +h 2 ) artifacts equal to zero and setting the taste splittings to zero in δf P,logs . The physical point is defined by setting the input meson  [126] 0.48 [128] 1.31 [134] f 0.59 [112] 0.41 [123] 0.88 [128] f ⊥ 0.64 [110] 0.32 [111] 0.66 [113] f + 0.59 [106] 0.29 [109] 0.60 [111] masses equal to their physical values, given in Table X. By construction, all quantities involving χ 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 X are chosen accordingly, following the prescription in Ref. [36]. The systematic uncertainty with the isospin-symmetric approximation is discussed below in Sec. VII B.
The results for the D → π 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 ∈ {0.9, 1.0, 1.1}. In all cases, the nearly coincident data around the physical charm mass (m h ≈ 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 → K and D s → K are quite similar and given in Appendix D 2. Table XI summarizes the fit quality for the preferred fits.
C. Alternative constructions of f + and f 0 Our default construction for f 0 is given by Eq. (2.7) and obtained in the preceding section.
In an analogous way we construct f + from the continuum-limit results for f 0 and f ⊥ in the preceding section following Eq. (2.8). As discussed in Sec. II, the PCVC relation in Eq. (2.9) provides the alternative constructions given in Eqs. (2.10) and (2.11). Additional freedom exists in whether the linear combinations in Eqs. (2.8), (2.10) and (2.11) are taken before or after the chiral-continuum limit. A comparison of the different constructions is given in Fig. 15 for D → π and D s → K (D → 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 ⊥ )+CL(f 0 ) versus CL(f ⊥ + f 0 ), respectively]. In all cases, f and f ⊥ are directly related to vector matrix elements via Eqs. (2.5) and (2.6). 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. (2.7). We construct f + via Eq. (2.8), using results for f ⊥ and f 0 given   Figure 15. Comparison of the form factors coming from different continuum-limit constructions for the decays D → π, and D s → 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 → K.

D. 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 → π and D → 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 ), where t + = (M D + M L ) 2 denotes the start of the multiparticle cut, L ∈ {π, K}, and t 0 ∈ [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   (5.24). The closest pole and the start of the cut are the same for both D → π and D s → K, since they both involve the same c → d quark-level transition. . Below, we take t 0 = 0, so Because −z max ≈ 0.332, 0.190, and 0.192, for D → π, D → K, and D s → 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, In these expressions, M J P refers to a possible sub-threshold (M 2 J P < t + ) pole, which requires explicit removal. For the scalar or vector form factors, the pole corresponds to any subthreshold 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 → K and D s → K with J P = 1 − . No sub-threshold poles are present for D → π, but the fits are more stable if the nearby poles are nevertheless included, as shown previously for D → π [57].
Since the input data from the continuum results in Sec. V C 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] × q 2 max . We then perform a joint correlated fit of these synthetic data to Eqs. (5.23) and (5.24), imposing the kinematic constraint f + (0) = f 0 (0) by taking a common coefficient for n = 0: a 0 ≡ b 0 . The pole masses entering Eqs. (5.23) and (5.24) are given in Table XII. Table XIII 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 × 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 XIII, 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] f with P (q 2 ) = 1 for D → π or z(q 2 , M 2 D * s ) for D → K and outer function φ(q 2 ) given by where t 0 = t + (1 − 1 − t − /t + ) and m c = 1.25 GeV. For comparison with the experimental determination of the shapes in Sec. VII C below, we use Eq. (5.25) together with the refitting procedure described in Ref. [62]. Table XIII. Correlated posterior values for a n and b m for the coefficients of the z expansion for the decays D → π, D → K, and D s → K. The simultaneous fit to Eqs. (5.23) and (5.24) constrains a 0 ≡ b 0 . The pole masses used in the fits are given in Table XII. The full correlation matrix is given in the supplementary material. (51)

E. Spectator dependence
From the hadronic perspective, the decay channels D → π and D s → 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 2% throughout the full kinematic range of the D s → K decay. The first experimental measurement of the decay D s → 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.   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 → π and D → K) and by HPQCD (for D → K) [62,63]. The more recent HPQCD calculation [63] includes the same set of D → K correlators as the earlier one [62], but they are analyzed to- Table XIV. Final results for f + (0) = f 0 (0), f + (q 2 max ), and f 0 (q 2 max ) for the decays D → π, D → K, and D s → 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 XV. Not included are additional systematic uncertainties associated with QED, isospin, and electroweak corrections (these effects are estimated in Sec. VII B). 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 max . process collaboration f 0 (0) D→π ν FNAL-MILC (Present work) ETMC 17 Figure 19. Comparison of our results for the D → π form factors and semimuonic differential decay rate (dΓ/dq 2 )(24π 3 /G 2 F )/|V cd | 2 with published results from ETMC [57]. No QED or electroweak corrections [cf. η EW in Eq. (2.1)] 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 B below.
gether with tensor-current three-point functions, data for heavier-than-charm quark masses, and D s → η s ν 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 → π 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Γ/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 ≈ 1σ. Similarly, our D → K results are compared with those of ETMC and HPQCD in Fig. 20. Mild tension, at the level of ≈ 2σ, 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 max ), and f + (q 2 max ) are summarized in Table XIV

VI. SYSTEMATIC ERROR ANALYSIS
The fits to the z expansion described in Sec. V D and given in Table XIII 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 XV for f + (0) = f 0 (0), f + (q 2 max ), and f 0 (q 2 max ) for all decay modes. As discussed in Secs. VI C and VI D, the very small corrections for the leading finite-volume shifts ( 0.01%) and the effect of nonequilibrated topological charge (relevant for a ≈ 0.042 fm only) have been applied to the form factors prior to fitting and thus are not included as separate errors. . Comparison of our results for the D → K form factors and semimunoic differential decay rate (dΓ/dq 2 )(24π 3 /G 2 F )/|V cs | 2 with published results from ETMC [57] and HPQCD [62]. No QED or electroweak corrections (cf. η EW in Eq. (2.1)) 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 B below. Table XV. Complete statistical and systematic error budget for the vector and scalar form factors at q 2 = 0 and q 2 max for the decays D → π, D → K, and D s → K. All values are given in percent. The breakdown of the chiral-continuum fit errors is discussed in Sec. VI B. Corrections for finitevolume and topological-charge effects, discussed in Secs. VI C and VI D, 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.
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 B.

A. 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 δf P,logs in Eq. (5.1). 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. (5.1) to include all the N 3 LO terms (i.e., terms cubic in the χ , χ H , and χ E ). Second, we consider variations to the model for discretization effects as given in Eqs. (5.16) and (5.17). 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 (λ = 0.1). Finally, the choice of data used in the fits is varied, rerunning after dropping the coarsest ensemble (a ≈ 0.12 fm) and after dropping the finest ensemble (a ≈ 0.042 fm).
As Fig. 21 shows, for D → π, that all variations are statistically consistent with the preferred fit at the level of one standard deviation. Stability plots for D → K and D s → K are similar and given in Figs. 47 and 48 in Appendix D 2.
The discussion in Sec. V C 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-χPT fit functions (without taste splittings in the chiral logarithms), give consistent results.

B. 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 π in Table X). 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 δf (a 2 +h 2 ) artifacts from c a 2 and c h 2 . The combined uncertainty from all other fit parameters in Eq. (5.1) is defined to be the error in the fit function and chiral interpolation. This error includes the uncertainty from the DD * π 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 XV for q 2 = 0 and q 2 max , and Figs. 22, 49 and 50 show the error budgets through the full kinematic range for D → π, D → K, and D s → 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 max , since this kinematic region involves an extrapolation (p → 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 ⊥ term vanishes in Eq. (2.8) at q 2 = 0, the contributions from statistical errors in f ⊥ 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.

C. 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 where I FV 1 (M ) is the finite-volume correction that vanishes exponentially for large volumes. The correction has the explicit form with the sum running over all nonzero lattice vectors n ∈ 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. (5.1). To quantify the overall size of the finite-volume effect, it is useful to compute the dimensionless ratio: As shown in Table XVI, the finite-volume corrections amount to 2% shifts in I 1 (m). In the chiral-continuum fits to Eq. (5.1), 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, O(0.01)%. Since the leading correction to the chiral logarithm has already been included in our fits to Eq. (5.1), 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.

D. 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 Q 2 and the simulation Q 2 sample in the extraction of heavy-light form factors by applying a correction factor ∆ Q f P , independent on q 2 , valid for all form factors considered in this work, and given by with f P,sample the simulation value of a given form factor at any value of q 2 , θ the vacuum angle, χ 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 χPT with θ = 0 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 → π(K) ν, m s for D s → K ν. The value of Q 2 sample is understood to be the measured value from the simulation. For the chiral susceptibility, we take the prediction from leading-order staggered χPT [117], where 1/M 2 = 2/M 2 ll,I + 1/M 2 ss,I involves the taste-singlet non-Goldstone states. At leading order, the only change from the familiar result [118] is the replacement of M 2 π by M 2 . The masses of the taste-singlet mesons are calculated using Eq. (5.10). Of the ensembles considered in this work (cf . Table I), the effects of nonequilibrated topological charge are relevant only for the finest ensemble (a ≈ 0.042 fm), for which Q 2 sample = 27.59 [116]. The resulting corrections, (∆ Q f P )/f P 0.0003, are applied to the form factor data on the a ≈ 0.042 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 → π, D → K, and D s → 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 A) explore several options for combining the experimental results with the lattice-QCD form factors and next (Sec. VII B) estimate QED and strong isopin-breaking effects. We are then in a position to determine via Eq. (2.1) the CKM matrix elements |V cd | and |V cs | with a full error budget (Sec. VII C) and to carry out tests of CKM unitarity (Sec. VII D). We also compute the Standard Model predictions for the LFU ratios R µ/e .

A. Experimental measurements
The differential decay rates dΓ/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 [109,[124][125][126][127]. Table XVII summarizes the published measurements according to decay channel. 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 → πµν [127] and D → Kµν [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→π/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 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.1) with a factor of m 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 µ /m e ) 2 ≈ 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 1% both from experiment and from the  [125] total rate only D + s → K 0 e + ν BES III 2019 [109] results of this paper, determinations of |V cd | and |V cs | require the inclusion of both terms in Eq. (2.1).

B. Systematic uncertainty from strong isospin effects and QED
The form-factor results reported in Table XIII 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 ↔ m d ), the leading contributions from SIB in the sea are of order (m d − m u ) 2 . This behavior appears in the χ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 X,  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 = 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. 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 → π − e + νγ) 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 (δ EM in Eq. (2.1)) have not been calculated for the decays D (s) → K/π ν. However, the analogous corrections to the decay amplitudes for K → π ν have been computed in the framework of χPT [133,134] and more recently in a hybrid framework combining χ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 XVIII, is that final states with a charged hadron (e.g., π − e + ) tend to have shifts of δ EM ≈ 1-1.5%, while the shifts for final states with a neutral hadron (e.g., π 0 e + ) are roughly a factor of 3-4 smaller. Differences between the decays with an e + or a µ + 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 δ EM . Instead, using the results for K → π ν as a rough guide, we include an additional systematic  Table X). 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 µ/e ends up being small due to phase-space suppression. Table XVIII. Long-distance electromagnetic corrections for the K 3 decay amplitude, taken from Ref. [41-43, 133, 134]. Since the shifts are computed for the amplitude, the factor of two is necessary for use with the decay rate. Entries correspond to 1 2 δ EM in %. Decay Cirigliano et al. [133,134] Seng et al. [41][42][43] K 0 → π − e + ν 0.50 ± 0.11 0.580 ± 0.016 K 0 → π − µ + ν 0.70 ± 0.11 0.77 ± 0.04 K + → π 0 e + ν 0.05 ± 0.13 0.105 ± 0.024 K + → π 0 µ + ν 0.08 ± 0.13 0.25 ± 0.05 uncertainty. For extractions of the CKM matrix elements, we add a conservative error of ±1% to the final value |V cd | or |V cs |. In all cases, the uncertainty is inflated without shifting the central values.
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 XVIII show, QED corrections for K + decays are 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 + πα/β, where β = 1 − M 2 L m 2 /(p L · p ) 2 is relative velocity between the charged final-state particles [133,[135][136][137][138]. For the decays considered here, the kinematics are such that β ≈ 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 → π ν and Table XVIII for guidance, the correlated [43] differences δ EM (µ) − δ 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 + απ), 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.

C. 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Γ/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.1) using the four-parameter z expansions for both f + (q 2 ) and f 0 (q 2 ) via Eq. (5.24) and Eq. (5.23). The CKM matrix element |V cx | 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 XIII, evaluated for q 2 at [0.1, 0.3, 0.63, 0.9]×q 2 max . The locations of these points are the same as the synthetic points used in Sec. V D. Since this method works directly with the full expression for the differential decay rate, Eq. (2.1), it makes no assumptions about the relative size of the vector and scalar contributions.
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 → π and D → K) and Fig. 26 (D s → K). For completeness, the fit posteriors for the z-expansion coefficients are given in Appendix 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 → π − e + ν, D + → π 0 e + ν, D 0 → π − µ + ν, and D + → π 0 µ + ν). We have experimented with different models for the missing off-diagonal blocks of the full correlation matrix, ranging from zero correlation to 100% correlation. Our results for |V cd | and |V cs | are extremely insensitive to the precise treatment of these off-diagonal correlations and give statistically indistinguishable results. We therefore report values from our preferred analysis, which uses a simple model for the correlations in which the off-diagonal blocks are taken to be constant, with correlation coefficient equal to the mean of the corresponding diagonal blocks. 10 . Regarding the measurements of D 0 → K − e + ν coming from BaBar [121], our fits drop the largest q 2 bin, since it is constructed by a normalization constraint (one minus the sum of the other bins). Especially 10 We thank the BES III collaboration for providing us with the correlations for the differential rate dΓ/dq 2 for the decays D + s → K 0 e + ν, D 0 → π − µ + ν, and D + → π 0 µ + ν as well as for information and guidance regarding the treatment of off-diagonal correlations between different decays (Lei Li, private communication, 22   The differential decay rates for D → π (top row) and D → K (bottom row) in the semielectronic (left) and semimuonic (right) channels. The blue curves shows the result of evaluating Eq. (2.1) using our lattice-QCD form factors, normalized by |V cx | 2 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 − /π − + ) final states, while the green points indicate experimental measurements for neutral-hadron (K 0 /π 0 + ) final states. In the top row, D → π 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. when fitting dΓ/dq 2 for semimuonic channels, including the scalar form factor is essential to achieving a good description of the data at low q 2 . However, the higher parameters b 1 and b 2 associated with the scalar form factor are constrained entirely by our precise synthetic data for f 0 . The influence of neglecting f 0 is considered below.
The plots in Figs. 25 and 26 also include comparisons with the shape obtained from our form factors, given by the parameters in Table XIII (lattice-QCD-only results), normalized by |V cx | 2 joint . Those correspond to the blue curves in the plots. In addition, fits are also conducted to the experimental data alone, although we do not use these results to extract the CKM matrix elements. All those z-expansion fits, as well as the joint fit, enforce the kinematic identity f + (0) = f 0 (0) by imposing a 0 = b 0 [cf. Eq. (5.23) and Eq. (5.24)]. The best-fist posterior values for those z-expansion fits are also given in Appendix D 5.
Another visualization of the form factors' shapes, which is independent of the overall normalization, comes from comparing the ratios r 1 ≡ a 1 /a 0 and r 2 ≡ a 2 /a 0 of the z-expansion coefficients from Eq. (5.25) after applying the refitting procedure of Ref. [62]. These ratios are displayed for D → π and D → K in Fig. 27. Our results, given by the black ellipses, show good agreement with the experimental shapes. For D → K, we also find good agreement with the lattice QCD calculation from HPQCD [62]. For D → π, we find r 1 = −2.009 (55) and r 2 = 0.14(36), with a correlation of ρ 12 = −0.58. For D → K, we find r 1 = −2.05 (11) and −0.59 (52), with a correlation of ρ 12 = −0.29.
The second method we use to extract CKM matrix elements, the binned method, combines lattice-QCD results with experimental data for the rate dΓ/dq 2 to give a binwise estimate of the CKM matrix element: where the quantity in the denominator is understood to be the binwise average (i.e., integrated over the bin) of the lattice-QCD form factors together with the appropriate kinematic factors appearing in Eq. (2.1), This expression depends on the lepton mass via ≡ m 2 /q 2 as well as the experimentally measured hadron masses M H and M L for each mode (e.g., D 0 and π − or D + and π 0 for D → π). A weighted, correlated average (i.e., a fit to a constant) then gives |V cx | Binwise . The binned method is entirely general and makes no assumptions about the relative size of the vector and scalar contributions. Results for |V cd | from D → π and |V cs | from D → K for each q 2 bin and experiment, as well as the correlated average over bins including only semielectronic (blue lines) or only semimuonic data (red lines), are shown in Fig. 28. The semimuonic results lie roughly 1σ below the semielectronic results for both |V cd | and |V cs |, so below we report the values in each channel as well as the combined results. Those combined extractions, including all leptonic channels, lie between the two bands in Fig. 28, and are statistically consistent with the individual determinations, as shown in Fig. 29. For |V cd | from D s → K, results are shown in Fig. 26. As argued above, with present statistical precision, the presence of the scalar form factor is quantitatively important for the differential rate dΓ/dq 2 , especially for semimuonic channels. Figure 30 shows the effect of dropping the contribution from f 0 for D → Kµν. Values for |V cs | binned are observed to shift by a few percent and, when considered as a function of q 2 , become statistically inconsistent with a constant. Similar few-percent shifts occur for D → πµν. |V cs | |V cs | binned , e + |V cs | binned , µ + BESIII 2019, K − µ + Figure 28. The binwise estimate of the CKM matrix element |V cx (q 2 i )| Binned from the decays D → π (top rows) and D → K (bottom row) in the semielectronic (left) and semimuonic (right) channels. The horizontal bands show the resulting values for |V cx | Binned from correlated fits to a constant in each channel. The result for |V cd | D→πµ + ν Binned (red) lies slightly below |V cd | D→πe + ν Binned (blue). The combined extraction using both channels lies between the two bands and is statistically consistent with each. A comparison of the different extractions of |V cx | is given in Fig. 29. For D → π, experimental data are taken from BaBar [122], CLEO [123], and BES III [124,126,127]. For D → K, experimental data are taken from BaBar [121], CLEO [123], and BES III [124,126,128]. Although all the correlated fits have good quality (χ 2 /DOF ≈ 1, p 0.05), the residuals for D → K are visually larger near q 2 max .  Figure 29.
Determinations of |V cd | and |V cs | using experimental measurements of the decays D → π and D → K. The outer and inner error bars and bands show the results with and without QED uncertainties, respectively.
Because the joint-fit and binned methods explicitly account for (potentially) percent-level contributions from the scalar form factor, they constitute our main extractions for |V cx |. For continuity with previous studies, we also consider the endpoint method, in which |V cx | is defined according to The experimental values are taken from the HFLAV world averages: |V cd |η EW f D→π + (0) = 0.1426 (18), |V cs |η EW f D→K + (0) = 0.7180(33) [44]. The resulting values for [|V cx |] Endpoint are shown in Fig. 29 and given in Table XIX. Although these endpoint results give a statistical precision comparable to our preferred extractions, it's worth emphasizing that our precise values for f + (0) were made possible by leveraging information about the form factor across the full kinematic range of the decays. The final errors can potentially be much larger in a simulation that works directly at the endpoint (q 2 = 0). For example, preliminary work by our collaboration has focused on q 2 ≈ 0 on many of the same ensembles and with comparable statistics [60]. Using the preliminary values of f + (0) from these proceedings gives values for |V cd | and |V cs | with errors that are roughly 2.5 to 3.5 larger than the final errors in the present work.
The results for |V cd | and |V cs | from the different methods described above and for different leptons in the final states are summarized in Fig. 29 and Table XIX. Since, as shown in the plot, the binned and joint-fit extractions give statistically consistent values well within 1σ, we take the joint-fit extractions to define our preferred results: where the first error comes from the experimental differential decay rate uncertainty, the second error comes from our form factor calculation (see Table XV), the third error shows the uncertainty in η EW , and the fourth and fifth from our estimate of SIB and long-distance QED corrections described in Sec. VII B. The errors in these expressions combine in quadrature to give the total errors in Table XIX. Since our preferred extraction of |V cs | includes both e + and µ + final states, the experimental contribution to the error is smaller by roughly a factor of two than in Ref. [62]. We also repeated our analysis separating charged-hadron and neutral-final states (e.g., π − e + versus π 0 e + ). No statistically significant difference was observed within the uncertainties, consistent with what was observed in Ref. [62].
For the first time, our calculation provides a value of |V cd | from D → π for which lattice QCD errors are at the same level as the experimental errors, ∼ 0.5% each. This represents an improvement by roughly a factor of six from the existing state of the art [57,139]. For |V cd | Ds→K , experimental errors dominate and are substantially larger than for D → π. Since the theoretical uncertainty is actually the smallest for D s → K, additional experimental measurements of this channel would be particularly welcome. On the other hand, theoretical error exceeds the experimental error by roughly a factor of two in the extraction of |V cs | from D → K, leaving room for improvements in the theory side. Experimental errors also dominate the CKM extractions from the semimuonic channels, where we have only included recent results from BES-III. Another key ingredient for improved semileptonic extractions of |V cd | and |V cs | would be the calculation of long-distance structure-dependent EM corrections or a more robust estimate of their effect on these decays, since our lack of knowledge of these corrections currently dominates the uncertainty of the most precise determinations.
Our correlated results for |V cd | and |V cs | also yield the ratio, where the correlation coefficient between |V cd | and |V cs |, neglecting QED, is 0.18. As described in Sec. VII B, we have taken a conservative 0.5% systematic uncertainty for QED effects in the ratio.
Using the latest measurements f D + |V cd | and f Ds |V cs | reported by HFLAV [44] and the ratio of decay constants f Ds /f D + computed by our collaboration in a similar set of ensembles and with the same action in Ref. [32], one finds [|V cd |/|V cs |] leptonic = 0.2212 (58), where the error is dominated by the experimental uncertainty. Both values are plotted in Fig. 32 together with previous leptonic [31,32,[140][141][142][143] semileptonic [48,49,57,62,139] determinations combined in averages by FLAG [36] and the result from the PDG global unitarity fit [25] (the global-fit methodologies of CKMfitter [144] and UTfit [145] give very similar results). The leptonic extraction above agrees with our semileptonic result within roughly 2σ, although, as plotted in Fig. 32, leptonic determinations tend to give smaller values of the ratio. The error in our result is more than a factor of two smaller than the leptonic one, with similar uncertainties from lattice QCD and experiment. Results for |V cd |/|V cs | from the PDG global fit assuming unitarity and from the ratio |V us |/|V ud | (see Sec. VII D below for more details) are also shown in Fig. 32 FNAL-MILC (Present work) FLAG Semileptonic HFLAV+FLAG Leptonic PDG Figure 31. Comparison of our preferred determinations of |V cd | D→π and |V cs | D→K (blue bands) with existing results in the literature. The outer and inner error bands show our preferred result with and without QED uncertainties, respectively. The world's first determination |V cd | Ds→K is also given. Results from FLAG are taken from Ref [36]. Results from the PDG appear in Ref. [25]. We emphasize that FLAG uses slightly different conventions for the semileptonic extraction of |V cd(cs) | as we used here; for instance they do not include short-distance electroweak corrections to G F or an error from QED. For the leptonic results, we combine the latest experimental averages reported in HFLAV [44] with the FLAG averages for f D and f Ds [36]. "CKM unitarity" denotes the global fit result reported by the PDG, which includes all available measurements (for all nine matrix elements) imposing three-generation unitarity.

D. Tests of CKM unitarity
Our results for |V cd | and |V cs | enable a test of unitary in the second row of the CKM matrix, including theoretical correlations between |V cd | and |V cs |. Using our preferred extractions in Eq. (7.4) and Eq. (7.6), and |V cb | incl+excl = (40.8 ± 1.4) × 10 −3 from a combined average of inclusive and exclusive semileptonic B-decays [25] 11 yields the following result for the deviation from unitarity in the second row:  uncertainty from QED in our extractions of |V cd | and |V cs |. We show the constraints on |V cd | and |V cs | from our calculation in Fig. 33, together with constraints coming from leptonic decays [32,44] and second-row unitarity. The leptonic inputs used for the green ellipse are summarized in Table XX. As the figure shows, semileptonic tests of second-row CKM unitarity are now slighty more precise than leptonic tests. The leptonic and semileptonic results are consistent at the level of roughly 1-2 standard deviations. One can perform further tests of the unitarity of the CKM matrix using the fact that in the Standard Model, |V cd | = |V us | + O(A 2 λ 5 ) and |V cs | = |V ud | + O(A 2 λ 4 ). Including the dominant corrections [146] with the Wolfenstein parameters taken from global unitarity fits by CKMFitter [144] (using values from the January 2022 update) gives |V cs | = 0.97282(32) from |V ud | 0 + →0 + , (7.9) |V cd | = 0.22317(53) from |V us | K 3 , (7.10) |V cd |/|V cs | = 0.22941(55) from |V us | K 3 /|V ud | 0 + →0 + , (7.11) using |V ud | = 0.97367(32) from superallowed 0 + → 0 + nuclear β decays [17,37] [25]. In all cases, the ellipses shows the correlated 1σ (68%) confidence intervals. The inner blue ellipse shows our result without the QED uncertainty. Table XX. Leptonic inputs used for comparison in Fig. 33. HFLAV reports the product η EW |V cx |f D (s) [44]. Following the prescription of the PDG [25], we include an EW+QED error of 0.7% for the product |V cx |f D (s) .

E. Lepton flavor universality
For a given semileptonic decay H → L ν, the lepton flavor universality (LFU) ratio R µ/e is defined as the ratio of branching fractions into muon versus electron final states where the total rates to each final state are defined in the usual way, In the SM, the LFU ratios are close but not identically equal to unity. This difference from unity arises from at least three effects. First, the lower boundary of the integration region in Eq. (7.14) depends on the lepton mass. Second, the differential decay rate in Eq. (2.1) itself depends on the lepton mass, with the scalar form factor contributing more for larger masses. Finally, QED corrections depend in principle on both the charges of the final state and the lepton mass. The coefficients G 2 F 24π 3 (η EW |V cx |) 2 are independent of q 2 and cancel in the ratio, meaning that predictions for R µ/e are entirely calculable using our lattice-QCD form factors, up to corrections from QED and SIB. The rates dΓ/dq 2 for the decay D → π, using as inputs our form factors f 0 (q 2 ) and f + (q 2 ) together with the estimates of systematic uncertainties from QED and SIB (see Sec. VII B), are shown in Fig. 34 for both semielectronic and semimuonic final states. When computing the rates, the meson masses were taken to be the average of the experimentally measured masses for the charged and neutral states (e.g., D 0 and D + or π 0 and π + ). The final results for the SM predictions of the ratios R µ/e are The dominant error is the systematic uncertainty from QED corrections, which we conservatively take to be 0.5%, as described in Sec. VII B. Our prediction for R D→K µ/e is in good agreement with a recent calculation by HPQCD, which found R D→K µ/e = 0.97594(19) QCD [500] QED and used the same estimate of the QED uncertainty [62]. 12 We also find good agreement with previous lattice QCD results by ETMC and experimental measurements of R D→π µ/e and R D→K µ/e , as shown in Fig. 35. The measurement of R D→π µ/e by BES III for the channel D 0 → π − lies below our result but is consistent at the 2σ level. Because the QED error is dominant for the lattice-QCD predictions of the LFU ratios, the insets in Fig. 35 compare the lattice-QCD results with the QED uncertainty removed.  with experimental HFLAV averages [44], which are dominated by measurements from BES III [125,127,128], and other SM predictions from lattice QCD [62,139]. In the main body of both figures, all lattice QCD results are presented with a QED uncertainty of 0.5%. The results from ETMC 18 were reported in the isospin-symmetric limit of QCD, without including QED or SIB uncertainties [139], so we have added the QED uncertainty for a like-to-like comparison. The insets compare lattice QCD results when QED uncertainty is removed.

VIII. CONCLUSIONS
We have calculated the hadronic form factors f + (q 2 ) and f 0 (q 2 ) relevant for the semileptonic decays D → π ν, D → K ν, and D s → K ν using lattice QCD. These decays occur at tree level in the SM and are important channels for determining the CKM matrix elements |V cd | and |V cs |. Our calculation uses N f = 2 + 1 + 1 flavors of dynamical staggered quarks and includes several ensembles with all quarks near their physical masses. The use of the HISQ action permits all the quarks to be treated with the same relativistic light-quark action and allows for nonperturbative renormalization using PCVC, Eq. (2.9). Our results improve significantly on the previous precision for the form factors for D → π and D s → K, and have precision comparable to that of recent N f = 2 + 1 + 1 calculations by HPQCD for D → K [62,63]. We agree well with HPQCD's D → K form factors over the entire kinematic range, especially with their latest results in Ref. [63], while for both D → π and D → K, our form factors are significantly larger near q 2 max than the N f = 2 + 1 + 1 results of ETMC [57]. Table XIII shows the z-expansion parameters from which our final results for the form factors, computed in isospin-symmetric QCD where m u = m d , can be reconstructed, while a complete error budget, including all statistical and systematic uncertainties, is given in Table XV for the edges of the kinematic range.
Our results suggest a very mild spectator dependence for D → π and D s → K, with close agreement at 2% level throughout the kinematic range between the respective form factors (cf. Fig. 18). This picture was also recently confirmed, within experimental uncertainty, by the first measurement of the decay D s → K by BES III [109].
When combined with the available experimental data for the corresponding decay rates, summarized in Table XVII, our form factors enable the extraction of the CKM matrix elements |V cd | and |V cs | with percent-level uncertainties. These extractions include correlations between all the lattice form factors and between the different experimental channels. 13 The values obtained from our preferred extractions are For |V cd | we obtain the most precise determination to date, with lattice-QCD form factors errors that, for the first time in a semileptonic extraction, are commensurate with experimental uncertainties. The improved determination of D → π form factors, together with the fact that we account for theoretical correlations among channels, also allows us to provide the most precise determination of the ratio |V cd |/|V cs |, around a factor of two more precise than the leptonic determination. The rate for D s → K was only recently measured for the first time by BES III [109], and our calculation delivers the first extraction of |V cd | Ds→K . Although this determination is not yet competitive with the one from D → π, the error is dominated by the statistics-limited experimental uncertainty. Our result for |V cd | Ds→K lies roughly 2σ above |V cd | D→π , albeit with large uncertainty. Experimental improvements for this Cabibbo-suppressed decay would immediately give improved precision for |V cd | Ds→K and help clarify the situation.
Our determinations of |V cd | and |V cs |, combined with the value of |V cb | from Ref. [25], give a precise test of second-row CKM unitarity. We find consistency with unitarity at the level of roughly 2% and one standard deviation, with an uncertainty dominated by the systematic effect of QED. As shown in Fig. 33, the precision of the semileptonic constraint is now slightly better than the corresponding leptonic one.
After demonstrating consistency between the form factor shapes from our calculations and those measured in experiments, we computed the SM prediction for the lepton flavor universality ratios R µ/e with sub-percent precision for all three decays: These results agree with previous N f = 2 + 1 + 1 lattice calculations, considerably improving the precision for D → π, and with experimental measurements within 2σ, for D → π and D → K.
With the total precision for |V cd | and |V cs | approaching the subpercent level, the effects of the scalar form factor in the differential rate, Eq. (2.1), become quantitatively important. For semielectronic decays, contributions from f 0 enter at roughly the 1% level in the lowest q 2 bin. For semimuonic decays, the effect is much larger, a roughly 10% effect in the lowest q 2 bin and a few-percent effect throughout the rest of the kinematic range. Figure 30 showed that naively neglecting contributions from f 0 can shift values for |V cs | by a few percent in the case of D → Kµν (similar results hold for D → πµν).
Future progress in the precision of |V cd |, |V cs |, and the LFU ratios will depend crucially on improved understanding of QED corrections to these decays, which are already the dominant source of uncertainty. The one exception is the decay |V cs | Ds→K + ν , for which the experimental error is still large. One avenue for improvement is through EFT calculations in the spirit of those for K → π ν [41-43, 133, 134], which were used in Sec. VII B to estimate our systematic uncertainties (cf. Sec. VII B). As usual, the intermediate mass of the charm quark (which is simultaneously too heavy for χPT to apply and too light for reliable application of HQET) may present a challenge for robust treatment with EFT. Another possibility is carrying out lattice simulations to compute the structure-dependent QED corrections to the semileptonic decay amplitudes. Such calculations have not yet reached a mature state, but the field is progressing rapidly, particularly for the QED corrections to leptonic decays [39,40,[147][148][149][150][151].
Regarding the pure QCD calculation, it should be straightforward to improve the precision of our form factor results. A leading contribution to the error budget is statistics (cf .  Table XV and Figs. 22, 49 and 50), for which the physical mass ensembles at a ≈ 0.06 fm and 0.09 fm play the largest role. As part of our ongoing work toward B-meson semileptonic decays, we are simulating on a finer physical-mass ensemble with a ≈ 0.04 fm. We expect that new data from this ensemble will reduce the uncertainties both from statistics and from the continuum extrapolation. Future calculations will also benefit from ongoing work in the community to improve scale-setting measurements (e.g., w 0 or the Ω-baryon mass) on the HISQ ensembles used in this work.

Appendix A: Analysis of staggered correlation functions
As demonstrated in Ref. [154], averaging over adjacent time slices can dramatically suppress contributions from oscillating states. Consider a two-point correlation function C 2 (t) Let E denote the ground state energy. Then the averaged two-point function is C 2 (t): where O is an interpolating operator as given in Table III and |∅ is the QCD vacuum. Similarly, consider a three-point correlation function C 3 (t, T ) with ground states E L and E H at the source and sink, respectively, and connected by the current J. The averaged three-point function is C 3 (t, T ): These averaged two-and three-point functions are used in Eqs. (4.10)-(4.12).
generic (heavy or light) quark x. Arguments from leading-order HQET [32] show that the conventional normalization can be restored, at leading order, by multiplying matrix elements containing the bilinear by the factor Z HQET, LO hx Ch(am 1 ) = cosh(am 1 ) 1 − 1 2 N(am 1 ) sinh 2 (am 1 ) , Residual discretization effects from next-to-leading HQET appear at order x 4 h and α s x 2 h , where x h is the parameter linear in the heavy quark mass defined in Eq. (5.14).
This result has important consequence for the near-singular covariance matrices encountered in practical problems. Let diag(λ) = U T ΣU denote the spectral decomposition of the "true" population covariance matrix of a statistical distribution, where U contains the eigenvectors and λ are the eigenvalues. The corresponding sample covariance matrix has decomposition diag(λ n ) = U T n S n U n . As usual, S n is an unbiased estimator of Σ. Therefore, U T S n U is also an unbiased estimator of diag λ. Unfortunately, one does not typically have access to the population eigenvectors of U and is instead obliged to work with the sample estimates of U n . As the preceding lemma makes clear, the sample eigenvalues λ n will be more widely dispersed than those of the population λ. Indeed, λ n is not an unbiased estimator of U T ΣU due to correlations between the eigenvectors in U n and the eigenvalues in λ n . The general idea behind shrinkage estimators is to apply some function which decreases the dispersion of the sample eigenvalues λ n to better approximate the population λ.
The remainder of this appendix is organized as follows. Appendix C 1 describes linear shrinkage, which is used in the chiral-continuum analysis (cf. Sec. V B). Appendix C 2 describes nonlinear shrinkage, which is used in the correlator fits (cf. Sec. IV).

Linear shrinkage
Linear shrinkage was introduced by Ledoit and Wolf in Ref. [157]. There seems to be some knowledge of this technique in the recent lattice-gauge-theory literature [158]. Because lattice data often vary over many orders of magnitude, it is common to invert the correlation matrix instead of the covariance matrix, with shrinkage techniques being applied to them instead.
The linear shrinkage estimatorĈ n is defined as the convex sum of two matrices: with λ ∈ [0, 1]. As the parameter λ is varied, the shrinkage estimator smoothly interpolates between the sample correlation matrix C n and the target matrix C target . Many options are possible for C target . Examples in the literature [157,158] advocate using the identity matrix as the shrinkage target. The idea is that suppressing the correlations by a small amount (say, λ = 0.05 or 0.1) will correct the small eigenvalues while preserving the rest of the correlated structure to the data. 14 Besides using the identity matrix, our analysis also experimented with block-diagonal matrices (e.g., to retain the full correlations between different momenta at fixed valence mass). The more complicated choices did not improve fit results compared with the simpler choice of the identity matrix. The preferred chiral-continuum analysis of Sec. V therefore uses only the identity matrix. Once a shrinkage estimator for the correlation matrix has been chosen, the corresponding covariance matrix follows in the usual way, where σ is a vector containing the standard deviations. The shrinkage estimator, which enjoys a smaller condition number and approximates the population covariance matrix better than the sample estimate, is then inverted to giveŜ −1 n , which is used in our fits.   Figure 36. Comparison of eigenvalue spectra resulting before and after shrinkage or an SVD cut, for the correlation matrices for C P Ds (t) (left) and C P K (t, 0) (right) on the physical-mass a ≈ 0.06 fm ensemble. Linear shrinkage was applied with λ = 0.1. An SVD cut of 10 −3 was chosen to have an effect on the spectrum similar to shrinkage.

Numerical examples of shrinkage
In this section, we present representative examples of correlation matrices appearing in our analysis. For concreteness, we consider the correlation matrices for the two-point functions C P Ds (t) and C P K (t, 0) (cf. Eqs. (4.1) and (4.3)) associated with the D s and K mesons on the physical-mass a ≈ 0.06 fm ensemble. The eigenvalue spectra associated with the correlation matrices are shown in Fig. 36, for raw data, nonlinear shrinkage, linear shrinkage with λ = 0.1, and an SVD cut of 10 −3 The raw spectra, shown in blue, span a range of roughly eight orders of magnitude. (In fact, not displayed are the last few eigenvalues, which are consistent with zero at double precision). For the given parameter choices, linear shrinkage and the SVD cut give similar results. With nonlinear shrinkage, the shape of the small-eigenvalue region of the spectrum retains some of its original curvature. In the case of the kaon (left in Fig. 36), the small eigenvalues from nonlinear shrinkage vary by approximately an order of magnitude over the region where they are roughly constant for linear shrinkage and SVD cut.
These methods all alter the covariance between pairs of data. Figures 37 and 38 show heat maps for the corresponding correlation matrices. As with the eigenvalue spectra in Fig. 36, the results for linear shrinkage and SVD cut are qualitatively similar. Compared with the other methods, nonlinear shrinkage tends to smooth the far off-diagonal correlation coefficients. All three correction methods suppress the near-diagonal correlations which are nearly unity in the raw data.                Table I were included in the fit.  Table I were included in the fit. after the fit to the z expansion. Contributions less than 0.01% are not shown.

z-expansion fits: Joint fits to lattice-QCD form factors and experimental data
Tables XXI to XXIII compare the results of the z-expansion fits for the decays D → π, D → K, and D s → K. The fits enforce the kinematic identity f + (0) = f 0 (0) by imposing a 0 = b 0 [cf. Eqs. (5.23) and (5.24)]. For the scalar form factor, the higher parameters b 1 , b 2 , and b 3 are unconstrained by the fits including experimental data. In the joint fit, the lattice QCD form factors include a systematic from SIB, as described in Sec. VII B. No uncertainty from QED is included in the fit, since this is applied directly to |V cx | as a final 1% uncertainty.