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. The possible detection of lepton number violating processes such as neutrinoless double beta decay (0νββ) is one of the current main goals for particle and nuclear physics research. In this process, an atomic nucleus decays into its neighbor with two neutron less and two proton more emitting only two electrons. Fundamental questions about the nature of the neutrino such as its Dirac or Majorana character, its absolute mass scale as well as its mass hierarchy can be determined if this process is eventually measured [1]. On the one hand, searching for 0νββ decays represents an extremely difficult experimental task because an ultra low background is required to distinguish the predicted scarce events from the noise. Recently, the controversial claim of detection in 76 Ge by the Heidelberg-Moscow (HdM) collaboration [2] has been overruled by the latest data released by . Nevertheless, these results are challenging the experiments that are already running or in an advanced stage of development to detect directly this process [3,[6][7][8][9][10][11][12][13][14]. On the other hand, in the most probable electroweak mechanism to produce 0νββ, namely, the exchange of light Majorana neutrinos [1,15], the half-life of this process is inversely proportional to the effective Majorana neutrino mass m ν , a kinematic phase space factor G 01 and the nuclear matrix elements M 0ν (NME):where m e is the electron mass and m ν = | k U 2 ek m k | is the combination of the neutrino masses m k provided by the neutrino mixing matrix U . The kinematic phase space factor can be determined precisely from the charge, mass and the energy available in the decay [16] while the nuclear matrix elements must be calculated using nuclear structure methods. The most commonly used ones are the quasiparticle random phase approximation [17][18][19][20][21] (QRPA), large scale shell model [22][23][24] (LSSM), interacting boson model [25,26] (IBM), projected HartreeFock-Bogoliubov [27] (PHFB) and energy density functional [28][29][30] (EDF). In recent years, most of the basic nuclear structure aspects of the NMEs have been understood within these different frameworks. In particular, the decay is favored when the initial and final nuclear states have similar intrinsic deformation [28,30,31]. Indications [18,21,23,28,30] about the strong sensitivity of the transition operator to pairing correlations suggest that fluctuations in this degree ...