After the introduction of the ionization-injection scheme in Laser Wake Field Acceleration and of related high-quality electron beam generation methods as two-color or the Resonant Multi Pulse Ionization injection, the theory of thermal emittance by C. Schroeder et al, has been used to predict the beam normalised emittance obtainable with those schemes. In this manuscript we recast and extend such a theory, including both higher order terms in the polinomial laser field expansion and non polinomial corrections due to the onset of saturation effects in a single cycle. Also, a very accurate model for predicting the cycle-averaged $3D$ momentum distribution of the extracted electrons, including saturation and multi-process events, is proposed and tested. We show that our theory is very accurate for the selected processes of Kr$^{8^+rightarrow10^+}$ and Ar$^{8^+rightarrow10^+}$, resulting in a a maximum error below $1%$ even in deep saturation regime. This highly accurate prediction of the beam phase-space can be implemented e.g., in laser-envelope Particle in Cell (PIC) or hybrid PIC-fluid codes, to correctly mimic the cycle-averaged momentum distribution without the need of resolving the intra-cycle dynamics. Finally, we introduce further spatial averaging with Gaussian longitudinal and transverse laser profiles, obtaining expressions for the whole-beam emittance that fits with Monte Carlo simulations in a saturated regime, too.