Massive blackhole binaries will be the loudest sources detectable by LISA. These systems are predicted to form during the hierarchical assembly of cosmic structures and coalesce by interacting with the surrounding environment. The hardening phase of their orbit is driven by either stars or gas and encodes distinctive features into the binary black holes that can potentially be reconstructed with gravitationalwave observations. We present a Bayesian framework to assess the likelihood of massive mergers being hardened by either gaseous or stellar interactions. We use stateoftheart astrophysical models tracking the cosmological evolution of massive blackhole binaries and construct a large number of simulated catalogs of sources detectable by LISA. From these we select a representative catalog and run both parameter estimation assuming a realistic LISA response as well model comparison capturing selection effects. Our results suggest that at least within the context of the adopted models future LISA observations can confidently constrain whether stars or gas are responsible for the binary hardening. We stress that accurate astrophysical modeling of the blackhole spins and the inclusion of subdominant emission modes in the adopted signal might be crucial to avoid systematic biases.
['Astrophysics - High Energy Astrophysical Phenomena', 'Astrophysics - Astrophysics of Galaxies', 'General Relativity and Quantum Cosmology']
Stars or gas Constraining the hardening processes of massive blackhole binaries with LISA
{'Stars or gas? Constraining the hardening processes of massive black-hole binaries with LISA': "Alice Spadaro , Riccardo Buscicchio , David Izquierdo-Villalba , 1, 2 3 3  \n1, 2, ∗ 1, 2, 3 1, 2  \nDavide Gerosa , Antoine Klein , and Geraint Pratten  \n1 Dipartimento di Fisica 'G. Occhialini', Universit'a degli Studi di Milano-Bicocca, Piazza della Scienza 3, 20126 Milano, Italy  \n2 INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, 20126 Milano, Italy  \n3 Institute for Gravitational Wave Astronomy & School of Physics and  \nAstronomy, University of Birmingham, Birmingham, B15 2TT, UK  \n(Dated: September 23, 2024)  \nMassive black-hole binaries will be the loudest sources detectable by LISA. These systems are predicted to form during the hierarchical assembly of cosmic structures and coalesce by interacting with the surrounding environment. The hardening phase of their orbit is driven by either stars or gas and encodes distinctive features into the binary black holes that can potentially be reconstructed with gravitational-wave observations. We present a Bayesian framework to assess the likelihood of massive mergers being hardened by either gaseous or stellar interactions. We use state-of-theart astrophysical models tracking the cosmological evolution of massive black-hole binaries and construct a large number of simulated catalogs of sources detectable by LISA. From these, we select a representative catalog and run both parameter estimation assuming a realistic LISA response as well model comparison capturing selection effects. Our results suggest that, at least within the context of the adopted models, future LISA observations can confidently constrain whether stars or gas are responsible for the binary hardening. We stress that accurate astrophysical modeling of the black-hole spins and the inclusion of subdominant emission modes in the adopted signal might be crucial to avoid systematic biases.", 'I. INTRODUCTION': "According to the current ΛCDM model of structure formation, galaxies form hierarchically [1]. In this context, the formation of massive black-hole (BH) binaries is a direct consequences of galaxy mergers [2]. Present evidence for massive BH binaries relies on observational signatures of active galactic nuclei, including offset or double-peaked spectral lines and periodic variations in light curves [3, 4]. While electromagnetic observations lack the capability to directly resolve massive BH binaries below the parsec scale, the recent evidence for a nanoHertz stochastic gravitational-wave (GW) signal in pulsar timing array datasets points to the existence of a large cosmic population of high mass ( > 10 7 M ⊙ ) merging BH binaries [5-8]. The detection of individual mergers of 10 4 -10 7 M ⊙ BHs is a cornerstone in the science case of future space-based GW detectors, notably the Laser Interferometer Space Antenna (LISA) [9].  \nThe astrophysical environments surrounding massive BH binaries are crucial in shaping the binary dynamics and promoting the final coalescence [10]. In particular, mergers can only be explained if dissipation mechanisms of astrophysical nature are at play in the early inspiral, before binaries enters their GW-driven regime. Modeling such processes is an active field of research, as different assumptions lead to different predictions for both properties and rates of merging massive BHs. Broadly speaking, hardening mechanisms can be divided into 'gas hardening'  \nand 'stellar hardening' [11]. The former, which dominates in gas-rich environments, is driven by the interaction of the massive BH binary with a gaseous circumbinary disk [12, 13]. Conversely, the latter takes place primarily in gas-poor environments, where the binary orbital separation shrinks due to three-body interactions with individual stars [14]. Both gas and stellar hardening affect the merger rate [15, 16] and can potentially be reconstructed with LISA observations [17].  \nWith this in mind, we present a Bayesian framework to compare single detections of massive BH binaries from LISA-simulated catalogs against state-of-theart astrophysical simulations. We use the L-Galaxies semi-analytic model [18] applied to the high-resolution Millenium-II dark-matter (DM) merger trees. We specifically investigate whether it will be possible to infer that putative massive BH binaries observed by LISA have evolved in either gas or stellar environments. We construct realistic LISA catalogs, perform parameter estimation including the full LISA response with the Balrog [19-21] code, and implement a Bayesian model comparison between the two hardening mechanisms while considering selection effects [22].  \nThis paper is organized as follows. In Sec. II, we illustrate the adopted astrophysical models. In Sec. III, we present the construction of our LISA catalogs. In Sec. IV, we describe our parameter-estimation and model-selection strategies. In Sec. V, we present and discuss our findings. Finally, in Sec. VI we summarize our results and outline possible future developments. Throughout this paper, we use units where c = 1.", 'A. Galaxy evolution model': 'We use the L-Galaxies semi-analytical model as presented in Refs. [18, 23, 24]. In particular, the model is built on top of merger trees from the MilleniumII [25] cosmological simulation, which follows the evolution of 2160 3 dark matter particles with a mass of 6 . 885 × 10 6 M ⊙ /h within a comoving box of side 100 Mpc /h , later rescaled [26] to match the Planck cosmology [27] (here h = 10 -2 H 0 Mpc s/km and H 0 is the Hubble constant).  \nL-Galaxies associates a fraction of baryonic matter to each newly resolved DM halo in the form of a diffuse, spherical, and quasi-static hot gas atmosphere. As the gas cools down, it settles into a disk [28]; this facilitates episodes of star formation, which in turn results in the assembly of a stellar disk component. Galaxies are also allowed to assemble an overdensity of stars in their central regions (i.e. their bulges). This is triggered by either internal processes, namely non-axisymmetric instabilities that redistribute the stellar matter, or galaxy mergers.', 'B. Black-hole growth': 'Each newly resolved DM halo is initialized with a massive BH seed spanning a mass range from ∼ 10 2 to ∼ 10 5 M ⊙ [24]. The initial dimensionless spin parameter χ is set uniformly in [0 , 0 . 998] [29].  \nThe evolution of masses and spins is influenced by the accretion of surrounding gas into the BHs as well as BH-binary coalescences. In particular, gas accretion is triggered by both galaxy mergers and disk instabilities [3032], and predominantly drives the mass growth [33]. The cold gas available for accretion settles in a reservoir of mass M res around the BH, which is progressively consumed through a sequence of transient accretion disks [23, 34]. As for the BH spin, its value is set by the frequency of accretion events consumed in prograde or retrograde orbits, which in L-Galaxies is linked to the coherence of the bulge kinematics [35, 36]. The contribution of binary coealescences to the BH mass and spin evolution is important at low redshifts ( z < 2), where galaxies have mostly depleted their gas reservoirs through star formation and feedback processes.  \nThe post-merger BH masses and spins are computed using fits to numerical-relativity simulations [37, 38]. During the BH evolution, L-Galaxies tracks the evolution of the spin magnitude χ but not the spin direction, which is however important for LISA observations. In this work, we construct this property for each binary component using the gas fraction of the environment f gas = M res / ( M res + M ), where M is the total mass of the BH binary. If f gas > 0 . 5, coalescence occurs in a gas-rich environment and the two spins are assumed to  \nFIG. 1. Density distribution of the total mass M and redshift z of massive BH binary mergers. Our fiducial population is shown in gray, where the contour encloses 99% of the estimated source density. The red (blue) distribution corresponds to the sub-population of mergers hardened by gas (stars); contours contain 50%, 90%, and 99% of the estimated source density. Thick contours refer to additional realizations to capture the statistical fluctuations in the underlying astrophysical model. Circles denote the 44 sources from TheCatalog used in our analysis, with colors indicating their respective f gas values. Note the two pairs of overlapping sources, c.f. Fig. 4.  \n<!-- image -->  \n/circledot  \nalign to the binary orbital angular momentum via the Bardeen-Petterson effect [39-41]. If instead f gas < 0 . 5, the environment is gas-poor and the BH spins are assumed to be isotropically distributed [33, 42].', 'C. Binary hardening': 'The dynamical pathway of massive BH binaries can be divided into three phases: pairing, hardening, and GW inspiral [2]. Following a galaxy merger, the pairing phase reduces the BH-BH separation from ∼ kpc to ∼ pc through dynamical friction. This is implemented in LGalaxies as in Ref. [43]. As the binary forms, different physical processes contribute to the shrinking of the binary semi-major axis depending on the gas fraction f gas . For f gas > 0 . 5, gas hardening dominates and the decay of the binary orbit is set by accretion torques from the circumbinary disk. In gas-poor environments with f gas < 0 . 5, instead, the binary evolution is driven by capture and ejection of stars.  \nGravitational-radiation reaction becomes dominant at roughly milliparsec scales and drives the system to the merger. For details on the L-Galaxies implementation see Ref. [23].  \nFIG. 2. Distribution of the gas fraction f gas of the environment across the fiducial population. Sources in blue (red) with f gas < 0 . 5 ( f gas > 0 . 5) are classified as mergers occurring in a gas-poor (gas-rich) environment. The two sub-population distributions are jointly normalized to the full population.  \n<!-- image -->', 'D. Fiducial population': "The model described in the previous subsections provides the astrophysical population of BH mergers we refer to as 'fiducial', which is depicted in grey in Fig. 1 and presented in Ref. [44]. This population comprises two distinct sub-populations with peaks at f gas ∼ 0 (stellar hardening) and f gas ∼ 1 (gas hardening) parameter (Fig. 2). The latter dominates the overall merger rate by a factor of ∼ 5. Figure 1 shows the respective mass and redshift distributions for each sub-populations. Systems evolving in gaseous environments merge on average at higher redshifts z ≲ 8 and have lower masses M ∼ 10 4 -10 6 M ⊙ compared to binaries hardened by stellar interactions which instead are more prominent at z ≲ 3 and have M ∼ 10 5 -10 9 M ⊙ . Nonetheless, the two populations overlap significantly in the M -z parameter space, making their distinguishability with LISA non trivial.", 'III. MOCK LISA CATALOGS': "We build mock LISA observation catalogs on top of the L-Galaxies simulations. The latter provides the number density of massive BH binary mergers, d 2 N/ d z d V c , per unit redshift z and comoving volume V c . We proceed as follows:  \n- (i) We compute the total merger rate of massive BH mergers predicted by the fiducial model [45]  \nd N d t obs = ∫ d N d z w ( z ) d z, (1)  \nwhere  \nw ( z ) = 4 π [ d L (1 + z ) ] 2 d z d V c . (2)  \nIn the equations above, d L is the luminosity distance to the source and t obs is the time measured at the detector.  \n- (ii) We draw the number of sources N astro according to a Poisson distribution with mean value λ = T × d N/ d t obs where T is the duration of the observing period. Massive BH binary signals in LISA typically last hours to months and gradually accumulate signal-to-noise ratio (SNR) over their inspiral. We conservatively set T = 5yr in the calculation above, which accommodates for signals merging after the nominal duration T LISA = 4yr [9].\n- (iii) We extract N astro sources from the astrophysical population with relative weights given by w i = w ( z i ). Each source is characterized by the binary component masses m 1 , 2 , mass ratio q = m 2 /m 1 , the aligned dimensionless spins χ 1 , 2 , and the redshift z .\n- (iv) From these, we further select sources with a cutoff frequency f cut > 1 × 10 -4 Hz and q > 0 . 01. We set f cut = 5 f ISCO where f ISCO is the detector-frame GW frequency at innermost stable circular orbit. This is to ensure that the dominant mode is in the LISA sensitivity band and that the adopted waveform model is reliable.\n- (v) We use the IMRPhenomXHM [46] waveform approximant which captures the full coalescence of quasi-circular, non-precessing BH binaries. The harmonics h lm are calibrated to numerical relativity and include the ( l, | m | ) = { (2 , 2) , (2 , 1) , (3 , 3) , (3 , 2) , (4 , 4) } multipoles. The implementation of the LISA response to such GW signal in the Balrog code has been presented in Ref. [20].\n- (vi) For each source, we draw its extrinsic parameters by sampling uniformly the time to merger t m ∈ [0 , 5yr], the initial phase ϕ 0 ∈ [0 , 2 π ], the polarization angle ψ ∈ [ -π/ 2 , π/ 2], the ecliptic sine latitude sin β ∈ [ -1 , 1] , the ecliptic longitude λ ∈ [0 , 2 π ], and the cosine inclination cos ι ∈ [ -1 , 1].\n- (vii) The SNR of each source, given its parameters θ , is evaluated using the three noise-orthogonal TDI channels h = { h k ; k = A,E,T } as follows  \nSNR 2 = ∑ k 〈 h k ( θ ) ∣ ∣ h k ( θ ) 〉 k . (3)  \nWe assume constant and equal-armlength approximation [47]. The inner product is given by  \n〈 a | b 〉 k = 2 ∫ f max f min ˜ a ( f ) ˜ b ∗ ( f ) + ˜ a ∗ ( f ) ˜ b ( f ) S k ( f ) d f, (4)  \n<!-- image -->  \n<!-- image -->  \nGas hardening  \nFIG. 3. Distribution of massive BH binary mergers across the generated catalogs. The left (gray), middle (blue), and right (red) panel shows the full population of BHs, the sub-population of stellar-hardened sources, and the sub-population of gas-hardened sources, respectively. Light histograms show the astrophysical distributions [steps (i) -(iii) in Sec. III] while dark histograms show the distributions of merger events detectable by LISA during an observation time T = 5yr [steps (iv)-(vii) in Sec. III]. The dashed distributions are a variation of the latter where we include higher-order harmonics and remove the condition on f cut [step (iv)]. The vertical dotted line denotes the selected catalog used to perform the statistical analysis of Sec. IV.  \n<!-- image -->  \nwhere ˜ a ( f ) denotes the Fourier transform of the time series a ( t ) and S k ( f ) the noise power spectral density of the k -th TDI channel. For the latter, we use the semi-analytical expression of Ref. [48] which models the superposition of LISA stationary instrumental noise and astrophysical confusion noise from unresolved Galactic binaries. We use a low-frequency cut-off f min = 0 . 1 mHz, as set in the ESA definition study report [9], and integrate up to f max = 1Hz, which is well above the maximum frequency of all generated GW signals. We consider sources detectable if they exceed the threshold SNR = 10 [9].  \nWe iterate the above procedure 1000 times, resulting in a distribution of LISA catalogs. The corresponding distributions for the number of merger events is shown in Fig. 3. Specifically, steps (i) -(iii) generate the light-gray histogram in the left panel, depicting the distribution of astrophysical mergers predicted by the fiducial model. The dark-gray histogram is obtained by applying steps (iv) -(vii) to the astrophysical catalogs and shows the distribution of massive BH merger events detectable by LISA assuming the dominant l = 2, | m | = 2 emission mode. The blue (red) histograms in the middle (right) panel show the subpopulations of mergers that evolved via stellar (gas)-hardening. Overall, the LISA catalogs have about one order of magnitude fewer events compared to the astrophysical catalogs, which is due to the detection criteria described above. The observable catalogs predominantly feature massive BH binary mergers occurring in gas-rich environments (46 on average) compared to systems evolving in gas-poor environments (1 on average). The averaged count of detectable sources mildly increases to 50 when we relax the assumption on the mass ratio, allowing q to be as small as 0.001. These sources, and those at even lower mass ratio, border the extreme mass-ratio inspiral regime and will need to be treated differently; we  \nleave this to future work.  \nThe dashed distributions in all the three panels of Fig. 3 show the detectable events modeled using the additional modes ( l, | m | ) = { (2 , 1) , (3 , 3) , (3 , 2)(4 , 4) } in the waveform approximant and removing the condition on f cut [step (iv)]. We find higher-order modes marginally increase the number of detectable sources (see Ref. [49]). This is more significant for massive binaries with M ≳ 10 8 M ⊙ , where the (2,2) mode frequency falls out of band. As shown in Fig. 3, this preferentially impacts binaries hardened via stellar processes: including higher-order modes in LISA analyses might have important repercussions on the astrophysical interpretation of the data.  \nIn the following, we study in detail one specific realization of our simulated LISA catalogs, which we refer to as TheCatalog . This specific catalog is close to the medians of the distributions in Fig. 3 and can thus be taken as representative. The 45 detectable sources of TheCatalog are shown as colored circles in the right panel of Fig. 4 together with the astrophysical population from which they are extracted. The left panel puts these sources into context by showing the common 'waterfall plot', i.e. the averaged SNR in the mass-redshift space [9]. LISA is expected to detect massive BH binary systems in the mass range of M tot ∼ 10 4 -10 7 M ⊙ and out to z ∼ 7, covering an SNR range from ∼ 10 to ∼ 3000. The 13 sources marked with white triangles in Fig. 4 are undetectable by LISA. Of these, 9 sources with M tot > 10 4 M ⊙ merge after the LISA mission nominal duration, thereby limiting their SNR growth, cf. step (ii) in Sec. III. Conversely, the 4 sources with M tot < 10 4 M ⊙ merge within the mission lifetime T LISA but remain with SNR < 10 nonetheless. These undetectable sources are excluded from the analysis presented in Sec. IV.  \nStellar hardening  \nFIG. 4. Left panel. Contour lines of constant SNR from massive BH binaries ( q = 0 . 5, χ 1 , 2 = 0 . 2) detectable by LISA, adapted from Ref. [9]. The SNR shown on the color scale has been averaged over sky location, polarization and inclination assuming the dominant quadrupole emission mode. Right panel. The gray contour encompasses 99% of the estimated astrophysical BH population simulated with L-Galaxies . Circles indicate the 45 detectable sources of TheCatalog , colored according to their SNR [see step (vii) in Sec. III]. Of note, two pairs of sources overlap in this plot: one, located at z ∼ 4 . 5 with M ∼ 6 × 10 3 M ⊙ , have the same weights w i [i.e. identical source parameters, see step (iii) in Sec. III]; the other pair, located at z ∼ 2 with M ∼ 4 × 10 4 M ⊙ , are characterized by very similar parameter values. White triangles indicate sources with SNRs below the detection threshold. Note that undetectable sources with high masses ( M > 10 4 M ⊙ ) are affected by the cut-off on the time to merger [step (ii) in Sec. III]. The dashed contour lines of constant SNR are the same of the left panel.  \n<!-- image -->  \n/circledot  \n/circledot", 'IV. STATISTICAL INFERENCE': 'We compare a single GW events against distributions of simulated sources using the Bayesian formalism spelled out in Ref. [22]; see also Refs. [50, 51]. Specifically, we quantify the relative degree of consistency between individual detections from TheCatalog and the two subpopulation models of binary hardening: gas (G) and stellar (S).  \nGiven the posterior distribution p ( θ | d, U) from each individual GW event obtained with some uninformative priors p ( θ | U), we compute the Bayes factor between models G and S as follows:  \nB G / S = B G / U B S / U = ∫ p ( θ | d, U) p ( θ | G) p ( θ | U) d θ ∫ p ( θ | d, U) p ( θ | S) p ( θ | U) d θ . (5)  \nThe astrophysical probability densities p ( θ | G) and p ( θ | S) act as new informative priors on the targeted parameters which we reconstruct from the simulated sources using Gaussian kernel density estimates (KDEs) [52].  \nParameter estimation is performed through Balrog on θ = { θ α , θ ζ } where θ α = {M c , δµ, χ 1 , χ 2 , d L } represent the quantities used to perform the astrophysical model selection [Eq. (6)] and θ ζ = { t m , ϕ 0 , ψ, sin β, λ, cos ι } are the additional extrinsic parameters characterizing the  \nmodeled signals. In particular, we use the redshifted chirp mass M c and the dimensionless mass difference δµ = ( m 1 -m 2 ) / ( m 1 + m 2 ) as they are the two mass parameters that enter the Post-Newtonian (PN) evolution at the leading- and next-to-leading order, respectively. We run full Bayesian inference on simulated data using the nested sampling algorithm [53] as implemented in Nessai [54]. All the injections are in zero noise and we choose uniform priors p ( θ | d, U) on each parameter over either its entire definition domain or a range that is sufficiently large to enclose the entire posterior. In particular, we classify sources into three SNR intervals: low (10 -100), moderate (100 -1000), and high ( > 1000). For each of these, we estimate prior limits by calibrating against a fiducial source selected at the lower boundary of that range. Therefore, the terms p ( θ α | U), p ( θ ζ | U), p ( θ ζ | G), p ( θ ζ | S) [see step (vi) of Sec. III] are constant, we can factor them out the integrals of Eq. (5) and evaluate the KDEs exclusively on θ α . This approach significantly improves the computational efficiency. Equation (5) thus simplifies to  \nB G/S = ∫ p ( θ α | d, U) p ( θ α | G)d θ α ∫ p ( θ α | d, U) p ( θ α | S)d θ α . (6)  \nNote that our detection statistics is approximate as GW detectability depends on the data realization (for detailed  \nFIG. 5. Logarithmic Bayes factor between the gas- and stellarhardening hypotheses as a function of the source SNR. Circles represent the 44 sources from TheCatalog used in the population analysis, with colors indicating their respective f gas values. Horizontal dashed lines denote the threshold values of the Bayes factor according to the Jeffrey scale. Black error bars refer to the statistical (but not systematic) fluctuations associated to the underlying astrophysical model.  \n<!-- image -->  \ndiscussions on this point see e.g. Refs. [55, 56]). We verified that all of our posterior samples have SNR > 10 and q > 0 . 01, which ensure consistency with the adopted detection criterion. Furthermore, it is important to note that posterior samples constrained to a parameter space θ α that is not supported by the gas [stellar] probability density, i.e. p ( θ α | G) = 0 [ p ( θ α | S) = 0], result in a formally null [infinite] value for the Bayes factor in Eq. (6). As a concrete example, see the discussion on the spin distributions in Sec.V.  \nCrucially, one needs to account for selection effects. The detectability-conditioned Bayes factor is given by [22]  \nD G / S = P (det | S) P (det | G) B G / S , (7)  \nwhere P (det | S) and P (det | G) are the fraction of sources from the sub-populations S and G, respectively, that can be detected a priori given adopted detection statistic SNR > 10 and q > 0 . 01. Our numerical implementation follows that described in Ref. [22].', 'A. Inference': "We apply the formalism from Sec. IV to the 45 mergers in TheCatalog , which LISA is expected to detect. Our parameter-estimation pipeline successfully recovers 44 sources, with posterior distributions on source parameters θ well confined within physically unbounded priors.  \nHowever, for one source with SNR ∼ 13, the stochastic algorithm fails to converge, leading to poor-quality posteriors. Therefore, we choose to exclude this source from the population analysis. Additionally, sources with SNRs in the range 10 -100 frequently exhibit non-linear correlations and multimodalities on θ ζ ; however both are sufficiently decoupled from the reconstructed θ α 's, which are the target parameters of our astrophysical model selection. In the following, we quote parameter estimates at 90% confidence interval.  \nAmong the target parameters, M c is the most precisely measured. Specifically, we constrain it with a relative precision ∆ M c / M c of ∼ 10 -6 -10 -4 for sources with moderate and high SNRs. For low SNRs, the relative precision is ∼ 10 -4 -10 -3 . The spin components χ 1 and χ 2 are measured with a relative precision of 0 . 006% -1% at high and moderate SNRs, and 1% to 65% at low SNRs.  \nSix systems in the low -SNR range exhibit mild biases on δµ due to the broad posteriors (∆ δµ/δµ ∼ 6% -160% at 90% credible level) and projection effects associated to non-linear parameter correlations. Of these, three systems have posterior distributions that include the injected values only in their 99 . 9% confidence inteval. For the other three sources, the true injected value lies within the 98 . 5% credible interval. Finally, we emphasize that LISA will be able to measure the source luminosity distance with a relative precision ∆ d L /d L of ∼ 0 . 3% -2% for high SNRs, ∼ 1% -45% for moderate SNRs, and ∼ 10% -115% for low SNRs.  \nUncertainties related to sky localization and host galaxy identification, which are crucial for multimessenger detections of gas-rich mergers [57], will be addressed in future work.", 'B. Model selection': "We use the posteriors described above to evaluate the detection-conditioned Bayes factors D G / S of Eq. (7). Results are presented in Fig. 5, where sources are ordered by their SNR and colored according to their true astrophysical sub-population. For all the sources we find 'decisive' evidence [58] in favor of the gas (stellar) sub-population model, with log 10 D G / S > 2 (log 10 D G / S < -2). Additionally, no significant correlation is observed between D G / S and the SNR of each source. This is largely due to the highly precise measurement of the target parameters, as described above. Even the wider posterior distributions are well localized within the sub-population parameter space, so increasing the SNR does not lead to improved discrimination between the two astrophysical models. We observe that most of the Bayes factors fall within the range 10 ≲ | log 10 D G / S | ≲ 25, with three exceptions that exhibit lower values. Among theese, two sources with SNR ∼ 2000 and 7 ≲ | log 10 D G / S | ≲ 10 evolved trough a gas-hardening phase ( f gas ∼ 0 . 72). Conversely, the lowest Bayes factor ( | log 10 D G / S | = 5 . 63) is associated with  \nFIG. 6. Marginal distribution of the aligned spin component of the primary BH χ 1 (left panel) and the secondary BH χ 2 (right panel) for mergers in gaseous (red) and stellar (blue) environments. Thick histograms indicate our fiducial model; thin histograms indicate additional realizations. Blue dots denote the median value of source s 1 from TheCatalog as well as the two additional sources s 2 , 3 , c.f. Sec. V. The error bars indicate the 90% confidence intervals. The two sub-population distributions are jointly normalized to the full population.  \n<!-- image -->  \nthe only source in TheCatalog that hardened through stellar interactions ( f gas ∼ 0 . 14).  \nWe further test the analysis pipeline by generating ten additional realizations of our astropyhysical model to evaluate the impact of simulation uncertainties on the results. The error bars in Fig. 5 indicate the resulting variability associated with each log-Bayes factor. Even when accounting for these statistical fluctuations, model G is confidently favored for all the G sources, thus validating the previous results. For the one S source, model S remains favored, though to a lesser extent. Crucially, these error bars capture the statistical fluctuations in the underlying astrophysical model and not the systematic effects due to the many assumptions entering the model itself.", 'C. Source parameters': 'Based on these results, we now investigate how the astrophysical properties of massive BH binaries might influence the distinguishability between models G and S. In Fig. 1, we show the 44 analyzed systems from TheCatalog together with the astrophysical sub-population distributions marginalized on the total mass and redshift. Notably, we observe that the gas-hardened sources are predominantly supported by the gas-hardening subpopulation, with most falling within (outside) the 90% contour level of the gas (stellar) distribution. Conversely, the stellar-hardened source, hereafter referred to as s 1 , falls outside the 90% contour level of the stellar-hardening sub-population but is located close to the peak (i.e.  \nwithin the 50% contour level) of the gas-hardening subpopulation. This suggest that additional parameters may be impacting its log 10 D G / S value, favoring the stellarhardening model.  \nTo investigate further, we analyze the one-dimensional distribution of the aligned component spins χ 1 and χ 2 , see Fig. 6. For the source s 1 considered so far, the primary component spin χ 1 = 0 . 49 +0 . 08 -0 . 08 lies below the 92 nd and 1 . 5 th percentiles of the stellar- and gas-hardening distributions, respectively, thus placing it confidently in the former. Additionally, the secondary component spin χ 2 = 0 . 39 +0 . 13 -0 . 13 is well below the 90 th percentile of the stellar-hardening distribution and lies within only the 0 . 3 rd percentile of the gas-hardening distribution. We find that the recovered δµ parameter does not provide significant support for model S, with the full posterior information falling below the 14 th percentile of the stellarhardening distribution but within the 50 th percentile of the gas-hardening sub-population.', 'D. Importance of the BH spins': 'From the discussion above, we conclude that, at least for the astrophysical models considered here, the component spins are the crucial parameters for inferring the hardening processes of massive BH binaries. We stress that the employed models describe a simplified scenario. Massive BH binaries evolving in gas-rich environments experience accretion torques that might align their component spins with the orbital angular momentum of the binary. The models used here assume that the spin alignment pro-  \ncess occurs rapidly, typically within the time-scale of the merger, thus leading to high and positive spin configuration at merger [42] ( χ 1 , 2 ∼ 0 . 98 on average, as indicated by the red distributions in Fig. 6). In contrast, binaries that evolve in gas-poor environments do not experience significant gas accretion and likely enter the GW dominated phase with roughly isotropic spin orientations, as reflected by the blue distributions in Fig. 6. Nevertheless, as pointed out in the literature, the overall picture is more complex due to the internal properties of the disk. Indeed, the effectiveness of disk alignment may be influenced by the gas temperature [40], accretion rate [59], as well as the disk viscosity [60] and its initial angle of misalignment, which can lead to critical disk-breaking configurations [61]. Accurate astrophysical models accounting for these properties are crucial for addressing such scenarios in future study.  \nFinally, we extend the analysis from TheCatalog by including two additional sources, s 2 and s 3 , selected from the LISA stellar-hardening distribution shown in Fig. 3. In order to provide a representative sample of LISA observations, we specifically select sources with SNR in the range 100 -1000, which is most densely populated interval for detectable massive BH binaries in LISA (see Fig. 4). This region is populated by sources with M ∼ 10 4 -10 6 M ⊙ at redshift z ∼ 2 -4. Additionally, we select the two sources such that their primary spin components fall in the tails of the stellar-hardening sub-population from Fig. 6. This allow us to further investigate the role of spins in the model comparison.  \nFirst, as shown in Fig. 6, the primary component spin χ 1 = -0 . 8 +0 . 1 -0 . 1 for source s 2 is entirely constrained within a region where p ( χ 1 | G) = 0, thus making the KDE evaluation unnecessary, as pointed out in Sec. IV. Therefore, we can definitely favor model S over model G (log 10 D G / S ≪-2) without even considering the influence of the other binary properties. For source s 3 , the analysis indicates a preference for the stellar-hardening scenario with log 10 D G / S = -4 . 65. This case is particularly challenging because the gas-hardening model has also strongly support over the source posterior. Specifically, s 3 has a mass M ∼ 7 × 10 4 M ⊙ and is located at a redshift z ∼ 1 . 6, positioning it at the peak of the M -z distribution for the gas-hardening sub-population (see Fig. 1). Additionally, it has a measured primary spin component of χ 1 = 0 . 98 +0 . 01 -0 . 03 . This value is notably high and exceeds the 99 . 5 th percentile of the stellar-hardening distribution and the 7 th percentile of the gas-hardening distribution. As illustrated in Fig. 6, s 3 is situated in a high-density region for gas-hardened sources but in a relative sparse region for stellar-hardened sources. Furthermore, the recovered value of δµ falls below the 9 th (38 th ) percentile for the stellar (gas) sub-population, providing no further support for model S with respect model G. On the other hand, the secondary component spin strongly favors model S. In fact, the recovered value χ 2 = 0 . 79 +0 . 04 -0 . 01 falls below the 95 th (1 . 2 nd ) percentile for the stellar- (gas-) hardening  \nsub-population.  \nThis further emphasizes that, at least in the astrophysical setup considered here, component spins play a crucial role in constraining the hardening mechanisms responsible for the final coalescence of massive BH binaries, even in ambiguous but rare cases (less than 1% of the stellar-hardening population sources detectable by LISA). However, simulation-induced fluctuations reduce the logBayes factor of source s 3 to as low as 0.35, indicating only a weak preference for model S. Incorporating spin precession into the population analysis may be crucial to enhance model discrimination and gain further insights into the galactic environments where massive BH binary mergers occur. This is left to future work.', 'VI. CONCLUSIONS': "In this work, we studied the astrophysical origin of merging massive BH binaries detectable by LISA, with a specific focus on the hardening mechanisms that drive the binary toward the final coalescence.  \nCombining results from large-scale cosmological simulations and state-of-the-art models for massive BH binary evolution, we constructed a large set of realistic mock catalogs for LISA. Our analysis showed that most detectable sources are expected to evolve in gas-rich environments. We then performed Bayesian inference and model selection on a representative set of LISA detections and quantifies the statistical evidence in favor of either the gas-hardening or stellar-hardening channel.  \nOur results demonstrated that component spins might play a crucial role in distinguishing between these two astrophysical sub-populations and LISA will be able to provide clear insight given the expected precision on measured parameters. This is a rather strong conclusion, which might however depend on the specific astrophysical model adopted in this paper. Reversing the argument, our findings imply that accurate modeling of the BH spins, both magnitudes and directions, is crucial to avoid major systematics in LISA's astrophysical inference.  \nIndividual source inference with a richer waveform phenomenology -including precession effects, smaller massratios, and potential eccentricity- as well as more complete astrophysical simulations will yield a more realistic estimate of the LISA potential to uncover the astrophysics of massive BHs. For instance, in this paper we find that including higher-order modes in the adopted GW signal preferentially impacts sources that evolved through specific astrophysical formation pathways. While we leave further explorations to future work, we stress that the methodology presented here is general and can be readily applied to more sophisticated models.", 'ACKNOWLEDGMENTS': "Data supporting this paper (including all mock LISA catalogs and posterior samples) are publicly released at Ref. [62]. We thank Matthew Mould, Monica Colpi, Alberto Sesana, Christopher Moore, Nathan Steinle, and Jonathan Gair for discussions. A.S. and D.G. are supported by ERC Starting Grant No. 945155-GWmining, Cariplo Foundation Grant No. 2021-0555, MUR PRIN Grant No. 2022-Z9X4XS. A.S., R.B. D.I.-V., and D.G. are supported by MUR Grant 'Progetto Dipartimenti di Eccellenza 2023-2027' (BiCoQ), and the ICSC National Research Centre funded by NextGenerationEU. R.B. is supported by Italian Space Agency Grant 'Phase  \n- [1] S. D. M. White and M. J. Rees, Mon. Not. R. Astron. Soc. 183 , 341 (1978).\n- [2] M. C. Begelman, R. D. Blandford, and M. J. Rees, Nature 287 , 307 (1980).\n- [3] A. De Rosa et al. , New Astron. Rev. 86 , 101525 (2019), arXiv:2001.06293 [astro-ph.GA].\n- [4] T. Bogdanovi'c, M. C. Miller, and L. Blecha, Living Rev. Relativ. 25 , 3 (2022), arXiv:2109.03262 [astro-ph.HE].\n- [5] G. Agazie et al. , Astrophys. J. Lett. 951 , L8 (2023), arXiv:2306.16213 [astro-ph.HE].\n- [6] J. Antoniadis et al. , Astron. Astrophys. 678 , A50 (2023), arXiv:2306.16214 [astro-ph.HE].\n- [7] D. J. Reardon et al. , Astrophys. J. Lett. 951 , L6 (2023), arXiv:2306.16215 [astro-ph.HE].\n- [8] H. Xu et al. , Research in Astronomy and Astrophysics 23 , 075024 (2023), arXiv:2306.16216 [astro-ph.HE].\n- [9] M. Colpi et al. , (2024), arXiv:2402.07571 [astro-ph.CO].\n- [10] P. Amaro-Seoane et al. , Living Rev. Relativ. 26 , 2 (2023), arXiv:2203.06016 [gr-qc].\n- [11] M. Dotti, A. Sesana, and R. Decarli, Advances in Astronomy 2012 , 940568 (2012), arXiv:1111.0664 [astro-ph.CO].\n- [12] A. I. MacFadyen and M. Milosavljevi'c, Astrophys. J. 672 , 83 (2008), arXiv:astro-ph/0607467 [astro-ph].\n- [13] J. Cuadra, P. J. Armitage, R. D. Alexander, and M. C. Begelman, Mon. Not. R. Astron. Soc. 393 , 1423 (2009), arXiv:0809.0311 [astro-ph].\n- [14] G. D. Quinlan, New Astron. 1 , 35 (1996), arXiv:astroph/9601092 [astro-ph].\n- [15] A. Sesana, Class. Quantum Grav. 30 , 244009 (2013), arXiv:1307.4086 [astro-ph.CO].\n- [16] A. Klein, E. Barausse, A. Sesana, A. Petiteau, E. Berti, S. Babak, J. Gair, S. Aoudia, I. Hinder, F. Ohme, and B. Wardell, Phys. Rev. D 93 , 024003 (2016), arXiv:1511.05581 [gr-qc].\n- [17] A. Toubiana, K. W. K. Wong, S. Babak, E. Barausse, E. Berti, J. R. Gair, S. Marsat, and S. R. Taylor, Phys. Rev. D 104 , 083027 (2021), arXiv:2106.13819 [gr-qc].\n- [18] B. M. B. Henriques, S. D. M. White, P. A. Thomas, R. Angulo, Q. Guo, G. Lemson, V. Springel, and R. Overzier, Mon. Not. R. Astron. Soc. 451 , 2663 (2015), arXiv:1410.0365 [astro-ph.GA].\n- [19] G. Pratten, P. Schmidt, H. Middleton, and A. Vecchio, Phys. Rev. D 108 , 124045 (2023), arXiv:2307.13026 [grqc].  \nA activities for the LISA mission' No. 2017-29-H.0. G.P. is supported by the Royal Society University Research Fellowship URF \\ R1 \\ 221500 and RF \\ ERE \\ 221015. D.G. is supported by MSCA Fellowships No. 101064542StochRewind and No. 101149270-ProtoBH. Computational work was performed at University of Birmingham BlueBEAR High Performance Computing facility, at CINECA with allocations through EuroHPC Benchmark access call grant EUHPC-B03-24, INFN, and Bicocca, and at Google Cloud through award No. GCP19980904.  \nSoftware : We acknowledge usage of of the following Python [63] packages for modeling, analysis, postprocessing, and production of results throughout: Nessai [54], matplotlib [64], numpy [65], scipy [66].  \n- [20] G. Pratten, A. Klein, C. J. Moore, H. Middleton, N. Steinle, P. Schmidt, and A. Vecchio, Phys. Rev. D 107 , 123026 (2023), arXiv:2212.02572 [gr-qc].\n- [21] A. Spadaro, R. Buscicchio, D. Vetrugno, A. Klein, and et al., Phys. Rev. D 108 , 123029 (2023), arXiv:2306.03923 [gr-qc].\n- [22] M. Mould, D. Gerosa, M. Dall'Amico, and M. Mapelli, Mon. Not. R. Astron. Soc. 525 , 3986 (2023), arXiv:2305.18539 [astro-ph.HE].\n- [23] D. Izquierdo-Villalba, A. Sesana, S. Bonoli, and M. Colpi, Mon. Not. R. Astron. Soc. 509 , 3488 (2022), arXiv:2108.11671 [astro-ph.GA].\n- [24] D. Spinoso, S. Bonoli, R. Valiante, R. Schneider, and D. Izquierdo-Villalba, Mon. Not. R. Astron. Soc. 518 , 4672 (2023), arXiv:2203.13846 [astro-ph.GA].\n- [25] M. Boylan-Kolchin, V. Springel, S. D. M. White, A. Jenkins, and G. Lemson, Mon. Not. R. Astron. Soc. 398 , 1150 (2009), arXiv:0903.3041 [astro-ph.CO].\n- [26] R. E. Angulo and S. D. M. White, Mon. Not. R. Astron. Soc. 405 , 143 (2010), arXiv:0912.4277 [astro-ph.CO].\n- [27] Planck Collaboration, P. A. R. Ade, N. Aghanim, C. Armitage-Caplan, M. Arnaud, M. Ashdown, F. AtrioBarandela, J. Aumont, C. Baccigalupi, A. J. Banday, and et al., Astron. Astrophys. 571 , A16 (2014), arXiv:1303.5076 [astro-ph.CO].\n- [28] R. S. Sutherland and M. A. Dopita, Astrophys. J. Supp. S. 88 , 253 (1993).\n- [29] K. S. Thorne, Astrophys. J. 191 , 507 (1974).\n- [30] E. Berti and M. Volonteri, Astrophys. J. 684 , 822 (2008), arXiv:0802.0025 [astro-ph].\n- [31] N. Fanidakis, C. M. Baugh, A. J. Benson, R. G. Bower, S. Cole, C. Done, and C. S. Frenk, Mon. Not. R. Astron. Soc. 410 , 53 (2011), arXiv:0911.1128 [astro-ph.CO].\n- [32] E. Barausse, Mon. Not. R. Astron. Soc. 423 , 2533 (2012), arXiv:1201.5888 [astro-ph.CO].\n- [33] Y. Dubois, M. Volonteri, and J. Silk, Mon. Not. R. Astron. Soc. 440 , 1590 (2014), arXiv:1304.4583 [astro-ph.CO].\n- [34] D. Izquierdo-Villalba, S. Bonoli, M. Dotti, A. Sesana, Y. Rosas-Guevara, and D. Spinoso, Mon. Not. R. Astron. Soc. 495 , 4681 (2020), arXiv:2001.10548 [astro-ph.GA].\n- [35] M. Dotti, M. Colpi, S. Pallini, A. Perego, and M. Volonteri, Astrophys. J. 762 , 68 (2013), arXiv:1211.4871 [astroph.CO].\n- [36] A. Sesana, E. Barausse, M. Dotti, and E. M. Rossi, Astro-\n- phys. J. 794 , 104 (2014), arXiv:1402.7088 [astro-ph.CO].\n- [37] W. Tichy and P. Marronetti, Phys. Rev. D 78 , 081501 (2008), arXiv:0807.2985 [gr-qc].\n- [38] E. Barausse and L. Rezzolla, Astrophys. J. Lett. 704 , L40 (2009), arXiv:0904.2577 [gr-qc].\n- [39] J. M. Bardeen and J. A. Petterson, Astrophys. J. Lett. 195 , L65 (1975).\n- [40] M. Dotti, M. Volonteri, A. Perego, M. Colpi, M. Ruszkowski, and F. Haardt, Mon. Not. R. Astron. Soc. 402 , 682 (2010), arXiv:0910.5729 [astro-ph.HE].\n- [41] D. Gerosa, G. Rosotti, and R. Barbieri, Mon. Not. R. Astron. Soc. 496 , 3060 (2020), arXiv:2004.02894 [astroph.GA].\n- [42] T. Bogdanovi'c, C. S. Reynolds, and M. C. Miller, Astrophys. J. Lett. 661 , L147 (2007), arXiv:astro-ph/0703054 [astro-ph].\n- [43] J. Binney and S. Tremaine, Galactic Dynamics: Second Edition (2008).\n- [44] D. Izquierdo-Villalba, A. Sesana, M. Colpi, D. Spinoso, M. Bonetti, S. Bonoli, and R. Valiante, Astron. Astrophys. 686 , A183 (2024), arXiv:2401.10983 [astro-ph.GA].\n- [45] E. Barausse, I. Dvorkin, M. Tremmel, M. Volonteri, and M. Bonetti, Astrophys. J. 904 , 16 (2020), arXiv:2006.03065 [astro-ph.GA].\n- [46] C. Garc'ıa-Quir'os, M. Colleoni, S. Husa, H. Estell'es, G. Pratten, A. Ramos-Buades, M. Mateu-Lucena, and R. Jaume, Phys. Rev. D 102 , 064002 (2020), arXiv:2001.10914 [gr-qc].\n- [47] M. Tinto and S. V. Dhurandhar, Living Rev. Relativ. 24 , 1 (2021), arXiv:gr-qc/0409034 [gr-qc].\n- [48] S. Babak, J. Gair, A. Sesana, E. Barausse, C. F. Sopuerta, C. P. L. Berry, E. Berti, P. Amaro-Seoane, A. Petiteau, and A. Klein, Phys. Rev. D 95 , 103012 (2017), arXiv:1703.09722 [gr-qc].\n- [49] N. Steinle, H. Middleton, C. J. Moore, S. Chen, A. Klein, G. Pratten, R. Buscicchio, E. Finch, and A. Vecchio, Mon. Not. R. Astron. Soc. 525 , 2851 (2023), arXiv:2305.05955 [astro-ph.HE].\n- [50] S. Vitale, D. Gerosa, W. M. Farr, and S. R. Taylor, in Handbook of Gravitational Wave Astronomy (2022) p. 45.\n- [51] I. Mandel, W. M. Farr, and J. R. Gair, Mon. Not. R. Astron. Soc. 486 , 1086 (2019), arXiv:1809.02063 [physics.data-an].\n- [52] D. W. Scott, Multivariate Density Estimation: Theory, Practice, and Visualization (2015).\n- [53] J. Skilling, Bayesian Analysis 1 , 833 (2006).\n- [54] M. J. Williams, J. Veitch, and C. Messenger, Phys. Rev. D 103 , 103006 (2021), arXiv:2102.11056 [gr-qc].\n- [55] R. Essick and M. Fishbach, Astrophys. J. 962 , 169 (2024), arXiv:2310.02017 [gr-qc].\n- [56] M. Mould, C. J. Moore, and D. Gerosa, Phys. Rev. D 109 , 063013 (2024), arXiv:2311.12117 [gr-qc].\n- [57] A. Mangiagli, C. Caprini, M. Volonteri, S. Marsat, S. Vergani, N. Tamanini, and H. Inchausp'e, Phys. Rev. D 106 , 103017 (2022), arXiv:2207.10678 [astro-ph.HE].\n- [58] R. E. Kass and A. E. Raftery, J. Am. Stat. Assoc. 90 , 773 (1995).\n- [59] D. Gerosa, B. Veronesi, G. Lodato, and G. Rosotti, Mon. Not. R. Astron. Soc. 451 , 3941 (2015), arXiv:1503.06807 [astro-ph.GA].\n- [60] G. Lodato and D. Gerosa, Mon. Not. R. Astron. Soc. 429 , L30 (2013), arXiv:1211.0284 [astro-ph.CO].\n- [61] C. Nixon, A. King, and D. Price, Mon. Not. R. Astron. Soc. 434 , 1946 (2013), arXiv:1307.0010 [astro-ph.HE].\n- [62] A. Spadaro and R. Buscicchio, github.com/RiccardoBuscicchio/lisa-mbhb-cats-andsamps,\n- doi.org/10.5281/zenodo.13787675 (2024).\n- [63] G. V. Rossum and F. L. Drake, Python 3 Reference Manual (CreateSpace, Scotts Valley, 2009) iSBN: 9781441412691.\n- [64] J. D. Hunter, Computing in Science & Engineering 9 , 90 (2007).\n- [65] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, and et al., Nature 585 , 357 (2020), arXiv:2006.10256 [cs.MS].\n- [66] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, and et al., Nature Methods 17 , 261 (2020), arXiv:1907.10121 [cs.MS]."}