  Neutrinoless double-beta ($0\nu\beta\beta$) decay is a heretofore unobserved
process which, if observed, would imply that neutrinos are Majorana particles.
Interpretations of the stringent experimental constraints on
$0\nu\beta\beta$-decay half-lives require calculations of nuclear matrix
elements. This work presents the first lattice quantum-chromodynamics (LQCD)
calculation of the matrix element for $0\nu\beta\beta$ decay in a multi-nucleon
system, specifically the $nn \rightarrow pp ee$ transition, mediated by a light
left-handed Majorana neutrino propagating over nuclear-scale distances. This
calculation is performed with quark masses corresponding to a pion mass of
$m_\pi = 806$ MeV at a single lattice spacing and volume. The statistically
cleaner $\Sigma^- \rightarrow \Sigma^+ ee$ transition is also computed in order
to investigate various systematic uncertainties. The prospects for matching the
results of LQCD calculations onto a nuclear effective field theory to determine
a leading-order low-energy constant relevant for $0\nu\beta\beta$ decay with a
light Majorana neutrino are investigated. This work, therefore, sets the stage
for future calculations at physical values of the quark masses that, combined
with effective field theory and nuclear many-body studies, will provide
controlled theoretical inputs to experimental searches of $0\nu\beta\beta$
decay.
