  Observation of neutrinoless double beta decay, a lepton number violating
process that has been proposed to clarify the nature of neutrino masses, has
spawned an enormous world-wide experimental effort. Relating nuclear decay
rates to high-energy, beyond the Standard Model (BSM) physics requires detailed
knowledge of non-perturbative QCD effects. Using lattice QCD, we compute the
necessary matrix elements of short-range operators, which arise due to heavy
BSM mediators, that contribute to this decay via the leading order $\pi^- \to
\pi^+$ exchange diagrams. Utilizing our result and taking advantage of
effective field theory methods will allow for model-independent calculations of
the relevant two-nucleon decay, which may then be used as input for nuclear
many-body calculations of the relevant experimental decays. Contributions from
short-range operators may prove to be equally important to, or even more
important than, those from long-range Majorana neutrino exchange.
