  Nuclear matrix elements (NME) for the most promising candidates to detect
neutrinoless double beta decay have been computed with energy density
functional methods including deformation and pairing fluctuations explicitly on
the same footing. The method preserves particle number and angular momentum
symmetries and can be applied to any decay without additional fine tunings. The
finite range density dependent Gogny force is used in the calculations. An
increase of $10\%-40\%$ in the NME with respect to the ones found without the
inclusion of pairing fluctuations is obtained, reducing the predicted
half-lives of these isotopes.
