  Starting from a manifestly Lorentz-invariant chiral Lagrangian, we present a
model-free prediction for the transition amplitude of the process
$nn\rightarrow pp e^-e^-$ induced by light Majorana neutrinos, which is a key
process of the neutrinoless double beta decay ($0\nu\beta\beta$) in heavy
nuclei employed in large-scale searches. Contrary to the nonrelativistic case,
we show that the transition amplitude can be renormalized at leading order
without any uncertain contact operators. The predicted amplitude defines a
stringent benchmark for the previous estimation with model-dependent inputs,
and greatly reduces the uncertainty of $0\nu\beta\beta$ transition operator in
the calculations of nuclear matrix elements. Generalizations of the present
framework could also help to address the uncertainties in $0\nu\beta\beta$
decay induced by other mechanisms. In addition, the present work motivates a
relativistic {\it ab initio} calculation of $0\nu\beta\beta$ decay in light and
medium-mass nuclei.
