  Neutrinoless double beta decay (\( 0 \nu \beta \beta \)) is a hypothetical
nuclear decay mode with important implications. In particular, observation of
this decay would demonstrate that the neutrino is a Majorana particle and that
lepton number conservation is violated in nature. Relating experimental
constraints on \(0 \nu \beta \beta\) decay rates to the neutrino masses
requires theoretical input in the form of non-perturbative nuclear matrix
elements which remain difficult to calculate reliably. This work marks a first
step toward providing a general lattice QCD framework for computing
long-distance \(0 \nu \beta \beta\) matrix elements in the case where the decay
is mediated by a light Majorana neutrino. The relevant formalism is developed
and then tested by computing the simplest such matrix element describing an
unphysical \( \pi^{-} \rightarrow \pi^{+} e^{-} e^{-} \) transition on a series
of domain wall fermion ensembles. The resulting lattice data is then fit to
next-to-leading-order chiral perturbation theory, allowing a fully-controlled
extraction of the low energy constant governing the transition rate,
\(g_{\nu}^{\pi \pi}(\mu = 770 \,\, \mathrm{MeV}) = -10.78(12)_{\rm
stat}(51)_{\rm sys}\). Finally, future prospects for calculations of more
complicated processes, such as the phenomenologically important \(n^{0} n^{0}
\rightarrow p^{+} p^{+} e^{-} e^{-}\) decay, are discussed.
