  Nuclear matrix elements (NMEs) for neutrinoless double-beta
($0\nu\beta\beta$) decay in candidate nuclei play a crucial role in
interpreting results from current experiments and in designing future ones.
Accurate NME values serve as important nuclear inputs for constraining
parameters in new physics, such as neutrino mass and the Wilson coefficients of
lepton-number-violating (LNV) operators. In this study, we present a
comprehensive calculation of NMEs for $0\nu\beta\beta$ decay in $^{76}$Ge,
$^{82}$Se, $^{100}$Mo, $^{130}$Te, and $^{136}$Xe, using nuclear wave functions
obtained from multi-reference covariant density functional theory (MR-CDFT). We
employ three types of transition potentials at the leading order in chiral
effective field theory. Our results, along with recent data, are utilized to
constrain the coefficients of LNV operators. The results demonstrate that the
combined NMEs based on the Feynman diagrams at the hadronic scale for the
nonstandard mechanisms lead to uncertainty by different nuclear models
comparable to that for the standard mechanism. The use of NMEs from various
nuclear models does not dramatically change the parameter space intervals for
the coefficients, although MR-CDFT yields the most stringent constraint.
Furthermore, our NMEs can also be used to perform a more comprehensive analysis
with multiple isotopes.
