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.