The idea is the following:

• take the integrated energy density in the gravitational strain, as expressed in units of the closure density <math tex>\Omega_{GW}</math> and given, for instance, by Barack and Cutler (it is necessary to remove the all-sky response averaging of 2/5 from their formulas, as well as the LISA-LIGO factor of 3/4)
• attribute this energy density to sources in the galaxy, proportionally to the binary density <math tex>\rho_0 \exp(-R/R_0) {\rm sech}^2(-z/z_0)</math> used by Nelemans et al.
• create a pixelization of the sky using Healpix, and use it to divide the energy incoherently among pseudo-random GW sources (two for each pixel, for the two incoherent polarizations)
• the pseudo-random GW sources are given by white noise, filtered to achieve the overall background strain spectrum of <math tex>f^{-7/3}</math>; the time-domain FIR (finite impulse response) filter used for this has 64 points, and was designed using MATLAB
• currently I follow this procedure to compute the power that must be given to each pixel; I generate a random R, z, and azimuthal angle according to the distribution given above (for R and z I use a transformation method from a uniform variate); I compute Galactic-center-centered Galactic coordinates for this “binary”; I transform to Earth/Sun-centered Galactic coordinates; I transform to Sun-Centered ecliptic coordinates; I attribute <math tex>1/d^2</math> in power, where d is the distance to the binary, to the Healpix pixel that contains the Ecliptic latitude and longitude; I repeat this, Monte Carlo fashion, for many million sources; I obtain the fractional power distribution between all the Healpix pixels in the sky at a given Healpix refinement

Discussion:

• this scheme may sound complicated, but the alternative within the simulators would be to attribute the GWDB signal directly to the phasemeter measurements, or to the TDI observables, which raises the problem of how to correctly model the correlations between these

Open points:

• the low-frequency behavior of this background is (correctly) to become buried in the noise because of the diminishing LISA response. The high-frequency behavior, instead, is to follow the <math tex>f^{-7/3}</math> slope until the shot-noise floor is reached. One may want to cut off the background at intermediate frequencies (by modified the filter) to model the individual resolution of binaries...