  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.
