Stochastic gravitational wave background from stellar origin binary black holes in LISA

We use the latest constraints on the population of stellar origin binary black holes (SOBBH) from LIGO/Virgo/KAGRA (LVK) observations, to estimate the stochastic gravitational wave background (SGWB) they generate in the frequency band of LISA. In order to account for the faint and distant binaries, which contribute the most to the SGWB, we extend the merger rate at high redshift assuming that it tracks the star formation rate. We adopt different methods to compute the SGWB signal: we perform an analytical evaluation, we use Monte Carlo sums over the SOBBH population realisations, and we account for the role of the detector by simulating LISA data and iteratively removing the resolvable signals until only the confusion noise is left. The last method allows the extraction of both the expected SGWB and the number of resolvable SOBBHs. Since the latter are few for signal-to-noise ratio thresholds larger than five, we confirm that the spectral shape of the SGWB in the LISA band agrees with the analytical prediction of a single power law. We infer the probability distribution of the SGWB amplitude from the LVK GWTC-3 posterior of the binary population model: at the reference frequency of 0.003 Hz it has an interquartile range of h2ΩGW(f = 3 × 10-3 Hz) ∈ [5.65, 11.5] × 10-13, in agreement with most previous estimates. We then perform a MC analysis to assess LISA's capability to detect and characterise this signal. Accounting for both the instrumental noise and the galactic binaries foreground, with four years of data, LISA will be able to detect the SOBBH SGWB with percent accuracy, narrowing down the uncertainty on the amplitude by one order of magnitude with respect to the range of possible amplitudes inferred from the population model. A measurement of this signal by LISA will help to break the degeneracy among some of the population parameters, and provide interesting constraints, in particular on the redshift evolution of the SOBBH merger rate.


Introduction
Stellar-origin binary black holes (SOBBHs) are among the targets of the Laser Interferometer Space Antenna (LISA) [1]. The emission of gravitational waves (GWs) from these binaries crosses the mHz frequency band, probed by LISA, while they are still far from coalescence. Given the recent constraints from the LIGO/Virgo/KAGRA (LVK) collaboration after three observation runs, we expect a large population of such systems contributing to the LISA data stream [2]. At least a few of these binaries will be individually detected [3][4][5][6][7][8], while the bulk of them will form a stochastic GW background (SGWB), as they are too faint/distant and/or because they produce long-lived overlapping time-domain signals. The characterization of both resolved and unresolved SOBBH sources is compelling since they are a source of confusion for other detectable sources in the LISA band. For example, the SOBBH SGWB contribution will act as a foreground for the detection of a possible signal of cosmological origin [9], see e.g. [10][11][12][13] for prospects about the detectability of a cosmological SGWB in the presence of a SGWB of astrophysical origin.
Previous papers have estimated the expected level of the SOBBH background. This can be achieved via direct extrapolation of the LVK observed merger rate, supplemented by a simple modelling of the black-hole (BH) population and of the time delay between the binary formation and the merger taken from [14][15][16], as done e.g. in [3,10,17,18]. In particular, [10,18] use the LVK observed event rate from the Gravitational Wave Transient Catalog 2 (GWTC-2) [19]. Alternatively, one can input more refined scenarios of BH formation from the evolution of different populations of stars, accounting for the cosmic chemical evolution, optical depth to reionisation, and metallicity of the interstellar medium, to evaluate the mass distribution of merging SOBBH and in turn the expected SGWB, as done e.g. in [20][21][22][23][24][25]. An estimate of the number of resolvable SOBBH in LISA using the GWTC-2 rate has been done e.g. in [6].
In this paper, we employ several methods to estimate the SGWB in the LISA band, using the most recent population constraints from the Gravitational Wave Transient Catalog 3 (GWTC-3) [2]. We evaluate the impact, on the SGWB amplitude, of the observational uncertainty on the population parameters, taken from the posterior parameter sample of GWTC-3: we find that the SGWB amplitude can vary by as much as a factor of five. When considered independently, we show that the parameter whose marginalised 2σ error influences the most the SGWB level is the power-law index of the redshift dependence of the merger rate. We also assess LISA's capability to detect and characterise the predicted SOBBH SGWB via a Monte Carlo (MC) analysis of simulated data, including the SGWB, the galactic binary (GB) foreground component, and the instrumental noise. The maximal marginalised error on the SGWB amplitude by LISA is ∼ 5%, i.e. much smaller than the variation due to the present (GWTC-3) observational uncertainty on the population parameters: this hints to the conclusion that LISA will have a role to play in constraining SOBBH population parameters via the SGWB measurement. Though future Earth-based GW detectors observations will improve on the GWTC-3 constraints by the time LISA flies, we expect that LISA will maintain an impactful constraining power, since the SGWB amplitude in the LISA band is influenced by the high-redshift behaviour of the merger rate, complementary to what will be accessible to ground-based detectors in the near future.
The paper is organised as follows. In Section 2 we describe the population model that we use, and the assumptions we consider, to construct the SOBBH catalogues. In particular, we disregard eccentricity in the waveform, as well as any redshift dependence of the population parameters, and we adopt a uniform distribution for the time-to-coalescence in the detector frame (Section 2.1). Faint and distant SOBBH contribute to the SGWB signal: we, therefore, need to complete the GWTC-3 merger rate, limited to low redshift, with a model for the star formation and evolution at higher redshift. As explained in Section 2.2, we assume that the merger rate tracks the cosmic star formation rate (SFR) up to high redshift [26]. We evaluate the impact of a time-delay between the binary star formation and the BBH merger on the SGWB amplitude in Section 4.2. In Section 2.3 we describe the other population parameters: for the mass distribution we adopt the Power Law + Peak model, and for the spin amplitudes a positive-exponents Beta distribution [19]; as for the remaining parameters, some of them are randomly generated (i.e. time-to-coalescence, initial phase, position in the sky, inclination, and polarization), whereas others are derived analytically (e.g. the initial frequency of the generated events, their distance...). We have also produced ten SOBBH catalogues at a benchmark fixed point in the population parameter space, that we use for consistency studies; their characteristics are presented in Section 2.4.
In Section 3 we present the four methods we have used to compute the expected SOBBH background signal. In order of sophistication: (i) the first procedure is based on an analytic evaluation of the characteristic strain as an integral over the number density of inspirals, as first proposed by [27] (Section 3.1); (ii) we then substitute the integral over the number density by an MC sum over a realisation of a population, (iia) first as a time-to-coalescenceaveraged sum, (iib) and then taking into account the time-to-coalescence of individual events and binning them according to their corresponding emission frequencies (Section 3.2); (iii) finally, in order to account for the actual detection process of the SOBBHs by LISA, we apply the iterative-subtraction method developed in [28], for which at each step we compare the signal-to-noise ratio (SNR) of each source i (ρ i ) to an SNR threshold (ρ 0 ), and if ρ i > ρ 0 , the source is classified as resolvable and is subtracted from the data. The iterative subtraction is performed on realistic LISA data-streams produced by injecting the time domain, spinning wave-form signals of the events, one by one. The latter procedure, despite being computationally expensive, yields a very accurate representation of the LISA data and allows for the evaluation of both the residual SGWB level and the subtracted sources (which we analyse in a companion paper [29]).
In Section 4 we present our results. We first check that the four methods give comparable SGWB levels (Section 4.1): since the number of subtracted sources is small [22,29], there is overall very good agreement. Method (i), i.e. the analytic integration of the background, while not capturing some detailed features of the signal, can safely be used to estimate the expected SGWB in the LISA band, for all points in the posterior parameter sample of the fiducial FidLVK model: the results are given in Section 4.2. In Section 4.3 we present the results of the MC analysis of simulated LISA data including the SOBBH SGWB, the galactic binary foreground, and the instrumental noise. We show that, also in presence of the GB foreground, with four years of data, LISA will be able to detect the SOBBH signal and reconstruct its amplitude and spectral index. Then, accounting for the estimated SOBBH signal and the GB foreground as extra noise contributions, in Section 4.4 we build the LISA power-law sensitivity (PLS) [10,[30][31][32]. Finally, in Section 4.5, we analyse how the precise measurement of the SOBBH SGWB by LISA would impact the inference on population parameters, as put forward in [33,34]. We find that the effect is most promising for the merger rate parameters, i.e. amplitude, and power-law index. We conclude in Section 5.

SOBBH population model
LISA is sensitive to the GW emission by SOBBHs in the inspiral phase. Within the timescale of the mission, which we assume of 4 years (i.e. 4.5 years with 89% duty cycle), the GW frequency emitted by most SOBBHs will slowly increase within the LISA frequency band, f ∈ 10 −4 , 0.1 Hz. A minority of SOBBHs will chirp (i.e. their GW emission will rapidly increase in frequency) and move throughout the band. Among the chirping SOBBHs, a fraction will be close to coalescence, so that the frequency of their GW emission will exit the LISA band and, shortly after, enter the ground-based detectors band, where they will merge. This opens up the possibility of multi-band observations and/or of archival analyses (see e.g. [3,[35][36][37]). On the other hand, no SOBBH entering the LISA band during the lifetime of the mission is statistically expected, as SOBBH with frequencies of the order of 10 −4 Hz are practically monochromatic during the lifetime of the experiment. 1 1 To give an example, by integrating the Newtonian relation dfGW/dt = 96/5 π 8/3 (GM/c 3 ) 5/3 f 11/3 GW one obtains that it takes about 10 8 years to shift the GW emission of an SOBBH with chirp mass M = 50M from 2 · 10 −5 Hz to 10 −4 Hz, where it will still be about 10 6 years away from the merger; while the same We aim at estimating the SGWB due to unresolved SOBBHs in LISA, accounting for the most recent population constraints from GWTC-3 [2]. For this aim, we generate catalogues of SOBBHs emitting in the LISA band, making some simplifying assumptions.
First of all, for simplicity, we neglect eccentricity in our analysis. LVK measurements poorly constrain the SOBBH eccentricities, but the eccentricity in the LISA band could be significant depending on the binary formation [38,39]. In addition, we neglect a possible dependence on redshift of the population parameters, since there are no strong constraints on how the SOBBH parameter distributions should vary with redshift, and state-of-the-art studies based on observations have not found conclusive evidence on the presence of any redshift dependence [34,40,41] (this possibility has been explored e.g. in [42]). Our methodology can incorporate a redshift dependence into the catalogue generation (albeit at a higher computational cost), if this will be constrained by future data.
Furthermore, we assume that the residual time to coalescence τ c (i.e. the amount of time that an observer in the source frame must wait in order to see the binary merge) is statistically uniformly distributed across the SOBBH population. This amounts to assuming that the formation, and therefore the coalescence rates, are in a steady state. Indeed, any change in the demographics of the binaries happens on a cosmic time-scale of O(10 9 ) yrs, i.e. much longer than the LISA observation time, which is the typical time over which our catalogues are representative. Furthermore, the maximal τ c that we consider in this analysis c,max ∼ O(10 4 ) yrs in the detector frame, also much smaller than the timescale over which the cosmic coalescence rate varies. We also neglect the possibility that the SOBBHs form on such a tight orbit that their GW emission at formation is already within the LISA band; this would indeed also break the uniform distribution hypothesis for τ c .
The above assumptions allow us to model the SOBBH population as follows. We consider the binaries emitting in the LISA band and observed by the detector at a given instant t, i.e. the time at which LISA switches on. Note that for the sake of the argument, we take this time in the source frame. Among the intrinsic parameters (masses, spins, phase, polarisation. . . ) and extrinsic ones (sky position, inclination. . . ) of each SOBBH, we single out the time-to-coalescence τ c = t c − t in the source frame (where t c denotes the time of coalescence of a given SOBBH) and the redshift of the source z, while ξ represents the remaining parameters. The population model parameterized in terms of some hyper-parameters θ provides the statistical distributions p( ξ| θ) of ξ (for simplicity we omit the vector symbol on ξ and θ from now on). The number of SOBBHs with given z, τ c , ξ, whose signals reach the interferometer at time t, is where V c is the universe's comoving volume, 2 and Within our assumptions, all values of τ c are equiprobable at any z: the rate density satisfies therefore R(z, τ c ) = R(z, τ c = 0), i.e. the one of merging SOBBHs. Moreover, in our binary will shift from 0.1 Hz to 1 Hz (where it will be about 16 minutes away from the merger) in about 5 days. 2 The redshift derivative of the comoving volume in Equation (2.1) accounts for the fact that spherical shells further from us enclose increasing amounts of volume and thus larger numbers of events for a given R(z). population, the number of SOBBHs with given z and ξ that are received by the interferometer at the times t and t + dt, are precisely N (z, τ c = 0, ξ) and N (z, τ c = 0 + dt, ξ), so one can equivalently interchange dτ c ↔ dt in Equation (2.2). All together, it follows R(z, τ c ) = d 2 N (z, τ c = 0, ξ, θ)/(dV c dt), which is precisely the merger rate density in the form that LVK is constraining [41]. Hereafter, we drop the τ c dependence in R(z, τ c ) as irrelevant.
Given our assumption that the merger rate R(z) is in a steady state, we can readily apply LVK findings to it. In the next section, we explain how we use Equation (2.1) to generate SOBBH catalogues compatible with the latest population constraints from LVK GWTC-3. However, the merger-events-based LVK constraints on R(z) are limited to small redshift, while we need to model sources also at high redshift, since they have a significant contribution to the SGWB. In order to simulate the high-redshift part of the SOBBH population, we therefore need to incorporate knowledge of the star formation and evolution at high redshift, as we will see below.

Implementing GWTC-3 posterior for the SOBBH population parameters
The SOBBH population model is determined by the merger rate R(z) and the distribution function p(ξ|θ) in Equation (2.1). In [41], the LVK collaboration has analysed a series of population models and produced inference on their parameters, finding that the most promising one to explain the SOBBH events gathered in GWTC-3 [2] is characterised by: (a) a powerlaw dependence of the merger rate with redshift, R(z) = R(0)(1 + z) κ ; (b) a population mass model, known as Power Law + Peak mass model, combining an inverse power-law dependence on the largest BH mass, with a Gaussian peak at approximately 30-40 M , and a power-law distribution for the mass ratio of the binary; (c) a population spin model in which the amplitudes are independent and follow positive-exponent Beta distributions favouring intermediate-valued spins, and whose tilt distribution is a mixture of an isotropic distribution and a truncated Gaussian. The distributions for masses and spins are explained in more detail in Appendix A.2.
For the sake of convenience, we will call this combination FidLVK, and we will use it as the fiducial model in our analysis. We provide population-averaged predictions of the SGWB, based on the publicly-available population parameter posterior distribution [43] for the FidLVK model conditioned to the SOBBH of GWTC-3 [41] (excluding low-masssecondary GW190814 and likely-NSBHB GW190917, as per the fiducial approach by LVK). The parameters θ of the mass and spin distribution p(ξ|θ) are imported directly from the LVK results. On the other hand, the parameters of the merger rate R(z) require a different treatment. Because of the interferometers frequency band, LVK probe the SOBBH population only at relatively low redshift, so the redshift dependence of the merger rate in LVK analyses is modelled as a power-law. Indeed, the GWTC-3 inferred merger rate posterior constrains the pivot rate R(0) and the power-law exponent κ only for z 0.5 [41]. In order to produce, from this posterior, SOBBH SGWB estimates valid in the LISA band, we need to extend the merger rate model towards higher redshift, since high redshift SOBBHs contribute significantly to the background.
For this purpose, we adopt a phenomenological approach and assume that the merger rate tracks the Madau-Fragos SFR [26], neglecting the presence of a time delay between the binary formation and merger. While in Section 4.2 we discuss the impact on the SGWB amplitude of including time delays, in the rest of the paper we parameterize the merger rate as where C ensures R 0 ≡ R(z = 0). The analysis of [26] finds the following best fit values for the SFR parameters: κ = 2.6, r = 3.6 and z peak = 2.04 (note the difference in the definition of z peak with respect to [34]). At redshift z 1, this behaves similarly to the R(z) = R 0 (1 + z) κ power law constrained by LVK, which finds a best fit κ ≈ 2.7. Motivated by this agreement, we incorporate the LVK GWTC-3 posterior into Equation (2.3) by matching R 0 and κ for each point in the population parameter sample with the fixed fiducial values r = 3.6 and z peak = 2.04 from [26]. The resulting posterior for the merger rate is by construction fully compatible with that of the low-redshift merger rate of LVK (see Appendix A.1 for further discussion and in particular Figure 14). 3 Finally, in order to keep consistency with LVK [41], we adopt the ΛCDM cosmological model with parameters fixed accordingly to the "Planck 2015 + external" data combination [44]. 4 These correspond to H 0 = 67.9 km/(s Mpc) for the local Hubble rate, and Ω m ≈ 0.3 and Ω Λ ≈ 0.7 for the matter and cosmological-constant energy densities. The cosmological model enters in the differential comoving volume per unit redshift, dV c (z)/dz, of Equation (2.1), and in the computation of the cosmological distances needed for the integration in Section 3.

SOBBH population synthesis
In Section 3 we propose four different methods to compute the SGWB due to unresolved SOBBHs. The first method consists of an integration of the number density in Equation (2.1) [27]. The other three are based on the superposition of the GW signals from SOBBHs populations, with different levels of sophistication. The latter methods provide a more refined evaluation of the SGWB and of its spectral shape and are also important to assess the size of statistical effects (e.g. the uncertainty due to the population realisation) and the consequences of other choices inherent to the catalogue simulation, such as the value of the maximal time-to-coalescence τ (det) c,max , see below. We thus need fast and reliable SOBBH population synthesis. We have written two independent population synthesis codes, which can be found in the following repositories: [45] and [46]. These have been also compared with [47]. In these implementations, the masses and spins are drawn from the LVK GWTC-3 distributions, briefly revised in Appendix A.2.
The redshift of the binaries is generated independently as an inhomogeneous Poisson point process, according to the z-dependent terms in Equation (2.1), between z min = 10 −5 (≈ 45 kpc of comoving distance, in order to exclude binaries within the Milky Way), and z max = 5, which is sufficient for an accurate SGWB computation, as we demonstrate in Section 3.1. Note that we will limit z max = 1 in the analyses based on catalogues whenever using a larger z max would prove too costly from the computational point of view, c.f. Section 3.3.
3 In [34], the LVK Collaboration also considers a similar high-redshift extension and finds mild constraints on r and z peak (using a different definition of the latter, see Appendix A.1) by combining the population parameters inferred from GWTC-2 resolved mergers with the upper limits imposed by the non-detection of the SGWB. We verify a posteriori the compatibility of our results with the upper limits on the SGWB amplitude presented in [34]. 4 Note that the ΛCDM parameter values used in [41] correspond to the incomplete Planck 2015 data combination plikHM_TE (high-T×E spectrum data only) instead of the fiducial plikHM_TTTEEE+lowTEB, which includes temperature-only and polarized data for the whole Planck multipole range. The difference is anyway negligible for the purposes of this paper.  Table 1: Priors for the parameters of individual SOBBHs. The uniform prior for the time-tocoalescence, which is source-dependent, is justified in Section 2.1. The priors on the ecliptic coordinates and the inclination impose statistical isotropy in the positions and orientation of the binaries.
The rest of the individual SOBBH parameters are generated from the priors presented in Table 1, based on physical considerations: isotropy for the sky position, inclination, and polarization; and uniform time-to-coalescence in the detector frame, as discussed in Section 2.1. From the randomly-sampled parameters, we compute the derived quantities necessary for the problem at hand, such as the frequency at the start of the LISA runtime, the LISA in-band time, cosmological distances, and so on.
The upper limit for the population synthesis time-to-coalescence τ (det) c,max needs to be high enough to give a faithful representation of the SOBBH SGWB signal in the LISA band (at least where it is the dominant contribution to the astrophysical-origin SGWBs), and at the same time it is conditioned by computational limitations. As discussed in Appendix A.3, τ (det) c,max = 10 4 yrs provides a good balance between these requirements.

Benchmark fixed-point catalogues for consistency studies
In addition to probabilistic GWTC-3-posterior forecasts, we also single out a fixed point in the population parameter space, which we use as a benchmark to compare different SGWB computation methods and assess the size of statistical and numerical effects. For this fixed point, we use values close to the median of the GWTC-3 FidLVK model posterior, indicated in Table 2, with an important modification.

Rate of events R(z)
Mass distribution Spin distribution  The determination of the mass range population parameters m min and m max in the LVK study [41] is sensitive to whether certain events from GWTC-3 (the extreme mass ratio binary GW190814 and the likely-NSBHB binary GW190917) are considered as outliers and excluded from the analysis. The inclusion of GW190814 suffices to push the lower mass bound down to m min ≈ 2.5 M .
Motivated by the possibility that such outliers may appear in future data, we enlarge the mass range for the benchmark fixed-point catalogues to [2.5, 100] M , also raising the upper bound up to the original prior boundary of the m max population parameter, previous to GTWC-3 constraints. The modification in particular of the lower mass boundary necessitates a further increase in the width of the low-mass smoothing function (see Appendix A.2), in order not to deviate strongly from the mean GWTC-3 mass probability density at masses m m min . This is achieved by increasing the value of the δ min population parameter (nevertheless, we have found that the SGWB calculation is not very sensitive to this choice).
Though this fixed point in the population parameter space does not pertain to the GWTC-3 posterior, due to the modifications to the mass distribution, it leads to an SGWB in the LISA band which is compatible with our posterior-based evaluations. It is therefore useful as a benchmark to gauge the sensitivity of the SGWB predictions to different assumptions on the population model.
We have generated a sample of 10 catalogues with parameters set to this benchmark fixed-point. Since they will be used only to compare and validate population-based SGWB computation methods, we have limited their redshift range to z max = 1, to reduce computational cost. As stated before, we limit the time-to-coalescence in these catalogues to τ

Computation of the SOBBH signal in the LISA band
We adopt four different methods to evaluate the SOBBH SGWB which allow us, by their different nature, to capture different features of the signal. In the following sections, we describe them.

Method (i): analytical evaluation
In this section we provide a brief description of the formalism employed for the analytic evaluation of the SOBBH SGWB, following [27]. The normalised SGWB energy density spectrum per logarithmic unit of frequency Ω GW (f ) can be defined from the total GW energy density present in the universe and emitted by the whole SOBBH population, expressed in the detector frame. Recalling Equation (2.1) this reads where ρ c = 3H 2 0 c 2 /(8πG) is the Universe's critical energy density and ρ (event) GW = t 00 denotes the energy density associated to a single SOBBH event, at the detector. Using [48], one can where all quantities are at the source: M is the SOBBH chirp mass, a · r its physical distance in the local wave zone, ι its orientation with respect to the detector, and ω s = πf s its orbital frequency. The second equality in Equation (3.2) has been obtained under the approximation of quasi circular motion for the binaryḟ s f 2 s , and we have averaged over the waveform phase. Substituting Equation (3.2) in Equation (3.1), expressing the differential comoving volume as dV c = c d 2 andΩ is the solid angle [49], and noting that a 2 r 2 = d 2 M /(1 + z) 2 , one gets, after integration over the solid angle (giving a factor 16π/5), where we have also used definition (2.1). One can change the integration variable from τ c to f s using the relation df s /dτ c = 96/5 π 8/3 (GM/c 3 ) 5/3 f 11/3 s , valid for quasi circular binaries in the Newtonian approximation, then change to the frequency at the detector f = f s /(1 + z), and equate the integrands in Equation (3.1), to obtain the SGWB energy density power spectrum: Among the set of binary parameters ξ, only the chirp mass is relevant within the Newtonian approximation. One can therefore express the SOBBH SGWB today as with f * an arbitrary pivot frequency, and In Figure 1 we plot the amplitude of the expected background at the reference frequency f * = 0.003 Hz, close to LISA's peak sensitivity [50], evaluated from the integral in Eq. (3.6). For the merger rate R(z), we adopt the phenomenological parameterization described in Section 2.2. While in Figure 1 we consider z max = 30, as we do not expect active sources at higher redshifts, in Figure 2, in contrast, we analyse the relative difference of considering smaller values for z max . The only other population parameters that enter the analytic evaluation are the masses m 1 , m 2 of the two compact objects, expressed in terms of  Table 2. The red arrow in κ indicates that the amplitude of the signal for κ max (h 2 Ω GW (f * ) 2.8 × 10 −12 ) is beyond the range of the plot. The signal amplitude grows with decreasing m min , contrary to all other parameters. the chirp mass: as previously stated, for their probability distribution p (m 1 , m 2 ), we adopt the Power Law + Peak model. Naturally, the amplitude of the background depends on the choice of the parameters in R(z) and p(m 1 , m 2 ): respectively, (R 02 ≡ R(z = 0.2), κ), and (α, δ min , m min , m max , λ peak , µ, σ, β), defined in Appendices A.1 and A.2. Following Section 2.4, we plot as a horizontal solid line the SGWB amplitude for the values indicated in Table 2 for each of the population parameters ; we also show the range of SGWB amplitudes (grey bars) obtained when each of the parameters is varied within its 5-95 percentile range according to the GWTC-3 posterior (see Appendices A.1 and A.2), while the rest of the parameters stay fixed to their values of Table 2. This shows how the different parameters in the model influence the SGWB amplitude when varied individually. Larger ranges for the SGWB amplitude translate into stronger constraining power from the measurement of the individual parameter; however, since this approach neglects degeneracies, large ranges for multiple parameters do not mean that these can be simultaneously constrained (the issue of SGWB-derived population parameter constraints will be further discussed in Section 4.5). Note that the signal amplitude grows with decreasing m min , contrary to its response to all other parameters. The population parameter with the largest impact on the SGWB amplitude is the power-law index of the merger rate κ, since its value controls the merger rate growth at intermediate redshift 1 < z < z peak , which strongly influences the outcome of the redshift integration in Equation (3.6). The red arrow in Figure 1 indicates that the SGWB amplitude obtained for κ max , h 2 Ω GW (f * ) 2.8 × 10 −12 , is beyond the range of the plot.
In Figure 2 we plot the relative percentage change of the SGWB amplitude when varying z max in Equation (3.6). Since the merger rate in Equation (2.3) decays at high redshift, the SGWB grows asymptotically towards a constant amplitude as we integrate over larger and larger redshift ranges. Taking z max = 30 as a reference, we plot for different values of z max . The figure indicates that integrating up to z max = 5 already allows to obtain ∼ 1% accuracy in the calculation of the SGWB amplitude. This is sufficient for the scope of this paper given that, as presented in Section 4.3, the typical error on the SGWB measurement by LISA is larger than that. The SGWB amplitude in Figure 2 has been evaluated adopting the parameter values corresponding to the benchmark fixed-point described in Section 2.4; the convergence trend is very similar when using the minimum or maximum values of a given parameter, keeping the others fixed to their central values.

Methods (iia) and (iib): Monte Carlo sum
An alternative method to compute the SGWB is to sum the GW signals, emitted by individual SOBBHs, over a realisation of the population drawn from the distribution represented by the number density. The simplest implementation consists in factoring out from Equation (3.4) the population number density of Equation (2.1), averaged over time-to-coalescence and skyposition, to obtain where f is the observed frequency, and M i , z i and d M,i are the chirp mass in the source frame, redshift and metric distance of the individual GW sources. The factor 1/τ c,max = (1+z)/τ (det) c,max comes from the time-averaging of the number density of Equation (2.1).
The SGWB amplitude resulting from the sum over a realisation of the SOBBH population is obviously realisation-dependent. We can assess its concordance, within the variance . The two methods are equivalent within the range of the population variance (∼ 0.2%), which is in turn much smaller than the integration error due to the z max = 5 cut (∼ 1%, see Figure 2). Such small effects will not be observable by LISA (see Section 4.3).
due to the population draws, with the analytical computation of Equation (3.6) by evaluating Equation (3.7) for a large number of realisations. The result is shown in Figure 3, assuming population parameters fixed to the benchmark values described in Section 2.4, and setting z max = 5. The population variance in terms of the ratio of the interquartile range to the mean of the realisations' amplitudes amounts to 0.2% only. 5 The difference between the SGWB amplitude obtained by averaging the realisations, and the one obtained by the numerical integration of Equation (3.6), is much smaller than the population variance, highlighting the equivalence between the two methods. Furthermore, the population variance uncertainty is much smaller than both the expected integration error due to fixing z max = 5, and the forecasted precision of the LISA measurement (see Section 4.3). A more refined approach to evaluating the SGWB can be obtained by summing the contribution of each SOBBH in the population, accounting for the actual frequency of emission of each source (while in Equation (3.7) only the chirp mass and the distance -equivalently redshift -pertain to the individual events). In order to do this, we rewrite the SGWB energy density starting from Equation (3.1), but re-expressing the number density as the number of events per unit of emission frequency f s using the relation df s /dτ c for quasi-circular Newtonian binaries, then changing the integration variable to the observed frequency, and equating the integrands in Equation (3.1) to single out the SGWB power spectrum: We can now express the integral in the above equation as an MC sum, as done previously in Equation (3.7), but this time computing the sum of the GW energy density emitted by every SOBBH per (detector-frame) frequency bin, where the latter is defined using some frequency sampling δ f as [(j − 1)δ f , jδ f ], N j being the subset of a population with emission frequencies (in detector frame) in bin j where the different powers of the per-source quantities with respect to Equation (3.7) can be explained by the frequency dependence of the number of sources in each bin. This assumes monochromatic sources, ignoring frequency drifting during the life of the mission. 6 The largest contribution to the background is produced by sources with f ∈ (10 −3 , 10 −2 ), whose frequency drifting is small; we can therefore choose e.g. the frequency with which they enter the LISA band (see Appendix A.3). We will show the result of both MC integrations, Eqs. (3.7) and (3.9), in Section 4.1.

Method (iii): iterative subtraction
The methods presented above are based on summing the signals of the SOBBH in the population, without accounting for the actual detection process, apart from restricting the maximal time-to-coalescence τ (det) c,max to a computationally manageable and detector-compatible value (for methods (iia) and (iib) of Section 3.2). However, we are ultimately interested in the SGWB signal in LISA, and the detector sensitivity can influence the SGWB spectral shape/amplitude. In order to consider such aspects, we also evaluate the SGWB following the methodology developed in [28], using ideas first presented in [52][53][54]. The procedure is based on generating LISA data-streams, by computing the waveform signals of all the events within the simulated population. Depending on the adopted waveform model, this can yield a very accurate representation of the LISA data, as far as SOBBHs are concerned. However, simulating millions of sources is computationally expensive, thus one has to allocate a considerable amount of computational resources to this task.
The procedure begins by fixing the mission duration T obs , here set to 4 years, and generating the signal to be measured by LISA. We compute the h + and h × waveforms for each source of the simulated catalogue, and then we project them onto the LISA arms. We use the IMRPhenomHM model [55], which describes spinning, non-precessing binaries. It is based on the IMRPhenomD [56,57] model, but it includes higher order modes. We use the lisabeta software [58,59] for our computations. When generating each waveform, we also compute their SNR in isolation, ρ iso i , with respect to the instrumental noise only, which will be used to reduce the computational requirements of the procedure, as explained below.
Next, we estimate the total power spectral density (PSD), S n, k , summing all the GW sources plus the instrumental noise. The index k refers to the iterative step. Since this PSD is very noisy, we compute its running median to produce a smoother version of it. We then evaluate the SNR ρ i of each source i using the smoothed S n, k as the total "noise" PSD. Note Frequency [Hz] that, to speed up the computation, this is performed only on the subset of sources with sizable pre-computed SNR in isolation ρ iso i (see [28]). The SNR ρ i are then compared to a threshold SNR ρ 0 : if ρ i > ρ 0 , the source is classified as resolvable, and is subtracted from the data. The smoothed residual PSD S n, k+1 is then re-evaluated after re-iterating through the catalogue of sources and subtracting the loud ones, and the procedure is repeated until the algorithm converges. Convergence is reached when all the sources are subtracted given the ρ 0 threshold, or if S n, k+1 and S n, k are practically identical at all frequencies considered. At the end of the procedure, we compute the final SNR of the recovered sources, with respect to the final estimate of S n, k final . Thus, as final products, we get both the SGWB due to the sources signal confusion, as well as the properties of the recovered sources (their number, waveform parameters, and final SNR).
Different realisations of the same population (with the same number density parameters) should yield different, though statistically compatible, sets of subtracted events, but a similar SGWB after smoothing. We have verified this statement by evaluating the SGWB from the 10 benchmark catalogues presented in Section 2.4; the result is shown in Figure 4a.
The crucial parameter of the iterative method is ρ 0 , the minimum SNR above which events are considered resolvable and thus subtracted from the total signal. We consider ρ 0 = 8 an appropriate choice [60,61], assuming that stochastic methods to sample the sources parameter space, more efficient than grid-based methods [62], can be used to analyse the LISA data streams. Archival searches will allow to further reduce the SNR threshold down to ρ 0 = 5 [37,61]. As can be appreciated from Figure 4b, as long as ρ 0 5, the number of detectable sources is too small to alter the shape and amplitude of the residual SGWB spectrum [3,22] (see also [29]). Our prediction for the SGWB level is therefore robust with respect to our choice of setting ρ 0 = 8. On the other hand, if values of ρ 0 4 will be justified in the context of future improvements in data analysis methods, or of archival searches using future groundbased detector data [63], the residual SGWB spectral shape must be adapted: as can be seen in Figure 4b, it no longer follows the analytical estimation of Section 3.1, which does not account the presence of the detector, but a dip on its amplitude appears at high frequencies.
Note that we have assumed uninterrupted measurement over the time frame T obs , and the instrumental noise, taken from [50,64], is assumed to be ideal, i.e. Gaussian and stationary. We also subtract each resolvable source from the data at its injection parameters, meaning that we generate "perfect residuals", or in other words, we neglect the uncertainty on the source parameters, which inevitably arises within the parameter estimation procedure. We, therefore, simulate an optimal case of the global fit scheme for the LISA SOBBHs. The above assumptions, while not totally realistic, allow us to simplify the analysis.

Comparison between SGWB computation methods in the LISA band
In this section we show the effect of fixing a maximal time-to-coalescence for the simulated populations on the SGWB spectral shape, and compare the SGWB signals resulting from the four methods described in Sections 3.1 to 3.3. As a benchmark, we use one of the fixedpoint catalogues presented in Section 2.4. The redshift range is limited to z ∈ [0, 1] (comoving distance up to ≈ 3 GPc) to guarantee the computational feasibility of the iterative-subtraction method. The amplitude of the SGWB signals shown in this section is therefore reduced (cf. Section 4.2), but this plays no role in the purpose of the tests performed here.
As discussed in Sections 2.3 and 2.4, in order to limit computational costs, synthetic populations are generated including events up to a maximum time-to-coalescence, that we fix to τ (det) c,max = 10 4 yrs in the detector frame. In order to investigate the effect of this assumption, one of the catalogues among the benchmark ones has been generated with τ (det) c,max = 1.5 × 10 4 yrs, and from it we have produced two sub-catalogues with τ (det) c,max = 1.0 × 10 4 and 5 × 10 3 yrs. The SGWBs inferred from these catalogues via the iterative-subtraction method are shown in Figure 5: excluding all sources beyond a given τ (det) c,max (appropriately redshifted in the source frame), results in a non-physical bending of the SGWB at low frequencies, depending on the maximal time-to-coalescence (in agreement with [28], see also Appendix A.3). It is therefore important to pick a value for τ (det) c,max ensuring a minimal loss of information while keeping the computational cost of generating the SGWB manageable: as discussed in Appendix A.3, we consider τ (det) cmax = 10 4 to be a good compromise. In Figure 6 we show the SGWBs computed using the three methods based on population synthesis, presented in Sections 3.2 and 3.3. The results are in very good agreement, for both the SGWB amplitude and spectral shape. In particular, those of the frequency-binned MC sum (method (iib)) and of the iterative subtraction (method (iii)) also follow the single powerlaw behaviour f 2/3 predicted by the analytical evaluation (Equation (3.4)), and taken over by the averaged power-law-like MC sum (Equation (3.7)). As far as the frequency-binned MC sum is concerned, this shows that our population catalogues are complete. As far as the iterative-subtraction method is concerned, instead, this is a consequence of the simulated detection process: the instrument sensitivity is such that the number of resolvable sources is c,max values on the resulting stochastic signal, computed using the iterative subtraction (method (iii)). The red dashed line represents the analytical result (method (i)). Imposing a maximum time-to-coalescence in generating the synthetic populations suppresses early-phase inspirals, producing a cutoff in the SGWB at low frequencies. This is not a physical effect, but a limitation of the population synthesis: the spectrum tends towards the expected power law as the upper limit in time-to-coalescence grows. , and the light orange one the frequency-binned sum (iib). The blue curves show the SGWB evaluated with the iterative subtraction (iii), Section 3.3, for two different data smoothing methods: in s1 (light blue) we have performed a running median over the PSD data using a rolling window of 1000 points, whereas in s2 (dark blue) we apply an additional Gaussian filter. The signals from the frequency-binned MC sum and from the iterative subtraction share some features, especially at low frequencies, where the differences due to neglecting the drifting and using simplified waveforms are less important. Both follow closely the SGWB of the averaged, power-law-like MC sum. Probability density LVK GWTC-3  too small, even at high frequencies, to alter the SGWB spectral shape, as already pointed out in [22] (see also [29]). The signals from the frequency-binned MC sum and from the iterative subtraction share some features, especially at low frequencies, despite the fact that the former uses simplified waveform and does not account for frequency drifts. Both approaches also follow closely the averaged power-law-like MC sum, which is distributed around the analytical calculation of the background, from Equation (3.6) (see Figure 3).

Expected SOBBH signal in the LISA band from GWTC-3
Having established the consistency of the four methods, we turn to the actual computation of the expected SGWB in the LISA band, based on the present knowledge about the SOBBH population. To this purpose, we rely on Equation (3.5) and evaluate the SGWB amplitude by integrating Equation (3.6) for all points in the LVK posterior parameter sample that is publicly available [43] for the FidLVK model [41], following the prescriptions described in Section 2.2. The distribution of the SGWB amplitude at the reference frequency f = 3 × 10 −3 Hz is shown in Figure 7 (blue solid line). On a logarithmic scale, it follows a lightly-right-skewed distribution with median h 2 Ω GW (f = 3 × 10 −3 Hz) = 7.87 × 10 −13 , and has an interquartile range of h 2 Ω GW (f = 3 × 10 −3 Hz) ∈ [5.65, 11.5] × 10 −13 .
The computation of the SGWB amplitude has been performed under the assumption that the merger rate inherits the functional redshift dependence of the SFR, Equation (2.3). As discussed in Section 2.1 and Appendix A.1, the agreement between the values of the merger rate parameters inferred from GW observations [41] with those of the SFR inferred from electromagnetic observations [26] supports this assumption at low redshift z 1.5. However, the merger rate remains untested at higher redshifts, and it is therefore important to investigate how much this assumption influences the final SGWB result. We do so by analysing one example of a more refined model for the merger rate, introducing a time delay t d between the formation of star binaries and their evolution into BBH systems. The merger rate is then given by the convolution of the SFR with the probability distribution of the time delay [20,[65][66][67]:  Figure 8) [41]. Looking at Equation (3.6), we see that (for redshiftindependent mass models) the redshift-dependent contribution to the background amplitude can be factored out. We can thus easily compute, for a given mass model, the ratio f between a time delayed model and our fiducial SRF case. For t d,min = 50 to 500 Myr, we find that they agree within 40%: accounting for the time delay, therefore, provides SGWB amplitudes close to the P 5 percentile of the median fiducial (SFR-extended) case (see Figure 7). The level of agreement drops to 36% for t d,min = 1 Gyr; however, from Figure 8, we can appreciate that the corresponding merger rate is rather in tension with LVK constraints. Our results, in terms of translating the GWTC-3 population constraints into a forecast for the SOBBH SGWB in the LISA band, appear to be robust within one order of magnitude: the highest contribution to the background comes in fact from the SOBBH population at z 1.5, for which the merger rate is well constrained by LVK GWTC-3. Note that all derived SGWB amplitudes fall well within LISA's detection capabilities (see Section 4.3). A more thorough study of the dependency of the SGWB amplitude on physically-motivated models for the merger rate and mass distribution can be found in [68].
Our results are also compatible with the latest constraints on the SGWB amplitude by LVK [34]. The upper bound on a power-law SGWB with spectral index 2/3 at f = 25 Hz is, at 95% credible level, 3.4 × 10 −9 (1.2 × 10 −8 ), when using a log-uniform (uniform) prior, which becomes in the LISA band h 2 Ω GW (f = 3 × 10 −3 Hz) < 3.8 × 10 −12 (1.3 × 10 −11 ). This upper bound applies to the total background, which contains other contributions together with the SOBBH confusion signal (for example the one from neutron star binaries). The actual limit on the SOBBH SGWB is therefore expected to be smaller. Nevertheless, the SGWB amplitude that we forecast remains compatible, being smaller than the LVK upper limit at 99% probability, with median amplitude being smaller by a factor of five (see Figure 7).
We also compare our results to a few other predictions for the SOBBH SGWB in the LISA band given in the literature, see Figure 9. In [69], taking into account early LVK constraints (from the first 6 events) for the merger rate, a time delay distribution p(t d ) ∝ 1/t d with t d,min = 50 Myr, and a different fiducial model for the mass distribution from the one used   Figure 9: Comparison between the SGWB amplitude posterior from this work (blue shaded area, also shown in Figure 7), with the median SGWB value evaluated accounting for timedelays (blue crosses, left-to-right in decreasing value of t d,min , cf. Figure 8), and with other recent results from the literature (black lines and crosses). The grey band to the right represents the LVK upper bound, derived in [34] for a power-law SGWB with index α = 2/3.
here, it was found that h 2 Ω GW (f = 3×10 −3 Hz) = 1.25 +1. 3 −0.7 ×10 −12 (90% credible level), which lies in the upper-half of our distribution (see Figure 7). In [23], the authors compute both the isotropic SOBBH SGWB component and its anisotropy, and find a lower prediction than in our analysis: h 2 Ω GW (f = 3 × 10 −3 Hz) ∈ [1.0, 1.9] × 10 −13 , for a number of astrophysicsmotivated models for the merger rate, adjusted to LVK GWTC-1 constraints. The latest LVK forecast [34], using the merger rate and the mass distribution inferred from GWTC-2, and the usual time-delay distribution, results in h 2 Ω GW (f = 3 × 10 −3 Hz) = 5.6 +1.9 −1.6 × 10 −13 (90% credible level), which is consistent with our results both when including and not including time delays. The analysis of [10], also based on the LVK GWTC-2 population model, uses powerlaw mass functions and the conventional time-delay distribution, and obtains h 2 Ω GW (f = 3 × 10 −3 Hz) = 2.9 +1.7 −1.4 × 10 −13 (approx 95% credible level): this prediction is compatible with our results, but towards the low side of the distribution in Figure 7. In [22], the authors use the population code Star-Track to model the binary formation, treating separately the contributions from population I/II and population III stars. The SGWB amplitude from SOBBHs formed by population I/II star is h 2 Ω GW (f = 3 × 10 −3 Hz) = 1.2 × 10 −12 , which lays in the upper part of our probability distribution. Population III stars contribute an additional 2%, h 2 Ω GW (f = 3 × 10 −3 Hz) = 2.25 × 10 −13 : since this is significantly larger than the expected uncertainty in LISA's measurement of the background (see Section 4.3), the presence of population III stars could be discriminated, provided that the population is known with sufficient certainty. Finally, in [70] it is found that the contribution of SOBBH to the SGWB is even lower than what found in [23], and subdominant in the LISA band with respect to the one from primordial BHs: h 2 Ω GW (f = 3 × 10 −3 Hz) 4.5 × 10 −14 .

SGWB Parameter Estimation
In this section we assess LISA's capability to detect and characterise the SOBBH SGWB. We perform an MC analysis of simulated data containing the instrumental noise, the stochastic foreground from binaries in the Galaxy, and different levels of the SOBBH SGWB, corresponding to the percentiles presented in Figure 7. The SOBBH SGWB is modelled following Equation (3.5), but both the amplitude and the spectral tilt are left as free parameters in the analysis: We apply a pre-processing procedure similar to the one employed in [32,71], which we briefly summarize here: assuming a mission duration of 4 years, we chunk the data stream into N c segments of 11.5 days each (corresponding to a frequency resolution ∆f 10 −6 Hz); we generate data in the frequency domain for each segment, including the instrumental noise, the GB foreground, and the SOBBH SGWB, and we average over these segments to get the simulated measured spectrum. Using the noise as an estimate for the variance, we define a likelihood consisting of a sum of Gaussian and log-normal components (the latter accounting for the skewness of the exact likelihood), as discussed in [32]. For the sake of speed and without loss of precision, this likelihood is applied to a coarse-grained version of the spectrum obtained by inverse variance weighting, the final data in frequency space being defined as where f k ij and D k ij are the coarse-grained frequencies and data respectively. Ω GW represents both the SOBBH component, with spectral shape defined by Equation (4.2), and the GB foreground component, based on the model from [28]. θ s is the vector of parameters of the signal: amplitude and spectral tilt of the SOBBH SGWB, while we reconstruct only the amplitude h 2 Ω Gal of the GB foreground. Ω n is the instrumental noise in omega units. We adopt a two-parameter noise model as typically done for LISA: the noise is characterized at low frequency by the acceleration component, parameterised by A, and at high frequency by the interferometric component, parameterised by P [71]. The two noise parameters form the vector θ n , and vary freely in our analysis. We sample over the joint ( θ s , θ n ) = (log 10 [h 2 Ω GW (f = 3 × 10 −3 Hz)], α, log 10 [h 2 Ω Gal ], A, P ) parameter space using the Nested Sampler Polychord [72,73] via its interface with Cobaya [74].
In Figure 10 we show the MC contours (2-σ contours) on the SOBBH signal parameters (log 10 (h 2 Ω GW (f * )), α), together with the parameters of the GB foreground and the noise (log 10 (h 2 Ω Gal ), A, P ), obtained by injecting each of the SOBBH SGWB posterior percentiles shown in Figure 7. For all the injected SGWB amplitudes, the reconstruction of both the signals and the noise is accurate, with all parameters consistent with the injected values at 2σ. In particular, the simultaneous reconstruction of the GB and SOBBH SGWB is achievable even when the amplitude of the latter is small, due to their different spectral shapes.
Rather than sampling over the tilt α, as we did in the present background-detection study, in a realistic data analysis pipeline searching for the SOBBH SGWB, the tilt would be fixed to α = 2/3. Thus, LISA's determination of the background amplitude could reveal more accurate, with respect to the tilt-marginalised errors presented here. On the other hand, realistic data would contain the contribution from all the other GW sources in the LISA band, which need to be extracted simultaneously to the SGWBs signals, possibly affecting the error on the SGWB amplitude compared to the simple MC evaluation performed here (see e.g. [75]).
As a sanity check, for the lowest value of the background amplitude, we have also performed a Fisher parameter estimation. In Figure 11 we present the comparison between the Fisher analysis and the corresponding MC result, showing that the two procedures are consistent in the reconstruction accuracy of the signal and noise parameters.
The results of this section show that LISA will be able to narrow down by one order of magnitude the current uncertainty on the SGWB amplitude due to the SOBBH population uncertainty inferred from GWTC-3 (see Figure 7). Moreover, we demonstrated that a clear detection of the SGWB is guaranteed, if the true signal falls within this uncertainty range. On the other hand, the lack of detection, or the detection of an SGWB outside the posterior prediction (likely lower), would indicate either that the population model needs to be changed, for example modifying the merger rate behaviour at high redshift, as discussed in Section 4.2, or possibly introducing a redshift-dependence in the mass probability density function; or, it could indicate that the nature of the SOBBHs is different from what assumed in this work, for example, they could have highly eccentric orbits. By the time LISA will perform the SGWB measurement (or constraint), the SGWB amplitude posterior predicted from ground-based observations will probably have narrowed, if not a detection be made by either 2G or (more likely) 3G detectors. Nevertheless, the LISA measurement will provide further insight into the population of inspiralling SOBBH, by probing the population properties at high redshift and with low masses, and by testing the SGWB signal in a different frequency window.

Impact on the Power-Law Sensitivity
The PLS represents the standard tool to estimate the observability of a given power-law SGWB. The PLS is normally defined assuming that the only stochastic component affecting the SGWB measurement is the instrumental noise [30][31][32]. In Figure 12 we present an improved version of the LISA PLS including the confusion noises generated by GBs and by SOBBHs. For the GBs we adopt the analytical template of [28] with all the parameters taken at their reference value; the SOBBH amplitude on the other hand is fixed to the median value evaluated in this analysis, see Section 4.2, and the tilt to 2/3. The GB contribution mainly affects the low-frequency range, while the SOBBH contribution is relevant at higher frequencies: this effect is reflected in the PLS. The inclusion of the GB confusion noise slightly modifies the PLS at low frequencies, while the impact of the SOBBHs is nearly negligible. Note that Figure 12 corresponds to Figure 2 of [10], while Figure 3 in the same reference is relative to a different treatment, meant to account for the effect of the SGWB amplitude uncertainty, evaluated from the GWTC-2 uncertainty on the merger rate at z = 0.

SGWB detection and the SOBBH population parameters
Intuitively, one might expect the constraining power of a measurement of the SGWB on the SOBBH population model to be very limited, regardless of its precision, since it would reduce the dimensionality of the population parameter space at most by one, leading to a highly-degenerate posterior. On the other hand, this can still have an important impact if the degeneracy associated with the SGWB measurement does not align with the correlations in the population parameter posterior associated with the detection of individual events, the misalignment being due to the fact that the population parameters influence the SGWB amplitude differently from how they influence the characteristics of the population of individually resolvable events. Indeed, it has been demonstrated that a SGWB measurement (or even upper limit) by LVK, in combination with resolved merger events, can constrain the redshift evolution of their merger rate [33,34] and possibly their mass distribution [76]. The high precision with which LISA is expected to measure the SOBBH background, as shown in Section 4.3, should render LISA especially suited to this task. In order to illustrate its potential constraining power, in Figure 13 we plot the GWTC-3 population parameters posterior sample as a scatter plot, highlighting the points compatible with a SGWB amplitude within the LISA 1-and 2-σ credible intervals, relative to a detection by LISA of a SGWB with amplitude corresponding to the median predicted SGWB level P 50 (see Figure 10). One can appreciate that the two-dimensional posterior shrinks significantly, depending on the  Figure 13: Impact, on the population parameters posterior inferred from GWTC-3, of the measurement by LISA of a SGWB with amplitude corresponding to the median value P 50 of the amplitude distribution given in Figure 7. The points coloured from blue to yellow represent the GWTC-3 posterior, and the color scale represents the corresponding SGWB amplitude. The points highlighted in yellow (red) represent the parameter values providing SGWB amplitudes within the 1-σ (2-σ) confidence region of the LISA measurement. Left panel: the initial merger rate R 0 versus its tilt κ. Right panel: the tilts (α, β) of the power-law distributions of m 1 and q = m 2 /m 1 respectively. A measurement of the SOBBH SGWB would break the degeneracy coming from constraints based on individual mergers, and the credible intervals would shrink correspondingly, especially for the merger rate parameters. Had we not fixed the high-redshift behaviour of the merger rate, but treated it probabilistically, the improvement with respect to the GWTC-3 constraints would be smaller, but still significant.
combination of population parameters.
In the left panel of Figure 13 we show the local merger rate R 0 versus its low-redshift tilt κ: the GWTC-3 posterior (points coloured in blue to yellow (for increasing SOBBH SGWB amplitude) presents a degeneracy due to the merger rate being best determined around z ≈ 0.2. Since the value of the low-redshift tilt κ has a strong impact on the SGWB amplitude, the latter varies considerably along this degeneracy (colour scale from blue to yellow). Thus, a precise SGWB measurement, as performed by LISA, would break this degeneracy by leading to a posterior, in the (R 0 , κ) parameter plane, almost perpendicular to the one inferred from the detection of individual SOBBH merger events by ground-based observatories.
The posterior distribution of the mass tilts (α, β), shown in the right panel, would also be significantly reduced. 7 Note that this could be further exploited by a measurement of the anisotropic component of the SOBBH background [77], since the relative amplitude of the anisotropic to the isotropic components appears to be correlated with the tilt of the mass distribution and with the maximal allowed mass [76].
The above results are valid within the assumptions of our analysis, in particular, that the merger rate at high redshift is fixed to the SFR as given in Equation (2.3), and that the LISA uncertainty on the SGWB amplitude is inferred from the MC analysis of a simulated 7 Figure 1 highlighted the effect on the SGWB amplitude of (λ, σ), the parameters of the peak in the Power Law + Peak m1 distribution, rather than the one of (α, β), the tilts of the power laws in m1 and q = m2/m1. The analysis performed in this section shows the converse, hence the choice of parameters in the right panel of Figure 13, apparently contradicting Figure 1. However, in Section 3.1 the impact of the population parameters uncertainty on the SGWB amplitude was considered for each parameter individually so that the influence of the posterior degeneracies was disregarded. data set containing exclusively the SOBBH SGWB, the GBs, and the instrumental noise. Allowing for variations in the high-redshift model of the merger rate, and/or performing a more realistic data analysis procedure accounting for the overlap of several categories of LISA sources, would likely reduce the potential of the SGWB measurement to shrink the population parameter posterior. However, these effects are not expected to alter the misalignment of the correlations in the posterior parameter space inferred from the measurements of individual events and from the measurement of the SGWB. Consequently, the latter will anyway retain, to some degree, its constraining power.

Conclusions
We have evaluated the SGWB expected in the LISA frequency band from SOBBHs, incorporating the most recent information on their mass function, spin distribution, and merger rate coming from LVK observations, in particular from the GWTC-3 posterior on the population parameters of the FidLVK model.
The LVK observations only probe the SOBBH population at low redshift, while faint and distant SOBBHs contribute the most to the background signal. In order to properly evaluate the SGWB, we have therefore extended the GWTC-3 power-law merger rate by assuming that it follows the SFR [26], since the low-redshift expansion of the latter is coherent with the GWTC-3 constraints. With the aim of assessing the impact that this assumption has on the SGWB amplitude, we have also added a time delay in the SOBBHs formation and found that (under a simple model for the time delay distribution, and reasonable values for the minimal time delay) this reduces the SGWB amplitude by at most 40%, remaining within the uncertainty inherent to the GWTC-3 posterior. Though the current precision of the model is not sufficient, future ground-based observations of individual merger events, together with a detection of the SGWB by LISA, will allow to constrain the merger rate and possibly time delays in the future.
We have used four methods to estimate the SGWB. The first method is based on analytic considerations and consists of the integration, over the number density of binaries, of their GW emission in the quasi-circular Newtonian approximation, resulting in a power-law SGWB with slope f 2/3 [27]. The analytical approach has been used to evaluate the impact of each population parameter on the amplitude of the SOBBH background, accounting for its marginalized 95% confidence levels from the GWTC-3 posterior. The power-law index κ of the low-redshift expansion of the merger rate is the parameter influencing the most the SGWB amplitude. We have then calculated the relative percentage change induced in the latter by varying the redshift upper cutoff: we found that integrating up to z max = 5 is sufficient to obtain ∼ 1% accuracy in the evaluation of the background amplitude, also well within the uncertainty induced by the GWTC-3 posterior.
The other three methods employed for the SGWB estimation, gradually increasing in complexity and accuracy, rely on synthetic SOBBH populations, which we have constructed following the GWTC-3 FidLVK posterior distribution. The second method consists in replacing the integration of the analytical method with a MC sum over the masses and redshift of the SOBBHs in the synthetic population realisation, averaging over the time-to-coalescence and the sky-position; while in the third method, the MC sum is performed accounting for the time-to-coalescence of individual events and binning them according to their corresponding emission frequency. These methods allow establishing that the impact of the population variance over the SGWB amplitude is of the order of 0.2%, negligible with respect to the effect of the maximal redshift choice, which is in turn smaller than the uncertainty due to the GWTC-3 posterior.
The fourth method incorporates the actual LISA detection process and consists in simulating LISA data-streams containing the waveforms of all the SOBBHs within the simulated population, and iteratively subtracting the loudest ones until only the confusion noise remains [28]. The threshold SNR used to single out the resolvable GW sources is set to ρ 0 = 8, but we find that the saturation threshold, above which the SOBBH signal is less sensitive to the choice of the threshold itself, is situated at ρ 0 ≈ 5.
We have checked that the four methods provide consistent results for the SGWB amplitude: this is indeed the case at frequencies higher than about 2 mHz, this threshold being exclusively due to the computational limitation of our synthetic populations. The SNR threshold choice results in fact in a limited number of resolvable events, so that the SGWB in the LISA frequency band does not deviate from the analytical power law prediction, which is reproduced also by the three methods based on population synthesis. However, if sources with SNR lower than five will be resolvable in the future, thanks to improvements in data analysis techniques, or to archival searches using future ground-based detector observations, it will be necessary to take into account that the shape of the SGWB in the LISA band deviates from the power law behaviour. This clearly stresses the importance of a precise identification of the resolved sources and of their subtraction, which we present in a follow-up paper [29].
The distribution of the SGWB amplitude at the reference frequency of 3 mHz is evaluated using the analytical method, for all points in the GWTC-3 posterior parameter sample. The interquartile range of the distribution is h 2 Ω GW (f = 3 × 10 −3 Hz) ∈ [5.65, 11.5] × 10 −13 . Our findings are in broad agreement with previous evaluations of the SOBBH stochastic signal and appear therefore to be robust with respect to assumptions such as the high-redshift behaviour of the merger rate and the mass distribution.
We have then performed a MC analysis of simulated LISA data to infer the parameters (i.e., amplitude and spectral tilt) of the SOBBH SGWB in the presence of instrumental noise and of the stochastic signal from GBs. We have found that, with four years of data, the template-based reconstruction of the parameters of both signals and of the noise is accurate to the percent level, with all parameters consistent with their injected values at 2-σ. In this simplified setting where no other GW source is present, and the GB background is static, the SOBBH SGWB can therefore be distinguished from the GB one, despite their overlap at low frequency. We have also compared the MC analysis with a Fisher Information Matrix analysis, finding good agreement, and derived the PLS accounting for the SOBBH and GB backgrounds.
The precision with which LISA will measure the amplitude of the SOBBH SGWB goes from at best 1% (at 1-σ), for the amplitude value corresponding to the 95th percentile of its posterior distribution, up to 5% for the fifth percentile. LISA will therefore reduce by one order of magnitude the current uncertainty on the SGWB amplitude predicted from the GWTC-3 population model. The accuracy of this measurement opens interesting perspectives. We have shown that LISA has the potential to break the degeneracy between some population parameters, since the correlations in the posterior parameter space inferred from the measurements of individual events and of the SGWB, are almost orthogonal. In particular, we forecast an important impact on the merger rate parameters, since the SGWB detection by LISA probes the population of inspiralling SOBBHs at high redshift, fully complementary to actual ground-based observations of low-redshift mergers.
Several extensions of our work are possible, tackling some of its underlying assumptions.
First of all, including eccentricity and precession in the waveforms might have an important effect on the SGWB [78,79]. While we have shown the effect of introducing a time delay between the star formation and the BHs mergers, the impact of the metallicity on the BH mass function has been neglected, see e.g., [20,[65][66][67] for recent studies. A further layer of complexity can be added including the possibility of a redshift dependence of the mass function [41]. The LISA error on the SGWB parameters should be forecasted including other types of sources in the data stream, both resolved and of stochastic nature. Extra-galactic neutron star binaries, for example, generate a collective signal that, although lower in amplitude, is similar to the SOBBHs one, and likely not negligible. Extreme mass ratio inspirals [80] also produce a background at mHz frequencies, although its amplitude is currently poorly constrained and its frequency dependence might not follow a simple power-law in the LISA band [51]. Finally, the effect of the SGWB measurement by LISA on the SOBBH population parameters demonstrated in this work should be properly evaluated via a joint analysis of simulated data from LISA and ground-based observatories, possibly 3G detectors which might be operative by the time LISA flies [81]. Such a joint analysis may also reveal deviations from the expected SOBBH SGWB spectrum, which could point towards a different origin for the BBHs (see e.g. [82]).

A Further information on the SOBBH population model
In this Appendix, we provide more detail on the SOBBH population model: we describe the characteristics of the probability distributions inferred from GWTC-3 observations [41], and justify some of our choices for the catalogues generation, in particular regarding the merger rate behaviour with redshift and the maximal time-to-coalescence.

A.1 Redshift-dependent SOBBH rate
As discussed in Section 2.2, the GWTC-3 constraints on the SOBBH merger rate variation with redshift, assumed to be a power law R(z) = R(0)(1+z) κ , are weak at z 0.5. Therefore, in order to produce accurate SGWB estimates, we need an Ansatz that extends the power law assumption towards higher redshift. We require R(z) to follow the redshift profile of the Madau-Fragos SFR [26]: with r, κ > 0 and R 0 ≡ R(z = 0), implying C = 1 + (κ/r) (1 + z peak ) −(κ+r) . Thus, along the evolution of the universe, from high to low redshift, the SOBBH merger rate initially rises as z −r as more stars are available, and eventually decreases as z κ after the peak of stellar formation. Different from previous studies, e.g. Ref. [34], we introduce the extra factor κ/r in the denominator of Equation (2.3) to guarantee that the function peaks precisely at redshift z peak ; otherwise, the actual peak of the function would deviate from the value of the nominal z peak parameter whenever κ/r = 1. Following this notation, the updated best fit values found in [26] are κ = 2.6, r = 3.6, and z peak = 2.04.
In order for the merger rate R(z) of Equation (2.3) to work as a reasonable highredshift extension of the GWTC-3 low-redshift constraints, we require it to reproduce the profile that LVK obtains for the FidLVK fiducial model fitting the GWTC-3 data [41]. In that study, inference is performed on a low-redshift power law R(z) ∝ (1 + z) κ , resulting in 8 κ = 2.7 +1.8 −1.9 and a pivot rate of R 0 = 17.3 +10.3 −6.7 Gpc −3 yrs −1 at z = 0, or alternatively R 0.2 = 28.3 +12.9 −9.0 Gpc −3 yrs −1 at z = 0.2. These constraints are represented by the blue-shaded region in Figure 14. At low redshift, the median value for the spectral index κ coincides with that of the SFR [26]: in order to extend R(z) at high redshift, we can therefore match the LVK posterior values for R 0 , κ with some values for r, z peak . The latter could e.g. be drawn from some prior distribution; for the purposes of this paper (comparing LISA's sensitivity to SOBBH SGWB amplitudes approximately compatible with the GWTC-3 population inference), it is enough to fix r, z peak to the SFR best fit values mentioned above [26]. The resulting, GWTC-3-compatible, high-redshift merger rate is displayed in Figure 14 in orange shading.
The separability of this distribution facilitates population synthesis since the parameters in the different components can be simulated independently (e.g. using inverse transform sampling in the single-parameter cases). 9 The mass density function is usually given in terms of the mass of the heaviest binary, by convention m 1 , and the mass ratio q = m 2 /m 1 ≤ 1: p(m 1 , m2) = π 1 (m 1 |m min , m max , δ min , α, λ peak , µ m , σ m ) × π 2 (q|m 1 , m min , δ min , β q ) , (A.2) 8 All parameter ranges are given as median ± its respective differences with the percentiles 5 and 95, taken from the public population posterior sample of GTWC-3 for the fiducial FidLVK model. 9 The data in GWTC-3 suggest some correlations that would break this separability, such as that between high spin and mass asymmetry. For the level of the analysis in this paper, it is safe to ignore this finding.  where π 1 is a mixture density function, times a low-mass smoothing: Here P [m min ,mmax] is a truncated power-law distribution with negative spectral index −α, normalized within the [m min , m max ] range, G [m min ,mmax] is a similarly-truncated Gaussian density function representing a possible mass pile-up of BBHs before the SN pair-instability gap [84], C 1 is an overall normalization factor (made necessary by the presence of the smoothing function), and S is a smooth cutoff for low masses that interpolates between 0 and 1 in the interval The probability density function for the mass ratio π 2 in Equation (A.2) is π 2 (q|m 1 , m min , δ min , β q ) = C q q βq S(qm 1 |m min , δ min ) , (A. 6) where C q (m 1 , m min , δ min , β q ) is a normalization factor. The fact that C q depends on m 1 in particular makes the distribution in Equation (A.2) non-separable. It can be computed as

A.3 Time-to-coalescence and frequency of emission
Here we discuss the role of the residual time to coalescence for the population synthesis and the SGWB computation. A correct prediction of the SOBBH SGWB in the LISA band implies catalogues complete enough to adequately simulate the signal. On the other hand, the only observational knowledge we have on these sources comes from LVK observations, which probe the population of merging SOBBHs. In Section 2.1 we have shown that, under the hypothesis that the binaries formation, and therefore their coalescence rates, is in a steady state, the binary rate R(z, τ c ) in Equation (2.1) is indeed equivalent to the one of the merging binaries, constrained by LVK observations. This allows us to construct the catalogues and consequently the SGWB estimation based on the LVK GWTC-3 posterior.
The hypothesis that the binary formation is in a steady state implies that we sample the time-to-coalescence of the binaries in the catalogues uniformly in the range τ c,max = 10 4 yrs, much smaller than the typical time over which the SOBBH population is expected to change, O(10 9 ) yrs. However, is this good enough to account for all the binaries emitting in the LISA band for the entire mission duration? In other words, are the simulated catalogues representative enough of the SOBBH population relevant for LISA? In what follows we demonstrate that, while not complete, our catalogues do indeed provide all the information necessary for a good characterisation of the SOBBH SGWB, as far as our choices on the time-to-coalescence are concerned.
The time interval over which we need to integrate the merger rate in order to obtain the appropriate number of observed events is where T obs denotes the total detector observation time, while T maxBand is the maximum, over all the binaries in the universe, of the portion of each binary's lifetime (i.e. of τ c ) which is spent in the detector frequency range. While in the case of LVK T maxBand is less than seconds, so that T tot T obs , LISA probes the SOBBH population at a different stage, when they are still far away from merging. Inserting the minimal LISA frequency 2 · 10 −5 Hz and the minimal mass in the catalogues m min = 2.5 M (see Section 2.4) in the Newtonian relation for circular orbits (here expressed at the detector, so that M z is the redshifted chirp mass) [59] f = 1 8π c,max = 10 4 yrs appears inappropriate by as much as 10 6 orders of magnitude.
In reality, τ (det) c,max = 10 4 yrs is a pertinent choice that, while preserving computational feasibility, still provides all the relevant information for the SGWB evaluation. By cutting the time-to-coalescence sampling at τ  Figure 15. Figure 5 shows the aggregated effect of this suppression in the SGWB (note that this figure is produced setting z max = 1, as explained in Section 2.4): it is clear from this figure that the relevant spectral property of the SGWB signal, i.e. the power-law behaviour in frequency, is still well captured by the signal produced via the simulated catalogues. The bending at low frequency is nonphysical and therefore irrelevant: the SGWB is expected to simply continue with the same power-law behaviour at low frequencies for synthetic populations with much higher τ c,max . Furthermore, in Figure 12, we can see that the GB foreground overcomes the SOBBH SGWB below 2-3 mHz. It is thus unlikely that an increase beyond τ c,max = 10 4 yrs would produce a noticeable effect in any realistic study. Given the growing computational cost of generating (and computing the SGWB of) synthetic populations with larger τ c,max , we conclude that τ (det) c,max = 10 4 yrs is a good compromise for the purposes of this study.