Short-distance matrix elements for $D^0$-meson mixing for $N_f=2+1$ lattice QCD

TL;DR

This paper computes short-distance matrix elements for $D^0$-meson mixing using three-flavor lattice QCD, providing essential data for understanding Standard Model and beyond Standard Model contributions to charm mixing.

Contribution

It presents the first comprehensive lattice QCD calculation of all five $ riangle C=2$ four-fermion operator matrix elements for $D^0$ mixing in the $N_f=2+1$ flavor setup, with detailed error analysis and scheme conversions.

Findings

Matrix elements evaluated at 3 GeV in MS-bar scheme.

Lower bounds on new physics scale from $D^0$ mixing data.

Correlation matrices provided for phenomenological applications.

Abstract

We calculate in three-flavor lattice QCD the short-distance hadronic matrix elements of all five $\Delta C=2$ four-fermion operators that contribute to neutral $D$-meson mixing both in and beyond the Standard Model. We use the MILC Collaboration's $N_f = 2+1$ lattice gauge-field configurations generated with asqtad-improved staggered sea quarks. We also employ the asqtad action for the valence light quarks and use the clover action with the Fermilab interpretation for the charm quark. We analyze a large set of ensembles with pions as light as $M_\pi \approx 180$ MeV and lattice spacings as fine as $a\approx 0.045$ fm, thereby enabling good control over the extrapolation to the physical pion mass and continuum limit. We obtain for the matrix elements in the $\overline{\text{MS}}$-NDR scheme using the choice of evanescent operators proposed by Beneke \emph{et al.}, evaluated at 3 GeV, $\langle D^0|\mathcal{O}_i|\bar{D}^0 \rangle = \{0.0805(55)(16), -0.1561(70)(31), 0.0464(31)(9), 0.2747(129)(55), 0.1035(71)(21)\}~\text{GeV}^4$ ($i=1$--5). The errors shown are from statistics and lattice systematics, and the omission of charmed sea quarks, respectively. To illustrate the utility of our matrix-element results, we place bounds on the scale of CP-violating new physics in $D^0$~mixing, finding lower limits of about 10--50$\times 10^3$ TeV for couplings of $\mathrm{O}(1)$. To enable our results to be employed in more sophisticated or model-specific phenomenological studies, we provide the correlations among our matrix-element results. For convenience, we also present numerical results in the other commonly-used scheme of Buras, Misiak, and Urban.

Tables

Table 1. Table 4: Strong coupling in the V 𝑉 V  scheme at the scale μ = 2 / a 𝜇 2 𝑎 \mu=2/a , and heavy-heavy and light-light vector-current renormalization factors Z V c ​ c 4 subscript 𝑍 superscript subscript 𝑉 𝑐 𝑐 4 Z_{V_{cc}^{4}} and Z V l ​ l 4 subscript 𝑍 superscript subscript 𝑉 𝑙 𝑙 4 Z_{V_{ll}^{4}} with statistical errors for the listed simulation κ c ′ superscript subscript 𝜅 𝑐 ′ \kappa_{c}^{\prime} values. Also shown are the nonperturbatively determined critical hopping parameter κ crit subscript 𝜅 crit \kappa_{\text{crit}} and the tadpole-improvement factors u 0 ​ P subscript 𝑢 0 𝑃 u_{0P} and u 0 ​ L subscript 𝑢 0 𝐿 u_{0L} obtained from the fourth root of the plaquette and the link in Landau gauge, respectively.

$\approx a$ (fm)	$a{m}_l^{\prime }/a{m}_s^{\prime }$	${\kappa }_c^{\prime }$	${\kappa }_{\text{crit}}$	$u_{0P}$	$u_{0L}$	${\alpha }_V\left(2/a\right)$	$Z_{V_{cc}^4}$	$Z_{V_{ll}^4}$
0.12	0.02/0.05	0.1259	0.14073	0.8688	0.837	0.3047	0.2912(1)	1.734(3)
	0.01/0.05	0.1254	0.14091	0.8677	0.835	0.3108	0.2947(1)	1.729(3)
	0.01/0.05	0.1280	0.14091	0.8677	0.835	0.3108	0.2786(1)	1.729(3)
	0.007/0.05	0.1254	0.14095	0.8678	0.836	0.3102	0.2946(1)	1.730(3)
	0.005/0.05	0.1254	0.14096	0.8678	0.836	0.3102	0.2946(1)	1.729(3)
0.09	0.0124/0.031	0.1277	0.139052	0.8788	—	0.2582	0.2761(2)	1.768(4)
	0.0062/0.031	0.1276	0.139119	0.8782	0.854	0.2607	0.2769(2)	1.766(4)
	0.00465/0.031	0.1275	0.139134	0.8781	—	0.2611	0.2776(2)	1.766(4)
	0.0031/0.031	0.1275	0.139173	0.8779	—	0.2619	0.2777(2)	1.765(4)
	0.00155/0.031	0.1275	0.139190	0.877805	—	0.2623	0.2777(2)	1.765(4)
0.06	0.0072/0.018	0.1295	0.137582	0.8881	—	0.2238	0.2614(2)	1.798(5)
	0.0036/0.018	0.1296	0.137632	0.88788	—	0.2245	0.2611(2)	1.797(5)
	0.0025/0.018	0.1296	0.137667	0.88776	0.869	0.2249	0.2612(2)	1.797(5)
	0.0018/0.018	0.1296	0.137678	0.88764	0.869	0.2253	0.2610(2)	1.796(5)
0.045	0.0028/0.014	0.1310	0.136640	0.89511	0.8797	0.2013	0.2498(2)	1.818(8)

Table 2. Table 5: One-loop renormalization coefficients defined in Eqs. ( 40 )–( 41 ) for the heavy-light Δ ​ C = 2 Δ 𝐶 2 \Delta C=2 local four-fermion operators at the renormalization scale μ = 3 𝜇 3 \mu=3  GeV. The renormalization coefficients are defined in Eqs. ( 42 ) and ( 45 ) 45 (\ref{eq:TadPT_zeta}) and calculated for the BBGLN  Beneke  et al. ( 1999 , 2002 ) choice of evanescent operators; their relation to the BMU  Buras  et al. ( 2000 ) scheme is given in Eqs. ( 44 ). The last two columns show the differences between the one-loop coefficients obtained in perturbation theory and mNPR for two choices of the tadpole-improvement factor u 0 ​ P subscript 𝑢 0 𝑃 u_{0P} and u 0 ​ L subscript 𝑢 0 𝐿 u_{0L} . Renormalization coefficients not listed are vanishing.

$\approx a\text{(fm)}$	$\begin{array}{c}a{m}_l^{\prime }\\ & \\ a{m}_s^{\prime }\end{array}$	${\rho }_{11}^{\left[1\right]}$	${\rho }_{12}^{\left[1\right]}$	${\rho }_{22}^{\left[1\right]}$	${\rho }_{21}^{\left[1\right]}$	${\rho }_{33}^{\left[1\right]}$	${\rho }_{31}^{\left[1\right]}$	${\rho }_{44}^{\left[1\right]}$	${\rho }_{45}^{\left[1\right]}$	${\rho }_{55}^{\left[1\right]}$	${\rho }_{54}^{\left[1\right]}$	${\zeta }_{P,ii}^{\left[1\right]}-{\rho }_{ii}^{\left[1\right]}$	${\zeta }_{L,ii}^{\left[1\right]}-{\rho }_{ii}^{\left[1\right]}$
0.12	0.4	$-0.2447$	$-0.0887$	$0.0107$	$0.0403$	$0.4103$	$-0.0097$	$-0.1728$	$-0.2820$	$0.0097$	$-0.3134$	$0.1382$	$0.0893$
$\left({\kappa }_c^{\prime }=0.1254\right)$	0.2	$-0.2602$	$-0.0927$	$0.0353$	$0.0423$	$0.3927$	$-0.0025$	$-0.1344$	$-0.2837$	$0.0000$	$-0.2974$	$0.1447$	$0.0960$
$\left({\kappa }_c^{\prime }=0.1280\right)$	0.2	$-0.2341$	$-0.0766$	$0.0369$	$0.0406$	$0.4133$	$-0.0065$	$-0.1396$	$-0.2767$	$0.0205$	$-0.3047$	$0.1178$	$0.0682$
	0.14	$-0.2605$	$-0.0929$	$0.0353$	$0.0423$	$0.3924$	$-0.0025$	$-0.1343$	$-0.2837$	$-0.0003$	$-0.2975$	$0.1450$	$0.0963$
	0.1	$-0.2606$	$-0.0930$	$0.0353$	$0.0423$	$0.3925$	$-0.0025$	$-0.1341$	$-0.2838$	$-0.0004$	$-0.2973$	$0.1451$	$0.0963$
0.09	0.4	$-0.1077$	$-0.0689$	$-0.2540$	$0.0212$	$0.5724$	$-0.0822$	$-0.5806$	$-0.2723$	$0.0884$	$-0.4724$	$0.1018$	$0.0519$
	0.2	$-0.1156$	$-0.0698$	$-0.2385$	$0.0223$	$0.5629$	$-0.0780$	$-0.5571$	$-0.2728$	$0.0839$	$-0.4631$	$0.1038$	$0.0538$
	0.14	$-0.1183$	$-0.0706$	$-0.2350$	$0.0226$	$0.5601$	$-0.0768$	$-0.5510$	$-0.2732$	$0.0821$	$-0.4606$	$0.1051$	$0.0552$
	0.1	$-0.1202$	$-0.0708$	$-0.2310$	$0.0229$	$0.5577$	$-0.0758$	$-0.5450$	$-0.2733$	$0.0810$	$-0.4584$	$0.1055$	$0.0556$
	0.05	$-0.1218$	$-0.0709$	$-0.2274$	$0.0232$	$0.5558$	$-0.0748$	$-0.5392$	$-0.2733$	$0.0801$	$-0.4562$	$0.1057$	$0.0558$
0.06	0.4	$0.0445$	$-0.0499$	$-0.5407$	$-0.0003$	$0.7495$	$-0.1615$	$-1.0218$	$-0.2611$	$0.1752$	$-0.6433$	$0.0598$	$0.0089$
	0.2	$0.0425$	$-0.0496$	$-0.5335$	$0.0001$	$0.7464$	$-0.1599$	$-1.0111$	$-0.2609$	$0.1745$	$-0.6394$	$0.0590$	$0.0082$
	0.14	$0.0407$	$-0.0498$	$-0.5297$	$0.0003$	$0.7442$	$-0.1590$	$-1.0055$	$-0.2610$	$0.1734$	$-0.6373$	$0.0595$	$0.0086$
	0.1	$0.0392$	$-0.0498$	$-0.5261$	$0.0006$	$0.7422$	$-0.1579$	$-1.0000$	$-0.2610$	$0.1726$	$-0.6352$	$0.0597$	$0.0088$
0.045	0.2	$0.1771$	$-0.0355$	$-0.7880$	$-0.0192$	$0.9027$	$-0.2306$	$-1.4020$	$-0.2514$	$0.2504$	$-0.7902$	$0.0227$	$-0.0290$

Table 3. Table 6: Tuned κ c subscript 𝜅 𝑐 \kappa_{c} values from Ref.  Bailey  et al. ( 2014 ) , and differences Δ ​ ( 1 / ( a ​ m 2 ) ) Δ 1 𝑎 subscript 𝑚 2 \Delta(1/(am_{2})) between the simulated and physical inverse charm-quark kinetic masses on each ensemble. For κ c subscript 𝜅 𝑐 \kappa_{c} , the first error is from statistics and fitting, and the second is from the uncertainty in the lattice scale r 1 subscript 𝑟 1 r_{1} . For Δ ​ ( 1 / ( a ​ m 2 ) ) Δ 1 𝑎 subscript 𝑚 2 \Delta(1/(am_{2})) the error is from the uncertainty in the tuned κ c subscript 𝜅 𝑐 \kappa_{c} .

$\approx a$ (fm)	$a{m}_l^{\prime }/a{m}_s^{\prime }$	${\kappa }_c$	$\Delta \left(1/\left(a{m}_2\right)\right)$
0.12	0.02/0.05	0.12452(15)(16)	$-$0.212(31)
0.12	0.01/0.05	0.12423(15)(16)	$-$0.168(29)
0.12	0.007/0.05	0.12423(15)(16)	$-$0.167(29))
0.12	0.005/0.05	0.12423(15)(16)	$-$0.167(29)
0.09	0.0124/0.031	0.12737(9)(14)	$-$0.089(43)
0.09	0.0062/0.031	0.12722(9)(14)	$-$0.099(42)
0.09	0.00465/0.031	0.12718(9)(14)	$-$0.082(42)
0.09	0.0031/0.031	0.12714(9)(14)	$-$0.092(41)
0.09	0.00155/0.031	0.12710(9)(14)	$-$0.101(41)
0.06	0.0072/0.018	0.12964(4)(11)	0.075(64)
0.06	0.0036/0.018	0.12960(4)(11)	0.000(63)
0.06	0.0025/0.018	0.12957(4)(11)	$-$0.016(62)
0.06	0.0018/0.018	0.12955(4)(11)	$-$0.026(61)
0.045	0.0028/0.014	0.130921(16)(7)	$-$0.083(75)

Table 4. Table 7: Slopes μ i subscript 𝜇 𝑖 \mu_{i} defined in Eq. ( 49 ) for the D 𝐷 D -mixing matrix elements renormalized in the continuum MS ¯ ¯ MS \overline{\rm MS} -NDR scheme and the BBGLN  Beneke  et al. ( 1999 , 2002 ) and BMU  Buras  et al. ( 2000 ) choices of evanescent operators in r 1 subscript 𝑟 1 r_{1} units. The uncertainties of r 1 / a subscript 𝑟 1 𝑎 r_{1}/a are not included in this table, but instead propagated while performing the chiral-continuum extrapolation as discussed in Sec.  VII.2.2 . The r 1 / a subscript 𝑟 1 𝑎 r_{1}/a values are precise enough to not affect uncertainties of the slopes at the reported level of precision.

$i$	1	2		3		4	5
		BMU	BBGLN	BMU	BBGLN
${\mu }_i$	$-$0.248(90)	0.084(45)	0.073(45)	$-$0.017(22)	$-$0.011(22)	$-$0.12(21)	$-$0.141(90)

Table 5. Table 8: Relationships between the coefficients of the wrong-spin terms in Eq. ( 51 ). The relationship between the primed LECs are analogous to those between the unprimed LECs, but with all entries in the table β j → β j ′ → subscript 𝛽 𝑗 superscript subscript 𝛽 𝑗 ′ \beta_{j}\rightarrow\beta_{j}^{\prime} in β i ′ ⁣ ( ξ ) subscript superscript 𝛽 ′ 𝜉 𝑖 \beta^{\prime(\xi)}_{i} .

		$\xi$
		P	A	T	V	I
	$1$	${\beta }_1$	$-8{\beta }_2-8{\beta }_3$	$-6{\beta }_1$	$8{\beta }_2+8{\beta }_3$	${\beta }_1$
	$2$	${\beta }_2$	$-{\beta }_1$	$-2{\beta }_2+4{\beta }_3$	${\beta }_1$	${\beta }_2$
$i$	$3$	${\beta }_3$	$-{\beta }_1$	$-4{\beta }_2-2{\beta }_3$	${\beta }_1$	${\beta }_3$
	$4$	$-{\beta }_4$	$-2{\beta }_5$	$0$	$-2{\beta }_5$	${\beta }_4$
	$5$	$-{\beta }_5$	$-2{\beta }_4$	$0$	$-2{\beta }_4$	${\beta }_5$

Table 6. Table 9: Priors for heavy-quark discretization terms in Eqs. ( 60 )–( 62 ).

$z_B$	$z_3$	$z_E$	$z_X$	$z_Y$	$z_2$
$0\left(2\right)$	$0\left(\sqrt{5}\right)$	$0\left(2\sqrt{2}\right)$	$0\left(2\sqrt{2}\right)$	$0\left(2\right)$	0(2)

Table 7. Table 10: Taste splitting and leading-order LEC r 1 ​ B 0 subscript 𝑟 1 subscript 𝐵 0 r_{1}B_{0} for the lattice spacings analyzed in this work, and in the continuum  Bailey  et al. ( 2014 ) . The labels A , T , V , I 𝐴 𝑇 𝑉 𝐼 A,T,V,I denote the axial-vector, tensor, vector, and scalar tastes respectively.

$\approx a$ (fm)	$r_1^2{a}^2{\Delta }_A$	$r_1^2{a}^2{\Delta }_T$	$r_1^2{a}^2{\Delta }_V$	$r_1^2{a}^2{\Delta }_I$	$r_1{B}_0$
0.12	0.2270	0.3661	0.4803	0.6008	6.832
0.09	0.0747	0.1238	0.1593	0.2207	6.639
0.06	0.0263	0.0430	0.0574	0.0704	6.487
0.045	0.0104	0.0170	0.0227	0.0278	6.417
continuum	0	0	0	0	6.015

Table 8. Table 11: Breakdown of the chiral-continuum fit error. The labels and estimation procedure are described in the text. Entries are in percent.

	stat.	inputs	$\kappa$ tuning	matching	chiral	LQ disc	HQ disc	$r_1/a$	fit total
$⟨{𝒪}_1⟩$	3.5	0.6	1.5	3.8	1.3	0.6	3.1	0.4	6.4
$⟨{𝒪}_2⟩$	1.8	0.5	0.4	2.2	0.8	0.4	2.4	0.5	4.0
$⟨{𝒪}_3⟩$	3.1	0.3	0.6	3.8	1.3	0.5	3.6	0.4	6.3
$⟨{𝒪}_4⟩$	2.2	0.6	0.5	2.0	0.9	0.3	2.6	0.5	4.2
$⟨{𝒪}_5⟩$	3.0	0.7	0.5	4.1	1.5	0.5	3.5	0.3	6.5

Table 9. Table 12: Total error budgets for D 𝐷 D -mixing matrix elements. Entries are in percent. The uncertainty due to the omission of the charmed sea quark is listed separately because the estimated error is significantly less quantitative than that from the other contributions.

	Fit total	$r_1$	FV	EM	Total	Charm sea
$⟨{𝒪}_1⟩$	6.4	2.1	0.1	0.2	6.8	2.0
$⟨{𝒪}_2⟩$	4.0	2.1	0.3	0.2	4.5	2.0
$⟨{𝒪}_3⟩$	6.3	2.1	0.3	0.2	6.6	2.0
$⟨{𝒪}_4⟩$	4.2	2.1	0.2	0.2	4.7	2.0
$⟨{𝒪}_5⟩$	6.5	2.1	0.2	0.2	6.8	2.0

Table 10. Table 13: D 𝐷 D -meson mixing matrix elements ⟨ 𝒪 i ⟩ delimited-⟨⟩ subscript 𝒪 𝑖 \left<\mathcal{O}_{i}\right> in the MS ¯ − NDR ¯ MS NDR \overline{\text{MS}}-\text{NDR} scheme at μ = 3 ​ GeV 𝜇 3 GeV \mu=3~{}\text{GeV} with total uncertainties. The first uncertainty is the error labeled “Total” in Table  12 , while the second is the estimated error from quenching the charm sea quark. Results are shown for both the BBGLN  Beneke  et al. ( 1999 ) and BMU  Buras  et al. ( 2000 ) evanescent-operator choices. Entries are in GeV 4 .

BBGLN	0.0805(55)(16)	-0.1561(70)(31)	0.0464(31)(9)	0.2747(129)(55)	0.1035(71)(21)
BMU	0.0806(54)(16)	-0.1442(66)(29)	0.0452(30)(9)	0.2745(129)(55)	0.1035(71)(21)

Table 11. Table 14: Correlations between the relative scale r 1 / a subscript 𝑟 1 𝑎 r_{1}/a on the ensembles used in this work. The central values and errors are given in Table  1 . Entries are symmetric across the diagonal.

$\approx a$ (fm)		0.12				0.09					0.06				0.045	0
	$m_l^{\prime }/m_s^{\prime }$	0.10	0.14	0.20	0.40	0.05	0.10	0.14	0.20	0.40	0.10	0.14	0.20	0.40	0.20	0.40
0.12	0.10	$1.000$
	0.14	$1.000$	$1.000$
	0.20	$1.000$	$1.000$	$1.000$
	0.40	$0.978$	$0.978$	$0.978$	$1.000$
0.09	0.05	$-0.177$	$-0.177$	$-0.177$	$0.006$	$1.000$
	0.10	$-0.174$	$-0.174$	$-0.174$	$0.008$	$1.000$	$1.000$
	0.14	$-0.170$	$-0.170$	$-0.170$	$0.010$	$0.999$	$1.000$	$1.000$
	0.20	$-0.165$	$-0.166$	$-0.166$	$0.013$	$0.998$	$0.999$	$1.000$	$1.000$
	0.40	$-0.139$	$-0.139$	$-0.139$	$0.032$	$0.987$	$0.990$	$0.993$	$0.995$	$1.000$
0.06	0.10	$0.538$	$0.538$	$0.538$	$0.410$	$-0.441$	$-0.425$	$-0.409$	$-0.391$	$-0.306$	$1.000$
	0.14	$0.540$	$0.540$	$0.540$	$0.411$	$-0.446$	$-0.430$	$-0.414$	$-0.396$	$-0.312$	$1.000$	$1.000$
	0.20	$0.541$	$0.541$	$0.541$	$0.412$	$-0.451$	$-0.436$	$-0.419$	$-0.401$	$-0.317$	$1.000$	$1.000$	$1.000$
	0.40	$0.543$	$0.544$	$0.544$	$0.413$	$-0.460$	$-0.445$	$-0.429$	$-0.411$	$-0.327$	$1.000$	$1.000$	$1.000$	$1.000$
0.045	0.20	$0.572$	$0.572$	$0.572$	$0.428$	$-0.583$	$-0.570$	$-0.555$	$-0.540$	$-0.463$	$0.983$	$0.984$	$0.985$	$0.987$	$1.000$
0	0.40	$-0.155$	$-0.155$	$-0.155$	$0.029$	$0.994$	$0.994$	$0.995$	$0.995$	$0.988$	$-0.426$	$-0.431$	$-0.436$	$-0.445$	$-0.569$	$1.000$

Table 12. Table 15: Priors employed in correlator fits in lattice-spacing units. Each prior is a Gaussian distribution with a central value and a one-sigma width, given in parentheses next to each central value. The uppermost panel shows constraints on energies of the two- and three-point correlators. The next panel gives constraints on the amplitudes of the two-point correlators. The lower two panels show additional constraints on the amplitude of the three-point correlators. All priors are in lattice units.

$\approx a$ (fm)	$E_0$	$\exp \left({\Delta }_{1,0}\right)$	$\exp \left({\Delta }_{k,k-2}\right)$
$0.12$	$1.0\left(0.1\right)$	$-1.2\left(0.5\right)$	$-1.0\left(1.5\right)$
$0.09$	$0.75\left(0.05\right)$	$-1.6\left(0.5\right)$	$-1.2\left(2.0\right)$
$0.06$	$0.56\left(0.04\right)$	$-2.0\left(1.0\right)$	$-1.6\left(2.0\right)$
$0.045$	$0.42\left(0.04\right)$	$-2.3\left(1.0\right)$	$-1.9\left(2.0\right)$
$\approx a$ (fm)	$Z_0$	$Z_1$	$Z_n$
$0.12$	$1.45\left(0.08\right)$	$0.7\left(0.5\right)$	$0.7\left(1.0\right)$
$0.09$	$1.42\left(0.08\right)$	$0.7\left(0.5\right)$	$0.7\left(1.0\right)$
$0.06$	$1.37\left(0.1\right)$	$0.7\left(0.5\right)$	$0.7\left(1.0\right)$
$0.045$	$1.36\left(0.1\right)$	$0.7\left(0.5\right)$	$0.7\left(1.0\right)$
$\approx a$ (fm)	${𝒵}_{00}^{O_1}$	${𝒵}_{00}^{O_2}$	${𝒵}_{00}^{O_3}$
$0.12$	$0.025\left(1.0\right)$	$-0.055\left(0.02\right)$	$0.012\left(0.006\right)$
$0.09$	$0.009\left(0.003\right)$	$-0.02\left(0.006\right)$	$0.005\left(0.002\right)$
$0.06$	$0.003\left(0.0015\right)$	$-0.007\left(0.003\right)$	$0.0017\left(0.0005\right)$
$0.045$	$0.0013\left(0.0008\right)$	$-0.0035\left(0.0015\right)$	$0.0008\left(0.0003\right)$
$\approx a$ (fm)	${𝒵}_{00}^{O_4}$	${𝒵}_{00}^{O_5}$	${𝒵}_{nm}^{O_i}$
$0.12$	$0.1\left(0.03\right)$	$0.04\left(0.015\right)$	$0.0\left(0.1\right)$
$0.09$	$0.04\left(0.01\right)$	$0.016\left(0.005\right)$	$0.0\left(0.05\right)$
$0.06$	$0.015\left(0.005\right)$	$0.006\left(0.002\right)$	$0.0\left(0.02\right)$
$0.045$	$0.007\left(0.002\right)$	$0.0027\left(0.0008\right)$	$0.0\left(0.01\right)$

Table 13. Table 16: Correlations between the D 𝐷 D -meson mixing matrix elements in the MS ¯ ¯ MS \overline{\rm{MS}} -NDR-BBGLN scheme given in Table  13 ; entries are symmetric across the diagonal. Correlations for the BMU scheme differ by less than 4%. The correlations include contributions from statistics and all systematics except the “charm sea” error, which is less well quantified. We suggest that an error of 2% be taken on all sums or differences of matrix elements and 0.5% on all ratios.

	1.0
	$-$0.2323	1.0
	0.0864	$-$0.2513	1.0
	0.2153	$-$0.3334	0.2163	1.0
	0.1865	$-$0.2246	0.1384	0.2574	1.0

Equations

$$i\frac{\partial }{\partial t}\left(\begin{array}{c}{D}^0\\ {\bar{D}}^0\end{array}\right)=$$

$$=$$

$$\Delta M$$

$$\Delta \Gamma$$

$$\Delta M-\frac{i}{2}\Delta \Gamma =\pm 2Q$$

$$Q^2$$

$${\varphi }_{12}\equiv \arg \frac{M_{12}}{{\Gamma }_{12}}.$$

$$y$$

$$x$$

$${\varphi }_{12}=-0.17^{\circ }\pm 1.8^{\circ },$$

$$2M_D\left(M_{12}-\frac{i}{2}{\Gamma }_{12}\right)=⟨D^0|{\mathcal{H}}^{\Delta C=2}|{\bar{D}}^0⟩+\sum _n\frac{⟨D^0|{\mathcal{H}}^{\Delta C=1}|n⟩⟨n|{\mathcal{H}}^{\Delta C=1}|{\bar{D}}^0⟩}{M_D-E_n+i{0}^{+}},$$

$$y$$

$$x$$

$${\mathcal{H}}^{\Delta C=2}=\sum _{i=1}^5{C}_i{\mathcal{O}}_i+\sum _{i=1}^3{\tilde{C}}_i{\tilde{\mathcal{O}}}_i,$$

$${\mathcal{O}}_1=$$

$${\mathcal{O}}_2=$$

$${\mathcal{O}}_3=$$

$${\mathcal{O}}_4=$$

$${\mathcal{O}}_5=$$

$${\tilde{\mathcal{O}}}_1=$$

$${\tilde{\mathcal{O}}}_2=$$

$${\tilde{\mathcal{O}}}_3=$$

$$2M_D{M}_{12}^{\text{NP}}=\sum _{i=1}^5{C}_i^{\text{NP}}(\mu )⟨D^0|{\mathcal{O}}_i|{\bar{D}}^0⟩(\mu ),$$

$$r_1=0.3117(22)\text{fm},$$

$${\rm Y}(x)=\Omega (x){\underline{\chi }}_q(x),$$

$$\Psi (x)=\left[1+a{d}_1\mathbf{\gamma }\cdot \mathbf{D}\right]\psi (x).$$

$$D^{†}(\mathbf{x},t)=\sum _{{\mathbf{x}}^{\prime }}\overline{{\rm Y}}(\mathbf{x},t)S(\mathbf{x},{\mathbf{x}}^{\prime }){\gamma }_5{\Psi }_c({\mathbf{x}}^{\prime },t),$$

$$C_{O_i}(t_x,t_y)=\sum _{\mathbf{x},\mathbf{y}}⟨D^{†}(\mathbf{y},t_0+t_y)O_i(0,t_0)D^{†}(\mathbf{x},t_0+t_x)⟩,$$

$$C(t-t_0)=\sum _{\mathbf{x}}⟨D(\mathbf{x},t)D^{†}(0,t_0)⟩.$$

$$C(t)=\sum _{n=0}^{N^{\text{2pt}}-1}(-1{)}^{n(t+1)}|Z_n{|}^2\left(e^{-E_{n}t}+e^{-E_n(T-t)}\right),$$

$$C_{O_i}(t_x,t_y)=\sum _{m,n=0}^{N^{\text{3pt}}-1}(-1{)}^{n(t_x+1)}(-1{)}^{m(t_y+1)}{\mathcal{Z}}_{nm}^{O_i}\frac{Z_n{Z}_m^{†}}{2\sqrt{E_n{E}_m}}{e}^{-E_n|t_x|}{e}^{-E_m{t}_y},$$

$$⟨O_i⟩={\mathcal{Z}}_{00}^{O_i}{M}_D.$$

$$M_{\mathrm{eff}}(t)={\text{cosh}}^{-1}\left(\frac{C(t+1)+C(t-1)}{2C(t)}\right)\stackrel{\text{large t}}{\to }{E}_0.$$

$$A_{\text{eff}}(t)=C(t)e^{t{M}_{\text{eff}}}\stackrel{\text{large t}}{\to }\frac{|Z_0{|}^2}{2E_0},$$

$$E_1\equiv E_0+{\Delta }_{1,0}$$

$${\Delta }_{k,j}\equiv \ln \left[\exp (E_k-E_j)\right].$$

Full text

I Introduction

The mixing between neutral $K$, $D$, $B$, and $B_{s}$ mesons and their antiparticles is loop suppressed in the Standard Model and, therefore, provides a window into new physics. Both indirect $CP$ violation in neutral kaon system ($\epsilon_{K}$) and the $B^{0}_{d}$ and $B^{0}_{s}$-meson oscillation frequencies ($\Delta M_{d}$ and $\Delta M_{s}$) have been measured to the subpercent level Patrignani et al. (2016); Amhis et al. (2017). Although still not as precise as experiment, the Standard Model theory for kaon and $B_{(s)}$-meson mixing is also under good control, owing to recent lattice-QCD calculations of the relevant hadronic matrix elements for kaons Dürr et al. (2011); Blum et al. (2016); Carrasco et al. (2015a) and for neutral $B_{(s)}$ mesons Bazavov et al. (2016). Neutral $D^{0}$-meson mixing remains the least understood of the four mixing processes, both theoretically and experimentally, but progress is being made on both sides.

In the Standard Model, neutral $D$-meson mixing is mediated at leading order in the electroweak interactions by intermediate down-type quarks, as illustrated in Fig. 1.

Hence, it provides unique information on new-physics contributions to the down-quark sector that is complementary to that provided by kaons and $B_{(s)}$ mesons, in which mixing is mediated by up-type quarks. In particular, $D$-meson mixing does not receive any top-quark enhancements at leading order. Further, mixing via the bottom quark is Cabibbo suppressed by $|V_{ub}V^{*}_{cb}|^{2}/|V_{u(d,s)}V^{*}_{c(d,s)}|^{2}\approx 0.2^{8}\sim\text{few}~{}10^{-6}$ relative to mixing via down and strange quarks. $D$-meson oscillations are thus, to a good approximation, facilitated by only two generations of quarks, and any observation of CP violation in $D$-meson mixing would be evidence of physics beyond the Standard Model (BSM).

At energies below the bottom quark mass, the electroweak box diagrams in Fig. 1 give rise to short-distance contributions from $\Delta C=2$ interactions and long-distance contributions from two $\Delta C=1$ interactions. The hadronic matrix elements of the former can be calculated within lattice QCD using standard methods and are the focus of this work. QCD calculations of hadronic matrix elements of the latter must wait for the development of better tools; we comment on the prospects for such calculations in Sec. X. Even though these long-distance effects are a dominant contribution to neutral $D$-meson mixing in the Standard Model, knowledge of the matrix elements of all short-distance $\Delta C=2$ operators that arise in the Standard Model and beyond can provide useful BSM discrimination Golowich et al. (2009), as described in more detail in Sec. II.

In this paper, we provide a new calculation of the $\Delta C=2$ $D$-mixing matrix elements on the MILC Collaboration’s $N_{f}=2+1$ gauge-field ensembles, which employ the $a^{2}$ tadpole-improved (asqtad) staggered action for the light quarks. We analyze the same set of ensembles as in our previous calculation of the $B$-mixing matrix elements Bazavov et al. (2016), and also follow an almost identical analysis procedure. Our results agree with previous $N_{f}=2$ and $N_{f}=2+1+1$ lattice-QCD calculations from the European Twisted Mass (ETM) Collaboration using twisted-mass fermions Carrasco et al. (2014a, 2015a), and have uncertainties commensurate with the projected experimental error on the mass difference $\Delta M$ [defined in Eq. (2)] from the LHCb and Belle II experiments Schwartz (2017); Morello (2016); Bona et al. (2016).

This paper is organized as follows. First, in Sec. II, we define the $\Delta C=2$ $D$-mixing matrix elements, and provide other important theoretical and phenomenological background. Next we present the set up of our numerical lattice-QCD calculation in Sec. III, including the lattice actions, simulation parameters, and two- and three-point correlation functions. We describe our two- and three-point correlator fits used to obtain the bare lattice mixing matrix elements in Sec. IV, followed by how we match these results to a continuum renormalization scheme in Sec. V. We adjust our matrix-element results to account for the slight difference between the simulation and physical charm-quark mass in Sec. VI before extrapolating our results to the physical light-quark mass and continuum limit in Sec. VII. In Sec. VIII, we discuss all sources of uncertainty in our calculation and provide a complete systematic error budget. Finally, we present our final $D$-mixing matrix-element results and discuss their phenomenological implications in Sec. IX, and conclude with an outlook for the future in Sec. X. Several appendixes provide additional details. Appendix A lists the correlations between the ratio of scales $r_{1}/a$ used in this work. The priors used in the two-point and three-point correlator fits are given in Appendix B, and Appendix C lists the correlations between our matrix-element results. Final results are provided in double-precision as Supplementary Material Dmi (2018).

II Theoretical and Phenomenological Background

The time evolution of a neutral-meson system, such as $D^{0}$ and $\bar{D}^{0}$, can be described by a Schrödinger equation

[TABLE]

where the mass matrix $M$ and decay matrix $\Gamma$ are Hermitian. The off-diagonal term $M_{12}-{\textstyle\frac{i}{2}}\Gamma_{12}$ mixes the flavor eigenstates into two mass (and width) eigenstates $D^{0}_{1}$ and $D^{0}_{2}$. Experiments have shown that CP violation in $D$-meson mixing is small, so the mass eigenstates are close to being CP eigenstates; usually 111HFAG’s charm web page HFA (2016) interchanges the labels $1\leftrightarrow 2$ on the eigenstates but ends up with the same physical convention for $\Delta M$ and $\Delta\Gamma$ as in Refs. Patrignani et al. (2016); Amhis et al. (2017). $|D^{0}_{1}\rangle$ is identified as the one with $\langle D^{0}_{1}|\text{CP}|D^{0}_{1}\rangle\approx+1$. Then, by convention, the mass and width differences between the two eigenstates are defined as Patrignani et al. (2016); Amhis et al. (2017)

[TABLE]

The signs of $\Delta M$ and $\Delta\Gamma$ are determined from experiment. The eigenvalue problem leads to the relation

[TABLE]

for the mass and width differences, where the sign is the (near) CP of the heavier state, and

[TABLE]

with

[TABLE]

Given measurements of $\Delta M$, $\Delta\Gamma$, and CP asymmetries sensitive to $\phi_{12}$, these formulas determine $M_{12}$ and $\Gamma_{12}$ apart from a mutual, unphysical phase Raz (2002); Ciuchini et al. (2007); Ball (2007).

The current measurements can be summarized as Amhis et al. (2017)

[TABLE]

and asymmetries consistent with zero. A fit then yields Amhis et al. (2017)

[TABLE]

and values for $y_{12}\equiv|\Gamma_{12}|/\Gamma$ and $x_{12}\equiv 2|M_{12}|/\Gamma$ the same as those for $y$ and $x$. It is expected that future measurements by LHCb and Belle II will reduce the uncertainties on $\Delta M$ to $10\%$ or less Schwartz (2017); Morello (2016); Bona et al. (2016).

The interpretation of these results within the Standard Model starts with the leading electroweak contributions, shown in Fig. 1. The $W$-boson mass, $b$-quark mass, and the typical scale of QCD, $\Lambda_{\text{QCD}}$, satisfy $m_{W}\gg m_{b}\gg\Lambda_{\text{QCD}}$, so the mixing matrix in Eq. (1) can be expressed as (see, e.g., Ref. Artuso et al. (2008))

[TABLE]

where $\mathcal{H}^{\Delta C=1}$ ($\mathcal{H}^{\Delta C=2}$) is an effective Hamiltonian changing charm by 1 (2) unit(s), obtained by integrating out the $W$ boson and $b$ quark (and, in general, any other massive particles). The absorptive part of the second contribution is the origin of $\Gamma_{12}$; the first term and the dispersive part of the second both contribute to $M_{12}$. The first contribution is local, stemming from processes in which all particles in Fig. 1 (and any in diagrams from new physics) propagate distances of $m_{b}^{-1}$ or less. In the second contribution, however, intermediate states, for example $K^{+}\pi^{-}$, can propagate a distance of order $\Lambda_{\text{QCD}}^{-1}$. This second contribution is difficult to compute because several multi-hadron intermediate states enter the sum, but not so many that an appeal to quark-hadron duality is likely to be successful.

It is easy to see via the conventional parametrization of the CKM matrix that the Standard Model predicts the angle $\phi_{12}$ to be very small: the imaginary parts of $M_{12}$ and $\Gamma_{12}$, and hence their phases, come from parts of the box diagrams carrying one or two factors of $\sin\theta_{23}\sin\theta_{13}\sin\delta=1.4\times 10^{-4}$, while the corresponding CKM factor in the real parts is $\sin\theta_{12}=0.225$. Because, with these conventions, both phases are small, so is the convention-independent difference $\phi_{12}$. For the same reason, the Standard Model real part of $M_{12}$ stems mostly from the long-distance contribution, the sum in Eq. (10). An estimate from a dispersion relation based on heavy quark effective theory yields Falk et al. (2004)

[TABLE]

which are compatible with the measurements, Eqs. (7) and (8).

Physics beyond the Standard Model could change $M_{12}$, $\Gamma_{12}$, or both. Many Golowich et al. (2009) extensions of the Standard Model alter only $\mathcal{H}^{\Delta C=2}$ and, hence, the magnitude and phase of $M_{12}$, but not $\Gamma_{12}$. In general, the $\Delta C=2$ effective Hamiltonian can be written

[TABLE]

where $C_{i}$ are the Wilson coefficients and the $\mathcal{O}_{i}$ are four-quark operators, given below. The determinations of $x_{12}$ and $\phi_{12}$ therefore can constrain the Wilson coefficients once the hadronic matrix elements $\langle D^{0}|\mathcal{O}_{i}|\bar{D}^{0}\rangle$ have been computed in (nonperturbative) QCD. Unfortunately, in view of the large range in Eq. (12), the constraint from $x$ will remain loose until new techniques for the long-distance term have been developed.

The operators in the $\Delta C=2$ effective Hamiltonian are

[TABLE]

where $\bar{c}$ and $u$ are the anticharm- and up-quark fields, with left and right Dirac projection matrices $L=\frac{1}{2}\left(1-\gamma_{5}\right)$ and $R=\frac{1}{2}\left(1+\gamma_{5}\right)$. The color indices are denoted $\alpha$ and $\beta$, and the spin indices are implied. All other Lorentz invariant four-quark operators can be reduced to this set by Fierz rearrangement Bouchard (2011). Further, parity conservation in QCD implies $\langle D^{0}|\mathcal{O}_{i}|\bar{D}^{0}\rangle=\langle D^{0}|\tilde{\mathcal{O}}_{i}|\bar{D}^{0}\rangle$, $i=1,2,3$.

In summary, the five matrix elements $\langle D^{0}|\mathcal{O}_{i}|\bar{D}^{0}\rangle$ suffice to obtain

[TABLE]

where the $C_{i}^{\text{NP}}(\mu)$ are the Wilson coefficients of the new-physics model, subsuming $\tilde{C}_{i}^{\text{NP}}$, renormalized in the same scheme as the matrix elements. They can be calculated in lattice QCD with the same methods as for $B^{0}_{(s)}$-$\bar{B}^{0}_{(s)}$ mixing Bazavov et al. (2016).

Given these matrix elements, Eq. (22) can be used to test models in which new physics does not change the phase of $\Gamma_{12}$. 222If the phase of $\Gamma_{12}$ were to change significantly, it would no longer be acceptable to treat the relative phase $\phi_{12}$ as the phase of $M_{12}$. A convenient way to do so is illustrated in Fig. 2, which plots $|x_{12}|e^{i\phi_{12}}$ as a complex number.

The colored contours show the fit results to the experimental data, while the gray horizontal bar shows the range given in Eq. (12). The gray bar is close to $\sin\phi_{12}=0$ and extends well beyond the range displayed here. With the proviso that new physics does not change the phase of $\Gamma_{12}$, the new-physics calculation from Eq. (22) yields a complex number $x_{12}^{\text{NP}}\approx|x_{12}^{\text{NP}}|e^{i\phi_{12}^{\text{NP}}}$, which should be added to the gray bar. If the total $|x_{12}|e^{i\phi_{12}}$ lands entirely outside the contours, the new-physics model is ruled out.

III Lattice simulation

In this section, we provide details of the numerical lattice calculations of the $D$-mixing matrix elements. We begin in Sec. III.1 with an overview of the gauge-field configurations and valence light- and charm-quark actions, and then define the two- and three-point correlation functions calculated in Sec. III.2. We conclude in Sec. III.3 with a discussion of statistical errors and autocorrelations.

III.1 Gauge-field configurations, light-, and heavy-quark actions

We use isospin-symmetric gauge-field configurations generated by the MILC Collaboration Bernard et al. (2001); Aubin et al. (2004a); Bazavov et al. (2010) with $N_{f}=2+1$ dynamical quarks; the degenerate up and down sea-quark masses span a range of values from $(0.4-0.05)\times m_{s}$, permitting a controlled chiral extrapolation to the physical value, while the strange sea-quark mass is close to the physical value. These ensembles employ the asqtad-improved staggered quark action, and have light-quark discretization errors of $\text{O}(\alpha_{s}a^{2},a^{4})$ Blum et al. (1997); Orginos and Toussaint (1998); Lagaë and Sinclair (1998); Lepage (1999); Orginos et al. (1999); Bernard et al. (2000a). The MILC asqtad ensembles were generated using the fourth-root procedure to yield a theory with one taste; both theoretical and numerical evidence indicate that the continuum limit of the rooted, staggered theory is indeed QCD Dürr (2006); Sharpe (2006); Kronfeld (2007); Golterman (2008). The gluons are simulated with the tadpole-improved Lüscher-Weisz action and have discretization errors starting at $\text{O}(\alpha_{s}a^{2},a^{4})$ Weisz (1983); Weisz and Wohlert (1984); Lüscher and Weisz (1985a, b). The lightest simulated pion mass $M_{\pi}=177$ MeV is very close to the physical value; heavier pion masses up to $M_{\pi}=555$ MeV are also included in the analysis to help guide the chiral extrapolation. Four lattice spacings ranging from $a\approx(0.12-0.045)$ fm provide good control over the extrapolation to the continuum. Figure 3 visually summarizes the range of pion masses and lattice spacings analyzed in this work. All ensembles have sufficiently large spatial lattice volumes ($M_{\pi}L\gtrsim 3.8$) that finite-volume effects are expected to be at the subpercent level Dürr et al. (2008). All ensembles have at least 500 configurations, and many have over 2000 configurations, yielding small statistical uncertainties. The set of ensembles used and simulation parameters are shown in Table 1.

The lattice spacing can be converted to $r_{1}$ units by multiplying with appropriate powers of the mass-independent ratio of scales $r_{1}/a$ Bazavov et al. (2010) listed in Table 1. Note that the ratios depend only on $\beta$, which through dimensional transmutation is linked to the lattice spacing. The scale $r_{1}$ is defined via the force $F(r)$ between two static quark by the condition $r_{1}^{2}F(r_{1})=1.0$ Sommer (1994); Bernard et al. (2000b). The relative scale is determined on each ensemble from the heavy-quark potential and then fit to a smoothing function in order to reduce sensitivity to lattice spacing Bazavov et al. (2010), and can be obtained with tiny statistical errors. We convert our final matrix-element results to physical units using the absolute scale Bazavov et al. (2012a),

[TABLE]

based on calculations of light pseudoscalar-meson decay constants from MILC Bazavov et al. (2009) and HPQCD Davies et al. (2010). A detailed discussion of the smoothing procedure may be found in Sec. IV B of Ref. Bazavov et al. (2012a).

Although the neutral $D$-meson has a valence up quark, in our simulations we generate correlation functions with seven to eight different light valence-quark masses on each ensemble, using the same asqtad action as for the sea quarks. The valence-quark masses employed in our simulations are given in Table 2, and correspond to pion masses from $M_{\pi}\approx 180$–880 MeV, enabling good control over the chiral extrapolation guided by SU(3) chiral perturbation theory. In Table 2 and throughout this work we denote the simulated valence light quark as $q$ with mass $am_{q}$, and reserve $m_{u}$ for the mass of the physical up quark.

For the heavy-quark propagators, we use the Sheikholeslami-Wohlert action with the Fermilab interpretation Sheikholeslami and Wohlert (1985); El-Khadra et al. (1997). The couplings in the theory are smoothly bounded for arbitrary values of the heavy-quark mass $am_{Q}$. After tree-level matching to QCD through heavy-quark effective theory (HQET), discretization errors from the action are of $\text{O}(\alpha_{s}a\Lambda_{\text{QCD}},a^{2}\Lambda_{\text{QCD}}^{2})$. The bare charm-quark mass in the Fermilab action is parametrized by the hopping parameter $\kappa_{c}$; the values employed in our simulations are tabulated in Table 2.

III.2 Lattice correlation functions

We calculate the two- and three-point correlation functions needed to obtain the matrix elements for neutral $D$-meson mixing using the same method as for $B$-meson mixing in Ref. Bazavov et al. (2012b). In particular, we first construct a specific combination of a light-quark propagator, heavy-quark propagator, and $\gamma_{5}$ with free spin and color indices known as an “open-meson propagator”. We then obtain the $D$-meson two-point correlators from a single open-meson propagator contracted with $\gamma_{5}$, and obtain the three-point correlation functions for the five $\Delta C=2$ four-quark operators from combining two open-meson propagators contracted with the appropriate Dirac structures. Here we describe the light- and heavy-quark propagators used to construct the open-meson propagators.

The valence light-quark propagators are generated using the asqtad action. We then construct the naive field $\Upsilon(x)$ from the staggered field $\chi(x)$ following Refs. Wingate et al. (2003); Bernard (2013),

[TABLE]

where $\Omega(x)=\gamma_{1}^{x_{1}}\gamma_{2}^{x_{2}}\gamma_{3}^{x_{3}}\gamma_{4}^{x_{4}}$ is the Kawamoto-Smit transformation and $\underline{\chi}$ denotes a vector of the four staggered copies. We remove tree-level, $\text{O}(a)$ discretization errors from the four-fermion operator by rotating the heavy-quark field $\psi(x)$ following Ref. Bazavov et al. (2012b),

[TABLE]

The simulation values of the rotation parameter $d_{1}$ are given in Table 2.

We form the $D$-meson interpolating operator from the light naive field and rotated heavy field as

[TABLE]

where $S(\bm{x},\bm{x}^{\prime})$ is a spatial smearing function. To improve the overlap with the $D$-meson ground state, we employ for the smearing function the 1S wavefunction of the Richardson potential Richardson (1979); Menscher (2005).

We obtain the $D$-meson mixing matrix elements from the zero-momentum three-point correlation functions,

[TABLE]

where the index $i=$1–5 labels the four-quark operator. We obtain the lattice operators, $O_{i}$, from the expressions (14)–(18) for the continuum operators, $\mathcal{O}_{i}$, by replacing the $u$ field with $\Upsilon$ and the $\bar{c}$ field with $\bar{\Psi}$. As shown in Fig. 4, the four-quark operator location is fixed at time $t_{0}$, while the $D$- and $\bar{D}$-meson times $t_{y}$ and $t_{x}$ vary. The construction of the three-point correlation function, as a result of the staggered formulation, introduces mixing between wrong-spin taste-mixing terms as discussed in detail in Sec. III B of Ref. Bazavov et al. (2016). This mixing is a lattice-discretization effect of $\text{O}(a^{2})$. The “wrong-spin” contributions are included in the chiral-continuum fit function, and hence removed when we take the continuum limit.

In order to extract the hadronic matrix element from the amplitude of the three-point correlator, we normalize the three-point correlators using the overlap function determined from the two-point correlator,

[TABLE]

III.3 Statistics and autocorrelations

We take advantage of the large temporal extents of the MILC lattices by computing the two- and three-point correlation functions with four to eight evenly-spaced time sources per configuration. Prior to analysis, the correlators are shifted to a common $t_{0}=0$, then averaged. Because the correlators from different time sources are only weakly correlated, this reduces the statistical errors by approximately a factor of $\sqrt{N_{\rm src}}$.

Because the lattices were generated with periodic boundary conditions, we gain another approximate factor of two in statistics by folding the data along the temporal midpoint $T/2$ so as to include the backward propagating signal. For the two-point correlator, we identify $t$ with $T-t$ and use the range $0\leq t\leq T/2$ in our correlator fits. For the three-point correlator shown in Fig.  4, we identify $-|t_{x}|$ with $-T-|t_{x}|$ and $t_{y}$ with $T-t_{y}$, and restrict our fit region to the $|t_{x}|<T/2$ and $t_{y}<T/2$ quadrant of the $|t_{x}|-t_{y}$ plane.

For the three-point correlation functions, we also exploit the parity symmetry of QCD to further increase statistics by averaging the matrix elements of parity-equivalent operators. The operators $\mathcal{\tilde{O}}_{1,2,3}$ in Eqs. (14)–(18) differ from $\mathcal{O}_{1,2,3}$ by an interchange of $L\leftrightarrow R$, and thus transform into each other under parity inversion. Consequently, the lattice matrix elements of these operators are equal up to statistical fluctuations. We therefore generate data for both $\langle\mathcal{{O}}_{1,2,3}\rangle$ and $\langle\mathcal{\tilde{O}}_{1,2,3}\rangle$ and take their average in order to gain an approximate factor of two in statistics.

We reduce autocorrelations between measurements computed on configurations close in Monte-Carlo simulation time by translating each gauge-field configuration by a random spatial shift $\vec{x}$ before calculating the valence light- and charm-quark propagators. We do not observe any remaining autocorrelations in the two- and three-point correlator data after this procedure. Figure 5 shows the scaled $D$-meson two-point correlator versus bin size on the coarsest and finest ensembles with $m_{l}^{\prime}/m_{s}^{\prime}=1/5$. The central values and errors are stable with increasing bin size. We also compute the integrated autocorrelation time, defined in Eq. (4.1) of Ref. Bazavov et al. (2016), and find it to be less than 1 on all ensembles. Based on these studies, we do not bin the data in this work.

IV Correlator analysis

In this section, we describe how we extract the $D$-mixing matrix elements from the two- and three-point correlation functions given in Eqs. (27) and (28). All dimensionful quantities are given in lattice spacing units, where the factors of $a$ are suppressed for simplicity.

IV.1 Correlator fit functions

Our choice of using the naive field $\Upsilon(x)$, defined in Eq. (24), for the valence up quark in the two- and three-point correlation functions results in contributions from $D$-meson states with positive parity Wingate et al. (2003). Their effects are included in the functional forms used to describe the Euclidean time-dependence of the correlation functions. The two-point function takes the form

[TABLE]

where the $Z_{n}$ are the overlap factors of the interpolating operators with the energy eigenstate labeled by $n$, $T$ is the temporal extent of the lattice, and terms with even (odd) $n$ describe the effects of negative (positive) parity states. The second term on the right-hand-side of Eq. (29) arises with periodic boundary conditions in time.

The three-point functions are described by

[TABLE]

where the desired matrix element is related to the ground-state amplitude $\mathcal{Z}_{00}^{O_{i}}$:

[TABLE]

The factor of the $D$-meson mass $M_{D}$ is due to the nonrelativistic normalization of states in Eq. (30). Effects from periodic boundary conditions are negligible in our three-point data and are therefore not included in Eq. (30).

IV.2 Method

As shown in Eq. (31), the three-point correlation function contains the desired matrix elements, however it also receives contributions from excited states. Our procedure for extracting the ground state matrix element is designed to account for the presence of excited states in the correlation functions while controlling and including the systematic errors associated with their residual contributions. We use Bayesian constraints with Gaussian priors in our fit functions in order to obtain a robust estimate of the uncertainty arising from excited state contamination. Correlations between two- and three-point functions are accounted for by performing a simultaneous fit to both data sets.

We implement the Bayesian constraints by minimizing the augmented $\chi^{2}$ function defined in Eq. (B3) of Ref. Bazavov et al. (2016), while the $Q$ parameter defined in Eq. (B4) of Ref. Bazavov et al. (2016) is used to assess the quality of the fits. The selection of the Bayesian priors is described in Sec/ IV.3. In the fits to the correlation functions we use intervals $t_{\text{min}}\leq t\leq t_{\text{max}}$ that do not span the entire time range. The time region used for the three-point functions is further restricted, as described in Sec. IV.4. We limit the number of states used in Eqs. (30) and (29) to $N^{\text{3pt}}=N^{\text{2pt}}\leq 6$. These choices are designed to optimize the extraction of the desired matrix elements from the correlation function data. Section IV.5 describes fit variations with different $t_{\text{min}}$, $t_{\text{max}}$, and $N^{\text{3pt}}$ to ensure our fit choices do not bias the fit results for the ground state parameters.

IV.3 Prior selection

The ground state parameters in the fit functions are well determined by the correlation function data. We therefore make sure that our selection procedure for the associated priors imposes only very loose constraints. At the same time, the data often do not provide good constraints for the excited state parameters in the fit functions. The purpose of the priors for the excited state parameters is to stabilize the fits and allow us to include the uncertainty due to residual excited state effects in the fit error. As discussed in Sec. IV.5.2, our main analysis employs two-and three-point functions with smeared $D$-meson operators. Below we describe the selection procedure and resulting choices for the priors and widths for these correlation functions. We start with the two-point function parameters, first for the ground state, followed by the excited states, and then discuss the three-point function parameters.

The priors for the ground state energy $E_{0}$ are obtained by examining the effective mass at large times $t$,

[TABLE]

A typical effective mass plot is shown in the left panel of Fig. 6. There is a clear plateau at large $t$ before the signal is overwhelmed by noise. We consider the different plateau values of the effective mass for the different light valence- and sea-masses, and determine a single $E_{0}$ prior distribution for each lattice spacing. Specifically, the central value of $E_{0}$ is chosen to lie in the center of the different effective mass values, and the prior width spans approximately twice the range of effective mass plateau values. For the ground state $Z_{0}$ priors we examine the scaled two-point function,

[TABLE]

where $M_{\rm eff}$ is an estimate of $M_{\rm eff}(t)$ in the large $t$ limit. An example of the scaled two-point function is shown in the right panel of Fig. 6; the plateau region at large time separation is again used to construct the prior range. In order to constrain $Z_{0}$ to be positive definite, we parametrize the corresponding fit parameter as the square root of $Z_{0}$. As illustrated in the left (right) panel of Fig. 6, the prior widths for $E_{0}$ ($Z_{0}$) are typically more than 100 (20) times larger than the errors on the corresponding posteriors.

In our data, the first excited state corresponds to the opposite-parity (scalar) ground state. We parametrize its energy $E_{1}$ in terms of the splitting $\Delta_{1,0}$ from the ground state with

[TABLE]

and

[TABLE]

This ensures that the parameter space for energy splittings is positive and enforces the ordering $E_{k}>E_{j}$. The prior central value and width for $\Delta_{1,0}$ is guided by the experimentally measured $D^{*}_{0}(2400)^{0}-D^{0}$ mass difference Patrignani et al. (2016); we use $\Delta_{1,0}\approx 450^{+400}_{-100}$ MeV. As illustrated in both panels of Fig. 6, oscillating states are clearly present in the correlation function data. We therefore expect the overlap with the opposite parity state ($Z_{1}$) to be nonzero. We set the prior central value for $Z_{1}$ to half the value of $Z_{0}$ because the smeared operators suppress the overlap with excited states (see Sec. IV.5.2). The prior width of $Z_{1}$ is set to be about one-$\sigma$ away from zero, reflecting our expectation for nonzero overlap.

For the remaining higher-state energies, we construct two towers of prior central values (and widths) starting at the pseudoscalar and scalar ground states,

[TABLE]

The choice of prior for the splittings of the radial excitations is motivated from quark model estimates Ebert et al. (2010), $\Delta_{k,k-2}\approx 650^{+2000}_{-500}$ MeV. The central values of the priors for the excited state overlap factors $Z_{k}$ are set to approximately half the central values of $Z_{0}$, again based on the expectation that smearing suppresses overlap with the excited states. The prior width is chosen to encompass $Z_{k}=0$, allowing for the possibility of an absence of signal in the data.

For the ground-state amplitude $\mathcal{Z}_{00}^{O_{i}}$ we examine the scaled three-point function at large times $t_{x},t_{y}$,

[TABLE]

an example of which is shown in the two panels of Fig. 7. The left panel illustrates the plateau along the diagonal where $|t_{x}|=t_{y}$, while the right panel shows the off-diagonal $|t_{x}|=t_{y}+1$. As with the two-point correlators, we obtain a prior central value and width at each lattice spacing from the variation of the scaled three-point functions with light valence- and sea-quark mass. As illustrated in both panels of Fig. 7, the resulting prior widths are typically 10 times larger than the fit errors. We also increased the width of the ground state prior by up to a factor of 3 and obtained identical results for the ground state amplitudes. The exponential decay visible in both panels is the dominant signal for excited states. The priors for the excited state amplitudes $\mathcal{Z}_{mn}^{O_{i}}$ (defined in Eq. (30)) are constrained to be the same order of magnitude as the ground state amplitude. The complete set of priors for the correlator analysis is listed in Table 15 in Appendix B.

IV.4 Fit region

The time regions over which we fit the two- and three-point function data is illustrated by the blue-colored areas in Fig. 8, where the horizontal and vertical red lines identify the extremum of the fit regions. Because the statistical errors are small for most of the available time range, there is a large region in the $|t_{x}|-t_{y}$ plane that is useful, in principle, for the three-point function analysis. However, the number of configurations in the ensembles used in this analysis while large, is not enough to resolve the covariance matrix over the entire time region bounded by $|t_{x}|,t_{y}\in\{t_{\text{min}},\ldots,t_{\text{max}}\}$, and we therefore constrain the three-point function analysis to the data along a bi-diagonal in the $|t_{x}|-t_{y}$ plane. Figure 7 illustrates that information along the off-diagonal direction is needed to account for excited state contributions that would not be easily resolved if the data were limited to just the diagonal. Another choice is to reduce $t_{\text{max}}$ for the three-point function but analyze the data in the entire region bounded by $|t_{x}|,t_{y}\in\{t_{\text{min}},\ldots,t_{\text{max}}\}$, i.e. including all the off-diagonal points in the $|t_{x}|-t_{y}$ plane. However, this would discard data at large time separations, where the ground state contribution dominates. We also consider other data reduction procedures, such as randomly drawing a certain number of data points from the $|t_{x}|,t_{y}\in\{t_{\text{min}},\ldots,t_{\text{max}}\}$ region. We find that all the choices for reducing the number of data points that we have studied yield results for the ground state parameters that are consistent with those from our main analysis. Table 3 lists our choices for $t_{\rm min}$ and $t_{\rm max}$ for our preferred fits. Across all four lattice spacings, we set $t_{\text{min}}=0.72$ fm while varying $t_{\text{max}}$ smoothly from 3 fm on the $a\approx 0.12$ fm ensembles, 2.7 fm on the $a\approx 0.09$ fm ensembles, 2.5 fm on the $a\approx 0.06$ fm ensembles, to 2.3 fm for the $a\approx 0.045$ fm ensemble. We choose one set of priors and time ranges for all correlation function fits at a given lattice spacing.

IV.5 Fit stability

Our preferred correlator fits are performed with the priors listed in Table 15, the time ranges listed in Table 3, the time region as described in Fig. 8, and with $N_{\text{state}}=2+2$ where the notation is used to denote the correlator fit ansatz with two normal parity and two opposite parity states. Here we describe fit variations that we use to test for systematic effects due to excited states, and other fit choices. We examine the dependence of the fit results for the ground state parameters under varying fit range, number of states, and operator smearing. In our simultaneous fits we vary each of the parameters ($t_{\text{min}}$, $t_{\text{max}}$, and $N^{\text{2pt}}=N^{\text{3pt}}$) for the two- and three-point functions.

IV.5.1 Fit range and number of states

Figure 9 (left) shows a typical example of a $t_{\text{min}}$ stability plot for $\mathcal{Z}_{00}^{O_{1}}$ with the preferred fit displayed with a solid blue point. The fits that include only 1+1 states show a strong $t_{\text{min}}$ dependence before reaching a plateau at $t_{\text{min}}\geq 9$, while fits with 2+2 or 3+3 states reach a plateau for $t_{\text{min}}\geq 3$. In addition, the 1+1 fit results have significantly smaller error bars than the 2+2 and 3+3 fits, while the errors going from the 2+2 to the 3+3 fits are essentially unchanged. We conclude that 2+2 fits with $5\leq t_{\text{min}}\leq 9$ are necessary and sufficient to account for excited state effects.

A typical example of a $t_{\text{max}}$ stability plot is shown Fig. 9 (Right) for $\mathcal{Z}_{00}^{O_{i}}$ where the preferred fit is marked with a solid blue point. The fit results (central values and error bars) do not change as $t_{\text{max}}$ is increased, indicating that contributions to the three-point function from periodic boundary conditions are negligibly small. However, the drift in the central values and increase in the error bars as $t_{\text{max}}$ is decreased to $t_{\text{max}}\lesssim 26$ indicate that the correlation function data at large time separations still contribute to the ground state signal and help stabilize the ground state posteriors. This informs our choices for $t_{\text{max}}$ in Table 3.

IV.5.2 Operator smearing

Figure 10 (left) shows an example comparing fit results for the two-point function ground state overlap factor $Z_{0}^{1/2}$ with 1S-smeared (blue circles) and local (red squares) source and sink operators. The 1S-smeared operator has better overlap with the ground state as evidenced by the larger central values for the corresponding overlap factor. In conjunction, the 1S-smeared operator has smaller overlap with excited states than the local operator, since the latter yields fit results that exhibit significant dependence on $t_{\text{min}}$ not seen with the former. A sample comparison of the ground state energies $E_{0}$ from fits to the same two-point functions is shown in Fig. 10 (right). The results are similar. The ground state energies obtained from the 1S-smeared two-point function do not vary with $t_{\text{min}}$, and have errors that stabilize for $t_{\text{min}}\gtrsim 5$. By comparison, the ground state energies from the local two-point function show strong $t_{\text{min}}$ dependence and differ at small $t_{\text{min}}$ from the ground state energies with 1S-smearing before they become consistent with them. The larger range of stability in $t_{\text{min}}$ observed for fit results from correlation functions with 1S-smeared source and sink operators makes it easier to obtain consistent fit regions across all ensembles and valence masses. We therefore use the correlation functions with 1S smeared sources and sinks in our main analysis.

IV.6 Error propagation

We propagate the distribution of the matrix element $\mathcal{Z}_{00}^{O_{i}}$ by bootstrap resampling. The distributions of the priors are included by randomizing the prior central value over the prior width under the bootstrap sampling Lepage et al. (2002). We generate 2000 bootstrap samples for each ensemble included in our analysis. An example comparing the bootstrap, single elimination jackknife, and Hessian posterior distributions of the ground state amplitude $\mathcal{Z}_{00}^{O_{1}}$ is shown in Fig. 11. This example is representative of the consistency seen amongst the bootstrap, jackknife, and Hessian distributions of our fit results and demonstrates that the posterior distributions are approximately Gaussian.

V Renormalization

In lattice-QCD calculations, the nonzero lattice spacing provides an ultraviolet cutoff. As a result, the matrix elements computed at different lattice spacings are regulated at different energy scales. In order to take the continuum limit, the matrix elements must be run to the same energy scale, and in order to combine lattice matrix-element results with continuum Wilson coefficients, they must be matched to a continuum renormalization scheme. In this paper we convert the bare lattice operators $O_{i}(a)$ evaluated at scale $a$, to the renormalized operators $\bar{O}_{i}(\mu)$ in the continuum $\overline{\text{MS}}$-scheme evaluated at the scale $\mu=3~{}\text{GeV}$. In noninteger dimensions, the Dirac algebra is infinite dimensional and is fully defined only after choosing a basis of evanescent operators Dugan and Grinstein (1991). This choice, which here only affects the renormalization of operators $\mathcal{O}_{2}$ and $\mathcal{O}_{3}$, is not unique; here we consider the schemes of Beneke et al. (BBGLN) Beneke et al. (1999, 2002) and of Buras, Misiak, and Urban (BMU) Buras et al. (2000), reporting results for both.

The $\Delta C=2$ four-fermion operators mix under renormalization. At $\text{O}{(\alpha_{s})}$, the renormalized operators are given by Evans et al. (2009); El-Khadra et al. (2017)

[TABLE]

where the $P_{p}$ are dimension-seven operators that are not needed in this paper, because with the $O_{i}$ described after Eqs.(27), the coefficient matrix $b_{ip}$ starts at order $\alpha_{s}$. Neglecting this contribution leads to a discretization error of order $\alpha_{s}a\Lambda_{\text{QCD}}$, the same as that from our choice of $c_{\text{SW}}$, which is the coefficient of the clover term in the Sheikholeslami-Wohlert action Sheikholeslami and Wohlert (1985).

We calculate the renormalization coefficients $Z_{ij}(a\mu)$ using the “mostly nonperturbative matching” (mNPR) method introduced in Refs. El-Khadra et al. (2001); Harada et al. (2002). We factor out the flavor-conserving renormalization coefficients $Z^{4}_{V_{qq}}$ using

[TABLE]

The $Z_{V^{4}_{qq}}$ are then calculated nonperturbatively from equal-mass vector current correlation functions, as discussed in Ref. Bazavov et al. (2012a). The flavor-changing coefficients $\rho_{ij}$ are calculated to one-loop in lattice perturbation theory via

[TABLE]

where the coefficients $Z^{[1]}_{V^{4}_{ll(cc)}}$ are defined through

[TABLE]

and the tree-level matching factors are $\mathcal{C}_{ll}=2u_{0}$ and $\mathcal{C}_{cc}=2\kappa^{\prime}_{c}u_{0}(1+m_{0}/u_{0})$ for our conventions for the staggered and clover fermion fields, respectively. The one-loop heavy-light four-fermion-operator-renormalization correction factor $\rho_{ij}$ is close to unity for $i=j$ because wavefunction renormalization diagrams cancel in the ratio $Z_{ii}/(Z_{V^{4}_{cc}}Z_{V^{4}_{ll}})$ at this order.

Table 4 lists the renormalization coefficients for flavor-conserving vector currents. The same value for $Z_{V_{ll}^{4}}$ is used for all valence-quark masses on a given ensemble because the observed quark-mass dependence is smaller than the statistical errors.

We calculate the $\rho$ factors in tadpole-improved lattice perturbation theory taking $\alpha_{s}=\alpha_{V}(q^{*})$ in the “$V$ scheme” Lepage and Mackenzie (1993) obtained from the static-quark potential Mason et al. (2005). We fix the scale to be $q^{*}=2/a$, which is the typical scale of the gluon loop momenta. Table 5 lists the one-loop coefficients $\rho_{ij}^{[1]}$ for the BBGLN choice of evanescent operators. The matching coefficients for the BMU evanescent operator prescription can be obtained from the BBGLN coefficients via Monahan et al. (2014); Bećirević et al. (2002)

[TABLE]

We also calculate the heavy-light matching coefficients $Z_{ij}$ at one loop in tadpole-improved perturbation theory. The tadpole-improved coefficients $\zeta^{[1]}_{ij}$ are related to the coefficients $Z^{[1]}_{ij}$ in Eq. (42) via

[TABLE]

We compute the renormalization coefficients taking the tadpole-improvement factor $u_{0P}^{[1]}=-0.76708(2)$ from the fourth root of the plaquette or $u_{0L}^{[1]}=-0.750224(3)$ from the link in Landau gauge. The numerical differences between the diagonal coefficients $\zeta^{[1]}_{ii}-\rho^{[1]}_{ii}$ obtained with one-loop perturbation theory using $u_{0P}$ and $u_{0L}$ and mNPR are given in the last two columns in Table 5, respectively. This difference is the same for all operators $i=1$–5.

[FIGURE:]

VI Charm-quark mass correction

The mass of the charm quark is set by the hopping parameter $\kappa$ in the Fermilab action. We determine the appropriate $\kappa_{c}$ by requiring that the $D_{s}$-meson kinetic mass obtained in our lattice simulations agree with the PDG value as described in Ref. Bernard et al. (2011); Bailey et al. (2014). In practice, initial low-statistics runs with several values of $\kappa$ were performed and used to determine the simulation values $\kappa_{c}^{\prime}$ for high-statistics data-production runs. The tuned value of the hopping parameter corresponding to the physical-charm quark mass, $\kappa_{c}$, was then determined using the high-statistics data. The physical $\kappa_{c}$ values on each ensemble are listed in Table 6.

We account for the slight difference between our simulation $\kappa_{c}^{\prime}$ and the physical $\kappa_{c}$ by incorporating a charm-quark mass correction into the chiral-continuum fit. To estimate this correction term, we start with the quark kinetic mass $am_{2}$, which is related to the hopping parameter as El-Khadra et al. (1997)

[TABLE]

where $am_{0}$ is the tadpole-improved bare-quark mass,

[TABLE]

and $\kappa_{\text{crit}}$ corresponds to the value of $\kappa$ where the mass of the pseudoscalar-meson mass vanishes. The nonperturbatively determined values of $\kappa_{\rm crit}$ are listed in Table 4. From heavy-quark power-counting, we expect the matrix elements to depend upon the heavy-quark mass as $1/m_{Q}$, which is identified with the kinetic mass in the Fermilab interpretation. We therefore adjust the matrix elements using a function linear in the inverse kinetic quark mass. We first compute on each ensemble the difference between the simulated and tuned inverse kinetic mass

[TABLE]

these values are given in Table 6. We then determine the slope of the matrix elements with respect to $1/m_{2}$ using data with $\kappa_{c}^{\prime}=0.1254$ and 0.1280 on the 0.12 fm, $m_{l}/m_{s}=0.2$ ensemble with valence-quark masses $m_{q}=0.0100$ and 0.0349. Table 7 gives the slopes

[TABLE]

obtained for $i$=1–5 from an unconstrained linear fit of the renormalized matrix elements $\left<\mathcal{O}_{i}\right>/M_{D}$ in $1/m_{2}$, while here $\Delta(1/am_{2})$ is the difference between the two simulated inverse kinetic masses. Figure 12 shows an example fit for $\left<\mathcal{O}_{1}\right>$.

We add the charm-quark mass correction to the chiral-continuum fit as a constrained fit parameter. This allows us to propagate the error stemming from the uncertainties in $\Delta(1/(am_{2}))$ and $\mu_{i}$ into the matrix-element uncertainties reported by the fit. We introduce priors for $\Delta(1/(am_{2}))$ with the central values and widths given in Table 6, where the errors are obtained by propagating the error from $\kappa_{c}$ through Eq. (46) and (47). We take the priors for $\mu_{i}$ from the linear fits described above and illustrated in Fig. 12; the prior values employed for $\mu_{i}$ are tabulated in Table 7.

VII Chiral and continuum extrapolation

VII.1 Chiral fit function

We extrapolate our renormalized lattice matrix-element results to the physical light-quark mass and continuum limit using SU(3), partially-quenched, heavy-meson, rooted staggered chiral perturbation theory (HMrS$\chi$PT) Aubin and Bernard (2006); Bernard (2013). HMrS$\chi$PT describes the dependence of the matrix elements on the light-quark masses and on the lattice spacing from taste-symmetry breaking in the staggered action. To incorporate additional systematics into the chiral-continuum fit, we supplement the HMrS$\chi$PT expression with terms to parametrize discretization errors from the heavy-quark, light-quark, and gluon actions, the uncertainty in the adjustment from the simulation to physical charm-quark mass, and higher-order terms in the operator matching procedure. Schematically, the fit function takes the form

[TABLE]

Because the lattice spacings on all ensembles differ, we bring all lattice masses and matrix elements into $r_{1}$ units during the physical point extrapolation. We discuss each term in turn in the following subsections.

VII.1.1 Chiral logarithms

We work at next-to-leading order (NLO) in HMrS$\chi$PT. The one-loop chiral logarithms describe nonanalytic dependence on the light-quark masses and lattice spacing, and are

[TABLE]

where repeated indices $\xi$ are summed. The coefficients $\beta_{i}$ and $\beta_{i}^{\prime}$ are the leading-order low energy constants (LECs) for the mixing matrix elements $\left<\bar{D}|\mathcal{O}_{i}|D\right>$ and $\left<\bar{D}^{*}|\mathcal{O}_{i}|D^{*}\right>$, respectively. The terms, $\mathcal{W}$, $\mathcal{T}$, and $\mathcal{Q}$ are the one-loop wavefunction renormalization, tadpole, and sunset corrections; their explicit expressions are given in Eqs. (63), (82)–(83), and (89) of Ref. Bernard (2013), respectively. The implementation of the staggered action employed in our simulations introduces mixing between spin and taste degrees of freedom. To account for this, we include the wrong-spin taste-mixing terms $\tilde{\mathcal{T}}$ and $\tilde{\mathcal{Q}}$ in our chiral-continuum fit with coefficients $\beta_{i}^{(\xi)}$ and $\beta_{i}^{\prime(\xi)}$, where $\xi=\{\text{P, A, T, V, I}\}$ labels the different taste contributions. For the wrong-spin terms $\tilde{\mathcal{T}}$ and $\tilde{\mathcal{Q}}$, we follow a different notation than in Ref. Bernard (2013) in order to separate the LECs $\beta_{i}^{(\prime)}$ from the one-loop diagram functions. The relationships between $\tilde{\mathcal{T}}$ and $\tilde{\mathcal{Q}}$ in Eq. (51) and the chiral-logarithm functions in Ref. Bernard (2013) are given in Eqs. (C1a)–(C2e) of Ref. Bazavov et al. (2016). The coefficients $\beta^{(\xi)}_{i}$ and $\beta^{\prime(\xi)}_{i}$ are not all independent; for convenience, their relationships are given in Table 8. Inspection of Table 8 shows that the matrix elements $\{\left<\mathcal{O}_{1}\right>,\left<\mathcal{O}_{2}\right>,\left<\mathcal{O}_{3}\right>\}$ mix, as do $\{\left<\mathcal{O}_{4}\right>,\left<\mathcal{O}_{5}\right>\}$.

To account for the discrete momentum spectrum dictated by the finite lattice spatial size and periodic boundary conditions, we use the finite-volume expressions for the NLO chiral logarithms, which are obtained by replacing the integrals over loop momenta by discrete sums. The explicit expressions are given in Eqs. (63) and (64) of Ref. Bernard (2013).

We work at leading order in the heavy-meson expansion, but include the largest $1/M_{D}$ effects from the hyperfine splitting $\Delta^{*}\equiv M_{D^{*}}-M_{D}$ and the SU(3)-flavor splitting $\delta_{su}\equiv M_{D_{s}}-M_{D^{0}}$. The parameter that characterizes the flavor splitting in HMrS$\chi$PT is $\lambda_{1}\equiv\delta_{su}/\left(M^{2}_{\eta_{s}}-M^{2}_{\pi^{0}}\right)$, where $M_{\eta_{s}}$ is the mass of a theoretical $\bar{s}s$ bound state.

VII.1.2 Analytic terms in the chiral expansion

The analytic terms

[TABLE]

are simple polynomials in the light-quark masses and lattice spacing. At nonzero lattice spacing, taste-symmetry breaking splits the masses of pions with different taste representations $\xi=\{P,A,T,V,I\}$. The tree-level staggered $\chi$PT relationship between the mass $M^{2}_{ab,\xi}$ of a pion with taste $\xi$ and the constituent staggered quark masses $m_{a}$ and $m_{b}$ is

[TABLE]

where the $a^{2}\Delta_{\xi}$ are the taste splittings, and the leading-order LEC $B_{0}$ is related to the chiral condensate. Following Ref. Bazavov et al. (2012a), we construct the analytic terms in Eq. (52) using the dimensionless variables

[TABLE]

where $M_{qq}$ is the mass of the taste-pseudoscalar $\bar{q}q$ meson ($q=u,l,s$) and $f_{\pi}$ is the pion decay constant, and

[TABLE]

where $a^{2}\bar{\Delta}=\frac{1}{16}\sum_{\xi}w_{\xi}a^{2}\Delta_{\xi}$ is the average taste splitting, and the weight factors for $\xi=\{P,A,T,V,I\}$ are $w_{\xi}=\{1,4,6,4,1\}$, respectively. The coefficients of the analytic terms are expected to be of $\mathrm{O}(1)$ from chiral power counting when written in terms of $x_{q}$ and $x_{\bar{\Delta}}$.

The NLO, NNLO, and NNNLO analytic terms are obtained by forming all combinations of $x_{i}$ ($i=u,d,s,\bar{\Delta}$) linear, quadratic, and cubic in $x_{i}$, respectively. We include the full set of NLO and NNLO analytic terms in our base fit used to obtain our matrix-element central values. They are:

[TABLE]

where $c_{i}$ and $d_{i}$ are LECs of the theory. The inclusion of NLO terms is needed to absorb the dependence of the nonanalytic one-loop terms in Eq. (51) on the scale $\Lambda_{\chi}$ in the chiral logarithms, while the NNLO terms capture higher-order effects that might contaminate the lower-order LECs. Although the NNLO terms are not needed to obtain an acceptable fit, they improve the $\chi^{2}_{\rm aug}/{\rm dof}$. We also perform fits including NNNLO analytic terms to check for fit stability and look for truncation errors. The N3LO terms are:

[TABLE]

VII.1.3 Heavy quark discretization effects

We parametrize the leading heavy-quark discretization errors of $\text{O}(\alpha_{s}a,a^{2},a^{3})$ in our data by adding the terms $F_{i}^{\text{HQ disc.}}$ to our fit function:

[TABLE]

where

[TABLE]

and the $\beta_{i}$ are the same LECs as in Eq. (51). The fit parameters $z_{h}^{i}$ combined with powers of $a\Lambda_{\text{HQ}}$ represent the HQET matrix elements. The “mismatch functions” $f_{h}(m_{0}a)$ with $h\in\left\{B,3,E,X,Y,2\right\}$ are smoothly varying functions of the bare lattice heavy-quark mass that encapsulate the short-distance differences between the lattice and continuum-QCD action and operator descriptions.

The tree-level $\text{O}(a^{2})$ and $\text{O}(a^{3})$ mismatch functions of the action were calculated in Ref. Oktay and Kronfeld (2008) by performing a matching calculation between lattice HQET and continuum HQET. The tree-level $\text{O}(a^{2})$ mismatch functions of the spinor, and consequently, the 4-quark operators, were worked out in Ref. El-Khadra et al. (1997). For $\text{O}(\alpha_{s}a)$ errors, mismatches of the action and operator are modeled, using the Fermilab interpretation, by a smooth function that has the correct $am_{0}\rightarrow 0$ and $am_{0}\rightarrow\infty$ limits. A complete list of the functional forms of the $f_{h}(m_{0}a)$ employed here is given in Appendix A of Ref. Bazavov et al. (2012a).

We also consider generic $\text{O}{(\alpha_{s}a^{2})}$ discretization errors when studying the stability of our fits by including the term

[TABLE]

This term receives one-loop corrections from the light-quark and gluon actions, and from heavy-quark discretization errors of higher order than those in Eqs. (60)–(62).

VII.1.4 Heavy-quark mass adjustment

Following Sec. VI, we adjust the matrix elements for the slight difference between the simulation and physical charm-quark masses within the chiral-continuum fit by adding the correction term

[TABLE]

We propagate the uncertainties in the slopes $\mu_{i}$ and differences $\Delta\left(\frac{1}{r_{1}m_{2}}\right)$ to the final fit error by including them as constrained fit parameters with Gaussian priors.

VII.1.5 Renormalization errors

As discussed in Sec. V, we calculate the renormalization coefficients at one loop; therefore, truncation errors start at $\text{O}(\alpha_{s}^{2},\alpha_{s}\Lambda_{\text{QCD}}/m_{c})$. The $\text{O}{(\alpha_{s}^{2})}$ truncation errors are estimated in the chiral-continuum fit by adding the terms

[TABLE]

where the coefficients $\rho^{[2]}_{ii}$ and $\rho^{[2]}_{ij}$ are free parameters (with $i\neq j$), and the $\beta_{i}$ are the leading-order LECs for matrix elements $\left<\mathcal{O}_{i}\right>$ defined in Eq. (51). The second term in Eq.(65) parametrizes the mixing between operators $\mathcal{O}_{i}$ under renormalization. We evaluate the renormalized coupling $\alpha_{s}$ at $q^{*}=2/a$. In tests of fit stability, discussed below, we consider effects of $\text{O}{(\alpha_{s}^{3})}$ by adding the terms

[TABLE]

The coefficients are set to $\rho^{[3]}_{ii}=\rho^{[3]}_{ij}=0$ in our base fit.

Errors from omitted $\text{O}{(\alpha_{s}\Lambda_{\text{QCD}}/m_{c})}$ terms in the perturbative calculation are absorbed by the coefficients $z_{\{B,3\}}$ of the mismatch functions $f_{\{B,3\}}(m_{0}a)$, which also scale as $a\Lambda_{\text{QCD}}/m_{c}$

VII.2 Chiral-continuum fit parameters

In the following section, we discuss the priors chosen for the chiral-continuum fit parameters. Briefly, we employ loose constraints based on power counting for the coefficients associated with the chiral and perturbative expansions in $\alpha_{s}$, and for those of the heavy-quark discretization terms. We incorporate the parametric uncertainties from our fit inputs into the final fit error by including them as constrained fit parameters with Gaussian prior widths corresponding to the errors on the inputs. We fix the values of a small number of inputs for which the uncertainty contribution to the final fit is negligible.

VII.2.1 Loosely constrained fit parameters

The LECs of HMrS$\chi$PT in Eqs. (51) and Eqs. (56)–(58) are expected to be of $\text{O}{(1)}$. We constrain them only loosely to allow their values to be determined by the data; the priors improve fit stability when adding higher-order terms in the chiral expansion. The priors for the coefficients of the chiral-logarithms $\beta_{i}$ and $\beta^{\prime}_{i}$ are set to

[TABLE]

where the central values are rough guesses based on the correlator fits. We take very wide prior widths for the coefficients of the NLO analytic terms, which are well determined by the data:

[TABLE]

and use prior widths of $\text{O}{(1)}$ as motivated by chiral power counting for the coefficients of the NNLO (and N3LO) analytic terms:

[TABLE]

Recall that our base fit includes terms only through NNLO.

When we include the generic discretization term $F_{i}^{\alpha_{s}a^{2}\,\text{gen}}$ in our fit, we constrain its coefficient to be $h_{0,i}=0(1)$.

We choose priors for the heavy-quark discretization terms based on HQET power-counting. We expect the individual coefficients $z_{i}$ to be of $\text{O}(1)$. In some cases, however, more than one operator shares the same mismatch function. We therefore choose the width for each prior such that the width-squared equals the number of terms sharing the corresponding mismatch function. The priors for the $z_{i}$ are given in Table 9. The heavy-quark discretization terms also depend upon the scale $\Lambda_{\rm HQ}$, which is the cutoff of the effective theory. We use $\Lambda_{\text{HQ}}=500~{}\text{MeV}$ based on studying the lattice-spacing dependence of our full-QCD matrix-element data adjusted to the same sea-quark masses via the chiral-continuum fit. Our physical continuum-limit matrix-element results are insensitive to reasonable variations in $\Lambda_{\rm HQ}$.

The priors for the unknown two-loop coefficients $\rho^{[2]}_{ii}$ and $\rho^{[2]}_{ij}$ are set to

[TABLE]

When $\text{O}{(\alpha_{s}^{3})}$ terms are included, we use the same constraints for the higher-order coefficients $\rho^{[3]}_{ii}=0(1)$ and $\rho^{[3]}_{ij}=0(1)$. These values are consistent with the observation that the one-loop coefficients in Table 5 are at most of $\text{O}{(1)}$.

VII.2.2 Constrained fit parameters

As discussed in Sec. III.1, we bring our renormalized matrix-element results on all ensembles into the same units before the chiral-continuum extrapolation using the intermediate scale $r_{1}/a$. We also convert all fit inputs taken from experiment into $r_{1}$ units using the physical scale $r_{1}$ in fm. We propagate the uncertainties in $r_{1}/a$ values to our fit parameter by including them as constrained fit parameters with prior central values and widths given by their values and errors in Table 1. The values of $r_{1}/a$ are correlated between ensembles because they are obtained from a fit of the data on all MILC asqtad ensembles to a smooth function of the coupling $\beta$ Allton (1996); Bazavov et al. (2010) of Appendix A. We include the correlations between $r_{1}/a$ values in our fit; the correlation matrix is given in Table 14, while double-precision values for $r_{1}/a$ and its correlations are provided as supplementary material. For the physical scale $r_{1}$, we use the prior $r_{1}=0.3117(22)$ fm taken from Ref. Bazavov et al. (2010).

The coefficients of the one-loop chiral logarithms depend upon the light pseudoscalar-meson decay constant $f_{\pi}$ and the $D^{*}$-$D$-$\pi$ coupling $g_{D^{*}D\pi}$. We constrain $f_{\pi}$ to the PDG value of the pion decay constant Rosner et al. (2016)

[TABLE]

where the uncertainties on $f_{\pi}$ are from $\Gamma$, $|V_{ud}|$ and higher-order radiative corrections, respectively. The error on the fit input $r_{1}f_{\pi}$ includes the uncertainties from both $f_{\pi}^{\rm PDG}$ and $r_{1}$ added in quadrature. The coupling $g_{D^{*}D\pi}$ has been studied in unquenched lattice QCD with 2 flavors, yielding $g_{D^{*}D\pi}^{N_{f}=2}=0.53(3)(3)$ Bećirević and Sanfilippo (2013), and with three flavors, yielding $g_{D^{*}D\pi}^{N_{f}=2+1}=0.55(6)$ Can et al. (2013). Based on these results, we constrain the coupling in our fit to be

[TABLE]

which covers the 1$\sigma$ ranges of both.

The one-loop HMrS$\chi$PT chiral logarithms also depend upon the parameters $a^{2}\delta^{\prime}_{A,V}$, which multiply contributions from quark-disconnected “hairpin” diagrams. Because the hairpin contributions arise from taste-symmetry breaking, the coefficients $a^{2}\delta^{\prime}_{A,V}$ scale with lattice spacing approximately as $\alpha_{s}^{2}a^{2}$. We constrain their values in our fit at $a\approx 0.12$ fm to the determinations from the MILC Collaboration’s chiral-continuum fit of pion and kaon masses and decay constants in Ref. Aubin et al. (2004b):

[TABLE]

For the remaining lattice spacings, we scale these values by the weighted average of the taste splittings $\bar{\Delta}_{a}/\bar{\Delta}_{0.12~{}{\rm fm}}$.

We constrain the hyperfine and flavor splittings in our fit to the experimentally-measured values. For the $D$-meson system, $\Delta^{*}\equiv M_{D^{*}}-M_{D}=142.020(71)~{}\text{MeV}$ Patrignani et al. (2016), which corresponds to

[TABLE]

including the uncertainty on $r_{1}$. The SU(3) flavor splitting $\delta_{su}=M_{D_{s}}-M_{D_{\pm}}=98.69(13)$ MeV Patrignani et al. (2016). Combining this with $M_{\pi^{0}}=134.9766(6)$ MeV Patrignani et al. (2016) and $M_{\eta_{s}}=685.8(4.0)$ MeV Davies et al. (2010), we obtain

[TABLE]

We incorporate the uncertainty from the charm-quark-mass correction into the fit error by constraining the slopes $\mu_{i}$ and differences in inverse kinetic masses $\Delta(1/(r_{1}m_{2}))$ that enter the correction term $F_{i}^{\kappa}$, Eq. (64), with Gaussian priors. We fix the prior central values and errors to the results obtained from our fits of the charm-quark-mass dependence of the matrix elements in Sec. VI. The values are listed in Tables 6–7, and include the errors from statistics, fitting, and $r_{1}$.

We extrapolate the matrix elements to the physical light-quark masses given in Table VIII of Bailey et al. (2014),

[TABLE]

by evaluating the fit function with the valence-light quark mass fixed to $m_{u}$, and the light and heavy sea-quark masses fixed to $\hat{m}\equiv(m_{u}+m_{d})/2$ and $m_{s}$, respectively. The errors on the physical light-quark masses include statistics and the dominant systematic uncertainties from the chiral-continuum extrapolation $r_{1}$. Finally, we convert the chiral-continuum extrapolated matrix elements to the relativistic normalization by dividing by the experimental $D^{0}$-meson mass Patrignani et al. (2016)

[TABLE]

VII.2.3 Fixed inputs

We fix the pseudoscalar-meson taste splittings and the leading-order LEC $B_{0}$ in the tree-level expression for the squared pion mass, Eq. (53), in the fit because their uncertainties are negligible compared to other contributions to the error. We use the values given in Table 10, which were obtained from an analysis of the staggered light pseudoscalar-meson spectrum.

We fix the scale $\Lambda_{\chi}$ in the chiral logarithms to the experimental $\rho$-meson mass $M_{\rho}^{\rm PDG}=775~{}\text{MeV}$ Patrignani et al. (2016), which corresponds to a value in $r_{1}$ units of $\left(\Lambda_{\chi}/r_{1}\right)^{2}=1.5$. We have checked that our physical continuum-limit matrix-element results are insensitive to reasonable variations in $\Lambda_{\chi}$.

VII.3 Chiral-continuum fit results

Our base fit used to obtain our matrix-element central values is to the function

[TABLE]

and includes terms to account for errors from truncating the chiral expansion, discretization errors from taste-symmetry breaking, heavy-quark discretization errors, errors from omitted higher order terms in the renormalization factors, and errors in the charm-quark-mass correction factors. We fit the renormalized lattice matrix elements for all five operators simultaneously including statistical correlations between data on the same ensembles. This reduces the error on the physical continuum-limit matrix elements, which share common LECs in HMrS$\chi$PT and mix under renormalization. Figure 13 shows our preferred chiral-continuum extrapolation as a function of the squared valence-meson mass $M_{\bar{q}q}^{2}$, which is approximately linear in the valence light-quark mass $m_{q}$. We obtain a good fit with a correlated $\chi^{2}_{\text{aug}}/\text{dof}=122.4/510$, where the quantity $\chi^{2}_{\text{aug}}/\text{dof}$, which is suitable for assessing the quality of constrained fits, is defined in Eq. (7.27) of Ref. Bazavov et al. (2016).

VIII Systematic error analysis

We now discuss all sources of systematic error that contribute to our $D$-meson-mixing matrix-element uncertainties and provide complete error budgets. We begin, in Sec. VIII.1, with a discussion of errors that are included in the chiral-continuum fit. Next, in Sec. VIII.2, we estimate the remaining contributions that must be added to the fit error a posteriori to obtain the total error. The one exception is the error due to the omission of charmed sea quarks, which we present separately in Sec. VIII.3. Throughout these sections we refer to alternate fits that we used to test error saturation. Figures in Sec. VIII.4 show the results of these alternate fits and additional consistency checks of our fit results and error estimates. Complete error budgets for all five $D$-mixing matrix elements are listed in Table 12.

VIII.1 Base chiral-continuum fit errors

As described previously, we constrain the parameters in our chiral-continuum fit with Gaussian priors. This enables us to account for the uncertainties in our input parameters, as well as to include higher-order terms in the chiral and heavy-quark expansions thereby incorporating possible truncation errors. We use the dependence of the best-fit parameters on each piece of information, including correlations, to separate the total fit error into approximate suberrors. The approximate breakdown of the total fit error into the suberrors for each matrix element is shown in Table 11. The first column shows the statistical error, which is obtained from the quadrature sum of the errors from all data points. The other suberrors are discussed in the following Secs. VIII.1.1–VIII.1.7. The three dominant sources of error for all matrix elements are statistics, matching, and heavy quark discretization effects, each of which contribute a similar amount to the total errors. The errors from tuning the simulation $c$-quark masses and other inputs, from the extrapolation to the physical light-quark mass and the continuum, and from the relative lattice-spacing determination all contribute at the percent or subpercent level.

VIII.1.1 Parametric inputs

The “inputs” column of Table 11 is given by the quadrature sum of the error contributions from most of the input parameters (the error from $r_{1}/a$ is considered separately), which are constrained with Gaussian prior widths given by their estimated uncertainties. The largest contribution to the “inputs” error is from the uncertainty in the coupling $g_{D^{*}D\pi}$. The uncertainties from the parameters in the chiral logarithms $f_{\pi},\Delta^{*},\lambda_{1},\delta_{V}^{\prime}$, and $\delta_{A}^{\prime}$ are subdominant. The input errors also include the uncertainties from the physical $u$-quark and $D^{0}$-meson masses, which are used to fix the physical-point in the chiral extrapolation and to convert the matrix-element results to physical units.

The parametric error from the pion decay constant is already included in the “inputs” uncertainty. We also check for stability of the chiral-continuum extrapolation against reasonable changes in the decay constant, which provides a measure of the error due to truncating the chiral expansion. We replace $f_{\pi}$ in the coefficient of the chiral logarithms with the PDG value of $f_{K^{\pm}}=155.6(4)$ Rosner et al. (2016), which corresponds to

[TABLE]

This leads to only a tiny shift in the matrix elements, as shown by the fit variation labeled “$f_{K}$ vs. $f_{\pi}$” in Figs. 14 and 15, indicating that the fit error indeed encompasses the chiral truncation error.

VIII.1.2 Charm-quark mass uncertainty

We adjust the simulation charm-quark masses to the physical tuned values before the chiral-continuum fit as described in Sec. VI. We propagate the uncertainty in the charm-quark mass correction by including the matrix-element slopes $\mu_{i}$ ($i$=1–5) and the shift in the kinetic mass $\Delta(1/(r_{1}m_{2}))$ as constrained parameters with prior widths given by the uncertainties listed in Tables 6 and 7. The sum of uncertainty contributions from the fit parameters associated with the charm-quark mass adjustment are listed in the column “$\kappa$ tuning” in Table 11.

VIII.1.3 Renormalization and matching uncertainty

We include terms of $\text{O}{(\alpha_{s}^{2})}$ in our base chiral-continuum fit with unknown coefficients $\rho^{[2]}_{ij}$ constrained to be of $\text{O}{(1)}$ to incorporate the uncertainty due to omitted higher-order renormalization and matching terms. The sum of uncertainty contributions from the fit parameters $\rho^{[2]}_{ij}$ are listed in the “matching” column in Table 11. The renormalization and matching uncertainties estimated from fitting our lattice simulation data range from 2.0% to 4.1%. Their values are compatible with the naíve estimate obtained from taking $\alpha_{s}=0.2$ from our finest lattice spacing and $\rho^{[2]}_{ij}=1$, which yields 6.5% for all operators.

We check that the inclusion of generic $\mathrm{O}(\alpha_{s}^{2})$ terms captures the uncertainty from truncating the perturbative expansion in $\alpha_{s}$ by performing two alternate fits: one with only the known renormalization and matching terms of $\text{O}{(\alpha_{s})}$, and another with terms through $\text{O}{(\alpha_{s}^{3})}$ and coefficients constrained to be of $\text{O}{(1)}$. The results of these fits are labeled “mNPR” and “mNPR+$\alpha_{s}^{3}$”, respectively, in Figs. 14 and 15. In both cases, the shift in central value is small. Without the $\text{O}{(\alpha_{s}^{2})}$ terms, the errors on the matrix elements are underestimated. The errors do not increase from those of the base fit, however, with the inclusion of $\text{O}{(\alpha_{s}^{3})}$ terms. We therefore conclude that the base fit includes the uncertainty from omitted higher-order matching and renormalization terms.

We also compare the results of the base fit, in which the matrix elements are renormalized with the mNPR approach, to those from fits in which the matrix elements are renormalized using 1-loop tadpole-improved perturbation theory, and $\alpha_{s}$ is obtained either from the fourth-root of the plaquette or from the gauge-fixed Landau link. In these fits, labeled “$\text{PT}_{\rm P}+\alpha_{s}^{2}$” and “$\text{PT}_{\rm L}+\alpha_{s}^{2}$”, respectively, in Figs. 14 and 15, we include $\text{O}{(\alpha_{s}^{2})}$ terms constrained as in the base fit. Here we observe larger changes in the matrix-central values, which are still less than 2$\sigma_{\rm fit}$ away from the base-fit results. The “$\text{PT}_{\rm P}+\alpha_{s}^{2}$” and “$\text{PT}_{\rm L}+\alpha_{s}^{2}$” results themselves differ by almost 1$\sigma_{\rm fit}$, indicating a systematic uncertainty in the perturbative matching associated with the choice of tadpole-improvement factor. This supports our expectation that the mNPR matching approach is more reliable.

VIII.1.4 Truncation of the chiral and heavy-meson expansion

We estimate the uncertainty from truncating the chiral expansion by summing contributions to the matrix-element errors from the NLO LECs $\{\beta_{i},\beta_{i}^{\prime}\}$ and all of the analytic LECs $\{c_{n},d_{n}\}$ that do not depend on lattice spacing. The results are given in the column labeled “chiral” in Table 11.

We check the robustness of the fit error and look for residual truncation effects by considering two variations of the chiral fit function with different sets of analytic terms. In the first fit, we remove all NNLO analytic terms. For this fit, we also restrict the matrix-element data included to only those points for which the valence-quark mass $m_{q}<0.65m_{s}$, since we expect heavier-mass data to be outside the validity of NLO $\chi$PT. In the second fit, we add all possible analytic terms of O(N3LO). The results of these fits are labeled “NLO ($m_{q}<0.65m_{s}$)” and “N3LO”, respectively, in Figs. 14 and 15. In both cases, the central values for the matrix elements shift only slightly. The NLO $\chi$PT fit without NNLO analytic terms underestimates the errors on the matrix elements and also has a significantly larger $\chi^{2}_{\rm aug}/{\rm dof}$ than the other variations shown. The errors do not increase from those of the base fit, however, with the inclusion of N3LO analytic terms. We also consider a fit variation in which we set the hyperfine and flavor splittings in the chiral logarithms, which are the leading corrections in the $1/M_{B}$ expansion, to zero. The result is labeled “no splitting” in Figs. 14 and 15. Again, the changes in matrix-element central values and errors are small. All of these tests demonstrate that the base fit includes the uncertainty from higher-order terms in the chiral and heavy-meson expansions.

We also study the impact of our prior constraints on the $\chi$PT LECs, which are based on expectations from chiral power counting. We perform three fits in which we double the prior widths of (1) the LO LECs; (2) the NLO LECs; or (3) the NNLO LECs. The results are labeled “LO x 2, NLO x 2, NNLO x 2”, respectively, in Figs. 14 and 15. The errors on the matrix elements are stable against increasing the prior widths by a reasonable amount, indicating that the priors are sufficiently unconstraining to allow the data to determine the fit results.

VIII.1.5 Light quark discretization errors

We estimate the uncertainties from light-quark discretization via the uncertainty in the coefficients $\{c_{n},d_{n}\}$ of all analytic terms which depend on the lattice spacing. This error is labeled “LQ disc” in Table 11.

The base chiral-continuum fit function includes taste-symmetry breaking effects in the expressions for the chiral logarithms. Therefore corresponding analytic terms proportional to $a^{2}\bar{\Delta}$ are needed to maintain invariance under variation of the chiral scale. These terms scale as $\alpha_{s}^{2}a^{2}$. Generic one-loop contributions from improving the gluon and light-quark actions at tree-level, however, give rise to discretization terms of $\text{O}{(\alpha_{s}a^{2})}$. We therefore perform an alternate fit including such a term to account for generic light-quark and gluon discretization effects. The result is labeled “generic $\mathrm{O}(\alpha_{s}a^{2})$” in Figs. 14 and 15, and is indistinguishable from the base-fit result. This indicates that the terms already included in the base chiral-continuum fit function are sufficient to describe the lattice-spacing dependence of the data, and that the base-fit error properly includes light-quark discretization errors.

VIII.1.6 Heavy quark discretization errors

We estimate the uncertainties from heavy-quark discretization via the uncertainty in the coefficients $z_{i}$ of the $\text{O}{(\alpha_{s}a,a^{2},a^{3})}$ terms in Eqs. (60)–(62). This error is labeled “HQ disc” in Table 11. We also perform two fits, each of which includes fewer heavy-quark terms than in the base fit. These are labeled “HQ $\mathrm{O}(\alpha_{s}a)$ only” and “HQ $\mathrm{O}(\alpha_{s}a,a^{2})$ only” in Figs. 14 and 15, the label referring to the type of terms included in the fit. The tiny changes in matrix-element central values and errors confirm that the base-fit errors properly include the uncertainty from truncating the heavy-quark expansion.

As a consistency check, we can compare the heavy-quark discretization errors estimated from the data with those based on power counting. We evaluate Eqs. (60)–(62) taking $\Lambda_{\text{HQET}}=500$ MeV for the heavy-quark scale, and coefficients $z_{i}$ set by the combinatoric factors $z_{E}=2\sqrt{2},z_{B}=2,z_{3}=\sqrt{5},z_{X}=2\sqrt{2},z_{Y}=2,z_{2}=2$. Assuming that all contributions enter with the same sign, this leads to a conservative 5% estimate for all matrix elements. This is larger than the data-driven heavy-quark-discretization-error estimates in Table 11, which range from 2.4% to 3.6%, but the similar size suggests that errors obtained from the fit are reasonable.

VIII.1.7 Relative scale ($r_{1}/a$) uncertainty

The relative scale $r_{1}/a$ is used to convert lattice data on each ensemble to $r_{1}$ units before the chiral-continuum fit. We incorporate the uncertainty from $r_{1}/a$ through the use of prior constraints in the same manner as the parametric “input” errors. The relative scale errors are given in the column labeled “$r_{1}/a$” in Table 11.

VIII.2 Additional errors

VIII.2.1 Absolute scale ($r_{1}$) uncertainty

The absolute lattice scale $r_{1}$ in fm enters the matrix-element analysis during conversions between lattice-spacing and physical units. The scale $r_{1}$ is used to convert the PDG average meson masses and pion decay constant, which are parametric inputs to the chiral-continuum fit, from GeV to $r_{1}$ units. We account for the error on $r_{1}$ during the unit conversion. The scale $r_{1}$ is also needed to obtain the $D$-mixing matrix elements in GeV. Because the absolute scale does not affect the minimization of the $\chi^{2}$ statistic, we add the error from $r_{1}$ due to this final unit conversion in quadrature to the fit error a posteriori; the results are listed in column “$r_{1}$” of Table 12.

VIII.2.2 Finite-volume effects

We employ the finite-volume expressions for the 1-loop chiral logarithms in the base chiral-continuum fit. To estimate the systematic uncertainty from finite-volume effects, we perform a second fit using the infinite-volume expressions. The results are labeled “no FV” in Figs. 14 and 15. The observed shifts in the central value of the matrix elements are approximately half a percent. We take half the value of the matrix-element shifts as the error due to finite-volume effects, noting that this is conservative because we are in fact including NLO finite-volume corrections, and the omitted terms are of NNLO and hence even smaller. The finite-volume errors are listed in the “FV” column of Table 12.

VIII.2.3 Isospin breaking and electromagnetism

We obtain the $D$-meson matrix elements in the chiral-continuum limit by evaluating the valence light-quark mass at the physical $m_{u}$ and the light sea-quark mass at the isospin average $\hat{m}=(m_{u}+m_{d})/2$. This accounts for the dominant isospin-breaking effects from the valence sector, but isospin-breaking effects from the sea sector must be included as a systematic uncertainty. Following the analysis of Sec.  III.B.4 in Ref. Bazavov et al. (2016), we estimate isospin-breaking effects to enter at $\mathrm{O}((m_{d}^{\text{sea}}-m_{u}^{\text{sea}})^{2})$, leading to a negligible $\sim 0.01\%$ uncertainty. The MILC asqtad ensembles employed do not include electromagnetism. Contributions from dynamical photons would enter at the one-loop level, and we estimate their size to be of $\text{O}{(\alpha_{\text{EM}}/\pi)}\sim 0.2\%$, again following Sec. VIII.B.4 of Ref. Bazavov et al. (2016). We add this error from the omission of electromagnetism in quadrature to the fit error and list it in the “EM” column of Table 12.

VIII.3 Omission of the charmed sea quark

The MILC asqtad ensembles do not include dynamical charm quarks in the sea. The effects of ignoring the charmed sea quark are discussed in detail in Sec. VIII.C of Ref. Bazavov et al. (2016), and are estimated from power counting to be $\sim$1–2% for reasonable choices of $\alpha_{s}$ and $\Lambda_{\rm QCD}$. To be conservative, we take the upper end of the range, or 2%, for the “Charm sea” error in Table 12. We note, however, that for the decay constants $f_{\pi}$, $f_{K}$, and $f_{D_{(s)}}$ where both 3- and 4-flavor simulation results are available, the observed differences are consistent with zero within errors.

VIII.4 Other consistency checks and error summary

Finally, we perform fits over various subsets of our data to check for overall consistency and further verify that our base-fit results are reasonable. These are included in Figs. 14 and 15. First we perform fits omitting data from the largest or smallest lattice spacing, which are labeled “$\text{no }a\approx 0.12\text{~{}fm}$” and “$\text{no~{}}a\approx 0.045~{}\text{fm}$”, respectively. The resulting matrix-element central values agree with the base fit within 1$\sigma_{\rm fit}$, providing further evidence that our wide range of lattice spacings is sufficient to control discretization errors. The resulting matrix-element uncertainties are larger, however, which is to be expected because a smaller data set is employed.

Our base-fit results are obtained from a single chiral-continuum fit to the matrix elements of all five operators, including correlations, to optimally constrain the shared LECs and parametric inputs. To test the impact of the correlations we perform five separate fits of the individual matrix elements; the results are labeled as “individual”. We observe large shifts in the matrix-element central values, as much as 2$\sigma_{\rm fit}$ in some cases, which are covered by equally substantial increases in the uncertainties for all matrix elements except $\langle{\mathcal{O}}_{3}\rangle$. In this case, the individual-fit error does not quite overlap that of the combined fit. The large errors obtained from the individual fits are associated with the uncertainties on the NLO LECs $\beta^{(\prime)}_{j}$. The majority LO, NLO and even NNLO LECs are well constrained by the data, while the parametric inputs are tightly constrained by priors. The LECs $\beta^{(\prime)}_{j}$, however, cannot be well determined using data for only a single matrix element, and the errors obtained in the individual fits are governed by the loose prior widths. The base fit resolves the $\beta^{(\prime)}_{j}$ because these coefficients multiply terms that mix under operator renormalization.

After considering all possible significant sources of uncertainty in the previous sections, Table 12 presents complete systematic error budgets for all five $D$-meson matrix elements. The error due to the omission of the charmed sea quark is listed separately after the total because the estimation of this error is far more rough and less quantitative than all others considered. Further, if more reliable estimates of charmed-sea-quark effects become available in the future, this separation will enable the errors on our results to be easily adjusted a posteriori.

IX Results

Here we present our final results with complete error budgets including all sources of uncertainty considered in the previous section. We first give results for the local neutral $D$-meson mixing matrix elements in Sec. IX.1. In Sec. IX.2, we discuss how to obtain bounds on generic sources of new physics and also illustrate how to apply our results to a specific model. In addition, Appendix C provides tables of correlations between our matrix-element results, and between our bag-parameter results, so that they can be employed in future phenomenological studies.

IX.1 Matrix elements

Table 13 presents our final results for the relativistically normalized $D$-meson mixing hadronic matrix elements of operators $\mathcal{O}_{i}$ ($i=1$–5) in Eqs. (14)–(18) including statistical and all sources of systematic uncertainty. We give the matrix elements in the $\overline{\text{MS}}$-NDR scheme at the scale $\mu=3~{}\text{GeV}$ obtained with the BBGLN Beneke et al. (1999) and BMU Buras et al. (2000) choices of evanescent operators. Although the choice of evanescent operators affects only the renormalization of operators $\mathcal{O}_{2}$ and $\mathcal{O}_{3}$, the numerical results for the other three matrix elements $\left<\mathcal{O}_{1,4,5}\right>$ differ slightly for the BBGLN and BMU schemes, because they are all obtained simultaneously with correlations in the joint chiral-continuum fit. The correlation matrix between the five matrix elements is given in Appendix C, Table 16. We quote the uncertainty due to the omission of charm sea quarks separately, because the estimate is semi-quantitative. Final results for the five matrix elements are also provided as supplemental material Dmi (2018) with double-precision.

Figure 16 compares our $D$-mixing matrix-elements with the lattice-QCD results obtained by the ETM Collaboration using $N_{f}=2$ Carrasco et al. (2014a) and $N_{f}=2+1+1$ Carrasco et al. (2015a) twisted-mass fermions. ETM presents values for bag parameters, which we then convert to matrix elements using Eq. (3) of Ref. Carrasco et al. (2014a), with their own $N_{f}=2$ and $N_{f}=2+1+1$ calculations of the quark masses Blossier et al. (2010); Carrasco et al. (2014c) and decay constants Carrasco et al. (2014b, 2015b). In each panel, the red error bars are obtained by propagating the decay-constant and quark-mass uncertainties in quadrature, while the blue error bars omit those uncertainties. This second approach may provide a rough picture of the errors on the matrix elements, given typical correlations between the three- and two-point correlation functions in the numerator and denominator. Our results for $\langle\mathcal{O}_{1}\rangle$, $\langle\mathcal{O}_{2}\rangle$, $\langle\mathcal{O}_{3}\rangle$ agree with the matrix elements converted from ETM’s bag parameters to within about 1–2$\sigma$, but those for $\langle\mathcal{O}_{4}\rangle$ and $\langle\mathcal{O}_{5}\rangle$ differ by 1.7–3.3$\sigma$ (assuming ETM’s quark-mass and decay-constant errors to be negligible). We find, however, that different choices for the quark masses and decay constants yield considerable variations in the converted matrix elements. If we convert the ETM bag parameters using the average quark masses and decay constants from the PDG Rosner et al. (2016); Patrignani et al. (2016), we obtain matrix element results that agree with ours within 0.25–1.52$\sigma$ for all five operators. Use of the averages from the Flavor Lattice Averaging Group Aoki et al. (2017) instead yields matrix-element results that lie in between these two determinations. We are planning to perform a correlated, combined analysis of the $D$-mixing matrix elements from this work with our collaboration’s $D$-meson decay constants calculated using the same lattice ensembles and parameters Bazavov et al. (2014); Neil et al. (2015). We will present the resulting the $D$-meson bag parameters in our forthcoming decay-constant paper. This will enable a more direct comparison with ETM’s results.

Similar tensions to those shown in Fig. 16 have been observed for the analogous neutral-kaon-mixing bag parameters $B_{4}$ and $B_{5}$ Bertone et al. (2013); Carrasco et al. (2015a); Choi et al. (2016); Garron et al. (2016), which the authors of Ref. Garron et al. (2016) attribute to the choice of intermediate renormalization scheme. Results obtained in Ref. Garron et al. (2016) using the symmetric regularization-independent momentum-subtraction (RI-SMOM) scheme Sturm et al. (2009) are in good agreement with an independent calculation that employs one-loop mean-field improved lattice perturbation theory Choi et al. (2016), but differ substantially with results obtained with the RI-MOM Martinelli et al. (1995) renormalization scheme. This discrepancy is attributed to underestimated systematic errors present in the RI-MOM scheme Garron et al. (2016); Sturm et al. (2009). We note that operator $\mathcal{O}_{4}$ has the largest negative anomalous dimension of the five $\Delta C=2$ operators, by about 50%, and is therefore most sensitive to running of the renormalization scale. Because $\mathcal{O}_{4}$ and $\mathcal{O}_{5}$ mix under renormalization, this also affects the results for $\langle\mathcal{O}_{5}\rangle$. Thus, $\mathcal{O}_{4}$ and $\mathcal{O}_{5}$ are likely to be the most sensitive to difficulties in renormalization.

IX.2 Implications for new physics

As discussed in Sec. II, neutral $D$-meson mixing is a sensitive probe of local $\Delta C=2$ interactions from physics beyond the Standard Model that contribute to the quantity $M_{12}$. The phase $\phi_{12}$ is particularly sensitive to new physics due to the CKM suppression of its contribution from the Standard Model. Assuming that any new physics does not change the phase of $\Gamma_{12}$, which is true in many BSM models Golowich et al. (2009), then the imaginary part of $M_{12}$ gives the most sensitive constraint on new physics. Here we use our matrix-element results to bound the scale of new physics in a generic model that alters the high-scale Wilson coefficients and in a specific flavor-violating Higgs model with tree-level flavor-changing neutral currents.

A general new-physics model will give nonzero values for some of the Wilson coefficients $C_{i}$ at a high scale $\Lambda_{\text{NP}}$. To evaluate the contribution to $M_{12}$ using Eq. (22), we must run the Wilson coefficients down to the scale $\mu=3$ GeV at which our matrix elements are renormalized. We use the one-loop running derived in Ref. Golowich et al. (2009), which works in a different four-fermion operator basis than the one used here. Fierz identities relate the two bases Bouchard (2011):

[TABLE]

The remaining operators $Q_{6}$, $Q_{7}$ and $Q_{8}$ are parity conjugates of $Q_{1}$, $Q_{4}$ and $Q_{5}$. Applying this change of basis to the formulas given in Appendix A of Ref. Golowich et al. (2009) yields the renormalization equations in our basis, which we then use directly. The operator running depends on the value of the strong coupling $\alpha_{s}$ at intermediate scales. We use the RunDec Mathematica library Chetyrkin et al. (2000) to compute $\alpha_{s}$ with four-loop running, including quark decoupling effects.

To illustrate the use of our results in constraining physics beyond the Standard Model, we first consider a simple new-physics model which gives rise only to operator $\mathcal{O}_{5}$ at an ultraviolet scale $\Lambda_{\text{NP}}$, with purely imaginary (and hence CP-violating) Wilson coefficient

[TABLE]

Running down to 3 GeV, we find

[TABLE]

where

[TABLE]

and all other Wilson coefficients at 3 GeV remain zero. Thus, operators $\mathcal{O}_{4}$ and $\mathcal{O}_{5}$ contribute a purely imaginary term to

[TABLE]

from new physics. (Here, $x^{\text{NP}}_{12}$ is defined as in Sec. II, $x^{\text{NP}}_{12}=M^{\text{NP}}_{12}/\Gamma_{D}$.) We can now determine the constraint on $\Lambda_{\text{NP}}$ in this model by summing the Standard Model and new physics contributions, and then comparing to experimental bounds, as depicted in Fig. 17. The gray box shows the region in which the Standard Model value for $x_{12}$ is expected to lie. The gold box includes a contribution from our simple scenario with $\Lambda_{\text{NP}}=40\,000$ TeV; this level of new contribution is consistent with the experimental bounds at 1$\sigma$. Finally, the red box shows a new-physics contribution with $\Lambda_{\text{NP}}=18\,000$ TeV, which is ruled out by experiment with very high confidence. 333Note that if the Wilson coefficient $C_{5}^{\text{NP}}(\Lambda_{\text{NP}})$ were purely real instead of purely imaginary, it would have to be much larger in magnitude to give a SM + NP region that lies outside the experimental contour, leading to a much weaker bound (on the order of 2000 TeV in this simple scenario).

We now apply this technique to place bounds on the scale of generic new-physics contributions to each operator. In this case, the Wilson coefficients are of the form

[TABLE]

where $F_{i}$ and $L_{i}$ are flavor and loop-counting factors Bona et al. (2008). In order to obtain bounds on $\Lambda_{\text{NP}}$ individually by operator, we set one $C_{i}^{\text{NP}}(\Lambda_{i,\text{NP}})=F_{i}L_{i}/\Lambda_{i,\text{NP}}^{2}$ at a time, with all other Wilson coefficients set to zero at the high scale. We then determine the value of $\Lambda_{i,\text{NP}}$ for which the new-physics prediction for $|x_{12}|e^{i\phi_{12}}$ is inconsistent with the experimental bound at 95% confidence, including the uncertainty in the matrix elements for $\left<\mathcal{O}_{i}\right>$. Here we assume that $\mathop{\mathrm{Im}}F_{i}L_{i}>0$, but taking the other sign gives a nearly identical constraint. (The experimental bound is $\mathop{\mathrm{Im}}x_{12}<0.0289\%$ ($\mathop{\mathrm{Im}}x_{12}>-0.0285$) in the positive (negative) imaginary direction.)

Using our matrix elements in the BBGLN scheme, we obtain

[TABLE]

Note that these bounds are from the imaginary part of $x_{12}$ only; for new physics with $\mathop{\mathrm{Im}}F_{i}L_{i}\approx 0$, the constraint on $\Lambda_{\text{NP}}$ from the $\mathop{\mathrm{Re}}x_{12}$ will dominate.

Our bounds in Eqs. (92)–(96) are stronger than those quoted by the ETM Collaboration Carrasco et al. (2014a), in part because we use more recent, tighter experimental bounds. For operators $\mathcal{O}_{3}$ and $\mathcal{O}_{5}$ in particular, our constraints on $\Lambda_{\text{NP}}$ are much higher (by factors of roughly 7 and 3, respectively). These two operators mix strongly with $\mathcal{O}_{2}$ and $\mathcal{O}_{4}$, respectively, such that their bounds stem principally from $C_{2}(3~{}\text{GeV})\mathcal{O}_{2}(\text{3~{}GeV})$ and $C_{4}(3~{}\text{GeV})\mathcal{O}_{4}(\text{3~{}GeV})$. If we artificially set $C_{2}(3~{}\text{GeV})=0$, then we obtain much weaker bounds of $\Lambda_{3,\text{NP}}\gtrsim 3\,330$ TeV and $\Lambda_{5,\text{NP}}\gtrsim 9\,700$ TeV, respectively, which are close to the values quoted by ETM in Ref. Carrasco et al. (2014a).

To further illustrate the use of our results in constraining new physics, we examine a specific model in which the Standard Model Higgs boson has flavor-violating couplings to quarks and leptons Harnik et al. (2013). A Higgs coupling of the form $Y_{uc}h\bar{u}_{L}c_{R}+Y_{cu}h\bar{c}_{L}u_{R}$ will induce $\Delta C=2$ four-fermion interactions at low energy. After integrating out the Higgs boson $h$, the effective Hamiltonian is

[TABLE]

where the Wilson coefficients at the scale $m_{h}$ are given by

[TABLE]

We write the two Yukawa couplings as $Y_{uc}=|Y_{uc}|e^{i\phi_{uc}}$, and similarly for $Y_{cu}$. Taking the Wilson coefficients to be purely imaginary and comparing the resulting $\mathop{\mathrm{Im}}x_{12}$ to the experimental bounds, we find from the $\mathcal{O}_{2}$ and $\mathcal{O}_{4}$ terms separately the bounds

[TABLE]

To be somewhat more general, we can instead marginalize over the phases in the couplings (integrating from 0 to $\pi/2$) and obtain exclusion contours in the $|Y_{uc}|$-$|Y_{cu}|$ plane, again solely from the stronger experimental bound on $\mathop{\mathrm{Im}}x_{12}$. The contours are shown in Fig. 18. Our bounds on the combinations of Yukawa couplings in Eqs. (101)–(102) are tighter than those in Ref. Harnik et al. (2013) by over an order of magnitude. The improvement stems primarily from the use of newer experimental measurements than employed by the authors of Ref. Harnik et al. (2013), who relied on the same data as the original model-independent analysis of $D$-meson mixing Ciuchini et al. (2007). We also explicitly run the Wilson coefficients from $m_{h}$ down to 3 GeV in order to compare to experiment, rather than attempting to obtain “model-independent” bounds on the Wilson coefficients at a generic high scale.

X Summary and outlook

In this paper we have presented a three-flavor lattice-QCD calculation of the neutral $D$-meson mixing matrix elements of all five $\Delta C=2$ dimension six local four-fermion operators. We obtain uncertainties comparable to those from earlier $N_{f}=2$ and $2+1+1$ calculations by the European Twisted Mass Collaboration Carrasco et al. (2014a, 2015a). Our results for $\langle D^{0}|\mathcal{O}_{i}|\bar{D}^{0}\rangle$ ($i=1$–3) agree with those of ETM to within about 1–2 standard deviations, but those for $i=4$ and 5 differ more significantly. These short-distance matrix elements are needed to evaluate the first term on the right hand side of Eq. (10), and can be combined with experimental measurements of the $D$-mixing parameters to yield useful constraints on theories beyond the Standard Model with sizable CP violation because $\phi_{12}$ is very small in the Standard Model. To illustrate the utility of our matrix-element results, we place bounds on generic new high-scale physics that would give rise to local $\Delta C=2$ interactions, finding $\Lambda_{\text{NP}}{(\mathop{\mathrm{Im}}F_{i}L_{i})^{-1/2}}\gtrsim 10$–$50\times 10^{3}$ TeV for the five local $\Delta C=2$ operators. These results are more stringent than previous bounds, in part because we use the latest experimental measurements, and in part because of the way we introduce new physics at the high scale, $\Lambda_{\text{NP}}$, and run down to 3 GeV.

The long-distance contributions described by the second term on the right hand side of Eq. (10) are expected to dominate the Standard Model prediction of the $D^{{0}}$-meson mixing observables, $\Delta M$ and $\Delta\Gamma$, but their size is not well known because they are difficult to calculate. Thus, despite the current precision of experimental measurements—and forthcoming improvements from the BES III, LHCb, and Belle 2 experiments—our results and those of Refs. Carrasco et al. (2014a, 2015a) will suffice for phenomenology until better methods for the long-distance contributions become available.

Clearly, reliable methods to compute the long-distance parts of $M_{12}$ and $\Gamma_{12}$ are needed. These must account for contributions from intermediate multi-hadron states that can propagate over hadronic distances. Fortunately, the theoretical framework for obtaining transition amplitudes, scattering lengths, and phase shifts for two-hadron systems in a finite spatial volume is already well developed Lüscher (1986, 1991); Rummukainen and Gottlieb (1995); Lellouch and Lüscher (2001); Kim et al. (2005); Christ et al. (2005), and the number of lattice-QCD calculations of these systems is growing rapidly. The first lattice-QCD calculations for systems with many open channels have been recently performed by the Hadron Spectrum Collaboration Dudek et al. (2014); Wilson et al. (2015a); Dudek et al. (2016); Wilson et al. (2015b); Moir et al. (2016). These techniques have been extended to neutral kaon mixing, and the first lattice-QCD calculations have recently become available Bai et al. (2014); Christ et al. (2015). There, long-distance effects are easier to quantify, because phase space suppresses all but the two-pion intermediate states. Efforts are underway by several groups to develop an analogous theoretical framework for three-hadron systems He et al. (2005); Hansen and Sharpe (2012); Briceño and Davoudi (2013a); Polejaeva and Rusetsky (2012); Briceño and Davoudi (2013b); Briceño (2014); Briceño et al. (2015); Hansen and Sharpe (2014); Briceño and Hansen (2015); Hansen and Sharpe (2015, 2016, 2017); Briceño et al. (2017), which are relevant for $D^{0}$ mixing. More recently, Hansen, Meyer, and Robaina proposed a new method to extract the spectral function for multi-hadron transition rates from finite-volume four-point correlation functions Hansen et al. (2017). In the coming years, these new approaches may be implemented in numerical lattice-QCD calculations and lead to more quantitative estimates of long-distance contributions to neutral $D$-meson mixing.

Acknowledgements.

We thank Alan Schwartz, Vittorio Lubicz, Silvano Simula, Roni Harnik, and Joachim Kopp for useful correspondence. Computations for this work were carried out with resources provided by the USQCD Collaboration, the National Energy Research Scientific Computing Center and the Argonne Leadership Computing Facility, which are funded by the Office of Science of the U.S. Department of Energy; and with resources provided by the National Institute for Computational Science and the Texas Advanced Computing Center, which are funded through the National Science Foundation’s Teragrid/XSEDE Program. This work was supported in part by the U.S. Department of Energy under Grants No. DE-FG02-91ER40628 (C.B.) No. DE-FC02-06ER41446 (C.D.) No. DE-SC0010120 (S.G.), No. DE-SC0010005 (E.T.N.), No. DE-FG02-91ER40661 (S.G., R.Z.), No. DE-FG02-13ER42001 (C.C.C., D.D., A.X.K.), No. DE-SC0015655 (A.X.K.), No. DE-FG02-13ER41976 (D.T.); by the U.S. National Science Foundation under Grants No. PHY10-67881 and No. PHY14-14614 (C.D.), No. PHY14-17805 (D.D., J.L.), and No. PHY13-16748 and No. PHY16-20625 (R.S.); by the Fermilab Fellowship in Theoretical Physics (C.M.B., C.C.C.); by the URA Visiting Scholars’ program (C.M.B., C.C.C., D.D., A.X.K.); by the MINECO (Spain) under Grant No. FPA2013-47836-C-1-P (E.G.); by the Junta de Andalucía (Spain) under Grants No. FQM-101 and No. FQM-6552 (E.G.); by the European Commission (EC) under Grant No. PCIG10-GA-2011-303781 (E.G.); by the German Excellence Initiative and the European Union Seventh Framework Program under Grant Agreement No. 291763 as well as the European Union’s Marie Curie COFUND program (A.S.K.). Brookhaven National Laboratory is supported by the Department of Energy under Contract No. DE-SC0012704. Fermilab is operated by Fermi Research Alliance, LLC, under Contract No. DE-AC02-07CH11359 with the United States Department of Energy, Office of Science, Office of High Energy Physics. This document was prepared by the Fermilab Lattice and MILC Collaborations using the resources of the Fermi National Accelerator Laboratory (Fermilab), a U.S. Department of Energy, Office of Science, HEP User Facility.

Appendix A Correlations among $r_{1}/a$ data

Table 14 provides the correlations among values of the relative scale $r_{1}/a$ on the ensembles, which are needed to convert quantities from lattice units to $r_{1}$ units prior to performing the chiral-continuum extrapolation. Correlations between $r_{1}/a$ arise after performing a smoothing fit simultaneously to all ensembles.

Appendix B Priors for two-point and three-point fits

Table 15 provides the priors employed in the joint two- and three-point correlator fits discussed in Sec. IV. The parameters $E_{0}$ and $Z_{n}$ are defined in Eq. (29), $\Delta_{k,j}$ is defined in Eq. (35), while the parameters $\mathcal{Z}^{O_{i}}_{nm}$, $i=1$–5, are defined in Eqs. (30)–(31). For the energy splitting the prior is defined such that $\exp(\Delta_{k,j})$ is Gaussian distributed.

Appendix C Correlations among matrix-element results

Table 16 provides the correlations among the $D$-meson mixing matrix elements for all five operators, to enable their use in future phenomenological studies.