We develop a hybrid formalism suitable for modeling scalar field dark matter, in which the phase-space distribution associated to the real scalar field is modeled by statistical equal-time two-point functions and gravity is treated by two stochastic gravitational fields in the longitudinal gauge (in this work we neglect vector and tensor gravitational perturbations). Inspired by the commonly used Newtonian Vlasov-Poisson system, we firstly identify a suitable combination of equal-time two-point functions that defines the phase-space distribution associated to the scalar field and then derive both a kinetic equation that contains relativistic scalar matter corrections as well as linear gravitational scalar field equations whose sources can be expressed in terms of a momentum integral over the phase-space distribution function. Our treatment generalizes the commonly used classical scalar field formalism, in that it allows for modeling of (dynamically generated) vorticity and perturbations in anisotropic stresses of the scalar field. It also allows for a systematic inclusion of relativistic and higher order corrections that may be used to distinguish different dark matter scenarios. We also provide initial conditions for the statistical equal-time two-point functions of the matter scalar field in terms of gravitational potentials and the scale factor.

This version matches the one accepted for publication in Physical Review D. There are some minor changes with respect to the previous version