  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.
