The role of scatter and satellites in shaping the large-scale clustering of X-ray AGN as a function of host galaxy stellar mass

The co-evolution between central supermassive black holes (BH), their host galaxies, and dark matter halos is still a matter of intense debate. Present theoretical models suffer from large uncertainties and degeneracies, for example, between the fraction of accreting sources and their characteristic accretion rate. In recent work we showed that Active Galactic Nuclei (AGN) clustering represents a powerful tool to break degeneracies when analysed in terms of mean BH mass, and that AGN bias at fixed stellar mass is largely independent of most of the input parameters, such as the AGN duty cycle and the mean scaling between BH mass and host galaxy stellar mass. In this paper we take advantage of our improved semi-empirical methodology and recent clustering data derived from large AGN samples at $z \sim 1.2$, demonstrate that the AGN bias as a function of host galaxy stellar mass is a crucial diagnostic of the BH--galaxy connection, and is highly dependent on the scatter around the BH mass--galaxy mass scaling relation and on the relative fraction of satellite and central active BHs. Current data at $z \sim 1.2$ favour relatively high values of AGN in satellites, pointing to a major role of disc instabilities in triggering AGN, unless a high minimum host halo mass is assumed. The data are not decisive on the magnitude/covariance of the BH-galaxy scatter at $z \sim 1.2$ and intermediate host masses $M_\mathrm{star} \lesssim 10^{11} \,\mathrm{M}_\odot$. However, future surveys like Euclid/LSST will be pivotal in shedding light on the BH--galaxy co-evolution.


INTRODUCTION
The presence of supermassive black holes at the centers of virtually every massive galaxy is an accepted paradigm. The masses of these central black holes appear to correlate with the properties of their host galaxies (e.g. Ferrarese & Ford 2005;Kormendy & Ho 2013;Graham & Scott 2015;Reines & Volonteri 2015;Shankar et al. 2019Shankar et al. , 2020 with tentative evidence for a link even with their host dark matter haloes (e.g. Ferrarese 2002;Bandara et al. 2009). The very existence of such correlations, which are observed to hold even at higher redshifts (e.g. Shankar et al. 2009b;Cisternas et al. 2011b;Shen et al. 2015;Suh et al. 2020), point to a degree of co-evolution between the black holes and their hosts. Thus, unveiling the shape, dispersion and evolution of these correlations represents a crucial task in present-day extra-galactic astronomy to acquire a full picture of galaxy formation and evolution.
Despite decades of observational and theoretical work aimed at deciphering the nature of the black hole-galaxy scaling relations, the causal link between the two remains still unsolved. A strong release of energy/momentum from an accreting central supermassive black hole (BH) shining as an Active Galactic Nucleus (AGN) naturally predicts a tight and constant correlation with velocity dispersion (e.g. ★ E-mail: akke.viitanen@helsinki.fi Silk & Rees 1998;Granato et al. 2004), which has been recognized as one of the most fundamental property linked to black hole mass (e.g. Bernardi 2007;Shankar et al. 2019;Marsden et al. 2020). Mergerdriven models of black hole growth would instead predict a tighter correlation with the host galaxy stellar mass (e.g. Jahnke & Macciò 2011;Hirschmann et al. 2010).
Given the large number of input assumptions in traditional ab initio cosmological models, a more phenomenological complementary approach has been introduced in the last decade to constrain the properties of black holes in a cosmological context. The latter method relies on identifying the overall populations of active and inactive black holes at a given redshift on statistical grounds via, e.g. continuity equation techniques (e.g. Marconi et al. 2004;Shankar et al. 2004Shankar et al. , 2013Aversa et al. 2015) and/or semi-empirical mock catalogues tuned to reproduce stellar mass functions, AGN X-ray luminosity functions and/or AGN clustering properties (e.g. Georgakakis et al. 2019;Aird & Coil 2021;Allevato et al. 2021) Further constraints on the co-evolution scenario are indeed provided by studying the spatial clustering of AGN, especially X-ray selected (e.g. Coil et al. 2009;Krumpe et al. 2010;Allevato et al. 2011;Krumpe et al. 2012;Mountrichas & Georgakakis 2012;Mendez et al. 2016;Powell et al. 2018;Viitanen et al. 2019;Allevato et al. 2019;Powell et al. 2020). Clustering provides an independent way of connecting accreting BHs to their large-scale environments via the link with their underlying dark matter halo population. Several clustering studies have argued that AGN environment could play a significant role in triggering the AGN phase, highlighting the importance of host dark matter halos in nuclear activity (e.g. Hickox et al. 2009;Allevato et al. 2011;Fanidakis et al. 2012;Gatti et al. 2016).
However useful, interpretation of X-ray AGN clustering results have suffered from not properly being able to take into account the host galaxy properties. Indeed what is currently debated is to which degree AGN clustering can be understood solely in terms of the underlying host galaxy properties (e.g. color, stellar mass star , starformation rate), and AGN selection effects (e.g. Hickox et al. 2009;Mendez et al. 2016;Mountrichas et al. 2019;Viitanen et al. 2019;Powell et al. 2020). Thus of special interest is the host galaxy stellar mass star due to its intimate connection with the underlying dark matter halo (e.g. Moster et al. 2010;Leauthaud et al. 2012;Behroozi et al. 2019).
In terms of the AGN-hosting dark matter halo connection, despite their extreme flexibility, all phenomenological AGN models suffer from serious degeneracies that cannot be easily broken. Shankar et al. (2020) have shown that the large-scale clustering of AGN, e.g. the bias, as a function of BH mass represents a powerful framework to constrain the normalization and slope of the correlation between BH mass and host galaxy stellar mass. Recently, Allevato et al. (2021) have found by using AGN mock catalogs built on observationally derived galaxy -BH scaling relations, AGN duty cycles (i.e.the probability for a galaxy hosting an active black hole), and Eddington ratio distributions, that the AGN large-scale bias is a crucial diagnostic to break degeneracies in the input AGN model. Other groups have, instead, built the correlation between black hole mass and host halo mass by passing the scaling relation with host galaxy stellar mass or velocity dispersion and the AGN duty cycle (e.g. Georgakakis et al. 2019;Aird & Coil 2021). Besides not allowing for any quantitative assessment of the underlying fundamental black hole-galaxy scaling relations, limiting the efficacy of these models in shedding light on the processes controlling the co-evolution of BHs and their hosts, the latter approach is still not immune to degeneracies between e.g. duty cycles and Eddington ratio distributions.
In particular, Allevato et al. (2021) have shown at ∼ 0.1 that the clustering at fixed stellar mass is largely independent of the input duty cycle, Eddington ratio distribution, and the starhalo relation. This implies that the clustering of normal galaxies matches the AGN clustering at a given stellar mass, provided the AGN hosts are a random subsample of the underlying galaxy population of the same stellar mass. However, as we will show in this work, the AGN largescale bias as a function of stellar mass is mainly set by the scatternot so much by the shape -of the BH mass -stellar mass relation and by the relative fraction of AGN in central and satellite dark matter halos. The role of the scatter in the BH -host mass relation on the AGN clustering has been investigated in different studies at high ∼ 4 (White et al. 2008;Wyithe & Loeb 2009;Shankar et al. 2010b;Bonoli et al. 2010). In particular, White et al. (2008) have shown that the very high bias measured for SDSS quasars at ∼ 4 by Shen et al. (2007) is in favour of a small scatter in the BH -galaxy scaling relation. In fact, the key idea is that any scatter in the input scaling relation increases the contribution from lower mass and less biased halos, thus decreasing the observed AGN large-scale bias.
However, these studies focus on high redshifts, where the population of supermassive BHs that powers the quasars is growing rapidly and inhabits the rarest, most massive halos. Moreover, they do not investigate the effect of the scatter on the AGN clustering in different stellar mass bins, also bypassing additional parameters that might affect the AGN bias normalization.
Instead, in this work we study the concurring effect of the scatter around the BH mass-galaxy mass relation and of the relative fraction of satellite and central active BHs on the AGN large-scale clustering, focusing our attention on the AGN bias as a function of host galaxy stellar mass at intermediate redshift ∼ 1. We thus expand on Allevato et al. (2021) by following their methodology to create realistic mock catalogs of galaxies and active BHs and (i) we include on top of a regular Gaussian scatter a positive covariant correlation between the BH and the host galaxy stellar mass at fixed halo mass; (ii) we focus on redshift ∼1.2, i.e. close to the peak of the AGN activity (e.g. Shankar et al. 2009a;Madau & Dickinson 2014); and (iii) we allow the satellite AGN fraction to vary. Recently Farahi et al. (2019) applied a similar covariant method approach in order to reveal an anti-correlation between hot and cold baryonic mass in galaxy clusters. In the context of BH-galaxy co-evolution, it would expected that a degree of correlation exists between the properties of the central BH and the galaxy, especially in any scenario involving an AGN feedback.
The AGN bias as a function of host galaxy stellar mass is thus a crucial diagnostic of the black hole-galaxy connection which, together with the bias as a function of black hole mass , can provide tight constraints on the co-evolution of black holes and galaxies. It is the ideal time to set out the methodology behind the modelling of AGN bias as a function of stellar mass based on the large amount of current and future data. Recent multi-wavelength surveys such as XMM-COSMOS, Chandra COSMOS legacy survey, AEGIS have in fact made it possible to study the AGN environment in terms of host galaxy properties. In particular, several AGN clustering studies have specifically considered AGN clustering as a function of the host galaxy stellar mass (e.g. Georgakakis et al. 2014;Mountrichas et al. 2019;Viitanen et al. 2019;Allevato et al. 2019). Moreover, in the next years the synergy among J-PAS, eROSITA, 4MOST, LSST, Euclid, and JWST will allow us to measure the AGN clustering with high statistics, and to derive host galaxy stellar mass estimates of moderate-to-high luminosity AGN up to ∼ 2 to measure the large-scale bias as a function of host galaxy properties of millions of objects at different .
The paper is organized as follows. In Section 2 we will detail the methodology for estimating the covariant scatter and the semiempirical model to build our AGN mock catalogues. In Section 3 we outline the results, while in sections 4 and 5 we present the discussion and conclusions of this work.

METHODOLOGY
In this work we will investigate the large-scale clustering properties of X-ray selected AGN, and moreover the imprint of covariant scatter between the central black hole mass BH and the host galaxy stellar mass star . For this purpose, we create mock catalogs of active and non-active galaxies by using semi-empirical models (SEMs) based on large N-body simulations, and measure the large-scale bias as a function of stellar mass ( star ) as well as the projected two-point correlation function p ( p ) of mock AGNs. The full description of numerical routines to create mock catalogs of galaxies and their BHs using SEMs is given in Allevato et al. (2021) at ∼ 0.1.
We extend this work to higher redshift = 1.22, and investigate the role of covariant scatter in assigning star and BH to mock AGNs. We only describe the important steps in the generation of the AGN mock catalogs and we refer the reader to Allevato et al. (2021) for the details.
We note that our results are specific to the redshift of interest, start-  ing from the input DM halo catalog and input scaling relations which evolve with redshift (e.g. Grylls et al. 2019), and upon assigning large-scale bias to DM halos = ( halo , ), with bias increasing with at a given halo (Sheth et al. 2001; van den Bosch 2002).

Input Scaling Relations
For building mock catalogs of AGNs, we start from the dark matter (DM) halo catalog extracted from the MultiDark Planck 2 (MDPL2) simulation (Riebe et al. 2013;Klypin et al. 2016) at redshift = 1.22, close to the redshift of recent X-ray AGN clustering measurements (e.g. Viitanen et al. 2019;Allevato et al. 2019). The MDPL2 simulation is a cosmological DM only N-body simulation, with a box size of 1 Gpc/ℎ, 3 840 3 DM particles, and a mass resolution of 1.51 × 10 9 M /ℎ. The cosmology used in the simulation is flat ΛCDM with ℎ = 0.6777, 8 = 0.8228, and Ω = 0.307115. We use the ROCKSTAR (Behroozi et al. 2013) DM halo catalog.The catalogs contain both central/parent halos and satellite halos with unstripped mass at infall, and the 3-dimensional positions for estimating the projected two-point correlation function p .
In the semi-empirical approach, each halo is assigned a galaxy and central BH (active or not) following the most up-to-date empirical relations. This approach has the benefit of avoiding the need to model physical AGN processes, while reproducing the galaxy/BH statistical properties, such as stellar mass function and luminosity function. First, halos and subhalos are populated with stellar masses star according to a star = star ( halo , ) (we drop the explicit dependence and focus on = 1.22 from now on) relation using the Moster et al. (2010) formulae, and the recent updated parameters of Grylls et al. (2019).
Galaxies are then populated with BHs according to an input BH ( star ) relation as derived in equation 6 of Shankar et al. (2016), inclusive of the scatter. This relation is significantly lower in normalization and steeper than relations inferred for BHs with dynamically measured masses (e.g. Kormendy & Ho 2013). In this work we will also explore how different input scaling relations (including relations derived for early-type (ET) and late-type (LT) galaxies) affect AGN clustering measurements as a function of the host galaxy star . It is worth noticing that we assume that all BHs share the same probability of residing in ET and LT galaxies, as well as star-forming and quenched galaxies. In Sec. 3 we will show the effect of relaxing this assumption.
As we investigate in this work, the scatter in the input star ( halo ) and BH ( star ) relations prove to be important parameters in dictating the clustering properties of mock AGN. For this purpose, we will use two different prescriptions for the scatter, which we label This work, and Reference hereafter. The two methods differ in the magnitude of the scatter in the input star ( halo ) and BH ( star ) relations, and in the case of covariant scatter we also explicitly define the covariance of the scatter between star and BH for a given DM halo mass halo , as detailed further in Section 2.2. To each BH is then assigned an Eddington ratio (and then an Xray luminosity) following a probability distribution described by a Schechter function as suggested in Bongiorno   The local galaxy data and the best-fit scaling relations used to derive the covariance matrix of the scatter for 'This work'. Left and right panels show the velocity dispersion against black hole mass, and K-band luminosity, respectively. The dashed black lines correspond to 16% and 84% quantiles of the best-fit relations based on 20 000 Monte Carlo samples. Aird et al. (2017); Georgakakis et al. (2017). The Schechter function is defined by two free parameters, the knee of the function ★ , and the index which describes the power-law behaviour at Eddington ratio values below the knee.
Each BH can be active or not according to an observationally deduced duty cycle (Schulze et al. 2015), i.e. a probability for each BH of being active. Following Schulze et al. (2015), for the highend of the black hole mass function (above 10 6 solar masses), the AGN duty cycle decreases with the BH mass at small redshifts and becomes almost constant at > 1. Figure 3 (left panel) shows the X-ray luminosity functions of our mock AGNs compared to data as derived in Miyaji et al. (2015) at = 1.098. The lines mark the contribution of AGN in different bins of host galaxy star . Moreover, the AGN host galaxy stellar mass function in X-ray luminosity bins is shown in the right panel.
We find that a Schechter Eddington ratio distribution with input parameters log ★ = −0.6 and = +1.2 (both parameters are dimensionless) we are able to reproduce the AGN X-ray luminosity function from Miyaji et al. (2015). We have verified that our results are robust to small changes in the input Eddington ratio parameters. Further dependencies on the input Eddington ratio are discussed in more detail in Allevato et al. (2021).
Finally, we follow Allevato et al. (2019) and assign an obscuration value H to each BH based on AGN luminosity following Ueda et al. (2014).

Covariant Scatter
In this work we implement two different methods to include the scatter in the input scaling relations. In the Reference case, for a given halo, we first assign the stellar mass following the star − halo inclusive of the scatter. Then, we use the scattered galaxy stellar mass to assign the black hole mass given BH ( star ) log star = log star (log halo ) + log star (1) where log star is a Gaussian random variable with zero mean and a scatter of 0.11 (Grylls et al. 2019) and log BH is a Gaussian random variable with star dependent intrinsic scatter of 0.32−0.1× [log( star / ) − 12] as suggested by Shankar et al. (2016). However, in the second approach (covariant scatter method, This Work), for a given halo , we assign the stellar mass star and black hole mass BH using the mean scaling relation for both properties. The multivariate normal scatter in star and BH is then assigned according to a given covariance matrix: where N is the 2-dimensional Gaussian distribution with covariance matrix C (with covariances ) Since data are sparse in terms of the covariance between star and BH for a given DM halo mass halo for AGNs, we use a local galaxy sample to derive the full covariance matrix by using the galaxy velocity dispersion and -band luminosities as proxies of halo and star , respectively. We use the most updated local galaxy sample from de Nicola et al.
(2019) which includes 84 galaxies with black hole masses measured either from stellar dynamics, gas dynamics, or astrophysical masers, velocity dispersion and K-band luminosities. In this work our interest lies in the massive end of the DM halo mass function, and given that the AGN duty cycle is not well constrained for BHs with masses , we de-select potential lowmass systems by including only galaxies with velocity dispersion > 10 2.1 km/s, leaving 71 galaxies in our final sample, which we show in the left and right panels of Fig. 2.
Our approach is outlined as follows. We start by re-evaluating the scaling scatter and the full covariance matrix for the scaling relations of the local galaxy sample of de Nicola et al. (2019). We use the K-band luminosity K and velocity dispersion as proxies of star and halo , in log-linear scale. We then estimate the BH -and K -scaling relations, from which we calculate the scatter and the covariance of BH and star for a given halo .
For fitting the scaling relations, we use the Bayesian routine (Kelly 2007), assuming: where K = log K , BH = log BH , and = log . The observational errors, prefixed by Δ in the equation, are assumed to be drawn from a multivariate gaussian distribution with a negligible covariance between the errors (Saglia et al. 2016;van den Bosch 2016). The intrinsic scatters are assumed to be drawn from normal distributions with zero means and 2 variance. We assume uninformative priors and do not mask the data, and use 20 000 Monte Carlo (MC) samples. We then derive the scatter (denoted by a prefixed ) from the log-linear relation for each MC sample = 1 − 20 000 and for each observation = 1 − 71 from the log-linear relations: Finally, we estimate the covariance from the scatter for each of the sample as follows (Farahi et al. 2019): where is the number of observations For our local sample of 71 galaxies, we report 16%, 50%, and 84% percentiles of the 20 000 MC samples in Table 1, as well as Fig. 2. In detail, for the scaling relations we find BH = −4.19, BH = 5.45, K = 2.06 and K = 3.83. For the covariance matrix, we find the 50% percentiles (log K , log K ) = 0.21, (log BH , log BH ) = 0.19 and (log K , log BH ) = 0.07, which as discussed previously, are also the values we utilize as the baseline for the full covariance between star and BH in the covariant scatter method. However, the covariant scatter derived here may be considered as an upper limit of the scatter in the input scaling relations. In building our AGN mock catalogs, we fix an upper limit to the scatter in the BH -stellar mass relation by requiring that the high-end tail ( star ≥ 10 11.5 ) of the galaxy stellar mass function at = 1.22 is not overproduced. This corresponds to a reduction of 2 log star by 50% (see Fig. 3).

Parameter Q
We also investigate the effect of increasing the relative probability of satellite BHs of being active. Formally, we define = Here refers to the numbers of total/central/satellite BHs. In principle, the parameter may not strictly be a constant but a function of stellar mass and/or environment. However, in the spirit of keeping a flexible and transparent semi-empirical approach, we will for simplicity keep a constant Allevato et al. 2021).
It is important to stress that varying the value of only modifies the relative contributions of central and satellite BHs, but it does not alter the total duty cycle and thus it does not spoil the match to the AGN luminosity function. The parameter is related to the fraction of AGN in satellite halos AGN sat by the relation where AGN sat is defined via summation over all DM halos : and BH sat = sat sat + cen (14) is the total fraction of (active and non active) BHs in satellites. In this work we will study how the parameter affects AGN clustering estimates, and in particular the AGN ( star ).

Mock AGN samples
For a closer comparison between our semi-empirical models and observations, we match our AGN simulated mocks to the observed samples in terms of X and stellar mass distributions. We also add to the intrinsic scatter in the stellar mass -halo mass relation, a measurement error of 0.20 dex to better reproduce the observed scatter in the stellar mass estimates in the COSMOS AGN samples ( In the following, we also compare with the AGN semi-empirical model of Aird & Coil (2021), for which we apply an X-ray luminosity cut at X > 10 42 . Lastly, we compare our predictions with the clustering estimates of SDSS quasars at ∼ 1.4 measured by Richardson et al. (2012), and for this test we impose a limit on luminosity of X > 10 44 erg s −1 and on Hydrogen column density of H < 10 22 cm −2 to select only optical, nominally Type I AGN (e.g., Ricci et al. 2017). We list some of the main features of each of our reference AGN mock subsample in Table 2.

Two-point correlation function and large-scale bias
For quantifying the clustering of mock AGNs, we use two complementary approaches. Firstly, we compute the projected two-point correlation function using the three-dimensional positions of the AGN host DM halos. Secondly, we calculate the AGN large-scale bias by averaging over the biased parent DM halo population.
We estimate the projected two-point correlation function of mock AGNs by using the 3-d positions of the hosting DM halos, following Davis & Peebles (1983): where ( p , ) is the 2-dimensional Cartesian two-point correlation function (e.g. Peebles 1980). Physical distances p and correspond to the perpendicular and parallel (defined with respect to a far-away observer) separations, defined for each AGN pair , ( p , ) separately. Using periodic boundary conditions in the simulation box with volume = 10 9 (ℎ −1 Mpc) 3 , ( p , ) is estimated using where is the sum of all unique mock AGN pairs , within cylindrical volume element Δ defined as the volume enclosed by log p ± Δ log p /2 and ± Δ /2 weighted by the AGN duty cycle , and is the expected number of randomly distributed mock AGN pairs within the same volume (e.g. Alonso 2012; Sinha & Garrison 2020). For estimating the correlation function, we use p = 0.1 − 100 ℎ −1 Mpc, and = 0 − 40 ( max = 40 ℎ −1 Mpc) with bin sizes Δ log( p / ℎ −1 Mpc) = 0.25 and Δ = 1 ℎ −1 Mpc, respectively. We use the publicly available C F code (Sinha & Garrison 2020).
For the large-scale bias ( star ) we estimate the bias for each DM halo separately. Halos are labeled as central or satellite. For each central halo, we assign a value of the large-scale bias according to van den Bosch (2002); Sheth et al. (2001). For each satellite DM halo we assign a large-scale bias value based on the mass of its parent halo, as each satellite traces the dense environment of the parent halo 1 . Then, we follow the formalism of Shankar et al. (2020); Allevato et al. (2021) to derive the bias of mock AGN and normal galaxies as a function of galaxy stellar mass, by using the parameter. The bias of mock objects with stellar mass in the range star ± d star /2 1 We further verify that the DM halo bias of Sheth et al. (2001) When = 1 and sat = cen , then on average central and satellite BHs share equal probabilities of being active.

RESULTS
In this section we present our main results in terms of projected twopoint correlation function p ( p ), AGN large-scale bias as a function of galaxy stellar mass ( star ) and mean AGN halo occupation distribution (HOD) ( halo ) , i.e. the average number of AGN as a function of the halo mass. In what follows, we label as "This work" all results based on the covariant scatter method, while we use the "Reference" label for all models characterised by independent Gaussian scatters in all the input scaling relations.
The left panel of Fig. 4 shows the projected correlation function p ( p ) of mock AGNs created using the covariant scatter method (solid lines) and the reference model (dashed red lines), for different host galaxy stellar mass bins. The right panel of Fig. 4 plots instead the relative bias between the two different approaches. On large scales, the difference is in the highest stellar mass bin 11 < log ( star / ) < 12, where the clustering strength of the covariant scatter is a factor of ∼ 0.7 lower compared to the reference case. At smaller stellar masses, the two cases differ by a more moderate factor of 0.9. Meanwhile, at lower scales p , and for the lower stellar mass bins, we find an opposite trend, where the clustering strength of the covariant scatter case is higher by up to a factor of ∼ 1.2.
In Fig. 5 we show the bias in bins of stellar mass of mock AGNs for the covariant scatter case (solid blue line), and the reference case (dashed red line), and for X-ray selected AGN at ∼ 1 Mountrichas et al. 2019;Viitanen et al. 2019). For these comparisons, we match the mock AGN sample in terms of X-ray luminosity to the Chandra-COSMOS sample , as explained in Sec. 2.4. In the covariant scatter case, the predicted AGN bias is almost constant as a function of stellar mass, as the covariance by design tends to naturally increase the scatter in stellar mass at fixed DM halo mass. In fact, further decreasing the intrinsic scatter in stellar mass to 40% of the value derived in the covariant scatter method, leads to an AGN large-scale bias that is closer to the reference case without covariance (dashed red line). Moreover, we note that the typical star measurement error in COSMOS (see Sec. 2) added on top of the intrinsic scatter bridges some of the gap in the difference of the biases at star 10 11 , where the inclusion of this additional scatter lowers the bias at a given stellar mass principally in the reference case.
For completeness, we also show the large-scale bias as a function of stellar mass of VIPERS (Marulli et al. 2013) and COSMOS (Jullo et al. 2012) for all galaxies at a similar redshift, which is well reproduced by our mock galaxies (i.e., bias not weighted by the duty cycle and without cuts in the X-ray luminosity, black dashed line in Fig. 5). It is worth noticing that, in the reference case, normal and active galaxies with the same luminosity cut follow the same bias -stellar mass relation. A larger degree of discrepancy in the bias in bins of stellar mass between AGN and normal galaxies is instead predicted when assuming a covariant scatter, especially at log ( star / ) > 11 (see Sec. 4 for more discussion). The predicted bias as a function of stellar mass of mock AGNs is consistent with the bias estimates of X-ray selected AGN characterized by a large uncertainty, irrespective of the used method for the scatter. However, the bias of mock AGNs appears systematically below by at least ∼ 20% the measurements derived in Allevato et al. (2019) for Chandra-COSMOS AGN.
There could be different, concurrent causes that could explain the offset between COSMOS and the mocks. The COSMOS field is renown to be characterized by overdensities at < 1 that might increase the AGN large-scale bias by up to 26% (e.g. Gilli et al. 2009;Mendez et al. 2016), and suffers from cosmic variance due to the small volume. In particular, to estimate the effect of such a finite volume in our AGN mocks, we provide an estimate of the cosmic variance by sampling the full MDPL2 simulation box with side length of 1 Gpc/h. We then extract a total of 36 unique subvolumes with area 27 × 27 (Mpc/ℎ) 2 approximately 1.4 deg 2 at = 1, comparable to the COSMOS field (Scoville et al. 2007), and length 1 Gpc/h, and calculate the large-scale bias for each of the 36 sub-boxes individually. We then use the standard deviation of the large-scale bias from the sub-volumes as an estimate of the effect of the cosmic variance. The filled areas in Fig. 5 shows the 1 error of the bias, upon selecting a COSMOS-like volume from the full MDPL2 simulation box, which is of the order of ∼0.1 dex at almost all stellar masses. Our results would thus imply that, by itself, cosmic variance cannot account for the full offset between our mocks and the Allevato et al. (2019) bias measurements. We will discuss additional factors that could contribute to the discrepancy between models and data in what follows below. We will specifically investigate the dependence of the AGN large-scale bias (and related AGN HOD) on other input parameters, namely the BH -galaxy scaling relation and the AGN satellite fraction (i.e., the number of AGN in massive galaxy groups/clusters) as parametrized by .

Dependence on Input Scaling Relations
In the creation of AGN mock catalogs we have assigned black hole masses to galaxies according to their host galaxy stellar mass following Shankar et al. (2016). Here we explore the impact of varying this input relation. We adopt the relation derived by Savorgnan & Graham (2016) for a sample of galaxies with dynamically measured BH masses, as presented in Eq. 3 of Shankar et al. (2019), which predicts higher BH for a given star . These two chosen relations broadly bracket the existing systematic uncertainties in the black hole mass-stellar mass relation of dynamically measured BHs in the local Universe, with other scaling relations broadly falling in between them (e.g., Terrazas et al. 2016;Davis et al. 2019;Sahu et al. 2019).
We find that, in the reference case, the AGN large-scale bias as a function of stellar mass is independent of the particular input BH mass -stellar mass relation. When assuming the covariant scatter method instead, the AGN bias predicted by the Savorgnan & Graham (2016) Figure 6. Large-scale bias as a function of host galaxy stellar mass and its dependence on the input BH ( star ) relation. The shaded area around Shankar et al. (2016) relation shows the bias weighted additionally by the probability of a galaxy of being quenched as a function of the DM halo mass and redshift (assuming = 2.5) as given by Zanisi et al. (2021). The limits correspond to AGN host galaxies consisting solely of quenched (upper limit), and non-quenched (lower limit) galaxies. Otherwise the symbols and linestyles are the same as in Fig. 5.
In our model we also assume that AGNs reside in a mixed population of host galaxies, i.e. active BHs share the same probability of being active in star forming and quenched galaxies. We can relax this assumption and explore the effect on the AGN bias as a function of stellar mass of having all AGNs in quenched galaxies. To this purpose, we assign to each AGN a probability of being in a quenched host galaxy ℎ as a function of halo mass and redshift following Rodríguez-Puebla et al. (2015) and Zanisi et al. (2021). In detail, we use quench = 1/[ 0 + ( ( , ) × 10 12 / halo )], where ( , ) = 0 + (1 + ) , with halo given in units of solar masses, 0 = 0.68, ∼ 1, and ∼ 2.5 (see Zanisi et al. 2021 for the details). We note that variations in = 1 − 4 affect the large-scale bias at most at the ∼ 5% level as compared to = 2.5. We show the result in Fig. 6 as a shaded area, where the upper (lower) limit corresponds to the AGN bias in bins of stellar mass when assuming all AGNs in quenched (non-quenched) galaxies. In particular, assuming that all mock AGNs reside in quenched host galaxies, which live preferentially in more massive and biased halos, produces a large-scale bias as a function of stellar mass slightly higher (1.05 times) than when considering AGN in a mixed population.

Dependence on Q
In this section we study the effect of varying the relative number of satellite AGNs parameterized by . Fig. 7 shows that higher values of tend to increase the normalization of the bias at all considered stellar masses, both in This Work and in the Reference case. This can be understood as an increase in the fraction of AGNs in highly biased massive systems, which are likely to host satellite AGNs in the first place. Increasing the probability of a satellite BHs being active thus increases the bias. We show this effect for moderate values of = 1, 2 (i.e. satellite BHs are 1−2 times more likely to be active than centrals), but also more extreme values of = 4, 6. Bias estimates for X-ray selected COSMOS AGN seem to favour values of ∼ 4, while XMM-XXL AGN are more in agreement with models with = 1. We can convert any value of to an actual fraction of satellite active BHs above a given luminosity threshold. Values of = 1, 2, 4, 6 correspond to AGN satellite fractions AGN sat = 0.10, 0.18, 0.31, 0.40, respectively (see Table 2). In Fig. 9 we compare our predicted satellite fractions as a function of parent host DM halo mass for different values, with models of Gatti et al. (2016) and Aird & Coil (2021) for AGNs with X > 10 42 erg/ s and data from Pentericci et al. (2013) and Martini et al. (2009). Measurements of the AGN satellite fraction in groups and clusters (Pentericci et al. 2013;Martini et al. 2009) suggest small values of ∼ 1 − 2, as also assumed in the semi-empirical model presented in Aird & Coil (2021). It is worth noticing that independently of the parameter , sat is by construction increasing as a function of the parent DM halo mass (see Fig. 9). In particular, higher increases the AGN satellite fraction in less massive parent halos. Our results show that the AGN satellite fraction is thus an input parameter that controls the normalization of the AGN large-scale bias as a function of stellar mass. Additional measurements of sat in groups and clusters at this redshift would help in putting independent constraints on .

AGN HOD and the 1-halo term
Finally, we estimate the projected two-point correlation function (2pcf) and corresponding mean halo occupation distribution (HOD) of mock AGNs. In particular, we measured the mock AGN 2pcf over the full scale range p = 0.1 − 100 ℎ −1 Mpc including the small scales within the 1-halo (< 1ℎ −1 Mpc) which is due to the correlation of AGN within the same DM halo. This regime is especially sensitive to the fraction of AGNs in satellite galaxies of galaxy groups and clusters. The mock AGN HOD is calculated as the average number of mock AGNs in bins of DM halo mass, separating the contribution of mock satellite and central AGNs. Fig. 8 shows the 2pcf of mock AGNs with X > 10 44 erg/s and H < 10 22 cm −2 for = 2, 4 and 6, compared with the 2pcf estimated Richardson et al. (2012). The comparison with data suggests high values of (= 4), with a corresponding satellite AGN fraction of AGN sat 0.23 (see Table 2). The corresponding mean HODs of mock AGNs and SDSS quasars are shown if Fig. 10 (right panel). While the SDSS quasar HOD is characterized by a significant central occupation only for DM halos with halo > 10 13 with a steep increase of the satellite quasar HOD as a function of the hosting halo mass, the occupation of our central mock AGNs is significant in halos with mass down to halo ∼ 10 12.3 ℎ −1 . Moreover, the satellite fraction derived for SDSS quasars in Richardson et al. (2012) is ∼ 1000 times smaller than the one measured in our AGN mock catalog. This comparison shows that AGN small scale clustering can be modelled by HODs characterized by different satellite AGN fractions and minimum mass of the hosting halos.
A similar result is also shown in the left panel, for mock AGN with X > 10 42 erg/s, compared to the HOD derived in Richardson et al. (2013) for X-ray selected COSMOS AGN at ∼ 1.2. Also for moderate luminosity AGN, they found an upper limit of AGN sat = 0.1 and a significant central occupation for central halos with mass halo > 10 12 ℎ −1 , i.e. almost 3 times larger than what derived from our mock AGNs. For comparison, we also show the AGN HOD   Aird & Coil (2021) model. The semi-analytic model prediction of Gatti et al. (2016) based AGN triggering through disc instabilities (DI) is shown as a dark dotted line. Martini et al. (2009);Pentericci et al. (2013) show measurements of X-ray AGNs in galaxy groups/clusters. parent halo masses. Their model also finds a significant AGN HOD for hosting halos with mass down to halo > 10 11.5 ℎ −1 .

DISCUSSION
In this work we create realistic mock catalogs of active BHs and galaxies at ∼ 1.2 based on semi-empirical models, to study the impact on the AGN large-scale bias as a function of the host galaxy stellar mass and focus in particular on the impact of (1) including a covariant scatter in the BH mass -stellar mass relation at fixed halo mass; (2) changing the satellite AGN fraction. We demonstrate that the AGN bias at fixed stellar mass is mainly dependent on the scatter in the input scaling relations and on the relative fraction of satellite and central active BHs.

Covariant Scatter
The scatter in the input BH mass -stellar mass relation strongly affects the AGN large-scale bias as a function of stellar mass. In particular we found that the larger the scatter, the flatter is the resulting bias -stellar mass relation. This behaviour is expected as the end product of a covariant scatter is to generate a larger scatter in the input In detail, the covariant scatter model produces an AGN bias versus star up to twice as small at star ∼ 10 11.5 than the Reference case. A model with covariant scatter, i.e. with a (positive) correlation between BH mass and galaxy mass at fixed halo mass, would then imply a bias as a function of stellar mass substantially different between AGN and the overall population of galaxies at the same stellar mass, at least for M star > 10 11 , where unfortunately AGN bias estimates are not yet available and/or still have large associated uncertainties (Mountrichas et al. 2019;Viitanen et al. 2019;Allevato et al. 2019). At star < 10 11 , instead, the difference in the largescale bias of normal and active galaxies at the same stellar mass is not significant, as also found in Mendez et al. (2016) and Powell et al. (2018).
Our predicted AGN large-scale bias is roughly consistent with the measurements by (Mountrichas et al. 2019;Viitanen et al. 2019), which are however characterized by large uncertainties, but it always falls short by a systematic ∼ 20% in reproducing the measured bias of CCL AGNs , at least for models with = 1. The COSMOS field (Scoville et al. 2007) is known to be characterized by large overdensities at < 1 which could increase the large-scale bias up to ∼ 20 − 50% (Gilli et al. 2009;Mendez et al. 2016;Viitanen et al. 2019). In fact, Viitanen et al. (2019) also estimate the AGN bias when removing from the sample AGN in galaxy groups and measured a bias ∼0.4 dex smaller in the lower stellar mass bin. In addition, the COSMOS field is affected by cosmic variance due to the small volume of the survey which introduces an additional 1 error of ∼0.1 dex to the bias measurements (see Sec.3 for more details). All these effects make it difficult to constrain models with and without covariant scatter. In the near future, precise clustering measurements of AGN that extend up to star > 10 11 will allow us to put constraints on these models.
The study of the degeneracy between AGN clustering strength and scatter in the BH-host scaling relations has already been emphasized by a number of groups at different redshifts (e.g., Shen et al. 2007;White et al. 2008;Wyithe & Loeb 2009;Bonoli et al. 2010;Shankar et al. 2010b,a). One of the main findings emphasised by these studies was the need for a small scatter in the input scaling relation to boost the AGN clustering signal to match some of the data. Our current results point to similar trends, whilst emphasizing the vital role of additional parameters, such as the relative fraction of AGN satellites, in amplifying the clustering signal.

Satellite AGN fraction
In this work we have highlighted the pivotal role of the parameter, i.e. the ratio between active satellite and central BHs, in regulating the normalization of the bias-stellar mass relation (e.g., Figure 7). In addition, constraining the parameter can shed light on AGN triggering models. For example, a high relative fraction of satellite AGN would be in line with the evidence put forward several times that secular processes and bar instabilities, and not only mergers (e.g., Hopkins et al. 2008), are efficient in triggering moderate-to-luminous AGN (e.g. Georgakakis et al. 2009;Allevato et al. 2011;Cisternas et al. 2011a;Schawinski et al. 2011;Rosario et al. 2011).
The value of is yet not well constrained at > 1, whilst most studies suggest that 2 at ≥ 1. Allevato et al. (2019) suggested that the ∼ 1 large-scale bias of CCL Type 2 AGN as a function of the host galaxy stellar mass can be reproduced assuming ∼ 2, which corresponds to a satellite AGN fraction sat AGN ∼ 0.15, and applying a cut in the AGN host halo mass of halo 10 12 ℎ −1 . Estimates of the AGN content in groups of galaxies in the two GOODS fields at ∼ 1 (Pentericci et al. 2013) favour ∼ 1 − 2. At smaller redshifts, an AGN satellite fraction sat AGN ∼ 0.18 has been suggested by Leauthaud et al. (2015) for COSMOS AGN at < 1. Allevato et al. (2012) performed direct measurements of the HOD for COSMOS AGN at < 0.2 based on the mass function of galaxy groups hosting AGN, and found that the duty cycle of satellite AGN is comparable or slightly larger than that of central AGN, i.e. ≤ 2. Georgakakis et al. (2019) used semi-empirical models to populate DM halos with AGN assuming no distinction between centrals and satellites, i.e.
= 1 and found a fair agreement, within the error budget, between the observationally derived AGN HOD (e.g. Allevato et al. 2012;Miyaji et al. 2011;Shen et al. 2013) and their AGN mock predictions. Shankar et al. (2020) also showed that the large-scale bias as a function of BH mass of X-ray and optically selected AGN at ∼ 0.25 can be reproduced by assuming ≤ 2. Recently, Allevato et al. (2021) showed that the Q parameter strongly affects the large-scale AGN bias as a function of stellar mass, BH mass and luminosity, with ∼ 1 − 2 more favoured by the data at ∼ 0.1.
We also investigated how the parameter affects the AGN HOD and in turn the 1-halo term. The HOD approach has been used by different authors to interpret AGN and quasar clustering measurements (e.g. Allevato et al. 2011;Miyaji et al. 2011;Richardson et al. 2012Richardson et al. , 2013Kayo & Oguri 2012). In particular, Richardson et al. (2012) and Kayo & Oguri (2012) performed clustering measurements of optically luminous quasars for both small (< 1 ℎ −1 Mpc) and large physical scales and performed HOD modelling to infer the relation between quasars and their host dark matter halos (see Fig. 10). Richardson et al. (2012) found at ∼ 1.4 a small fraction of luminous SDSS quasars to be in satellites DM halos, with AGN sat ∼ 7.4 × 10 −4 and that the central (satellite) occupation becomes significant only at masses above halo ∼ 10 13 ℎ −1 (∼ 10 14 ℎ −1 ). Kayo & Oguri (2012) also modelled the clustering of SDSS quasars at ∼ 1.4 and reported a satellite fraction ∼100 times higher than Richardson et al. (2012), by using a different HOD parameterization. We also showed that by increasing the input parameter we were able to reproduce the full 2pcf of luminous SDSS quasars with resulting HODs very different from those derived by Richardson et al. (2012), further highlighting the degeneracies in HOD modelling. We note that such degeneracies have already been reported from HOD mod-

Alternatives to high Q values
The comparison of our mock predictions with the 2pcf at all scales of SDSS quasars suggest high values of (∼4), which are in contrast with previous measurements from HOD modelling in Richardson et al. (2012) and Kayo & Oguri (2012) and previous findings at lower (e.g. Leauthaud et al. 2012;Allevato et al. 2012). There could be various and concomitant causes that could determine this offset. An explanation could be that the MultiDark simulations used to create the AGN mock catalogs are characterized by missing satellite halos possibly due to the low mass resolution (1.5 × 10 9 ℎ −1 ) and/or stripped halos. Another possibility is that luminous quasars do not reside in central DM halos with mass halo < 10 13 ℎ −1 as suggested by the HOD modelling of Richardson et al. (2012). In fact, as shown in Fig. 8, the 2pcf (including the 1-halo term) of SDSS quasars can be easily reproduced by mock AGNs with = 1 and applying a cut in the minimum mass of central AGN hosting halos at halo ∼ 10 13 ℎ −1 . Similarly, one might expect that the bias estimates as a function of stellar mass for COSMOS AGN Viitanen et al. 2019) can also be reproduced assuming smaller (∼1) and a minimum mass in the central halos of halo ∼ 10 12 ℎ −1 (following Richardson et al. 2013). These results suggest that the minimum central AGN hosting halo mass and sat are degenerate. The smaller the AGN satellite fraction is, the higher the mass needed for a central halo to host an AGN above a given luminosity. Moreover, the smaller sat , the steeper is the increase of the satellite AGN occupation as a function of the parent halo mass. It is certain that a fraction of quasars must be satellites to produce the small-scale clustering of AGN. Currently available satellite AGN fraction estimates in groups and clusters of galaxies presented by Martini et al. (2009) andPentericci et al. (2013) at ≤ 1 suggest small values of . Additional independent measurements of the AGN satellite fraction in groups and clusters at ≥ 1 would help in breaking the degeneracy in these model parameters.

CONCLUSIONS
Numerous degeneracies in the input parameters of cosmological models still prevent solid progress on the open issue of the coevolution of supermassive black holes (BHs) and their host galaxies and dark matter halos. Building on previous work from our group and by making use of advanced and diverse semi-empirical routines, also inclusive of the covariance among some input parameters, we here show that: • The overall shape and normalization of the large-scale bias as a function of AGN host galaxy stellar mass ( star ), is largely independent of the input stellar mass -halo mass relation, duty cycle and Eddington ratio distribution, while it is mostly driven by the dispersion in -not so much by the shape of -the input stellar mass-black hole mass relation.
• A model with covariant scatter, i.e. with a (positive) correlation between BH mass and galaxy mass at fixed halo mass, predicts an AGN bias almost independent of the stellar mass and substantially different from the bias of the underlying population of galaxies of the same stellar mass, at least in the range star > 10 11 . Present AGN clustering estimates at ∼ 1.2 do not allow us to clearly distinguish between models with and without a covariant scatter.
• The other parameter controlling the normalization of the AGN bias ( star ) is , the relative fraction of AGN hosted in satellite and central BHs of a given mass. Increasing the probability of AGN to be hosted in satellites rather than in centrals of equal BH mass, naturally increases the large-scale clustering as the bias becomes more heavily weighted towards more massive host halos.
• The comparison with the large-scale bias of COSMOS AGN at ∼ 1.2 and with the two-point correlation function of SDSS quasars at ∼ 1.4, suggests ∼ 4 which corresponds to a relative fraction of AGN hosted in satellites AGN sat ∼ 0.2−0.3. However, the data are also reproduced by models that adopt ≤ 2, i.e. values more consistent with independent estimates at the AGN fraction at lower z, as long as a cut is applied in the minimum mass of central AGN hosting halos, as also suggested by the HOD modelling in some clustering studies. This result unveils a strong degeneracy between the AGN satellite fraction and the minimum halo mass hosting AGN above a given luminosity. Independent estimates of the fraction of active satellites in groups at ≥ 1 will help in breaking this degeneracy.
In the next years, current and imminent extragalactic surveys, such as Euclid, eROSITA and LSST will precisely measure the clustering AGN at different masses and redshifts allowing to set invaluable constraints on many important features of AGN demography, such as limits on the covariance between AGN and galaxies, on the minimum halo mass hosting AGN, on the relative fraction of AGN satellites, and several others which will in turn provide essential constraints on the still puzzling co-evolution of BHs and their host galaxies and dark matter halos.