We present new results for the phase space factors involved in double beta decay for β − β − transitions to ground states and excited 0 + 1 states, for isotopes of experimental interest. The Coulomb distortion of the electron wave functions is treated by solving numerically the Dirac equation with inclusion of the finite nuclear size and electron screening effects, and using a Coulomb potential derived from a realistic proton density distribution in the daughter nucleus. Our results are compared with other results from literature, obtained in different approximations, and possible causes that can give differences are discussed. Introduction Within the Standard Model (SM) the double beta decay (DBD) can occur through several decay modes, but the only measured at present is that with the emission of two electrons and two antineutrinos (2νβ − β − ) and which conserves the lepton number. However, beyond SM theories allows this process to occur without emission of neutrinos as well, and this possibility makes DBD a nuclear process of major interest for testing the lepton number conservation (LNC) and for understanding the neutrino properties. There are recent excellent reviews, containing also a comprehensive list of references [1] - [5], where the reader can find a complete information on this subject. The DBD lifetimes can be factorized, in a good approximation, as follows: