A complete calculation of phase space factors (PSFs) for Majoron-emitting 0νβ − β − decay modes is presented. The calculation makes use of exact Dirac wave functions with finite nuclear size and electron screening and includes lifetimes, single-electron spectra, summed electron spectra, and angular electron correlations. Combining these results with recent microscopic interacting boson model nuclear matrix elements (NMEs) we make half-life predictions for the ordinary Majoron decay (spectral index n=1 ). Furthermore, comparing theoretical predictions with the obtained experimental lower bounds for this decay mode we are able to set limits on the effective Majoron-neutrino coupling constant ⟨g M ee ⟩.