  The amplitude for the neutrinoless double $\beta$ ($0\nu\beta\beta$) decay of
the two-neutron system, $nn\to ppe^-e^-$, constitutes a key building block for
nuclear-structure calculations of heavy nuclei employed in large-scale
$0\nu\beta\beta$ searches. Assuming that the $0\nu\beta\beta$ process is
mediated by a light-Majorana-neutrino exchange, a systematic analysis in chiral
effective field theory shows that already at leading order a contact operator
is required to ensure renormalizability. In this work, we develop a method to
estimate the numerical value of its coefficient in analogy to the Cottingham
formula and validate the result by reproducing the charge-independence-breaking
contribution to the nucleon-nucleon scattering lengths. Our central result,
while derived in the $\overline{\text{MS}}$ scheme, is given in terms of the
renormalized amplitude $\mathcal{A}_\nu(|\mathbf{p}|,|\mathbf{p}^\prime|)$,
matching to which will allow one to determine the contact-term contribution in
regularization schemes employed in nuclear-structure calculations. Our results
thus greatly reduce a crucial uncertainty in the interpretation of searches for
$0\nu\beta\beta$ decay.
