Quasars at intermediate redshift are not special; but they are often satellites

Understanding the links between the activity of supermassive black holes (SMBH) at the centres of galaxies and their host dark matter haloes is a key question in modern astrophysics. The final data release of the SDSS-IV eBOSS provides the largest contemporary spectroscopic sample of galaxies and QSOs. Using this sample and covering the redshift interval $z=0.7-1.1$, we have measured the clustering properties of the eBOSS QSOs, Emission Line Galaxies (ELGs) and Luminous Red Galaxies (LRGs). We have also measured the fraction of QSOs as a function of the overdensity defined by the galaxy population. Using these measurements, we investigate how QSOs populate and sample the galaxy population, and how the host dark-matter haloes of QSOs sample the underlying halo distribution. We find that the probability of a galaxy hosting a QSO is independent of the host dark matter halo mass of the galaxy. We also find that about 60\% of eBOSS QSOs are hosted by LRGs and about 20-40\% of QSOs are hosted by satellite galaxies. We find a slight preference for QSOs to populate satellite galaxies over central galaxies. This is connected to the host halo mass distribution of different types of galaxies. Based on our analysis, QSOs should be hosted by a very broad distribution of haloes, and their occurrence should be modulated only by the efficiency of galaxy formation processes.


INTRODUCTION
Two outstanding questions in extragalactic astrophysics are the manner in which galaxies sample the dark matter (DM) halo mass function, and how active galactic nuclei (AGN) sample the galaxy population. These questions are central because it is now believed that the energy and kinematics associated with AGN are crucial in understanding how galaxies form and regulate their star formation (see reviews by e.g. Fabian 2012; Kormendy & Ho 2013;Madau & Dickinson 2014;King & Pounds 2015;Somerville & Davé 2015;Xue 2017;Padovani et al. 2017).
In the case of a Gaussian random field, the two-point correlation function statistic provides a full characterisation of large-scale structure (see e.g. Peebles 1980; Bardeen et al. 1986;Wang et al. 2013). A given LSS tracer will display biased clustering, with an amplitude that increases for objects associated with rare massive haloes. Thus, by measuring the 2PCF of QSOs and assuming an underlying form of the dark matter halo mass function, it is possible to associate the QSOs with a given DM halo mass. The bias parameter, , is typically determined on linear scales ∼ 5 − 30 ℎ −1 Mpc, and this measurement and comparison has been carried out for a range of QSO redshifts, luminosities and colours.
Results from the last 10-20 years have traditionally placed QSOs in a mean DM halo mass of a few ×10 12 (e.g. Croom et al. 2005;Coil et al. 2007;Myers et al. 2007;da Ângela et al. 2008; Ross et al. 2009). However, one critical drawback of these measurements is that often only a single 'typical' halo mass is reported; given the very large range of black hole masses in QSOs ( BH = 10 6 − 10 10 ℎ −1 ), this single effective halo mass is relatively uninformative. Ideally, one would like to understand the full distribution of halo masses associated with QSOs, and then relate this to the underlying dark matter halo mass functions, and indeed to the way in which galaxies in general populate haloes. Are QSOs a subset of the full galaxy population? Or is luminous AGN activity found preferentially in one type of galaxy? Or in a particular environment?
An exception to the 'single halo mass' approach is White et al. (2012) who assume that QSOs reside in haloes with a lognormal distribution of masses, centred on a characteristic mass that scales with luminosity. This model results in a range of masses for luminous QSOs at redshift ≈ 2.4 from 0.8 − 4 × 10 12 ℎ −1 , with a central value of 2 × 10 12 ℎ −1 . Miyaji et al. (2011) and Krumpe et al. (2012) also studied the cross-correlation of ROSAT AGN with SDSS galaxies and constrain the host halo mass distribution of the ROSAT AGN. They predict that the satellite fraction of AGN reduces with host halo mass in contrast to luminous galaxies.
Several groups and authors have used the 2PCF to infer the halo occupation distribution (HOD) of the QSO population, and the 'satellite fraction'. The HOD provides a complete description of the relation between QSOs and dark matter at the level of individual virialized haloes (Berlind & Weinberg 2002;Kravtsov et al. 2004;Zheng et al. 2005;Chatterjee et al. 2012;Richardson et al. 2012). The QSO HOD is defined by ( | ℎ ), the conditional probability that a halo of virial mass ℎ contains QSOs above some specified (luminosity) threshold. A DM halo may contain zero, one, or more than one QSO. If more than one, the most massive galaxy is deemed to be at the centre of the potential well, and the less massive QSO(s) are the 'satellites' in that halo (even though they may still be relatively massive, e.g. 10 11 themselves). Richardson et al. (2012) present estimates of the 2PCF for QSOs, and interpret them with the HOD framework. In order to explain the small-scale clustering, the HOD model requires that a small fraction, sat = (7.4 ± 1.4) × 10 −4 , of the QSOs be satellites in DM haloes at 1.4. The median masses of the host haloes of central and satellite QSOs at these redshifts are constrained to be cen = 4.1 +0.3 −0.4 × 10 12 ℎ −1 and sat = 3.6 +0.8 −1.0 × 10 14 ℎ −1 respectively. Note that even thought centrals are expected to be the most massive QSOs in a given halo the satellites resides only in the massive haloes, whereas centrals QSOs can be found in relatively low mass haloes. Therefore, the median mass of host haloes of satellites QSOs is found to be two order of magnitude larger than central QSOs. Shen et al. (2013) also present measurements of the 2PCF, this time via the cross-correlation of 8,200 SDSS QSOs and 350,000 massive red galaxies from the SDSS-III Baryonic Oscillation Spectroscopic Survey (BOSS) at 0.3 < < 0.9. They estimate a QSO linear bias of Q = 1.38 ± 0.10 at = 0.53 corresponding to a characteristic host halo mass of 4 × 10 12 ℎ −1 (compared with a characteristic host halo mass for galaxies of 10 13 ℎ −1 ). Comparing these measurements with HOD models suggests that QSOs reside in a broad range of host haloes. The host halo mass distributions significantly overlap with each other for QSOs at different luminosities, implying a poor correlation between halo mass and instantaneous QSO luminosity. Shen et al. (2013) also find that the QSO HOD parameterisation is largely degenerate such that different HODs can reproduce the cross-correlations equally well, but with different satellite fractions and host halo mass distributions. Georgakakis et al. (2019) study the distribution of AGN host haloes using semi-empirical modelling. This model was shown to be consistent with current clustering measurements and found that AGNs host halo mass distribution is broad. They also predict that the fraction of satellite AGN increases towards the massive haloes. The Sloan Digital Sky Survey (SDSS) has been the state-ofthe-art in spectroscopic QSO surveys for the last 15 years. The Extended Baryon Oscillation Spectroscopic Survey (eBOSS) is the culmination of the SDSS-I, -II, -III and -IV quasar programmes and has recently completed a spectroscopic survey of >500,000 QSOs over 6,000 square degrees, covering redshifts 0.7 < < 3.5 (Lyke et al. 2020). eBOSS is currently the premier dataset to measure QSO clustering (Ross et al. 2020). In near future eROSITA (Merloni et al. 2012) will provide the most promising X-ray AGN sample. Eftekharzadeh et al. (2019) measure the quasar clustering signal across four orders of magnitude in scale, (0.01 100 ℎ −1 Mpc) at 1.5 using data from eBOSS. Using the HOD prescription, these authors find a satellite fraction of sat = 0.071 +0.009 −0.004 and minimum mass of min = 2.31 +0.41 −0.38 × 10 12 ℎ −1 for the host dark matter haloes best describes the quasar clustering on all scales. Rodríguez-Torres et al. (2017) used a modified Sub-halo abundance matching method to model eBOSS QSOs showing the mean host halo mass of 5 × 10 12 ℎ −1 . In this paper we extend these measurements of the clustering of galaxies and QSOs in eBOSS in order to: (i) understand the relation of the active QSOs to the general galaxy population; and (ii) understand the relation of the large-scale structure traced by the QSOs to the underlying DM halo distribution.
We will make progress by employing and expanding on recent work using the 'Multi-Tracer HOD' model (MTHOD: Alam et al. 2019a, hereafter Paper I) and apply this method to the latest version (Data Release 16) of the SDSS-IV eBOSS data. Our goal is to investigate a series of halo occupation distribution (HOD) models that described how QSOs populate the distribution of dark matter haloes. We will use the luminous red galaxies (LRGs) and the emission line galaxies (ELGs) that inhabit the same cosmological volume as the QSOs to perform these tests and investigations. This will allow us to discriminate between the models. This paper is organized as follows. In Section 2, we present our data sample. In Section 3, we describe the techniques involved in our measuring the 2PCF, and note several effects that could give rise to systematics in the measurements. In Section 4 we describe our Multi-tracer Halo Occupation Distribution (MTHOD) model and our Galaxy-QSO Occupation Distribution (GQOD) model that we use to model the galaxy and AGN population. In Section 5, we present our clustering measurements and the derived parameters. In Section 6 we place our new results in a broad context and note our main findings. We conclude in Section 7. Appendix A gives technical details. We assumes a flat ΛCDM cosmology with Ω = 0.307, Ω = 0.048, ℎ = 0.67, = 0.96 and 8 = 0.82. Our assumed cosmology is close to the best fit parameters reported in Planck Collaboration (2018) and motivated by the fiducial cosmology assumed in the -body simulation Prada et al. (MPDL2;) that we employ in our HOD models.

DATA
In this section we describe the spectroscopic data from the SDSS eBOSS survey that we will use for our clustering measurements. We also will utilise new deep public imaging data from the Hyper Suprime-Cam (HSC) Subaru Strategic Programme. . Sky coverage and number density distribution of our sample used in this paper. The top panel shows a ∼ 4 × 2 deg 2 patch of the eBOSS sample between redshifts 0.82 and 0.88. The three tracers LRG, ELG and QSO are shown in red squares, blue stars and cyan circles respectively. The varying symbol sizes represent the absolute AB -magnitude of each individual object. This clearly shows that QSOs are rare and bright compared to LRG and ELG galaxies. The bottom panel shows the number density distribution of each tracer in our sample. The red solid line, blue dotted-dashed line and cyan dashed line represents LRG, ELG and QSO respectively. In this paper we apply a redshift cut at 0.7 and 1.1 for our analysis.

SDSS-IV: eBOSS
We use data obtained from the extended Baryon Oscillation Spectroscopic Survey (eBOSS: Dawson et al. 2013) . This is one of the programmes of the wider 5 year Sloan Digital Sky Survey-IV (SDSS-IV; Blanton et al. 2017) using the BOSS spectrograph (Smee et al. 2013) on the Sloan Telescope (Gunn et al. 2006).
We use a subset of the eBOSS samples covering two fields between redshifts 0.7 and 1.1 where LRGs, ELGs and QSOs overlap in volume. The QSOs and LRGs sample the same area of the sky, but the ELG sky coverage is smaller. However, the ELG volume lies mostly inside the QSO/LRG volume by construction. Using the cross-correlation clustering technique, this overlap region can be also used to study the environemental dependence of QSOs, resulting in a sampling of the underlying dark matter distribution. We now briefly describe the relevant aspects of eBOSS sample selection, with more details being given by Prakash et al. (2016) for LRGs, Raichoor et al. (2017) for ELGs and Myers et al. (2015) for QSOs. Table 1 gives the number of objects from each selection used in this study. The total volume where all three tracers are observed is 0.64 ℎ −3 Gpc 3 between redshifts of 0.7 and 1.1.

Luminous Red Galaxy (LRG) selection
The LRGs are the most luminous and reddest galaxies, residing in massive dark matter haloes with high bias. The eBOSS LRGs are selected from SDSS imaging data (Albareti et al. 2017) in combination with infrared photometry from WISE (Wright et al. 2010) using the following target selection rules: where , and are the 'model' magnitudes of the SDSS photometric bands and W1 refers to the WISE magnitude in the 3.4 m channel.
The selections in equations 1, 2 and 3 are designed to achieve the redshift range, reduce stellar contamination and reduce lowinterlopers, respectively. The details of how these rules were derived and additional considerations are discussed in Prakash et al. (2016).

Emission Line Galaxy (ELG) selection
The eBOSS ELGs are expected to be star-forming galaxies at high redshift, and are thus selected based on high OII flux. The ELG sample is selected from DECAM Legacy survey (DECaLS: Dey et al. 2019). The target selection of ELGs in the North Galactic Cap (NGC) and South Galactic Cap (SGC), are slightly different due to availability of deeper data in SGC. We only use the SGC part of the ELG sample due to overlap with other tracers hence we only describe the SGC selection here. The ELG selection has two parts; the first part is to select star-forming galaxies corresponding to OII emission lines using the band flux cuts 21.825 < < 22.825.
The second rule for ELG selection is to preferentially select galaxies around redshift 1.0 given by following equations: where , , are the observed magnitude in DECaLS , and photometric bands. More details of how these rules were derived and additional considerations are discussed in Raichoor et al. (2017). Myers et al. (2015) describe in detail the requirements and how the eBOSS QSO sample is selected. First a supersample of QSOs is selected from SDSS imaging with either < 22 or < 22 and > 17, where and are the PSF magnitude of SDSS photometric bands and is the FIBER2MAG. This supersample is passed through the XDQSOz algorithm (Bovy et al. 2012), which assigns a probability for each object of being a QSO in a given redshift range. The eBOSS sample uses a probabilistic cut of QSO ( > 0.9) > 0.2. An infrared cut in WISE imaging is also used to remove stellar contamination. The final QSOs sample with good redshifts is obtained using (IMATCH=1 or 2) along with a target completeness cut ( eBOSS > 0.5) and spectroscopic completeness cut ( > 0.5) as described in Ross et al. (2020).

Quasar brightest 50% sample
We would also like to investigate how quasars populate their DM haloes as a function of QSO luminosity. QSO luminosity should de- pend on SMBH mass and mass accretion rate. With the link between SMBH mass and bulge mass established at low ( 0.1) redshift, and with a connection between bulge and halo mass, one might suspect that more luminous quasars (with more massive SMBHs) might populate their DM haloes in a different manner.
We split the full QSO sample into the brightest 50% of objects, as given by their −band absolute PSF magnitude, noting that the observed band samples 3660-4530Å rest-frame wavelength at = 0.86, where there are no strong broad emission lines. This is 28,742 objects which we call this sample the 'Brightest 50%'.
The top panel of Figure 1 shows a ∼ 4 × 2 deg 2 patch of eBOSS between redshifts 0.82 and 0.88. The three tracers LRG, ELG and QSO are shown in red, blue and cyan coloured filled symbols, respectively. The size of the symbols represents the absolute -band (AB) magnitude of individual objects. This clearly shows that QSOs are rare and bright compared to the LRGs and ELGs. The bottom panel of Figure 1 shows the number density distribution of each tracer in our sample. In this paper we apply redshift cuts at = 0.7 and = 1.1 for our analysis. We note that around redshift = 0.86, the mean redshift of our measurements, the number density of LRG, ELG and QSO are 10 −4 , 4 × 10 −4 and 2 × 10 −5 ℎ −1 Mpc −3 respectively. The redshift distribution of the Brightest 50% QSO sample is also shown in the bottom panel. The redshift distribution of the full QSO and Brightest 50% QSO samples are the same, so a direct comparison between the two is reasonable.

Hyper Suprime-Cam imaging
Deep imaging data from Hyper Suprime-Cam (HSC; Miyazaki et al. 2018) exist for a portion of our spectroscopic data. Imaging data from the the Hyper Suprime-Cam Subaru Strategic Programme (HSC-SPP) cover part of the sky shown in Figure 1. Using the 8.2m Subaru Telescope (Iye et al. 2004), the HSC-SSP (Aihara et al. 2018(Aihara et al. , 2019 currently offers the best combination of depth and image quality for a ground based survey. The Wide Layer achieves depths of g = 26.6, r = 26.2, i = 26.2, y = 25.3 and z = 24.5 in the five broad-band filters. The seeing ranges from 0.58 to 0.77 (Aihara et al. 2019). All data products are available for the Second Data Release at https://hsc-release.mtk.nao.ac.jp/doc/.
The HSC-SSP data are of high enough quality to see galaxy groups out to ∼ 1 (e.g. Umetsu et al. 2020). Thus, we will use the HSC data to visually inspect the environments of the = 0.7 − 1 eBOSS QSOs.

Clustering and the 2-Point correlation function
Here we give a brief description of the 2PCF; for a more formal treatment the reader is referred to e.g. Peebles (1980). The 2PCF, ( ), is defined by the joint probability that two objects (e.g. galaxies) are found in the two volume elements 1 and 2 placed at separation , with being the object number density. To calculate ( ), points are given inside a window of observation, which is a 3D body of volume ( ). We calculate the position of each galaxy in 3dimensional space by converting the measured redshift to a lineof-sight distance using our fiducial cosmology. As usual, we also generate a catalogue of random points, with the same window function as the data, but without correlated positional information.
We then measure the galaxy auto-correlation function using the minimum variance Landay-Szalay estimator (Landy & Szalay 1993) given by: where , and are the number of galaxy-galaxy, galaxyrandom and random-random pairs as a function of vector separation in 3-dimensional space. The galaxy, QSO and random catalogues for the large-scale structure measurements for the SDSS DR16 sample are publicly available 1 . The cross-correlations are measured using the following estimator: where 1 2 and 1 2 are the number of galaxy-galaxy and galaxy-random pairs from different samples. As is conventional, we project the 3-dimensional space onto a 2-dimensional space that decomposes pair separation vectors along the line-of-sight ( ) and perpendicular to the line-of-sight ( ⊥ , or ). This gives us the 2-dimensional correlation function ( , ).
In practice, we measure the projected correlation ( ) by integrating the 2-dimensional correlation function along the line-ofsight between = −40 ℎ −1 Mpc to = +40 ℎ −1 Mpc and using 25 equally spaced bins in logarithmic scale for ⊥ between 0.1 ℎ −1 Mpc and 30 ℎ −1 Mpc. Typically, the projected correlation function is integrated to = 100 ℎ −1 Mpc or larger to avoid the need to model redshift space clustering. But our model is evaluated in redshift space with full non-linearity and hence we do not have any constraint on minimum for the projected correlation function. The projected correlation function ( ) helps us constrain the HOD parameters that govern the galaxy-halo connection. To estimate the errors of our measurements, we create 86 jackknife regions for our sample and calculate the jackknife covariance of the measurements. Note that our jackknife region in the overlap sky area corresponds to approximately 6 deg ×1.3 deg. The 1.3 deg at the mean redshift of the sample corresponds to roughly 40 ℎ −1 Mpc and hence large enough for our measurements.

Potential systematic errors
The clustering measurement is sensitive to the completeness of the observed galaxy sample for a given selection scheme. Therefore, it is important to account for variation in the number of detected galaxies as a function of their position in the sky, plus various selection biases due to the systematic errors introduced by instrumentation and measurement.
The number of detected galaxies and the spectroscopic success rate can be correlated with, for example, stellar density, extinction, sky brightness, airmass or position in the fibre plate (see e.g., Bolton et al. 2012). To remove these correlations, Ross et al. (2012) introduced the use of systematic weights, and the use of systematic weights for LRGs, QSOs and ELGs is investigated by Bautista et al. (2018), Gil-Marín et al. (2018) and Raichoor et al. (2021), respectively. We show in Paper I that measurement of at scales between 1-30 ℎ −1 Mpc is insensitive to the resulting corrections by any of the introduced systematic weights. We also note that due to fibrecollision the galaxy sample becomes highly incomplete below the fibre scale, which is approximately 65 on the sky, corresponding to a scale of 0.5ℎ −1 Mpc at = 0.86 (Anderson et al. 2012;Guo et al. 2012;Bianchi & Percival 2017). For the purpose of this paper we will not use the measurements or our models at scales smaller than the fibre collisions.

Environmental measures
We also measure galaxy environment following the method described in Alam et al. (2019b). Here we briefly summarise the method. We focus on the measurement of local overdensity of galaxies around our sample.
In order to measure galaxy overdensity, we first create the Voronoi tesselation of the survey volume using the method developed in Alam et al. (2019b). This partitions the volume into disjoint cells each containing only one galaxy. We use the random catalogue provided with the Large Scale Structure catalogue and count the number of randoms in each Voronoi cell ( rand cell ). We then estimate the density as the ratio of the mean density of randoms ( rand ) to the random counts for each cell, and associate this with the galaxy in that cell, cell : The measured density is then smoothed at this chosen scale to determine the smoothed density ( smth ) as with smth = 5 ℎ −1 Mpc. This is then converted to an overdensity, 5 : This allows us to assign a value of 5 for each QSO and galaxy in the observed sample living in the overlapped volume. We finally measure QSO fraction as the function of 5 using: where = log 10 (1 + 5 ) and tracer ( ; Δ) gives the weighted count of number of object of a particular tracer with between − Δ and + Δ. In these counts, we only consider objects with −1 < < 1 and divide this in five bins with Δ = 0.1. We use only five bins in order to keep the size of the covariance matrix small while still having the overall trend of QSO with . Each object in the sample is weighted such that the redshift distributions of all objects are the same.

MODELLING THE GALAXY AND QSO POPULATIONS
Our aim is to use the (cross-)clustering measurements of the LRG, ELG and QSO populations as constraints for our HOD models in order to understand the QSO population. We employ two HOD models: (i) the Multi-tracer Halo Occupation Distribution (MTHOD) model to model the overall galaxy population and (ii) the Galaxy-QSO Occupation Distribution (GQOD) model in order to model the statistical properties of QSOs as a distinct sub-population from the parent galaxy population. The anzatz we use to model the QSO population is that the active galactic nuclei (AGN) observed as QSOs are not special in their inherent host galaxy properties, but are a sub-sample of the global galaxy population.

MTHOD galaxy catalogues
To model the galaxy population we use the Multi-tracer Halo Occupation Distribution (MTHOD) model and catalogue of Alam et al. (2019a). The MTHOD model introduces a new approach to model multiple tracers in the same volume. In general each of the tracers can have its own occupation recipe for the central and satellite galaxies. At the same time, the MTHOD ensures that the joint occupation probabilities are well behaved by limiting the total probability of central galaxies in a halo to 1. It also forbids the non-physical situation of multiple types of galaxies at the centre of the same dark matter halo. The key parameters in MTHOD models are the separate parameters for the occupation probability of central and satellite galaxies for each tracer; these are given in the Appendix.
The DM haloes are then populated using the following equations for central and satellite galaxies as a function of halo mass, halo : where the sum is over all tracers in the list, TR = {LRG, QSO, ELG}. This equation requires a constraint of tot cen ≤ 1 for any halo mass. In paper I, all three tracers (i.e. LRGs, ELGs and QSOs) are modelled within the MTHOD framework. However, in this paper, we take a different approach. We only use the LRG and ELG galaxies, and do not use the QSOs from the default model. Using the MTHOD mock catalogue, we can measure the clustering and the central/satellite properties for the LRG and ELG populations.
We assume that the MTHOD galaxy catalogue models the complete set of galaxies that host eBOSS QSOs. The MTHOD model samples galaxies starting from a minimum halo mass of 2.1 × 10 11 . The mean halo mass of eBOSS QSOs is shown to be 5 × 10 12 (see Figure 4 of Alam et al. 2019a). The eBOSS galaxies thus cover the entire halo mass range needed to model the eBOSS QSO selection. Therefore, in the absence of any strong environmental effect or assembly bias, it is reasonable to assume that the MTHOD galaxy catalogue models the complete set of galaxies that host eBOSS QSOs. We do expect this assumption to fail in detail (as indicated by Figure 8 of Alam et al. 2019a) but this should be a second order effect given the current avaialble constraints on such effects

Galaxy QSO Occupation Distribution (GQOD) model
A second model, the Galaxy QSO Occupation Distribution (GQOD) model is employed to model the statistical properties of QSOs as a distinct sub-population from the parent MTHOD galaxy population. The probability of a galaxy being a QSO is given by: where on is the probability for the galaxy to turn on with a given host halo mass, ( type ) is the probability for a galaxy to be a QSO given its host galaxy type and ( pos ) is the probability for a galaxy to be a QSO given its host galaxy position (central/satellite). A summary of these key parameters of the GQOD model are given in Table 2 and we do consider the three probabilities to be independent of each other. Here, "turn on" is shorthand for the physical processes involved having sufficient mass accretion in the QSO central engine that the QSO becomes luminous enough to be detected in our survey volume. We assume on , the fraction of galaxies that will have a QSO turn on, is a function of DM halo mass, halo . This might have a wide distribution in halo mass (see e.g., the deduced wide range of halo mass for central galaxies in Richardson et al. 2012Richardson et al. , 2013. Alternatively, QSOs might reside in dark matter haloes of a certain particular mass range and hence the turn-on probability will have a narrow distribution with halo mass (see, e.g., Kayo & Oguri 2012;Eftekharzadeh et al. 2019). We model on ( halo ) with a 'linear spline sampling' of 8 halo masses between 10 11 ℎ −1 and 10 15 ℎ −1 . Moving to type -the probability that a host galaxy of a QSO may depend on host galaxy type -we introduce the parameter LRG to define the fraction of QSOs with LRG host galaxies in the sample. We define LRG = 1 to mean that only LRGs can host QSOs; consequently, LRG = 0 will mean that only ELGs can be QSO hosts. If the posterior likelihood of either of these extremes is zero, then we can rule out the QSOs being turned on in only one type of galaxy. In a model where the QSO probability is independent of host galaxy type, the fraction of QSOs with a given type of host galaxy should be the same as the fraction of the host galaxy type in the parent population. LRG (and thus ELG = 1 − LRG ) can be measured directly from data, as is done in Alam et al. (2019a, their Figure 7). However, the match between the measured and modelled LRG and ELG is very close, so we are able to just use the smoother models.
The probability a host galaxy is a QSO may also depend on 'position', that is, whether it is a central or satellite galaxy. Studies including Zheng et al. (2009);Richardson et al. (2012Richardson et al. ( , 2013 assumed different halo distributions for central and satellite galaxies, while Kayo & Oguri (2012) and Eftekharzadeh et al. (2019) assumed the halo to be indifferent in hosting a QSO as either a satellite or central galaxy. We include the parameter pos to encapsulate the positional information. We adopt two model flavours to study this aspect of QSO physics. In the first flavour, there is no dependence of QSO probability on whether the host galaxy is a central or satellite. Thus, the fraction of QSOs hosted by satellite galaxies will be equal to the fraction of satellite galaxies in the parent population. We call this an 'inherent' satellite fraction and label it ' sat (inherent)'. In the second flavour, the QSO probability can depend on whether the host galaxy is a satellite or a central and this is denoted by an additional parameter, sat . This allows an extra degree of freedom where the satellite fraction in the QSO population does not have to represent the satellite fraction of the parent galaxy population. We call this an 'enforced' satellite fraction, labelled ' sat (enforced)'.
The full list of parameters in our model is given in Table 2. We have a total of 9 parameters in the inherent sat model with 8 of them to model the halo mass dependence of on using linear spline with 8 knots listed in Table 2 and the 9th parameter is to model fraction of QSO with host galaxy as LRG ( LRG ). The enforced sat model) has an additional 10th parameter to model the fraction of QSOs that are satellite galaxies sat .
Once we have simulated the QSO catalogue, we can predict the projected auto-correlation of quasars as well as the cross-correlation of quasars with LRGs and ELGs. Given the denser galaxy population one can also measure the fraction of QSOs as a function of overdensity. We use these four measurements to constrain the parameters of the model.

Model constraints and parameter estimation
For any point in the parameter space one can evaluate the probability of a galaxy to turn on using equation 17 and hence obtain a sample of QSOs in the MTHOD catalogue. The LRG-auto, ELG-auto and LRG×ELG cross clustering measurements are used to constrain the MTHOD model. Since these are presented in Paper I, we do not report them again here. The QSO-auto, QSO×LRG and QSO×ELG cross-correlation, and QSO fraction are the (new) measurements used to constrain the GQOD.
We use emcee (Foreman-Mackey et al. 2013) to sample the parameters for each of the two models in this work via a Markov chain Monte Carlo (MCMC) process. We then evaluate the auto and cross-correlation of QSOs with LRG and ELG samples using methods described in Section 3. This measurement requires pair counting to be performed at each step of sampling. We use the publicly available code corrfunc (Sinha 2016) to evaluate pair-counts efficiently at each iteration. We also pre-compute the overdensity field for each galaxy in the MTHOD catalogue using the method described in Alam et al. (2019b) -which is then used to evaluate the fraction of QSOs as a function of envoronmental overdensity. We developed a python library to efficiently create QSO samples as a fuction of our model parameters 4 . In this process we use the number density of the QSO sample, treated as fixed at the mean redshift.

Parameters Description Used with
on ( halo ) This models the probability of a galaxy to turn on with a given host halo mass. This is modelled as a linear spline with 8 knots at locations Inherent & Enforced log 10 ( knots ) = [11. 2, 11.8, 12.15, 12.5, 12.75, 13.0, 13.5, 15.2] It also require the constraint of ∫ ∞ 0 on ( halo ) gal ( halo ) halo = QSO where gal ( halo ) and QSO are the number density of galaxies and QSO respectively.
Probability of galaxy to turn into QSO given its host galaxy type We then evaluate a model 2 using following equation: where represents sampling parameters and −1 is the inverse covariance matrix obtained from jackknife analysis including the errors in estimate of model prediction. The MCMC process then samples the parameters according to the 2 given by the model.

RESULTS
In this section we present our results. Our convention is to report data measurements with black data points. The two flavours of the GQOD model will be represented by red (solid) and blue (dashed) lines for the 'inherent' and 'enforced' model fits, respectively.

Clustering results
We analyse the eBOSS QSO auto-correlation and the QSO-LRG and QSO-ELG cross-correlations between redshifts of = 0.7−1.1. The measurement of these statistics along with the best-fit GQOD models is shown in Figure 2. The top-left, top-right and bottomleft panels of Figure 2 show the QSO auto-correlation, ( ), the QSO×ELG cross-correlation ( ), and the QSO×LRG cross-correlation ( ) functions, respectively. The measurement from eBOSS data is shown with black circle with error bars estimated from the jackknife sampling method. The vertical black dashed line shows the fibre collision scale at the mean redshift ( = 0.86) and the projected correlation function below this scale are shown with open circle and not used in the analysis. The dashed blue and solid red lines with shaded region shows the enforced and inherent sat models respectively with its 1 constraint in all panels. Note that different scales are correlated which is accounted for in our full covariance matrix.
The models describe the data very well in the range of scale considered in the current analysis. The smaller 1 regions in the QSO×ELG correlation function is due to the higher gal number density of the ELG sample. Since we do not sample the smaller ('1-halo' term) scales, there is no discernible difference in the models whether or not the satellite fraction is an additional free parameter.

Environment and QSO fractions
We measure the QSO fraction as a function of galaxy overdensity environment. Recall, QSO ( ) is measured straight from the data using the Voronoi cell estimation method (Section 3.3). This is presented in the bottom-right panel of Figure 2.
The QSO fraction is around 3% across the whole range of environments sampled (black points in Figure 2, bottom right). The best-fit models (red and blue lines with 1 shadings) again fit the data very well. The models rise monotonically from around 2.6% for the least-dense regions to around 3.4% for the most dense environments.
In Figure 2, we also show the LRG fraction, LRG , as the solid (black) line and the ELG fraction (where ELG = 1 − LRG ), as the dashed-dotted (black) line for comparison. The behaviour of the LRG and ELG fraction with overdensity is probably dominated by galaxy quenching being more efficient in the high halo mass regime (Alam et al. 2019a), hence leading to the high LRG fraction in the most overdense regions. We note that the QSO fraction has essentially a flat dependence on overdensity. This is in contrast to both the LRG or ELG fractions, which have an environmental density dependence. This implies that the QSO population must contain a mixture of LRGs and ELGs as host galaxies. Figure 3 shows the best-fit GQOD values derived for on (the fraction of galaxies that will have a QSO turned on), the fraction of QSOs that are satellites ( sat ) and the number density of QSOs, QSO , as a function of DM halo mass. As before, the best-fit sat (inherent) model is given by the solid (red) line, with the sat (Enforced) model given by the dashed (blue) line. 1 errors are given by the shaded regions.

Dependence on DM halo mass
In the top panel of Figure 3 we see the fraction of galaxies that have a QSO turned on is essentially independent of halo mass for both the Enforced or inherent models, i.e., the probability that a galaxy has quasar activity is independent of halo mass. This is a key result: it implies that the halo mass distribution of QSOs is very broad, despite the model having freedom to choose a narrow range of halo mass through our spline fit. In order to test that the flat nature of on is indeed a better model we tested a lognormal model for on . We looked at 2 value in a grid of mean halo mass and scatter for lognormal model finding the minimum 2 /dof = 72/53 for a model with mean halo mas of 10 13 ℎ −1 with width 1.05. We also observed that the smaller width is strongly ruled out by having large value of 2 whereas larger width increases 2 marginally. When this is compared with our default model giving best fit with 2 /dof = 42/45 then we can be more confident that the flat on describes the data very well.
The middle panel shows the satellite fraction as a function of halo mass. The red line shows the satellite fraction for the best-fit inherent sat model, while the blue line shows the satellite fraction for the best-fit enforced sat model. The fraction of satellites being QSOs rises from 0% at masses above 5×10 12 . Around 2×10 13 . Both models put the fraction of QSOs being in satellites galaxies at over 50%. For the dashed blue line that represents the satellite fraction in the enforced sat , the model slightly prefers the QSOs to turn on in satellite galaxies (more notable at lower halo mass) leading the blue dashed line to always be above the red dashed line (see section 6.2 for more discussion). Georgakakis et al. (2019) estimate the satellite fraction for AGN as the function of X-ray luminosity finding it to be around 10-20% consistent with our satellite fraction for eBOSS QSOs.
In this middle panel, we also show for comparison the fraction of satellites for LRGs and ELGs with solid black and dashed dotted lines. There are no satellite galaxies that are LRGs in haloes of halo < 10 13 (as found in previous LRG HOD studies: e.g. Zehavi et al. 2005;Ross et al. 2008;Reid & Spergel 2009;White et al. 2011). In haloes of 10 13 < halo < 10 14.5 , the fraction of satellites that are LRGs increases from 0 to just under 60%. Then, for haloes with halo > ∼ 10 14.5 the fraction of satellites that are LRGs jumps to near 80%. This is because a halo can only have one central but more massive haloes can have multiple objects defined as satellites (e.g. in groups and clusters); thus the satellite fraction can, and will, approach 100%. The reason sat levels at 60% for the LRGs is because we have many more satellites that are ELGs; in a model without ELGs, then by default the satellite fraction for LRGs would be 100% in the most massive haloes.
The bottom panel shows the number density of QSOs per unit logarithmic halo mass. The black solid line shows the same distribution for the parent population scaled down by a factor of 25. We note that the halo mass distribution of QSOs is very similar to that of the host population, which comes from the fact that the turn-on probability of QSOs is independent of halo mass. We note that this predicts a very broad distribution of QSO host DM halo mass. Overall, quasars inhabit dark matter haloes in essentially an Figure 3. The top, middle and bottom panels show on , sat and number density of QSOs per unit logarithm of halo mass as the function of halo mass respectively. In all panels the dashed blue and solid red lines show the two models with enforced sat and inherent sat respectively with shaded regions representing 1 errors. The on in the top panel presents the probability of a galaxy to turn on, which is also called the duty cycle of the QSO. The middle panel shows the fraction of QSOs living in satellite galaxies. For comparison, we also show in this panel the fraction of satellites for LRGs and ELGs, using black solid and dashed-dotted line respectively. The bottom panel shows number density of QSOs per unit log 10 ( ℎ ) and the black solid line shows the distribution for the parent galaxy population, scaled down by a factor of 25 for comparison. identical way to the full galaxy population, although they are not as common.

QSO dependence on luminosity
In Figure 3 we also show, with the cyan line and shading, the best-fit GQOD values derived for the enforced sat model, but only using the brightest 50% of the QSO data.
Using the Brightest 50% data, and refitting the best-fit models, we see from Figure 3 in the top panel that the probability for a galaxy to have quasar activity is 50% for the brightest half of the quasar population, but remains independent of halo mass. The amplitude of the model fit simply comes from the number density; indeed, Figure 4. This shows the two-dimensional posterior of the fraction of satellites in QSOs ( sat ) and the fraction of LRGs in QSOs ( LRG ) for the two models. We note that for the red contours is a derived parameter, whereas forthe blue contours it is a free parameter. This shows that QSOs cannot come entirely from either LRGs or ELGs but rather are a mixture with roughly equal proportions. We also note that QSOs could have significantly larger satellite fractions and in the blue contour where we allow satellites to have non-equal probabilities of converting to QSOs then we find that the data prefer even larger satellite fractions.
if one were to reduce the QSO population via a random sampling fraction , then the probability of a galaxy having quasar activity will also reduce by a factor . The key thing to note is this function remains flat, i.e. independent of halo mass.
From Figure 3 and the cyan line in the middle panel (which is very similar to, but under the blue line/shade), the fraction of satellite galaxies that are quasars is independent of quasar luminosity.
Finally, in the bottom panel of Figure 3, we see that overall, the 50% more luminous quasars inhabit dark matter haloes, as function of mass, in essentially an identical way to the full QSO population. As stated above, the probability of turn on for the smaller number of more luminous QSOs is reduced, though we have same number of galaxies in our parent sample. Consider the inherent model (red contours): since the models do not contain LRG = 1 or LRG = 0, this model rules out (at the > 3 level) the possibility that the QSO population comes entirely from either the LRG or ELG population (when we do not allow the extra degree of freedom in the inherent satellite fraction model). It also shows that roughly 20% of QSO host galaxies are satellite galaxies with about 60 − 70% of QSO hosts being an LRG galaxy.

QSO dependence on galaxy type and position
When we allow the additional freedom from satellites in the Figure 5. Cutout images from HSC data release 2 around six different QSO randomly selected in our sample. Each of the cutouts is centred around an eBOSS QSO, with a field of view about 1 across. This shows a sampling of the different environments that QSOs inhabit. The cutouts were created using the HSC public data access tools. enforced sat model (shown in blue), we are still able to rule out the possibility of the LRGs hosting the eBOSS QSO population at the 3 level; but QSO host galaxies being entirely ELGs is possible (at the 2 level). Allowing the additional freedom in the model, the enforced sat prefers a larger fraction, ∼40%, of QSO hosts to be satellite galaxies, but a lower fraction (20 − 60%) of LRG host galaxy (compared to the inherent model). This is due to both the LRGs and satellites residing in haloes with relatively higher halo mass. Therefore one can get the same clustering by either having a large LRG fraction with lower satellite fraction as in the inherent model, as by having a large satellite fraction with a smaller LRG fraction (as is the case for the enforced model).
It is also interesting to note that in either case the satellite fraction between 20% − 40% is much larger than host galaxy population and disagrees with other analysis of QSO satellite fraction (Starikova et al. 2011a;Richardson et al. 2012). Generally it is considered that large satellite fraction will mean enhanced clustering within the 1-halo term and hence our large satellite fraction might mean we overpredict the QSO auto-correlation at scales below 1 ℎ −1 Mpc. But this is not the case, as shown in the top left panel of Figure 2 where our model is consistent with the observed QSO clustering at these smaller scales. This is due to the QSO number density being very low; thus, even when we allow the model and QSOs to have a large satellite fraction, there will only be a few QSOs per halo and hence the central-satellite or satellite-satellite terms will not be as strong.

HSC Imaging of QSO Groups
Our results suggest that a large percentage of QSO host galaxies are satellites. We turn to the deep HSC imaging data to see if we can find direct examples of this.
A collage of HSC cutouts centred on six QSOs from the eBOSS sample is shown in Figure 5. The filters were mapped to a RGB colour scheme using the HSC public data access tools 5 .
The six QSOs are the ones shown in the top panel of Figure 1 with a declination cut > −1.0 • to restrict to the region with HSC imaging. The field of view is 30 on a side, corresponding to a scale length of 240 kpc at = 0.86 for our given cosmology.
In the collage we see several different example environments for the QSO. Note, the QSO appears in the centre of the image; this does not indicate at all that the QSO lies the centre of the group, or at the peak of the DM halo. From Figure  Without having spectroscopic redshifts for the objects in these cutouts, it is of course very tricky to confirm galaxy group membership. However, given the depth and seeing quality of these images, and the lack of obvious foreground galaxies, it is entirely reasonable to assume we are seeing the environments of the QSOs.

DISCUSSION
In this section, we discuss our key results and their implications. We first discuss how the QSO population traces the underlying dark matter halo mass function (Section 6.1). We then place our results for the QSO satellite fraction in context, comparing with previous studies in Section 6.2. In Section 6.3 we think about the inferences that may be made from our clustering results regarding the growth history of stellar and black hole mass. In Section 6.4 we note the immediate implications of this current work and point towards future investigations.

QSOs and the dark matter halo mass function
One key motivation for this study was to answer the question: how do QSOs populate their dark matter haloes as a function of mass? Figure 3 shows us the answer here: the parameter on is flat, and so the fraction of galaxies that have a QSO turned on is essentially independent of halo mass for both the enforced or inherent models. The probability that a galaxy has quasar activity is therefore independent of halo mass. The result in the bottom panel of Figure 3 immediately follows: the QSOs sample the same DM halo mass function as the parent galaxy sample. Conroy & White (2013) find the same result. In fact, their model is consistent with QSOs being equally likely to exist in galaxies, and therefore dark matter haloes, over a wide range in masses and suggests a single QSO duty cycle at redshift < 3.
In fact, this observation has been already recognised in the X-ray community (e.g. Leauthaud et al. 2015;Mendez et al. 2016;Powell et al. 2018;Plionis et al. 2018;Krishnan et al. 2020). All these studies find no significant differences in the clustering properties of X-ray AGN compared to a matched galaxy sample. That is, X-ray AGN inhabit DM haloes that are consistent on average with the overall inactive galaxy population. This result runs from the local Universe (0.01 < < 0.1; Powell et al. 2018) to high redshifts, 4.5 (Krishnan et al. 2020). Studying narrow-line AGN, and comparing to a matched control samples of inactive galaxies, Li et al. (2006) find that AGN have almost the same clustering amplitude as the control galaxies, on scales larger than a few Mpc. Here we are showing, for the first time, that the same is true for the blue optically selected broadline QSOs at modest redshift.
For the Brighter 50% QSO sample, we find that on remains flat, and thus this sample also directly maps to the same DM halo mass function as the parent galaxy sample. This is another key result as it immediately explains the lack of dependence on luminosity for QSO clustering (e.g. da Ângela et al. 2008;Shen et al. 2009Shen et al. , 2013Chehade et al. 2016;Powell et al. 2020). This is consistent with the result that the QSO luminosity possibly has large scatter at fixed black hole mass and hence do not particularly correlates with the host halo mass Hickox et al. (2014).

The fraction of QSOs that are satellite galaxies
We remind the reader that our estimate of the fraction of QSOs that are in satellite galaxies comes from our GQOD model and the MTHOD mock catalogue: we derive the QSO satellite fraction, sat , by comparing the model to the data. Our results suggest that above a halo mass of 10 12 a substantial number of QSOs can be found in satellite galaxies.
A range of previous studies give a very large range of measured satellite fractions for QSOs. Richardson et al. (2012) report a satellite fraction for 1.4 QSOs as sat = (7.4 ± 1.4) × 10 −4 . This is a much smaller figure than our large percentage, though it is not straightforward to compare since Richardson et al. (2012) quote a probability density function of the satellite fraction as given by all their HOD models, whereas we have the satellite fraction as a function of DM halo mass. Also, and using very similar data to Richardson et al. (2012), Kayo & Oguri (2012) find sat = 0.054 +0.017 −0.016 , i.e. nearly an order of magnitude smaller than our figure of >30% at halo 3×10 12 . This apparent discrepancy in satellite fraction can possibly be explained by realizing that these studies use smallscale pairs built from samples of binary quasars. If, for instance, binary quasars consisted of two QSOs that had very similar masses, then it is not clear that one being a satellite and the other being a central is a meaningful distinction. In this scenario, on average, a "non-central" member of each binary quasar could be interpreted as part of a small fraction of high-mass satellites. We do note that we do not sample the 1-halo term very well (due to fiber collisions), and so suggest that our sample samples quasars in e.g. groups where the QSO can be a satellite, and this occupation fraction might be (very) different from e.g. binary QSOs in the same halo but potentially different sub-haloes. Starikova et al. (2011b) present results showing that the Chan-dra/Boötes AGN are predominantly located at the centres of dark matter haloes and tend to avoid satellite galaxies in haloes of this or higher mass. However, also using moderate luminosity X-ray AGN (at < 1 from the COSMOS field), Leauthaud et al. (2015) report a mean satellite fraction of sat = 18 ± 2 %. Wang & Li (2019) is another key result here. Using a sample of 100,000 AGN from SDSS (with the AGN being classified via the BPT diagram, Baldwin et al. 1981, and lying mostly below = 0.3), these authors also perform a clustering measurement in order to investigate the DM halo properties of narrow-line AGN. Wang & Li (2019) investigate the central/satellite fraction for the AGN as a function of stellar mass, and host galaxy colour, so a direct comparison to our results is tricky. However, these authors do see a substantial fraction of AGN hosts being non-central galaxies, especially at lower, ★ < 10 10.5 stellar mass (their Figure 7). Jiang et al. (2016) study the AGN population at low redshift and found that Type I and Type II AGN resides in dark matter haloes with similar masses. But the satellite fraction of Type I AGN are smaller than the Type II AGN. They suggest that this has interesting implications in the QSO unified model as they detect environmental differences between Type I and Type II AGN. It will be interesting to see in the future if a similar difference can be observed at high redshift using our model.

Host galaxy type for eBOSS QSOs
In the local Universe we observe a relationship between SMBH and host-galaxy bulge velocity dispersion ( ) and luminosity (Bell 2008;Gültekin et al. 2009;Volonteri 2010;McConnell & Ma 2013;Kormendy & Ho 2013); supermassive Black holes are known to be related to their host bulges. Thus, it is not surprising that galaxies with generally massive bulges, such as LRGs, would be considered QSO host galaxies.
Returning to Wang & Li (2019), they find that AGN in galaxies with blue colours at all masses, or massive red galaxies with ★ 10 10.5 , show almost identical clustering amplitudes at all scales to control galaxies of the same mass, colour, and structural parameters.
As mentioned above, direct comparison is difficult since our results are reported as a function of DM halo mass, and Wang & Li (2019) report stellar mass. Moreover, there is definitely not a 1-to-1 mapping of the red/blue galaxy population in Wang & Li (2019) to the LRG/ELG population we study. That acknowledged, we can look at broad trends. Wang & Li (2019) find that there is slight preference for AGN to trigger in red satellite galaxies. This is broadly consistent with our findings. Wang & Li (2019) also find the blue AGN are less likely to be satellite then the general blue galaxy population. This is again consistent with our findings (our Figure  3, middle panel), which shows that if we allow sat to be a free parameter, then the QSO satellite fraction is higher at lower halo mass, but become less than the blue galaxies at higher halo mass. We find the transition halo mass is around halo ∼ 10 13 . Wang & Li (2019) model their results with a simple halo model -where central fraction of the AGN is the only free parameter -and place AGN preferentially in the dark matter halo centres, but requiring a mass-dependent central fraction. Their results suggest that the mass assembly history of dark haloes may play an additional role in the AGN activity in low-mass red galaxies.
Interestingly, Krishnan et al. (2020) note that the most important property in determining the AGN clustering signal is the fraction of AGN in passive host galaxies. This is true for our study as well using the inherent model. We note that our results are in contrast with those of Matsuoka et al. (2014), who studied high signal to noise spectra and inferred a very small red fraction for QSO hosts.

Quasars at intermediate redshift are not special, but they are often satellites
The results in our paper illuminate several outstanding issues in QSO physics. The mean halo mass of QSO has historically been measured to be a few ×10 12 , (e.g. Croom et al. 2005;Ross et al. 2009;White et al. 2012), essentially corresponding to a small group mass (cf. The Local Group having a total mass of 5.27 × 10 12 : Li & White 2008). Yoon et al. (2019) find that, on average, both massive quasars and massive galaxies reside in environments more than 2 times as dense as those of their less massive counterparts with log 10 ( BH / ) < 9.0. However, massive quasars reside in environments about ∼2 times less dense than inactive galaxies with log 10 ( BH / ) 9.4, and only about one third of massive quasars are found in galaxy clusters, while about two thirds of massive galaxies reside in such clusters. This indicates that massive galaxies are a much better signpost for galaxy clusters than massive quasars.
This is also what we are seeing. We are also finding QSOs to be hosted by galaxies with massive bulges, i.e. the LRGs. But the key thing is that neither the QSOs or the LRGs are necessarily central galaxies. This explains why the QSOs generally have a lower mean clustering amplitude than the massive galaxies, unless they were e.g. radio loud objects Shen et al. 2009), such as the radio loud giant ellipticals (e.g. M87 or similar). Thus, a key conclusion of our work is that quasars are not special, but they may well often be satellites.
Also at ∼ 4, recent studies of the quasar environment find no strong evidence of luminous quasars to reside in dense environments or be associated with proto-clusters (see, e.g., Uchiyama et al. (2018); also Overzier (2016) and references therein). They do find quasars to reside in small size haloes that that are much more in accord with typical average halo masses found at lower redshifts (Eftekharzadeh et al. 2015;White et al. 2012) than reported halo mass by Shen et al. (2007).
We are also finding that QSOs can inhabit smaller haloes, of the kind that are dominated by star-forming galaxies. Possibly these two classes of object have different triggering mechanisms, but as far as optical luminosity is concerned we can not differentiate the two cases.

CONCLUSIONS
In this paper we have used the final SDSS eBOSS DR16 spectroscopic dataset for luminous red galaxies, emission line galaxies and QSOs to make multi-tracer clustering measurements. The motivation is to (i) investigate how the QSO population samples the galaxy population; and (ii) to understand how the QSO host dark matter haloes sample the underlying dark matter halo distribution.
Our main conclusions are: • The probability that a galaxy has quasar activity is independent of dark matter halo mass.
• QSOs host galaxies have a large satellite fraction, probably due to their low number density (and this is possible, even without a large one-halo term).
• We infer the halo mass distribution of QSOs to be very broad, independent of assumptions in modelling about the parent population.
• All QSOs can not be in LRG host galaxies (at more than the 3 level).
• Likewise, all QSOs cannot correspond only to ELG host galaxies (at the ∼2 level).
• Given that the spline model works and that the parameter on is flat, the error function is generally a good model to describe QSO ( halo ).
The discussion of how environmental influences or assembly bias affects the QSO population is left for future studies. In the broadest sense, this self-sufficient study, in which the internally observed and measured correlation functions constrain the characteristics of the native halo catalogue, provides a much more self-consistent picture to the nuanced world of quasar occupation: (i) This likelihood-driven study grants a fair chance to multiple characteristic dependencies to arise from the measurement. The inferred indifference of the satellite fraction to the host halo mass in this picture, distinguishes itself from the previous studies and reported discrepancies where the partially outsourced measurement, halo model or adapted halo mass distribution, invite a host of added assumptions and justifications. See, e.g., Eftekharzadeh et al. (2019) and Kayo & Oguri (2012) for similarly adapted functional forms for the halo profile and halo mass distribution model, and yet significantly different satellite fraction that was attributed to a plausible luminosity dependency of the small scale clustering as opposed to lack thereof on large scales. (ii) This study highlights the notion that sampling quasars that belong to a group or are otherwise members of a close-pair system (a.k.a. 'binary'), leads to diverging conclusions on satellite occupation (see sec 6.2). This could be viewed as further evidence for hierarchical growth studies that have found earlier formation times for close halo pairs compared to their distant counterparts (Sheth & Tormen 1999;Harker et al. 2006).
(iii) Halo mass measurement for quasars in early eBOSS data inferred two plausible scenarios for the relatively constant characteristic halo mass between ∼ 1 and ∼ 2 ( Figure 10 in Laurent et al. (2017)) to be either due to massive haloes dominating the average or less luminous quasars inhabiting a wide range of halo masses and therefore putting the moderately luminous quasars in a different evolutionary path than the ∼ 4000 highly luminous quasars sampled by Shen et al. (2007) with a dramatically higher average halo mass. The inferred broad halo mass range in this study provides an elaborate case for the latter scenario using quasars in the same parent sample and luminosity class.
As the next generation surveys including eROSITA (Merloni et al. 2012 Takada et al. 2014) now come online, we will be able to test our findings and fully investigate the host galaxy population of luminous AGN across cosmic history.