In the case of a globally hyperbolic spacetime, what you want is a smooth Cauchy temporal function (the gradient is everywhere timelike, not just causal, and each level set is a Cauchy surface that is necessarily spacelike). That global hyperbolicity is also sufficient the the existence of a smooth temporal function was also shown by Bernal and Sanchez, in the followup [arXiv:gr-qc/0512095] to their original paper that you cite.
I'm not sure if you are also interested in foliations of non globally hyperbolic spacetimes. I'm not completely sure what the right conditions would be then.