Impact of isovector pairing fluctuation on neutrinoless double-beta decay in multi-reference covariant density functional theory

We extend the multi-reference covariant density functional theory (MR-CDFT) by including fluctuations in quadrupole deformations and average isovector pairing gaps simultaneously for the nuclear matrix elements (NMEs) of neutrinoless double-beta $(0\nu \beta \beta )$ decay in the candidate nuclei ${}^{76}$Ge, ${}^{82}$Se, ${}^{100}$Mo, ${}^{130}$Te, and ${}^{136}$Xe assuming the exchange of either light or heavy neutrinos. The results indicate a linear correlation between the predicted NMEs and the isovector pairing strengths, as well as the excitation energies of $2_1^{+}$ and $4_1^{+}$ states. By adjusting the pairing strengths based on the excitation energies of the $2_1^{+}$ states, we calculate the NMEs for $0\nu \beta \beta$ decay, which are reduced by approximately $12\mathrm{\%}$ to $62\mathrm{\%}$ compared to the results obtained in the previous studies by Song et al. [Phys. Rev. C95, 024305 (2017)]. Additionally, upon introducing the average isovector pairing gap as an additional generator coordinate in the calculation, the NMEs increase by a factor ranging from $56\mathrm{\%}$ to $218\mathrm{\%}$.