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 π − → π + 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.Introduction.-Neutrinoless double beta decay (0νββ) is a process that, if observed, would reveal violations of symmetries fundamental to the Standard Model, and would guarantee that neutrinos have nonzero Majorana mass [1, 2]. Such decays can probe physics beyond the electroweak scale and expose a source of leptonnumber (L) violation which may explain the observed matter-antimatter asymmetry in the universe [3,4]. Existing and planned experiments will constrain this novel nuclear decay [5][6][7][8][9][10][11][12][13][14][15][16], but the interpretation of the resulting decay rates or limits as constraints on new physics poses a tremendous theoretical challenge.The most widely discussed mechanism for 0νββ is that of a light Majorana neutrino, which can propagate a long distance within a nucleus. However, if the mechanism involves a heavy scale, Λ ββ , the resulting L-violating process can be short-ranged. While naïvely short-range operators are suppressed compared to long-range interactions due to the heavy mediator propagator, in the case of 0νββ, the long-range interaction requires a helicity flip and is proportional to the mass of the light neutrino. In a standard seesaw scenario [17][18][19][20][21], this light neutrino mass is similarly suppressed by the same large mass scale, so the relative importance of long-versus short-range contributions is dependent upon the particle physics model under consideration and in general cannot be determined until the nuclear matrix elements for both types of processes are computed.Both long-and short-range mechanisms present substantial theoretical challenges if we hope to connect high energy physics with experimentally observed decay rates. The former case is difficult because one must understand long-distance nuclear correlations. In the latter case the short-distance physics is masked by QCD effects, requiring non-perturbative methods to match few-nucleon matrix elements to Standard Model operators.Effective field theory (EFT) arguments show that at leading order (LO) in the Standard Model, there are nine local four-...