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

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.

We calculate in three-flavor lattice QCD the short-distance hadronic matrix elements of all five ∆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π ≈ 180 MeV and lattice spacings as fine as a ≈ 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 MS-NDR scheme using the choice of evanescent operators proposed by Beneke et al., evaluated at 3 GeV, D 0 |Oi|D 0 = {0.0805 (55) (16), −0.1561 (70) (31), 0.0464(31)(9), 0.2747(129) (55), 0.1035(71)(21)} 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×10 3 TeV for couplings of 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.

II. THEORETICAL AND PHENOMENOLOGICAL BACKGROUND
The time evolution of a neutral-meson system, such as D 0 andD 0 , can be described by a Schrödinger equation where the mass matrix M and decay matrix Γ are Hermitian. The off-diagonal term M 12 − i 2 Γ 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 1 |D 0 1 is identified as the one with D 0 1 |CP|D 0 1 ≈ +1. Then, by convention, the mass and width differences between the two eigenstates are defined as [1,2] (2. 3) The signs of ∆M and ∆Γ are determined from experiment. The eigenvalue problem leads to the relation for the mass and width differences, where the sign is the (near) CP of the heavier state, and Given measurements of ∆M , ∆Γ, and CP asymmetries sensitive to φ 12 , these formulas determine M 12 and Γ 12 apart from a mutual, unphysical phase [14][15][16].
The current measurements can be summarized as [2] y ≡ ∆Γ 2Γ = 0.61 ± 0.07%, (2.7) x ≡ ∆M Γ = 0.41 +0.14 −0.15 %, (2.8) and asymmetries consistent with zero. A fit then yields [2] φ 12 = −0.17 • ± 1.8 • , (2.9) and values for y 12 ≡ |Γ 12 |/Γ and x 12 ≡ 2|M 12 |/Γ the same as those for y and x. It is expected that future measurements by LHCb and Belle II will reduce the uncertainties on ∆M to 10% or less [9][10][11]. 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, Λ QCD , satisfy m W m b Λ QCD , so the mixing matrix in Eq. (2.1) can be expressed as (see, e.g., Ref. [17]) where H ∆C=1 (H ∆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 Γ 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 −1 b or less. In the second contribution, however, intermediate states, for example K + π − , can propagate a distance of order Λ −1 QCD . 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 φ 12 to be very small: the imaginary parts of M 12 and Γ 12 , and hence their phases, come from parts of the box diagrams carrying one or two factors of sin θ 23 sin θ 13 sin δ = 1.4 × 10 −4 , while the corresponding CKM factor in the real parts is sin θ 12 = 0.225. Because, with these conventions, both phases are small, so is the convention-independent difference φ 12 . For the same reason, the Standard Model real part of M 12 stems mostly from the long-distance contribution, the sum in Eq. (2.10). An estimate from a dispersion relation based on heavy quark effective theory yields [18] y ∼ 10 −2 , (2.11) x ∼ 10 −3 to 10 −2 , (2.12) which are compatible with the measurements, Eqs. (2.7) and (2.8).
Physics beyond the Standard Model could change M 12 , Γ 12 , or both. Many [7] extensions of the Standard Model alter only H ∆C=2 and, hence, the magnitude and phase of M 12 , but not Γ 12 . In general, the ∆C = 2 effective Hamiltonian can be written iÕi , (2.13) where C i are the Wilson coefficients and the O i are four-quark operators, given below. The determinations of x 12 and φ 12 therefore can constrain the Wilson coefficients once the hadronic matrix elements D 0 |O i |D 0 have been computed in (nonperturbative) QCD. Unfortunately, in view of the large range in Eq. (2.12), the constraint from x will remain loose until new techniques for the long-distance term have been developed. The operators in the ∆C = 2 effective Hamiltonian are O 3 =c α Ru βcβ Ru α , (2.21) wherec and u are the anticharm-and up-quark fields, with left and right Dirac projection matrices L = 1 2 (1 − γ 5 ) and R = 1 2 (1 + γ 5 ). The color indices are denoted α and β, and the spin indices are implied. All other Lorentz invariant four-quark operators can be reduced to this set by Fierz rearrangement [19]. Further, parity conservation in QCD implies D 0 |O i |D 0 = D 0 |Õ i |D 0 , i = 1, 2, 3.
In summary, the five matrix elements D 0 |O i |D 0 suffice to obtain (2.22) where the C NP i (µ) are the Wilson coefficients of the new-physics model, subsumingC NP i , 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) -B 0 (s) mixing [6]. Given these matrix elements, Eq. (2.22) can be used to test models in which new physics does not change the phase of Γ 12 . 2 A convenient way to do so is illustrated in Fig. 2, which plots |x 12 |e iφ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. (2.12). The gray bar is close to sin φ 12 = 0 and extends well beyond the range displayed here. With the proviso that new physics does not change the phase of Γ 12 , the new-physics calculation from Eq. (2.22) yields a complex number x NP 12 ≈ |x NP 12 |e iφ NP 12 , which should be added to the gray bar. If the total |x 12 |e iφ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 A 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 B. We conclude in Sec. III C with a discussion of statistical errors and autocorrelations.
A. Gauge-field configurations, light-, and heavy-quark actions We use isospin-symmetric gauge-field configurations generated by the MILC Collaboration [20][21][22] 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) × 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 O(α s a 2 , a 4 ) [23][24][25][26][27][28]. 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 [29][30][31][32]. The gluons are simulated with the tadpole-improved Lüscher-Weisz action and have discretization errors starting at O(α s a 2 , a 4 ) [33][34][35][36]. The lightest simulated pion mass M π = 177 MeV is very close to the physical value; heavier pion masses up to M π = 555 MeV are also included in the analysis to help guide the chiral extrapolation. Four lattice spacings ranging from a ≈ (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 π L 3.8) that finite-volume effects are expected to be at the subpercent level [37]. 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 I.
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 [22] listed in Table I. Note that the ratios depend only on β, 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 2 1 F (r 1 ) = 1.0 [38,39]. 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 [22], and can be obtained with  [22]. From left to right are shown the precise input value of β which is related to the gauge coupling β = 10/g 2 , the approximate lattice spacing a, the mass-independent ratio of the scale and lattice spacing r1/a, the simulated light-to-strange sea-quark mass ratio am l /am s , the lattice volume (L/a) 3 × (T /a), the taste-Goldstone pion mass Mπ and RMS sea-pion mass M RMS π , the spatial extent in units of the pion mass MπL, and the number of configurations N conf in each ensemble. The primes on the light-quark masses distinguish the simulation values from the physical (unprimed) values m l = (mu + m d )/2 and ms.
. tiny statistical errors. We convert our final matrix-element results to physical units using the absolute scale [40], based on calculations of light pseudoscalar-meson decay constants from MILC [41] and HPQCD [42]. A detailed discussion of the smoothing procedure may be found in Sec. IV B of Ref. [40].
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 II, and correspond to pion masses from M π ≈ 180-880 MeV, enabling good control over the chiral extrapolation guided by SU(3) chiral perturbation theory. In Table II 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 [43,44]. The couplings in the theory are smoothly bounded for arbitrary values of the heavy-quark mass am Q . After treelevel matching to QCD through heavy-quark effective theory (HQET), discretization errors from the action are of O(α s aΛ QCD , a 2 Λ 2 QCD ). The bare charm-quark mass in the Fermilab action is parametrized by the hopping parameter κ c ; the values employed in our simulations are tabulated in Table II.

B. Lattice correlation functions
We calculate the two-and three-point correlation functions needed to obtain the matrix elements for neutral Dmeson mixing using the same method as for B-meson mixing in Ref. [45]. In particular, we first construct a specific combination of a light-quark propagator, heavy-quark propagator, and γ 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 γ 5 , and obtain the three-point correlation functions for the five ∆C = 2 four-quark operators from combining two open-meson propagators contracted with the appropriate Dirac structures. Here we describe the lightand 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 Υ(x) from the staggered field χ(x) following Refs. [46,47], where Ω(x) = γ x1 1 γ x2 2 γ x3 3 γ x4 4 is the Kawamoto-Smit transformation and χ denotes a vector of the four staggered copies. We remove tree-level, O(a) discretization errors from the four-fermion operator by rotating the heavy-quark field ψ(x) TABLE II. Parameters of the valence-quark propagators used in this work. Each ensemble is labeled by the approximate lattice spacing a and the ratio of simulated light to strange sea-quark masses, am l /am s . The same valence light-quark masses amq are used on all ensembles with the same approximate lattice spacing except on the a ≈ 0.09 fm, am l /am s = 0.00155/0.031 ensemble. The next three columns list the parameters of the simulated charm-quark propagators: the clover coefficient cSW [43], the hopping parameter κ c , and the rotation coefficient d 1c . The primes on the hopping parameter and rotation coefficient distinguish the simulation values from the physical (unprimed) values. The last column shows the number of time sources per configuration Nsrc.  3. Left: MILC asqtad ensembles used in our analysis. The colors label the approximate lattice spacings a ≈ 0.12 fm (yellow), 0.09 fm (green), 0.06 fm (blue), and 0.045 fm (purple), while the symbol size is proportional to the number of configurations. The cyan star shows the physical point (continuum limit and physical pion mass.) Right: Valence pion masses used in our analysis. The diagonal line denotes full-QCD points with M val Lattice three-point correlation function CO i (tx, ty). The thick and thin lines denote the charm-and light-quark propagators, respectively. The D-meson is created at time t0 + tx < t0, while theD meson is later annihilated at time t0 + ty > t0. The four-quark operator Oi is fixed at time t0.
following Ref. [45], The simulation values of the rotation parameter d 1 are given in Table II. We form the D-meson interpolating operator from the light naive field and rotated heavy field as where S(x, x ) 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 [48,49]. We obtain the D-meson mixing matrix elements from the zero-momentum three-point correlation functions, where the index i =1-5 labels the four-quark operator. We obtain the lattice operators, O i , from the expressions (2.14)-(2.18) for the continuum operators, O i , by replacing the u field with Υ and thec field withΨ. As shown in Fig. 4, the four-quark operator location is fixed at time t 0 , while the D-andD-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. [6]. This mixing is a lattice-discretization effect of 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,

C. 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 √ N 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 ≤ t ≤ 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Õ 1,2,3 in Eqs. (2.14)-(2.18) differ from O 1,2,3 by an interchange of L ↔ 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 O 1,2,3 and Õ 1,2,3 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 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 /m s = 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. [6], 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. (3.5) and (3.6). All dimensionful quantities are given in lattice spacing units, where the factors of a are suppressed for simplicity.

A. Correlator fit functions
Our choice of using the naive field Υ(x), defined in Eq. (3.2), for the valence up quark in the two-and three-point correlation functions results in contributions from D-meson states with positive parity [46]. 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 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. (4.1) arises with periodic boundary conditions in time.
The three-point functions are described by where the desired matrix element is related to the ground-state amplitude Z Oi 00 : The factor of the D-meson mass M D is due to the nonrelativistic normalization of states in Eq. (4.2). Effects from periodic boundary conditions are negligible in our three-point data and are therefore not included in Eq. (4.2).

B. Method
As shown in Eq. (4.3), 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 χ 2 function defined in Eq. (B3) of Ref. [6], while the Q parameter defined in Eq. (B4) of Ref. [6] is used to assess the quality of the fits. The selection of the Bayesian priors is described in Sec/ IV C. In the fits to the correlation functions we use intervals t min ≤ t ≤ t 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 D. We limit the number of states used in Eqs. (4.2) and (4.1) to N 3pt = N 2pt ≤ 6. These choices are designed to optimize the extraction of the desired matrix elements from the correlation function data. Section IV E describes fit variations with different t min , t max , and N 3pt to ensure our fit choices do not bias the fit results for the ground state parameters.

C. 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 E 2, our main analysis employs two-and threepoint 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, 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 valenceand 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, where M eff is an estimate of M 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 ∆ 1,0 from the ground state with and 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 ∆ 1,0 is guided by the experimentally measured D * 0 (2400) 0 − D 0 mass difference [1]; we use ∆ 1,0 ≈ 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 E 2). The prior width of Z 1 is set to be about one-σ 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, The choice of prior for the splittings of the radial excitations is motivated from quark model estimates [50], ∆ k,k−2 ≈ 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 Z Oi 00 we examine the scaled three-point function at large times t x , t y , 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 Z Oi mn (defined in Eq. (4.2)) 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 XV in Appendix B.

D. 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 ∈ {t min , . . . , t 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 max for  the three-point function but analyze the data in the entire region bounded by |t x |, t y ∈ {t min , . . . , t 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 ∈ {t min , . . . , t 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 III lists our choices for t min and t max for our preferred fits. Across all four lattice spacings, we set t min = 0.72 fm while varying t max smoothly from 3 fm on the a ≈ 0.12 fm ensembles, 2.7 fm on the a ≈ 0.09 fm ensembles, 2.5 fm on the a ≈ 0.06 fm ensembles, to 2.3 fm for the a ≈ 0.045 fm ensemble. We choose one set of priors and time ranges for all correlation function fits at a given lattice spacing.

E. Fit stability
Our preferred correlator fits are performed with the priors listed in Table XV, the time ranges listed in Table III, the time region as described in Fig. 8, and with N 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 min , t max , and N 2pt = N 3pt ) for the two-and three-point functions. Figure 9 (left) shows a typical example of a t min stability plot for Z O1 00 with the preferred fit displayed with a solid blue point. The fits that include only 1+1 states show a strong t min dependence before reaching a plateau at t min ≥ 9, while fits with 2+2 or 3+3 states reach a plateau for t min ≥ 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 ≤ t min ≤ 9 are necessary and sufficient to account for excited state effects.

Fit range and number of states
A typical example of a t max stability plot is shown Fig. 9 (Right) for Z Oi 00 where the preferred fit is marked with a solid blue point. The fit results (central values and error bars) do not change as t 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 max is decreased to t max 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 max in Table III. 2. Operator smearing Figure 10 (left) shows an example comparing fit results for the two-point function ground state overlap factor Z 1/2 0 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 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 min , and have errors that stabilize for t min 5. By comparison, the ground state energies from the local two-point function show strong t min dependence and differ at small t min from the ground state energies with 1S-smearing before they become consistent with them. The larger range of stability in t 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.

F. Error propagation
We propagate the distribution of the matrix element Z Oi 00 by bootstrap resampling. The distributions of the priors are included by randomizing the prior central value over the prior width under the bootstrap sampling [51]. 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 Z O1 00 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Ō i (µ) in the continuum MS-scheme evaluated at the scale µ = 3 GeV. In noninteger dimensions, the Dirac algebra is infinite dimensional and is fully defined only after choosing a basis of evanescent operators [52]. This choice, which here only affects the renormalization of operators O 2 and O 3 , is not unique; here we consider the schemes of Beneke et al. (BBGLN) [53,54] and of Buras, Misiak, and Urban (BMU) [55], reporting results for both.
The ∆C = 2 four-fermion operators mix under renormalization. At O(α s ), the renormalized operators are given by [56,57] where the P p are dimension-seven operators that are not needed in this paper, because with the O i described after Eqs.(3.5), the coefficient matrix b ip starts at order α s . Neglecting this contribution leads to a discretization error of order α s aΛ QCD , the same as that from our choice of c SW , which is the coefficient of the clover term in the Sheikholeslami-Wohlert action [43]. We calculate the renormalization coefficients Z ij (aµ) using the "mostly nonperturbative matching" (mNPR) method introduced in Refs. [58,59]. We factor out the flavor-conserving renormalization coefficients Z 4 Vqq using The Z V 4 qq are then calculated nonperturbatively from equal-mass vector current correlation functions, as discussed in Ref. [40]. The flavor-changing coefficients ρ ij are calculated to one-loop in lattice perturbation theory via where the coefficients Z and the tree-level matching factors are C ll = 2u 0 and C cc = 2κ 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 ρ 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 IV lists the renormalization coefficients for flavor-conserving vector currents. The same value for Z V 4 ll 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 ρ factors in tadpole-improved lattice perturbation theory taking α s = α V (q * ) in the "V scheme" [60] obtained from the static-quark potential [61]. We fix the scale to be q * = 2/a, which is the typical scale of the gluon loop momenta. .
We also calculate the heavy-light matching coefficients Z ij at one loop in tadpole-improved perturbation theory. The tadpole-improved coefficients ζ [1] ij are related to the coefficients Z We compute the renormalization coefficients taking the tadpole-improvement factor u   ii − ρ [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 V, respectively. This difference is the same for all operators i = 1-5.

VI. CHARM-QUARK MASS CORRECTION
The mass of the charm quark is set by the hopping parameter κ in the Fermilab action. We determine the appropriate κ c by requiring that the D s -meson kinetic mass obtained in our lattice simulations agree with the PDG value as described in Ref. [64,65]. In practice, initial low-statistics runs with several values of κ were performed and used to determine the simulation values κ c for high-statistics data-production runs. The tuned value of the hopping parameter corresponding to the physical-charm quark mass, κ c , was then determined using the high-statistics data. The physical κ c values on each ensemble are listed in Table VI. We account for the slight difference between our simulation κ c and the physical κ 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 [44] 1 where am 0 is the tadpole-improved bare-quark mass, and κ crit corresponds to the value of κ where the mass of the pseudoscalar-meson mass vanishes. The nonperturbatively determined values of κ crit are listed in Table IV. 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 these values are given in Table VI. We then determine the slope of the matrix elements with respect to 1/m 2 using data with κ c = 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 VII gives the slopes obtained for i=1-5 from an unconstrained linear fit of the renormalized matrix elements O i /M D in 1/m 2 , while here ∆(1/am 2 ) is the difference between the two simulated inverse kinetic masses. Figure 12 shows an example fit for O 1 .
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 ∆(1/(am 2 )) and µ i into the matrix-element uncertainties reported by the fit. We introduce priors for ∆(1/(am 2 )) with the central values and widths given in Table VI, where the errors are obtained by propagating the error from κ c through Eq. (6.1) and (6.2). We take the priors for µ i from the linear fits described above and illustrated in Fig. 12; the prior values employed for µ i are tabulated in Table VII.
L,ii − ρ   [65], and differences ∆(1/(am2)) between the simulated and physical inverse charm-quark kinetic masses on each ensemble. For κc, the first error is from statistics and fitting, and the second is from the uncertainty in the lattice scale r1. For ∆(1/(am2)) the error is from the uncertainty in the tuned κc.

A. 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χPT) [47,66]. HMrSχ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χ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 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.  TABLE VIII. Relationships between the coefficients of the wrong-spin terms in Eq. (7.2). The relationship between the primed LECs are analogous to those between the unprimed LECs, but with all entries in the table βj → β j in β

Chiral logarithms
We work at next-to-leading order (NLO) in HMrSχPT. The one-loop chiral logarithms describe nonanalytic dependence on the light-quark masses and lattice spacing, and are For the wrong-spin termsT andQ, we follow a different notation than in Ref. [47] in order to separate the LECs β ( ) i from the one-loop diagram functions. The relationships betweenT andQ in Eq. (7.2) and the chiral-logarithm functions in Ref. [47] are given in Eqs. (C1a)-(C2e) of Ref. [6]. The coefficients β (ξ) i and β (ξ) i are not all independent; for convenience, their relationships are given in Table VIII. Inspection of Table VIII shows that the matrix elements 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. [47].
We work at leading order in the heavy-meson expansion, but include the largest 1/M D effects from the hyperfine splitting ∆ * ≡ M D * − M D and the SU(3)-flavor splitting δ su ≡ M Ds − M D 0 . The parameter that characterizes the flavor splitting in HMrSχPT is where M ηs is the mass of a theoreticalss bound state.

Analytic terms in the chiral expansion
The analytic terms where the a 2 ∆ ξ are the taste splittings, and the leading-order LEC B 0 is related to the chiral condensate. Following Ref. [40], we construct the analytic terms in Eq. (7.3) using the dimensionless variables where M qq is the mass of the taste-pseudoscalarqq meson (q = u, l, s) and f π is the pion decay constant, and where a 2∆ = 1 16 ξ w ξ a 2 ∆ ξ is the average taste splitting, and the weight factors for ξ = {P, A, T, V, I} are w ξ = {1, 4, 6, 4, 1}, respectively. The coefficients of the analytic terms are expected to be of O(1) from chiral power counting when written in terms of x q and x∆.
The NLO, NNLO, and NNNLO analytic terms are obtained by forming all combinations of x i (i = u, d, s,∆) 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: 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. (7.2) on the scale Λ χ 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 χ 2 aug /dof. We also perform fits including NNNLO analytic terms to check for fit stability and look for truncation errors. The N 3 LO terms are:

Heavy quark discretization effects
We parametrize the leading heavy-quark discretization errors of O(α s a, a 2 , a 3 ) in our data by adding the terms F HQ disc. i to our fit function: where 13) and the β i are the same LECs as in Eq. (7.2). The fit parameters z i h combined with powers of aΛ HQ represent the HQET matrix elements. The "mismatch functions" f h (m 0 a) with h ∈ {B, 3, E, X, Y, 2} 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 O(a 2 ) and O(a 3 ) mismatch functions of the action were calculated in Ref. [67] by performing a matching calculation between lattice HQET and continuum HQET. The tree-level O(a 2 ) mismatch functions of the spinor, and consequently, the 4-quark operators, were worked out in Ref. [44]. For O(α 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 → 0 and am 0 → ∞ limits. A complete list of the functional forms of the f h (m 0 a) employed here is given in Appendix A of Ref. [40].
We also consider generic O(α s a 2 ) discretization errors when studying the stability of our fits by including the term 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. (7.11)-(7.13).

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 We propagate the uncertainties in the slopes µ i and differences ∆ 1 r1m2 to the final fit error by including them as constrained fit parameters with Gaussian priors.

Renormalization errors
As discussed in Sec. V, we calculate the renormalization coefficients at one loop; therefore, truncation errors start at O(α 2 s , α s Λ QCD /m c ). The O(α 2 s ) truncation errors are estimated in the chiral-continuum fit by adding the terms F α 2 s renorm i = α 2 s (q * ) ρ [2] ii β i + ρ [2] ij β j (7. 16) where the coefficients ρ [2] ii and ρ [2] ij are free parameters (with i = j), and the β i are the leading-order LECs for matrix elements O i defined in Eq. (7.2). The second term in Eq.(7.16) parametrizes the mixing between operators O i under renormalization. We evaluate the renormalized coupling α s at q * = 2/a. In tests of fit stability, discussed below, we consider effects of O(α 3 s ) by adding the terms ii β i + ρ [3] ij β j .
The coefficients are set to ρ [3] ii = ρ 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 α 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.

Loosely constrained fit parameters
The LECs of HMrSχPT in Eqs. (7.2) and Eqs. (7.7)-(7.9) are expected to be of 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 β i and β i are set to β 1,3,4,5 = 1(1) , β 2 = −1(1) , β 2,3,4,5 = 0(1), (7.18) 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: Recall that our base fit includes terms only through NNLO. When we include the generic discretization term F αsa 2 gen i 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 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 IX. The heavy-quark discretization terms also depend upon the scale Λ HQ , which is the cutoff of the effective theory. We use Λ HQ = 500 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 Λ HQ .

Constrained fit parameters
As discussed in Sec. III A, 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 I. 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 β [22,68] of Appendix A. We include the correlations between r 1 /a values in our fit; the correlation matrix is given in Table XIV, 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. [22].
The coefficients of the one-loop chiral logarithms depend upon the light pseudoscalar-meson decay constant f π and the D * -D-π coupling g D * Dπ . We constrain f π to the PDG value of the pion decay constant [69] f π = 130.50(1)(3)(13) MeV , (r 1 f π ) 2 = 0.04249(61) (7.23) where the uncertainties on f π are from Γ, |V ud | and higher-order radiative corrections, respectively. The error on the fit input r 1 f π includes the uncertainties from both f PDG π and r 1 added in quadrature. The coupling g D * Dπ has been studied in unquenched lattice QCD with 2 flavors, yielding g [70], and with three flavors, yielding g N f =2+1 D * Dπ = 0.55 (6) [71]. Based on these results, we constrain the coupling in our fit to be g D * Dπ = 0.53(8) , (7.24) which covers the 1σ ranges of both. The one-loop HMrSχPT chiral logarithms also depend upon the parameters a 2 δ A,V , which multiply contributions from quark-disconnected "hairpin" diagrams. Because the hairpin contributions arise from taste-symmetry breaking, the coefficients a 2 δ A,V scale with lattice spacing approximately as α 2 s a 2 . We constrain their values in our fit at a ≈ 0.12 fm to the determinations from the MILC Collaboration's chiral-continuum fit of pion and kaon masses and decay constants in Ref. [72]: (7) , r 2 1 a 2 δ A = −0.28 (6) . We incorporate the uncertainty from the charm-quark-mass correction into the fit error by constraining the slopes µ i and differences in inverse kinetic masses ∆(1/(r 1 m 2 )) that enter the correction term F κ i , Eq. (7.15), 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 VI-VII, 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

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. (7.4), in the fit because their uncertainties are negligible compared to other contributions to the error. We use the values given in Table X, which were obtained from an analysis of the staggered light pseudoscalar-meson spectrum.
We fix the scale Λ χ in the chiral logarithms to the experimental ρ-meson mass M PDG ρ = 775 MeV [1], which corresponds to a value in r 1 units of (Λ χ /r 1 ) 2 = 1.5. We have checked that our physical continuum-limit matrixelement results are insensitive to reasonable variations in Λ χ .

C. Chiral-continuum fit results
Our base fit used to obtain our matrix-element central values is to the function  and includes terms to account for errors from truncating the chiral expansion, discretization errors from tastesymmetry 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χPT and mix under renormalization. Figure 13 shows our preferred chiral-continuum extrapolation as a function of the squared valence-meson mass M 2 qq , which is approximately linear in the valence light-quark mass m q . We obtain a good fit with a correlated χ 2 aug /dof = 122.4/510, where the quantity χ 2 aug /dof, which is suitable for assessing the quality of constrained fits, is defined in Eq. (7.27) of Ref. [6].

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 A, with a discussion of errors that are included in the chiral-continuum fit. Next, in Sec. VIII B, 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 C. Throughout these sections we refer to alternate fits that we used to test error saturation. Figures in Sec. VIII D 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 XII. A. 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 bestfit 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 XI. 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 A 1-VIII A 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.

Parametric inputs
The "inputs" column of Table XI 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π . The uncertainties from the parameters in the chiral logarithms f π , ∆ * , λ 1 , δ V , and δ A 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 π in the coefficient of the chiral logarithms with the PDG value of f K ± = 155.6(4) [69], which corresponds to This leads to only a tiny shift in the matrix elements, as shown by the fit variation labeled "f K vs. f π " in Figs. 14 and 15, indicating that the fit error indeed encompasses the chiral truncation error.

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 matrixelement slopes µ i (i=1-5) and the shift in the kinetic mass ∆(1/(r 1 m 2 )) as constrained parameters with prior widths given by the uncertainties listed in Tables VI and VII. The sum of uncertainty contributions from the fit parameters associated with the charm-quark mass adjustment are listed in the column "κ tuning" in Table XI.

Renormalization and matching uncertainty
We include terms of O(α 2 s ) in our base chiral-continuum fit with unknown coefficients ρ [2] ij constrained to be of O(1) to incorporate the uncertainty due to omitted higher-order renormalization and matching terms. The sum of uncertainty contributions from the fit parameters ρ [2] ij are listed in the "matching" column in Table XI. 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 α s = 0.2 from our finest lattice spacing and ρ [2] ij = 1, which yields 6.5% for all operators.
We check that the inclusion of generic O(α 2 s ) terms captures the uncertainty from truncating the perturbative expansion in α s by performing two alternate fits: one with only the known renormalization and matching terms of O(α s ), and another with terms through O(α 3 s ) and coefficients constrained to be of O(1). The results of these fits are labeled "mNPR" and "mNPR+α 3 s ", respectively, in Figs. 14 and 15. In both cases, the shift in central value is small. Without the O(α 2 s ) 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 O(α 3 s ) 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 α s is obtained either from the fourth-root of the plaquette or from the gauge-fixed Landau link. In these fits, labeled "PT P + α 2 s " and "PT L + α 2 s ", respectively, in Figs. 14 and 15, we include O(α 2 s ) terms constrained as in the base fit. Here we observe larger changes in the matrix-central values, which are still less than 2σ fit away from the base-fit results. The "PT P + α 2 s " and "PT L + α 2 s " results themselves differ by almost 1σ 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.

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 {β i , β i } 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 XI.
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 χPT. In the second fit, we add all possible analytic terms of O(N 3 LO). The results of these fits are labeled "NLO (m q < 0.65m s )" and "N 3 LO", respectively, in Figs. 14 and 15. In both cases, the central values for the matrix elements shift only slightly. The NLO χPT fit without NNLO analytic terms underestimates the errors on the matrix elements and also has a significantly larger χ 2 aug /dof than the other variations shown. The errors do not increase from those of the base fit, however, with the inclusion of N 3 LO 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 χ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.

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 XI.
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∆ are needed to maintain invariance under variation of the chiral scale. These terms scale as α 2 s a 2 . Generic one-loop contributions from improving the gluon and light-quark actions at tree-level, however, give rise to discretization terms of O(α 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 O(α 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.

Heavy quark discretization errors
We estimate the uncertainties from heavy-quark discretization via the uncertainty in the coefficients z i of the O(α s a, a 2 , a 3 ) terms in Eqs. (7.11)-(7.13). This error is labeled "HQ disc" in Table XI. We also perform two fits, each of which includes fewer heavy-quark terms than in the base fit. These are labeled "HQ O(α s a) only" and "HQ O(α 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. (7.11)-(7.13) taking Λ HQET = 500 MeV for the heavy-quark scale, and coefficients z i set by the combinatoric factors 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 XI, which range from 2.4% to 3.6%, but the similar size suggests that errors obtained from the fit are reasonable.

Relative scale (r1/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 XI. B. Additional errors

Absolute scale (r1) 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 χ 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 XII.

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 XII.

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 averagem = (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. [6], we estimate isospin-breaking effects to enter at O((m sea d − m sea u ) 2 ), leading to a negligible ∼ 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 O(α EM /π) ∼ 0.2%, again following Sec. VIII.B.4 of Ref. [6]. We add this error from the omission of electromagnetism in quadrature to the fit error and list it in the "EM" column of Table XII.

C. 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. [6], and are estimated from power counting to be ∼1-2% for reasonable choices of α s and Λ QCD . To be conservative, we take the upper end of the range, or 2%, for the "Charm sea" error in Table XII. We note, however, that for the decay constants f π , 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.

D. 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 "no a ≈ 0.12 fm" and "no a ≈ 0.045 fm", respectively. The resulting matrix-element central values agree with the base fit within 1σ 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σ fit in some cases, which are covered by equally substantial increases in the uncertainties for all matrix elements except O 3 . In this case, the individual-fit  After considering all possible significant sources of uncertainty in the previous sections, Table XII 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 A. In Sec. IX B, 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.

A. Matrix elements
Table XIII presents our final results for the relativistically normalized D-meson mixing hadronic matrix elements of operators O i (i = 1-5) in Eqs. (2.14)-(2.18) including statistical and all sources of systematic uncertainty. We give the matrix elements in the MS-NDR scheme at the scale µ = 3 GeV obtained with the BBGLN [53] and BMU [55] choices of evanescent operators. Although the choice of evanescent operators affects only the renormalization of operators O 2 and O 3 , the numerical results for the other three matrix elements O 1,4,5 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 XVI. 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 [12] 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 [8] and N f = 2 + 1 + 1 [5] twisted-mass fermions. ETM presents values for bag parameters, which we then convert to matrix elements using Eq. (3) of Ref. [8], with their own N f = 2 and N f = 2 + 1 + 1 calculations of the quark masses [73,75] and decay constants [74,76]. 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 O 1 , O 2 , O 3 agree with the matrix elements converted from ETM's bag parameters to within about 1-2σ, but those for O 4 and O 5 differ by 1.7-3.3σ (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 [1,69], we obtain matrix element results that agree with ours within 0.25-1.52σ for all five operators. Use of the averages from the Flavor Lattice Averaging Group [77] 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 [78,79]. 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 [5,[80][81][82], which the authors of Ref. [82] attribute to the choice of intermediate renormalization scheme. Results obtained in Ref. [82] using the symmetric regularization-independent momentum-subtraction (RI-SMOM) scheme [83] are in good agreement with an independent calculation that employs one-loop mean-field improved lattice perturbation theory [81], but differ substantially with results obtained with the RI-MOM [84] renormalization scheme. This discrepancy is attributed to underestimated systematic errors present in the RI-MOM scheme [82,83]. We note that operator O 4 has the largest negative anomalous dimension of the five ∆C = 2 operators, by about 50%, and is therefore most sensitive to running of the renormalization scale. Because O 4 and O 5 mix under renormalization, this also affects the results for O 5 . Thus, O 4 and O 5 are likely to be the most sensitive to difficulties in renormalization.
FIG. 16. Comparison of the three-flavor D-mixing matrix elements obtained in this work (filled symbols) with the two-and four-flavor results from the ETM Collaboration [5,8] (unfilled symbols). For the ETM results, we have converted their quoted bag parameters to matrix elements using their two-and four-flavor quark masses and decay constants from Refs. [73][74][75][76]. The total uncertainty quoted from this work does not include the error from quenching the charm sea quark. On the ETM points, the larger red error bars include the uncertainties from fD and mq in quadrature, while the smaller blue error bars omit those uncertainties. ETM's two-flavor results do not include an estimate of the error due to quenching the strange sea quark.

B. Implications for new physics
As discussed in Sec. II, neutral D-meson mixing is a sensitive probe of local ∆C = 2 interactions from physics beyond the Standard Model that contribute to the quantity M 12 . The phase φ 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 Γ 12 , which is true in many BSM models [7], 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 Λ NP . To evaluate the contribution to M 12 using Eq. (2.22), we must run the Wilson coefficients down to the scale µ = 3 GeV at which our matrix elements are renormalized. We use the one-loop running derived in Ref. [7], which works in a different four-fermion operator basis than the one used here. Fierz identities relate the two bases [19]: 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. [7] yields the renormalization equations in our basis, which we then use directly. The operator running depends on the value of the strong coupling α s at intermediate scales.
We use the RunDec Mathematica library [85] to compute α 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 O 5 at an ultraviolet scale Λ NP , with purely imaginary (and hence CP-violating) Wilson coefficient Running down to 3 GeV, we find Im C NP 5 (3 GeV) = and all other Wilson coefficients at 3 GeV remain zero. Thus, operators O 4 and O 5 contribute a purely imaginary term to from new physics. (Here, x NP 12 is defined as in Sec. II, x NP 12 = M NP 12 /Γ D .) We can now determine the constraint on Λ 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 Λ NP = 40 000 TeV; this level of new contribution is consistent with the experimental bounds at 1σ. Finally, the red box shows a new-physics contribution with Λ NP = 18 000 TeV, which is ruled out by experiment with very high confidence. 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 where F i and L i are flavor and loop-counting factors [86]. In order to obtain bounds on Λ NP individually by operator, we set one C NP i (Λ i,NP ) = F i L i /Λ 2 i,NP at a time, with all other Wilson coefficients set to zero at the high scale. We then determine the value of Λ i,NP for which the new-physics prediction for |x 12 |e iφ12 is inconsistent with the experimental bound at 95% confidence, including the uncertainty in the matrix elements for O i . Here we assume that Im F i L i > 0, but taking the other sign gives a nearly identical constraint. (The experimental bound is Im x 12 < 0.0289% (Im x 12 > −0.0285) in the positive (negative) imaginary direction.) Using our matrix elements in the BBGLN scheme, we obtain Λ 1,NP (Im F 1 L 1 ) 1/2 × 7 630 TeV, (9.12) Λ 2,NP (Im F 2 L 2 ) 1/2 × 24 100 TeV, (9.13) Λ 3,NP (Im F 3 L 3 ) 1/2 × 23 100 TeV, (9.14) Λ 4,NP (Im F 4 L 4 ) 1/2 × 48 500 TeV, (9.15) Λ 5,NP (Im F 5 L 5 ) 1/2 × 26 900 TeV, (9.16) Note that these bounds are from the imaginary part of x 12 only; for new physics with Im F i L i ≈ 0, the constraint on Λ NP from the Re x 12 will dominate. Our bounds in Eqs. (9.12)-(9.16) are stronger than those quoted by the ETM Collaboration [8], in part because we use more recent, tighter experimental bounds. For operators O 3 and O 5 in particular, our constraints on Λ NP are much higher (by factors of roughly 7 and 3, respectively). These two operators mix strongly with O 2 and O 4 , respectively, such that their bounds stem principally from C 2 (3 GeV)O 2 (3 GeV) and C 4 (3 GeV)O 4 (3 GeV). If we artificially set C 2 (3 GeV) = 0, then we obtain much weaker bounds of Λ 3,NP 3 330 TeV and Λ 5,NP 9 700 TeV, respectively, which are close to the values quoted by ETM in Ref. [8].
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 [87]. A Higgs coupling of the form Y uc hū L c R + Y cu hc L u R will induce ∆C = 2 four-fermion interactions at low energy. After integrating out the Higgs boson h, the effective Hamiltonian is We write the two Yukawa couplings as Y uc = |Y uc |e iφuc , and similarly for Y cu . Taking the Wilson coefficients to be purely imaginary and comparing the resulting Im x 12 to the experimental bounds, we find from the O 2 and O 4 terms separately the bounds |Y uc | 2 + |Y cu | 2 1.04 × 10 −10 , (9.21) |Y cu Y uc | 2.15 × 10 −11 . (9.22) To be somewhat more general, we can instead marginalize over the phases in the couplings (integrating from 0 to π/2) and obtain exclusion contours in the |Y uc |-|Y cu | plane, again solely from the stronger experimental bound on Im x 12 . The contours are shown in Fig. 18. Our bounds on the combinations of Yukawa couplings in Eqs. (9.21)-(9.22) are tighter than those in Ref. [87] by over an order of magnitude. The improvement stems primarily from the use of newer experimental measurements than employed by the authors of Ref. [87], who relied on the same data as the original model-independent analysis of D-meson mixing [15]. 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 ∆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 [5,8]. Our results for D 0 |O i |D 0 (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. (2.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 φ 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 ∆C = 2 interactions, finding Λ NP (Im F i L i ) −1/2 10-50 × 10 3 TeV for the five local ∆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, Λ NP , and run down to 3 GeV.
The long-distance contributions described by the second term on the right hand side of Eq. (2.10) are expected to dominate the Standard Model prediction of the D 0 -meson mixing observables, ∆M and ∆Γ, 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. [5,8] 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 Γ 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 [88][89][90][91][92][93], 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 [94][95][96][97][98]. These techniques have been extended to neutral kaon mixing, and the first lattice-QCD calculations have recently become available [99,100]. 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 [101][102][103][104][105][106][107][108][109][110][111][112][113], 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 [114]. 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.
Table XV 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. (4.1), ∆ k,j is defined in Eq. (4.7), while the parameters Z Oi nm , i = 1-5, are defined in Eqs. (4.2)-(4.3). For the energy splitting the prior is defined such that exp(∆ k,j ) is Gaussian distributed.

Appendix C: Correlations among matrix-element results
Table XVI provides the correlations among the D-meson mixing matrix elements for all five operators, to enable their use in future phenomenological studies.   Table XIII; 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.