Next-to-leading-order prediction for the neutrinoless double-beta decay

The neutrinoless double-beta decay ($0\nu\beta\beta$) of two neutrons$nn \rightarrow ppee$ is the elementary subprocess of $0\nu\beta\beta$ decay in nuclei. Accurate knowledge of the $nn \rightarrow ppee$ amplitude is required to pin down the short-range contributions in the nuclear matrix elements of the candidate nuclei for large-scale $0\nu\beta\beta$ searches. In this Letter, we report the first next-to-leading-order prediction of the nn \rightarrow ppee amplitude, with Bayesian uncertainty quantification. This is made possible by the development of the relativistic chiral effective field theory, in which no unknown contact term is required up to next-to-leading order. The theory is validated by reproducing in a parameter-free way the available data on the charge independence and charge symmetry breaking contributions in the two-nucleon scattering. The present work makes an essential step towards addressing the uncertainty in the theoretical calculations of the nuclear matrix elements relevant for $0\nu\beta\beta$ searches.