Photometric Objects Around Cosmic Webs (PAC). VI. High Satellite Fraction of Quasars

The Photometric objects Around Cosmic webs (PAC) approach developed in Xu et al. has the advantage of making full use of spectroscopic and deeper photometric surveys. With the merits of PAC, the excess surface density n¯2wp of neighboring galaxies can be measured down to stellar mass 1010.80 M ⊙ around quasars at redshift 0.8 < z s < 1.0, with the data from the Sloan Digital Sky Survey IV extended Baryon Oscillation Spectroscopic Survey and the Dark Energy Spectroscopic Instrument Legacy Imaging Surveys. We find that n¯2wp generally increases quite steeply with the decrease of the separation. Using the subhalo abundance-matching method, we can accurately model the n¯2wp both on small and large scales. We show that the steep increase of n¯2wp toward the quasars requires that a large fraction fsate=0.29−0.06+0.05 of quasars should be satellites in massive halos, and we find that this fraction measurement is insensitive to the assumptions of our modeling. This high satellite fraction indicates that the subhalos have nearly the same probability of hosting quasars as the halos for the same (infall) halo mass and that the large-scale environment has negligible effect on the quasar activity. We show that even with this high satellite fraction, each massive halo on average does not host more than one satellite quasar, due to the sparsity of quasars.


INTRODUCTION
The large-scale structure of the universe is primarily influenced by dark matter that only undergoes gravitational interactions.Various luminous tracers, such as luminous red galaxies (LRGs), quasars, and emission line galaxies, are employed to understand the distribution of the underlying dark matter.Quasars are the most luminous tracer among them and can reach the high redshift (z s1 > 2.0), offering invaluable information for investi-gating structure formation and cosmology in the early universe (Neveux et al. 2020;du Mas des Bourboux et al. 2020;Hou et al. 2021;Brieden et al. 2022).Moreover, quasars are believed to be powered by the mass accretion onto supermassive black holes (SMBHs) located at the center of each galaxy (Soltan 1982).The strong feedback from the SMBHs regulates the formation of their host galaxies, shaping a process of co-evolution (Heckman & Best 2014).This underscores the importance of understanding quasars in the broader context of galaxy formation.
In the realm of both cosmological and galaxy formation studies employing quasars, it is crucial to accurately delineate the connection between quasars and the underlying dark matter, which often involves precise measurements of the number density and multi-scale clustering of quasars.Since the pioneering work of Osmer (1981), extensive efforts have been dedicated to quantifying quasar clustering through correlation functions.Among these, the two-point correlation function (2PCF) stands out as the predominant method.The availability of datasets such as the 2dF Quasi-Stellar Object Redshift Survey (2QZ) (Croom et al. 2004), SDSS-I/II/III/BOSS, and SDSS-IV/eBOSS (Schneider et al. 2002(Schneider et al. , 2010;;Eisenstein et al. 2011;Lyke et al. 2020), along with the emergence of surveys like the Dark Energy Spectroscopic Instrument (DESI) (DESI Collaboration et al. 2016), Subaru Prime Focus Spectrograph (PFS) survey (Tamura et al. 2018), and Euclid (Laureijs et al. 2011), has enabled high-precision measurements of the 2PCF at large scales (> 1 h −1 Mpc).However, accurately measuring quasar clustering at small scales (< 1 h −1 Mpc) remains challenging due to the technical issue of fiber collisions and the low quasar number density.While several up-weighting schemes have been proposed to address the fiber collision problem (Jing et al. 1998;Guo et al. 2012;Reid et al. 2014;Hahn et al. 2017;Percival & Bianchi 2017;Yang et al. 2019;Sunayama et al. 2020), and additional observations have been conducted to increase the number of quasar pairs selected from 2QZ and SDSS images (Hennawi et al. 2006;Myers et al. 2008;Richardson et al. 2012;Eftekharzadeh et al. 2017), the inherent difficulties persist.Moreover, numerous studies seek to enhance small-scale measurements by cross-correlating quasars with other spectroscopic tracers, such as LRGs (Kauffmann & Haehnelt 2002;Padmanabhan et al. 2009;Miyaji et al. 2011;Zhang et al. 2013;Shen et al. 2013).Although smallscale measurements can be more precise than relying solely on quasar auto-correlation, the obstacles remain as the the other tracers are also sparse at high redshift and the fiber collision issue may persist.
On the modeling front, various studies aim to establish connections between observed quasars and the underlying dark matter.In the simplest scenario, the linear bias factor can be determined by examining the relative amplitude of dark matter and quasars in the 2PCF at large scales (Jing 1998).This linear bias factor, in turn, can be employed to provide a rough estimate for the typical mass of the dark matter halos that host quasars (Ross et al. 2009).The halo occupation distribution (HOD) method has been utilized to model the 2PCF for an extensive range.Numerous studies have consistently found that the typical mass of quasars falls within the range of 10 12 h −1 M ⊙ ∼ 10 13 h −1 M ⊙ , showing a remarkable degree of independence with respect to redshift or quasar luminosity (Porciani et al. 2004;Croom et al. 2005;Myers et al. 2006;Coil et al. 2007;Myers et al. 2007;Shen et al. 2007;da Ângela et al. 2008;Padmanabhan et al. 2009;Richardson et al. 2012;White et al. 2012;Shen et al. 2013;Mitra et al. 2018;Powell et al. 2020).
The satellite fraction of quasars can be derived using the HOD framework, and this measurement proves particularly sensitive to small scale clustering.As noted previously, accurately measuring the small-scale clustering of quasars poses a challenge, leading to discrepant results in the literature (Starikova et al. 2011;Kayo & Oguri 2012;Leauthaud et al. 2015;Eftekharzadeh et al. 2019;Georgakakis et al. 2019;Alam et al. 2021).For instance, the satellite fraction of quasars is determined as 0.068 +0.034 −0.023 and 0.099 +0.046 −0.036 through the crosscorrelation between quasars and LRGs under two distinct quasar HOD models at redshift 0.3−0.9(Shen et al. 2013).Additionally, it is measured to be 7.3 +0.6  −1.5 × 10 −4 using the auto-correlation of quasars at redshift 0.4−2.5 (Richardson et al. 2012).However, it is important to determine this fraction, as the information is crucial not only for determining the small scale clustering of quasars (e.g., the close pairs) but also for studying the halo environment effect on the formation of quasars.
In this study, to achieve precise measurements of quasar clustering across various scales and accurately constrain the quasar-halo connection and satellite fraction, we utilize the Photometric Objects Around Cosmic Webs (PAC) approach.On the basis of Wang et al. (2011), Xu et al. (2022b) developed the PAC method, which capitalizes on the benefits of leveraging both spectroscopic and deeper photometric surveys.PAC can measure the excess surface density n2 w p of photometric objects with specific physical properties around spectroscopic tracers.With the PAC method, Xu et al. (2022a) and Xu et al. (2023) precisely determined the stellar-halo mass relation (SHMR) and galaxy stellar mass function (GSMF) down to the stellar mass limit where the spectroscopic sample is already highly incomplete, emphasizing the effectiveness of the PAC method.Given the absence of fiber collision problems between photometric and spectroscopic surveys, PAC becomes a valuable tool for accurately measuring the environment on h −1 Mpc scale for quasars.Additionally, PAC enables the assessment of the spatial distribution of galaxies with diverse properties around quasars, offering extensive information to study the quasar-halo connection.We use the quasar sample from the Sixteenth Data Release (DR16) in SDSS-IV/eBOSS (Prakash et al. 2016;Ahumada et al. 2020;Lyke et al. 2020) and the photometric sample from the Ninth Data release (DR9) in DESI Legacy Imaging Surveys (Legacy Surveys, Dey et al. 2019).We concentrate on a limited redshift range of 0.8 < z s < 1.0 due to constraints posed by the sur-vey depth of the Legacy Surveys.We employ the subhalo abundance matching method (SHAM) to model the n2 w p measurements from PAC and simultaneously constrain the SHMR for galaxies and the quasar-halo connection.Our analysis of the quasar-halo connection reveals a substantial satellite fraction f sate = 0.29 +0.05 −0.06 .This satellite fraction implies that a subhalo has nearly the same chance to host a quasar as a halo for the same (infall) mass, and quasars are more likely triggered by the galactic internal process instead of the large scale (halo) environment.
The paper is structured as follows: Section 2 provides a description of PAC and its measurements.Section 3 covers the N-body simulation, the subhalo abundance matching method description, and the results of modeling, along with corresponding discussions.Finally, Section 4 presents our conclusions.Throughout the paper, we adopt a spatially flat ΛCDM cosmology with Ω m,0 = 0.268, Ω Λ,0 = 0.732, and H 0 = 100 hkms −1 Mpc −1 = 71 kms −1 Mpc −1 (Hinshaw et al. 2013).

OBSERVATIONS AND MEASUREMENTS
In this section, we provide a brief overview of the PAC concept, and detail the spectroscopic and photometric samples utilized in this study.We assess the completeness of the photometric sample in terms of stellar mass.Subsequently, we present the measurement of n2 w p .

Photometric objects Around Cosmic webs (PAC)
Suppose we have two populations of objects, one from a spectroscopic catalog and the other from a photometric catalog.PAC can measure the excess surface density n2 w p of photometric objects with certain physical properties around spectroscopic objects across a wide range of scales: where n2 is the mean number density of photometric objects and S2 is the mean angular surface density of photometric objects.r 1 is the comoving distance to the spectroscopic objects.w p (r p ) and ω 12,weight (θ) are the projected cross-correlation function (PCCF) and the weighted angular cross-correlation function (ACCF) between spectroscopic objects and photometric objects, with r p = r 1 θ.To account for the case that r p varies with r 1 at fixed θ caused by the redshift distribution of the spectroscopic objects, ω 12 (θ) is weighted by 1/r2 1 .With PAC, we can take advantage of the deep photometric surveys to statistically obtain its rest-frame physical properties without making use of photo-z.We list the key steps of PAC as follows and refer to Xu et al. (2022b) for a detailed aacount: 1. Divide spectroscopic objects into several subsamples with narrower redshift bins to limit the range of r 1 in each bin.
2. Assign the median redshift in each redshift bin to the entire photometric sample.The physical properties of photometric objects can be estimated through spectral energy distribution (SED) fitting with the assigned redshift.Consequently, in each redshift bin, there is a physical property catalog for the photometric sample.
3. In each redshift bin, select photometric objects with specific physical properties and calculate n2 w p (r p ) according to Equation 1. Foreground and background objects with incorrect redshifts are effectively eliminated through ACCF, leaving only the photometric objects around spectroscopic objects with correct physical properties.
4. Combine the results from different redshift bins by averaging with appropriate weights.

Spectroscopic and Photometric Samples
For the photometric sample, we utilize the DR9 catalog from the Legacy Imaging Surveys (Dey et al. 2019), which includes the Dark Energy Camera Legacy Survey (DECaLS), the Mayall z-band Legacy Survey (MzLS), and the Beijing-Arizona Sky Survey (BASS) 2 .They cover approximately 14000 deg 2 in the Northern and Southern Galactic caps (NGC and SGC) with 9000 deg 2 at decl.≤ 32 • and 5100 deg 2 at decl.> 32 • .The sources are extracted using Tractor (Lang et al. 2016) and subsequently modeled with profiles convolved with a specific point spread function (PSF).These profiles include a delta function for point sources, and exponential law, de Vaucouleurs law, and a Sérsic profile for extended sources.We use their best-fit model magnitudes throughout the paper with the Galactic extinction corrected with the maps of Schlegel et al. (1998).We merely include the footprints observed at least once in g, r and z bands.To remove bright stars and bad pixels, we use MASKBITS3 offered by the Legacy Surveys.

Completeness and Designs
As in Xu et al. (2022b), we employ the z-band 10σ point source depth to assess the mass completeness of photometric objects.This analysis is conducted within the matched footprint of the Legacy Surveys and the eBOSS quasar sample.As presented in Figure 2, 90% of the regions covered by the Legacy Surveys are deeper than 22.34 in z band.Therefore, we take z = 22.34 as the galaxy depth for the Legacy Surveys.
We calculate galaxy physical properties using the SED code CIGALE (Boquien et al. 2019) with grz band fluxes.We adopt the stellar spectral library provided by Bruzual & Charlot (2003) to build up stellar population synthesis models.The initial stellar mass function (IMF) in Chabrier (2003) is assumed.We take Z/Z ⊙ = 0.4, 1.0, 2.5 as three different metallicities in our model.A delayed star formation history (SFH) ϕ(t) ≃ t exp(−t/τ ) is taken, with the timescale τ varies from 10 7 to 1.258×10 10 yr with an equal logarithm space of 0.1 dex.We apply the starburst reddening law of Calzetti et al. (2000) to consider the dust attenuation, in which the color excess E(B − V ) changes from 0 to 0.5.Following Xu et al. (2022b), we select the deepest 50 deg 2 regions in the Legacy Surveys comprising 738 bricks, where the z-band 10σ point source depth exceeds 23.37, to study the stellar mass completeness C 95 (M * ) with photometric redshift (z p ) calculated by Zhou et al. (2021).C 95 (M * ) is defined that 95% of the galaxies are brighter than C 95 (M * ) in z band for a given stellar mass M * .In Figure 3, we show the stellar mass-z band magnitude distribution of the deep photometric sample at redshift 0.8 < z p < 1.0.We mark C 95 (M * ) with red line in Figure 3, from which the complete stellar mass is 10 10.80 M ⊙ at redshift 0.8 < z s < 1.0 for the depth of z band 22.34.Hence, we use the stellar mass of photometric objects at the interval of [10 10.80 , 10 11.80 ] M ⊙ at redshift 0.8 < z s < 1.0 separated by four equal red- .The distribution of stellar mass versus z-band magnitude for DECaLS DR9 galaxies with photo-z between 0.8 and 1.0.The red line represents the faint z band magnitude limit so that 95% of galaxies at each stellar mass are brighter than the limit.The horizontal gray dashed line shows the survey depth in z-band determined in Figure 2. The vertical gray dashed line is the completeness limit in stellar mass 10 10.8 M⊙ for the photometric galaxies at the redshift of our interest.
shift bins.Furthermore, to constrain the SHMR at lower masses, we incorporate the incomplete stellar mass bins within the range of [10 10.00 , 10 10.80 ] M ⊙ through appropriate modeling, as illustrated later.

Measurements
We follow the approach in Xu et al. (2022a) to average the result at different sky regions and redshift bins with proper weights to acquire final n2 w p (r p ).
For representation simplicity, let A ≡ n2 w p (r p ). Supposing A is calculated for N r redshift bins and N s sky regions, we further divide each region into N sub subregions for error estimation.In our specific scenario, N r is set to 4, corresponding to the redshift bins [0.80,0.85],[0.85,0.90],[0.90,0.95]and [0.95,1.00].Additionally, N s is defined as 3, representing BASS+MzLS NGC, DE-CaLS SGC and DECaLS NGC while N sub is established as 10.According to Equation 1, A i,j,k can be measured in the ith redshift bin, jth sky region and kth sub-sample through the Landy-Szalay estimator (Landy & Szalay 1993).Firstly, we average the measurements at different sky regions using the region area w s,j as a weighting factor to take the sky area difference into account: . (2) Secondly, we obtain the mean values and the error vector σ i of the mean values for each redshift bin from N sub sub-samples: Finally, results from different redshift bins are averaged according to the error vector σ i .Let (5) We measure n2 w p (r p ) and the auto-correlation w p of quasars in the radial range of 0.1 h −1 Mpc < r p < 15 h −1 Mpc.These scales are sufficiently broad to capture the distributions of both centrals and satellites.
The measurements of n2 w p (r p ) are depicted as data points with error bars in Figure 4.The error bars represent the vector σ.These measurements are overall good across all stellar mass bins within the entire radial range, except for the highest stellar mass bin [10 11.6 , 10 11.8 ] M ⊙ .Additionally, the auto-correlation w p of quasars, with the weights assigned in the same way as Hou et al. (2021), is presented in Figure 5 to enhance the precision in constraining the mean host halo mass of quasars.

SIMULATION AND RESULTS
In this section, we outline the N-body simulation and the subhalo abundance matching (SHMR) method employed in this study for modeling the measurements.We present the best-fit results for the SHMR within the redshift range of 0.8 < z s < 1.0 and explore the quasar-halo connection.Additionally, we highlight the high satellite fraction of quasars from our findings and assess the assumptions underlying the model.

N-body Simulation
We employ the CosmicGrowth simulation suite (Jing 2019), a grid of high-resolution N-body simulations utilizing the computationally-efficient adaptive parallel P 3 M method (Jing & Suto 2002;Xu & Jing 2021).We utilize the ΛCDM simulation with the following cosmological parameters: Ω m,0 = 0.268, Ω Λ,0 = 0.732, h = 0.71, n s = 0.968, and σ 8 = 0.83.This simulation comprises 3072 3 dark matter particles within a 600 h −1 Mpc box, with a softening length of η = 0.01 h −1 Mpc.Groups are characterized with the friends-of-friends algorithm (Davis et al. 1985), employing a linking length of 0.2 times the mean particle separation.The halos are then processed with HBT+ (Han et al. 2012(Han et al. , 2018) ) to identify subhalos and trace their evolution histories.Merger timescales of the subhalos with fewer than 20 particles are estimated using the fitting formula in Jiang et al. (2008).Subhalos that have already merged into central subhalos are disregarded.The adopted simulation has a mass resolution of m p = 5.54 × 10 8 h −1 M ⊙ , which proves to be sufficiently fine for our study.We utilize snapshot 76 at a redshift of approximately 0.92 to align with the observations.

Subhalo Abundance Matching
To associate galaxies with (sub)halos in the N-body simulation, we implement the SHAM method.We utilize the widely-used five-parameter formula for the SHMR, a double power law with a constant scatter (Wang & Jing 2010;Yang et al. 2012;Moster et al. 2013;Xu et al. 2023): Here, M acc is defined as the virial mass M vir of the (sub)halo when it was last to become a central dominant object, where M vir is defined through the fitting formula in Bryan & Norman (1998).The scatter in log(M * ) at a given M acc is assumed to follow a Gaussian distribution with a width of σ.The same set of parameters is applied for both centrals and satellites.(Wang & Jing 2010;Behroozi et al. 2019;Xu et al. 2023).In addition to the standard SHAM parameters described above, given the inclusion of measurements from incomplete stellar mass bins, we introduce four additional parameters k 1 , k 2 , k 3 , k 4 to account for the incompleteness of the four stellar mass bins.The modeled results from the simulation are obtained by multiplying the n2 w p (r p ) from the entire sample by the incompleteness, under the assumption that the incompleteness only affects the amplitude of the observed n2 w p (r p ).
To associate quasars with (sub)halos, we assume that the probability that a (sub)halo hosts a quasar follows a Gaussian distribution of logarithmic halo mass log 10 [M acc /(h −1 M ⊙ )] with a mean value µ and a dispersion σ q .We use the same µ and σ q for both halos and subhalos, Further we introduce an additional parameter B to adjust the overall probability for a subhalo to host a quasar relative to a halo at the same mass.This ensures that the subhalos have B times probability to host a quasar relative to a halo of the same mass.Under the above schedule, the satellite fraction f s of quasars is defined as: where N sate and N cen represent the numbers of subhalos and halos selected based on the same Gaussian distribution.
After assigning galaxies and quasars to (sub)halos, we compute the correlation functions using Corrfunc (Sinha & Garrison 2020) in the simulation.To compare with the measurements in observation, we define the χ 2 as: σ 2 (r j ) .
Here, N m and N r denote the numbers of stellar mass bins and of radial bins respectively.We explore the parameter space using the Markov Chain Monte Carlo (MCMC) sampler emcee (Foreman-Mackey et al. 2013), employing maximum likelihood analyses for the three sets of parameters {M 0 , α, β, k, σ}, {k 1 , k 2 , k 3 , k 4 }, and {µ, σ q , B}.

SHMR and quasar-halo connection
The best-fit results and errors for the parameters are presented in the first row of Table 1.Additionally, the joint posterior distributions of the parameters are illustrated in Figure A1 using corner (Foreman-Mackey 2016).The best-fit n2 w p (r p ) and w p results from the model are shown in Figure 4 and Figure 5 with solid lines.The fits are generally good for all stellar mass bins, both on the small and large scales.
For quasars, we find that µ, σ q , and B are 13.01 +0.20 −0.13 , 0.65 +0.11  −0.12 , and 1.24 +0.26 −0.22 , respectively.The results indicate a very high satellite fraction of quasars, with f sate = 0.29 +0.05 −0.06 .B is equal to 1 within the 1σ error, implying that subhalos have nearly the same chance to host quasars as the halos, and the quasar activity is not affected by the larger halo environment.After determining the host halo masses for subhalos and adjusting the number density of quasars to match the observed value of approximately 1.29 × 10 −5 h 3 Mpc −3 , we illustrate the HOD of our quasar sample in Figure 6.We obtain median host halo masses for central and satellite quasars of log 10 [M med h,cen /(h −1 M ⊙ )] = 12.05 +0.60 −0.60 and log 10 [M med h,sat /(h −1 M ⊙ )] = 12.90 +0.68 −0.72 .From the HOD, we observe that although the satellite fraction of quasars looks very high, each massive halo, on average, does not host more than one satellite quasar due to the low total number density of quasars.To more effectively depict the mass distributions, we also present the host halo mass distributions for central and satellite quasars in Figure B1.
For galaxies, we present the SHMR at 0.8 < z s < 1.0 in Figure 7  quasar samples and deeper photometric catalogs are required to more precisely constrain the low mass end (e.g., β in SHMR) and to explore higher redshifts.
Moreover, we determine the incompleteness of photometric samples for the four lowest stellar mass bins, yielding k 1 = 0.43 +0.03 −0.05 , k 2 = 0.62 +0.05 −0.05 , k 3 = 0.70 +0.04 −0.05 and k 4 = 0.83 +0.04 −0.05 .These results highlight the capability of the PAC method to constrain the stellar mass incompleteness of photometric samples and utilize the incomplete sample to explore lower mass regions.

Testing the robustness of the results
The satellite fraction of quasars in our results is f sate = 0.29 +0.05 −0.06 .This high satellite fraction suggests that subhalos have nearly the same probability (B = 1.24 +0.26  −0.22 ) as halos to host quasars for the same host (infall) mass, and the large-scale environment has little effect on quasar activity.To verify the robustness of such a high satellite fraction of quasars, we conducted three tests as follows.Since BASS is a bit shallower in g and r bands than the DECaLs, we measure n2 w p (r p ) in the DECaLS and BASS respectively, and a compari- son of the results shows that the measurement is robust (Figure C1).Initially, we observe a degeneracy between the parameters µ and σ q in Figure A1.To investigate the impact of this degeneracy on the satellite fraction of quasars, we fix σ q at 0.50, a value consistent with Richardson et al. (2012), and constrain the remaining parameters through MCMC analysis.The resulting best-fit parameters are presented in the second row of Table 1, and we find that a high satellite fraction is still necessary to align with the observations, and all parameters remain consistent with the previous results within the 1σ interval.For better clarity, we depict the n2 w p results under these best-fit parameters for two stellar mass bins (10 10.5 M ⊙ and 10 10.9 M ⊙ ) in Figure 8 using solid blue lines.We find that the measurements can also be well fitted by this model.This test confirms that the satellite fraction of quasars remains unaffected by the degeneracy between µ and σ q .
Secondly, to examine whether a low satellite fraction of quasars, as suggested by some previous studies (Richardson et al. 2012;Shen et al. 2013), is consistent with our measurements, we set the B parameter to be 0.16, corresponding to f sate = 0.05 according to Equation 8, while keeping the remaining parameters unchanged.The simulation results in two stellar mass bins (10 10.5 M ⊙ and 10 10.9 M ⊙ ) are presented in Figure 8 with solid orange lines.In this case, we observe that the model fails to match the steep increase of n2 w p at small scales (r p < 1 h −1 Mpc), suggesting that a higher Table 1.The best-fit parameter results and errors for the SHMR, incompleteness, and quasar-halo connection, assuming a Gaussian distribution for the quasar probability.The first row displays the fiducial results with all parameters free, while the second row presents the results for the remaining parameters after setting σq to 0.5.M0 is in unit of h −1 M⊙ and k is in M⊙. ) photo 10 10.5 M SHAM(gauss) q = 0.50 SHAM(gauss) B = 0.16 SHAM(erf) 10 1 10 0 10 1 r p (h 1 Mpc) photo 10 10.9 M Figure 8.The n2wp measurements in two stellar mass bins, 10 10.5 M⊙ and 10 10.9 M⊙.Dots with error bars represent the observational measurements.The orange solid lines depict the best-fit results with fsate set to 0.05 (equivalently B = 0.16), while keeping other parameters the same as in the first row of Table 1.The blue solid lines show the best-fit results with σq set to 0.50.The green solid lines represent the best-fit results assuming a probability of a (sub)halo hosting a quasar modeled by an error function.
satellite fraction is required to reproduce the observation.
Thirdly, considering our assumption of a Gaussian format for the quasar probability in our modeling, we aim to verify whether this assumption affects the satellite fraction of quasars.Therefore, we replace the Gaussian format with the error function format for both central and satellite quasars shown as follows: ) .
(10) where P (M h ) is the probability of a halo becoming a quasar.log 10 (M min ) and σ log 10 (M h ) are the characteristic mass scale and transition width of a softened step function.The best-fit parameters are presented in Table 2, where we still find a high satellite fraction of quasars indicated by the parameter B = 1.22 +0.26  −0.18 with f sate = 0.28 +0.05 −0.04 .For a clearer illustration, simulation results in two stellar mass bins (10 10.5 M ⊙ and 10 10.9 M ⊙ ) are shown in Figure 8 with solid green lines.We observe that the measurements can also be wellfitted by this model.To compare the difference between the Gaussian format and the error function, we show the HOD from the error function model in Figure 6.We observe that HODs for satellites are nearly identical for these two models.Although some discrepancy is found for centrals with M h > 10 13.0 M ⊙ , this does not significantly impact f sate due to the low number densities of these massive halos.This test confirms that our results are not sensitive to the assumptions in the quasar-halo connection.
Our higher satellite fraction of quasars than that determined by Shen et al. (2013) may come from two sources.One is the measurements of the cross correlation between galaxies and quasars.In this paper, we have used galaxies with stellar mass larger than 10 10.0 M ⊙ , compared with LRGs used by Shen et al. (2013).While we find the slopes of the cross correlation functions within r p = 1 h −1 Mpc are consistent for massive galaxies (i.e.> 10 11.2 M ⊙ ) between the two works, our measurement for less massive galaxies also requires more satellite quasars around small central galaxies (cf.left panel of Figure 8).The other source is different assumptions for the modeling.We use the abundance matching method, while Shen et al. (2013) used the HOD method.We note that the satellite fraction is quite sensitive to the HOD forms used for the central quasars in Shen et al. (2013).It changes from f sat = 0.068 +0.034  −0.023 when the error function form is used to f sat = 0.099 +0.046  −0.036 when the lognormal form is used instead.Furthermore they require that satellite quasars exist only in halos with M h > 10 13.0 h −1 M ⊙ (cf their Figure 14).In comparison, we allow many more small halos with M h < 10 13.0 h −1 M ⊙ (Figure B1) to host satellite quasars which are required to reproduce the cross correlation between quasars and galaxies with stellar mass less than 10 11.0 M ⊙ (e.g.left panel of Figure 8).Therefore, our higher satellite fraction mainly comes from those satellite quasars in halos with M h < 10 13.0 h −1 M ⊙ that were not probed by Shen et al. (2013).Because central halos as small as M h > 10 11.5 h −1 M ⊙ can host quasars in the both works (cf.their Figure 14 and our Figure 6), it is more reasonable that halos with M h < 10 13.0 h −1 M ⊙ can host satellite quasars, as there are subhalos that are massive enough.

CONCLUSIONS
In this paper, we utilize a spectroscopic quasar sample from SDSS-IV eBOSS DR16 and a photometric galaxy sample from the Legacy Surveys.We employ the PAC method to measure n2 w p (r p ) between quasars and galaxies, in conjunction with the quasar auto-correlation w p , at redshift 0.8 < z s < 1.0.Leveraging the advantages of PAC, we obtain reliable quasar clustering down to 0.1 h −1 Mpc.We model the measurements in N-body simulations using the SHAM approach and assume a Gaussian probability for a (sub)halo to host a quasar.We constrain the SHMR of galaxies and the quasar-halo connection.We verify the assumptions and confirm that a high satellite fraction of quasars is required to reproduce the observation.The main results are listed as follows: 1.Under the assumption that the probability of a halo becoming a quasar follows a Gaussian distribution of logarithmic halo mass log 10 [M h /(h −1 M ⊙ )], we find that the median host halo masses for central and satellite quasars are log 10 [M med h,cen /(h −1 M ⊙ )] = 12.05 +0.60 −0.60 and log 10 [M med h,sat /(h −1 M ⊙ )] = 12.90 +0.68 −0.72 at redshift 0.8 < z s < 1.0 and a high satellite fraction of quasars of f sate = 0.29 +0.05 −0.06 .This high satellite fraction of quasars indicates that subhalos have nearly the same probability to host quasars as the halos for the same host (infall) mass, and the largescale environment has little effect on the quasar activity.
2. We constrain the SHMR of galaxies at redshift 0.8 < z s < 1.0, which can be described by a double power law with log 10 [M 0 /(h −1 M ⊙ )] = 11.83

Figure 1 .
Figure 1.Spatial coverage of the two samples used in this work.The photometric sample is the Legacy Imaging Surveys DR9 where the DECaLS part is shown in blue color and the BASS+MzLS part is represented by dark orange color.The spectroscopic sample is SDSS-IV eBOSS DR16 quasars shown in blue color.

Figure 2 .
Figure 2. Cumulative distribution function of survey area to the z band 10 σ point source depth in the Legacy Surveys matched with eBOSS quasar footprint.
Figure3.The distribution of stellar mass versus z-band magnitude for DECaLS DR9 galaxies with photo-z between 0.8 and 1.0.The red line represents the faint z band magnitude limit so that 95% of galaxies at each stellar mass are brighter than the limit.The horizontal gray dashed line shows the survey depth in z-band determined in Figure2.The vertical gray dashed line is the completeness limit in stellar mass 10 10.8 M⊙ for the photometric galaxies at the redshift of our interest.

Figure 5 .
Figure 5.The measurement (data points with error bars) and the model fitting (line) of the projected auto-correlation wp of quasars.

Figure 6 .
Figure6.The mean HODs of quasars computed with the best-fit parameters.The results based on the Gaussian function for the quasar-halo relation are displayed in blue, with the parameters from the first row of Table1.For comparison, the results based on the error function are shown in orange, with the parameters from Table2.Central and satellite quasars are depicted by dotted and dot-dashed lines, respectively.

Figure 7 .
Figure 7.The stellar-halo mass relation at redshift 0.8 < zs < 1.0.The mean SHMR and 1σ errors are shown with blue line and shadow.

APPENDIXA.
POSTERIOR DISTRIBUTIONS OF THE PARAMETERSWe present the posterior PDFs of the 12 parameters in our fiducial model in this section, as shown in FigureA1.B. HOST HALO MASS DISTRIBUTIONS OF QUASARSWe present the mass distributions of host halos of central and satellite quasars in FigureB1.C. COMPARISON OF N2 W P (R P ) MEASUREMENTSWe compare the n2 w p (r p ) measurements with DR9 catalog from DECaLS and from BASS+MzLS that are shown in FigureC1.Gui et al.

Figure A1 .Figure B1 .
Figure A1.The posterior distributions of the 12 parameters in our model obtained through MCMC.The vertical green line indicates the median value of each parameter, and the dashed blue lines denote the 16th and 84th percentiles after marginalizing the parameters.

Table 2 .
The best-fit parameter results and errors for the SHMR, incompleteness, and quasar-halo connection, assuming an error function for the quasar probability.M0 and Mmin are in units of h −1 M⊙ and k is in M⊙.
No. B20019.We gratefully acknowledge the support of the Key Laboratory for Particle Physics, Astrophysics and Cosmology, Ministry of Education.This work made use of the Gravity Supercomputer at the Department of Astronomy, Shanghai Jiao Tong University.Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions.SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah.The SDSS website is www.sdss.org.SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional/MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.The Legacy Surveys consist of three individual and complementary projects: the Dark Energy Camera Legacy Survey (DECaLS; Proposal ID #2014B-0404; PIs: David Schlegel and Arjun Dey), the Beijing-Arizona Sky Survey (BASS; NOAO Prop.ID #2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band Legacy Survey (MzLS; Prop.ID #2016A-0453; PI: Arjun Dey).DECaLS, BASS and MzLS together include data obtained, respectively, at the Blanco telescope, Cerro Tololo Inter-American Observatory, NSF's NOIRLab; the Bok telescope, Steward Observatory, University of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOIRLab.The Legacy Surveys project is honored to be permitted to conduct astronomical research on Iolkam Du'ag (Kitt Peak), a mountain with particular significance to the Tohono O'odham Nation. project