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.