Photometric Objects around Cosmic Webs (PAC) Delineated in a Spectroscopic Survey. I. Methods

We provide a method for estimating the projected density distribution $\bar{n}_2w_p(r_p)$ of photometric objects around spectroscopic objects in a redshift survey. This quantity describes the distribution of Photometric sources with certain physical properties (e.g. luminosity, mass, color etc) Around Cosmic webs (PAC) traced by the spectroscopic objects. The method can make full use of current and future deep and wide photometric surveys to explore the formation of galaxies up to medium redshift ($z_s<2$), with the aid of cosmological redshift surveys that sample only a fairly limited species of objects (e.g. Emission Line Galaxies). As an example, we apply the PAC method to the CMASS spectroscopic and HSC-SSP PDR2 photometric samples to explore the distribution of galaxies for a wide range of stellar mass from $10^{9.0}{\rm M_\odot}$ to $10^{12.0}{\rm M_\odot}$ around massive ones at $z_s\approx 0.6$. Using the abundance matching method, we model $\bar{n}_2w_p(r_p)$ in N-body simulation using MCMC sampling, and accurately measure the stellar-halo mass relation (SHMR) and stellar mass function (SMF) for the whole mass range. We can also measure the conditional stellar mass function (CSMF) of satellites for central galaxies of different mass. The PAC method has many potential applications for studying the evolution of galaxies.


INTRODUCTION
Modeling galaxy formation in the cosmological context is one of the greatest challenges in astrophysics and cosmology today. In past few decades, the broad contours of galaxy formation physics are investigated and built (Mo et al. 2010;Silk & Mamon 2012;Somerville & Davé 2015). In theory, cosmic structures grow by gravitational instability from the initial tiny quantum fluctuations generated during the inflationary epoch. Dark matter (DM) halos, which are defined as dark matter objects with the density hundreds times the density of the Universe, are formed around the density peaks of the Universe. Halos grow by accreting surrounding smaller DM halos and diffuse matter including DM and gas (Zhao et al. 2009). The gas is heated by shocks when accreted into a halo (Dekel & Birnboim 2006). Through radiation, the hot gas loses its energy and cools, and the cold gas spirals into the halo center to form a galaxy, which is typically a spiral. Under the hierarchical growth of the DM halo, surrounding galaxies can also be accreted into the halo and become its satellites. Some satellite galaxies, especially relatively massive ones, spiral into the center and coalesce with the central galaxy due to dynamical friction (Lacey & Cole 1994;Jiang et al. 2008). A merger of a central disk galaxy with a satellite of comparable mass (say, a mass ratio ≥ 0.2) may significantly change its internal structure, either producing an elliptical galaxy or a disk galaxy with a significant bulge (Naab et al. 2006;Cox et al. 2006;Bekki & Couch 2011). Studies find that black hole (BH) mass is highly correlated with bulge mass (Kormendy & Ho 2013), so we expect there are supermassive black holes (SMBHs) within elliptical galaxies or bulges. The strong energy output and/or material outflow produced by an SMBH can heat and/or blow out the surrounding cold gas, thus suppressing the gas reservoir from which stars form and making the host galaxy look red (Fabian 2012). Therefore, the bimodal distribution of galaxy color is highly correlated with the galaxy morphology distribution, with ellipticals being red and spirals being blue.
To fully understand the galaxy properties and distributions, one needs to consider all the complicated physical process involved in galaxy formation and evolution mentioned above. Although there exists a broad contour of the galaxy formation process, a fully predictive framework from first principles has yet to be established (Naab & Ostriker 2017). Physical models such as semianalytic models (White & Frenk 1991;Kauffmann et al. 1993;Kang et al. 2005;Guo et al. 2013) and hydrodynamic simulations (Schaye et al. 2015;Pillepich et al. 2018;Davé et al. 2019) approximate physics below their respective resolution scales to simulate the effects of supernovae, radiation pressure, multi-phase gas, black hole accretion, active galactic nucleus (AGN) feedback and metallicity evolution in galaxy formation. However, different approximations lead to different galaxy properties (Somerville et al. 2008;Lu et al. 2014) and the uncertainty still remains.
Empirical modeling such as halo occupation distribution (HOD, Jing et al. (1998); Peacock & Smith (2000); Ma & Fry (2000); Seljak (2000); Berlind & Weinberg (2002); Yang et al. (2003); Zheng et al. (2005); Zu & Mandelbaum (2015)) and abundance matching (AM, Wechsler et al. (1998); Wang & Jing (2010); Moster et al. (2013); Behroozi et al. (2019)) uses significantly weaker prior and the physical constraints come almost entirely from observations. These models connect the average galaxy properties such as occupation numbers and stellar mass to halos as a function of halo mass or halo circular velocity, and determine the parameters by fitting the observed properties of galaxies. Although more complex models have been attempted recently to incorporate properties of gas (Popping et al. 2015), metallicity (Rodríguez-Puebla et al. 2016) and dust (Imara et al. 2018) to compare with observations, the stellarhalo mass relation (SHMR) is still one of the most commonly used relations to model galaxy-halo connection (Guo et al. 2010;Wang & Jing 2010;Moster et al. 2013), in which larger halos host larger galaxies with a relatively tight scatter. Recent studies found that galaxies with different properties, such as color, may have differ-ent SHMR, which may indicate the existence of so called galaxy assembly bias that causes the scatter in the average SHMR (Cooper et al. 2010;Wang et al. 2013;Zentner et al. 2014;Hearin et al. 2015;Mandelbaum et al. 2016;Cui et al. 2021). Similarly, for HOD, the galaxy occupation number may also depend on halo properties such as halo concentration and environment other than the halo mass (Hadzhiyska et al. 2020). The best secondary parameter for the halo characterization is still being sought in the development of HOD and AM .
Stellar mass function (SMF) and galaxy clustering (GC) are the two most commonly used properties to constrain the parameters in HOD and AM. The measurement of these quantities usually relies on spectroscopic surveys with redshift information. In the past two decades, there has been significant progress in large spectroscopic surveys (York et al. 2000;Colless et al. 2001;Steidel et al. 2003;Le Fèvre et al. 2005;Ahn et al. 2012;Bolton et al. 2012;Garilli et al. 2014;Takada et al. 2014;DESI Collaboration et al. 2016;Ahumada et al. 2020). In the local universe (z s ∼ 0), large spectroscopic surveys, in particular the Two Degree Field Galaxy spectroscopic survey (2dFGRS; Colless et al. (2001)) and the Sloan Digital Sky Survey (SDSS; York et al. (2000)), have been used to measure the SMF and GC down to 10 9.0 M (Cole et al. 2001;Norberg et al. 2002;Li et al. 2006;Li & White 2009;Peng et al. 2010), although the accuracy of the measurements, especially the GC, is still very limited by the survey volume for faint galaxies. At a higher redshift (z s = 0.5 ∼ 1.0), the DEEP2 Galaxy spectroscopic survey (Davis et al. 2003), the VIMOS-VLT Deep Survey (VVDS; Le Fèvre et al. (2005)) and the VIMOS Public Extragalactic spectroscopic survey (VIPERS; Garilli et al. (2014)) have been used to successfully measure SMF and GC for galaxies with M * > 10 10.0 M (Pozzetti et al. 2007;Meneux et al. 2008;Mostek et al. 2013;Marulli et al. 2013;Davidzon et al. 2013). However, despite the huge efforts, measurement for fainter objects is still very difficult, and stellar mass limited samples are usually very small at even higher redshift. Fortunately, huge next generation spectroscopic surveys are being constructed for cosmological studies at intermediate to high redshift (Laureijs et al. 2011;Takada et al. 2014;DESI Collaboration et al. 2016). Due to limited wavelength coverage and sensitivity of the spectrographs, different populations of galaxies, such as Emission Line Galaxy (ELG), QSOs and Lyman Break Galaxy (LBG), are targeted at different redshifts. These populations are all expected to trace large-scale structures or the Cosmic Webs, so they can be used to extract information for cosmological studies. However, it is difficult to use these surveys to con-struct stellar mass limited samples for the target selections used by the surveys.
Compared to spectroscopic surveys, photometric surveys, which take the images and obtain the photometric information, are usually deeper and more complete in terms of the stellar mass. However, without precise redshift measurement, the usefulness of photometric surveys is limited. People attempt to infer the photometric redshift from their broad band magnitudes in the photometric surveys (Ilbert et al. 2006(Ilbert et al. , 2009Sánchez et al. 2014;Nishizawa et al. 2020). Photo-z has been used to measure SMF (Fontana et al. 2006;Pérez-González et al. 2008;Ilbert et al. 2010;Bielby et al. 2012;McLeod et al. 2021) and GC (Crocce et al. 2016;Ishikawa et al. 2020;Wang et al. 2021b), though one has to be very careful about the error and systematics of photo-z especially for faint galaxies. There will be the next generation large and deep multi-band photometric surveys, such as Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST; Ivezić et al. (2019)) and Euclid (Laureijs et al. 2011), which are expected to fairly sample galaxies to an unprecedented faint limit.
As mentioned above, photometric and spectroscopic surveys both have their advantages and disadvantages. Future large spectroscopic surveys have the precise redshift information but only for bright objects and/or selected (or biased) populations, while photometric surveys are deeper and more complete but without accurate redshift measurement. So far, the measurement of SMF and GC is mainly from the spectroscopic surveys. Thus, studies based on SMF and GC, such as HOD and AM, are focused on the local universe or the massive end at higher redshift, resulting in a poor understanding of the faint end and the redshift evolution. To study small and faint objects, quite a few studies attempted to combine a spectroscopic survey with a photometric survey, which can reach several magnitudes fainter than pure spectroscopic surveys. With the Cosmic Webs traced by objects in spectroscopic surveys, properties and distribution of photometric objects around the Cosmic Web can be studied using the photometric surveys. Results can be achieved by stacking satellite and neighbor counts around a large sample of central galaxies, with foreground and background sources subtracted statistically (Phillipps & Shanks 1987;Lorrimer et al. 1994;Wang et al. 2011;Guo et al. 2012;Ménard et al. 2013;Newman et al. 2015;Lan et al. 2016). However, most previous studies only use luminosity and color properties of a photometric sample, which is not suitable for quantitative HOD and AM studies, which requires physical properties of galaxies such as stellar mass, rest-frame color and star formation rate (SFR).
In this paper, we provide a method to measure the distributions and properties (luminosity, mass, color, SFR, morphology et.al.) of Photometric objects Around Cosmic webs (PAC) represented by spectroscopic objects in a spectroscopic survey. The basic idea is that for spectroscopic source i at redshift z s,i , only those objects in the photometric sample at around z s,i are correlated the source i and share the similar redshift. Thus for the source i, we calculate the physical properties for all sources in the whole photometric sample by assuming they were all at the redshift z s,i . Through the crosscorrelation of the photometric and spectroscopic samples, foreground or background galaxies with wrong redshift information can be canceled out with the help of random samples, and the true distribution of photometric sources with specified properties around the spectroscopic sources can be obtained. Because both the spectroscopic and photometric samples are huge in the next generation surveys, we will develop a method to speed up the computation for the physical properties.
We introduce the details of PAC in Section 2. In Section 3, we apply PAC to observations. The measurement is modeled in N-body simulation using AM in Section 4. And a brief conclusion is given in Section 5. We adopt the cosmology with Ω m = 0.268, Ω Λ = 0.732 and H 0 = 71 km/s/Mpc through out the paper.

METHODOLOGY
In this section, we introduce a method for estimatingn 2 w p (r p ) from w 12 (ϑ), where w p (r p ) and w 12 (ϑ) are the projected cross-correlation function (PCCF) and the angular cross-correlation function (ACCF) between a given set of spectroscopically identified galaxies and a large sample of photometric galaxies, andn 2 is the mean number density of the photometric galaxies. The quantityn 2 w p (r p ) has clear physical meaning that it measures the true excess of the photometric objects around the spectroscopic objects on the sky projection. If we choose photometric galaxies at a given stellar mass, the measurement over a range of the stellar mass gives the information on the stellar mass function and the clustering as a function of the stellar mass, which are the key ingredients for understanding the connection of galaxies to dark matter halos. Therefore, we extend this method to statistically measuring the distribution of the photometric galaxies with specified physical properties (i.e stellar mass, SFR, color etc) around spectroscopically identified galaxies.
2.1. Estimatingn 2 w p (r p ) from w 12 (ϑ) Throughout this section, we call a spectroscopic sample population 1 and a photometric sample population 2.
Assuming an object in population 1 is at distance r 1 , the number of objects dN 2 in population 2 within a solid angle element dΩ 2 in the direction r 2 is: Where n 2 is the mean number density of population 2 at the distance r 2 . The expected number of population 2 objects around a population 1 object is: S 2 is the mean angular surface density of population 2, ξ 12 is the cross-correlation function (CCF) between the two populations and w p (r p = r 1 θ) ≡ ξ 12 ( r 2 p + π 2 )dπ is the projected cross-correlation function. The approximation holds if θ is small. Then we have: Here D 1 D 2 and D 1 R are the cross pair counts between population 1 and population 2, and between population 1 and the random sample of population 2. Hence, with our method,n 2 w p (r 1 θ)r 2 1 can be estimated fromS 2 w 12 (ϑ) with only the redshift of population 1. The physical meaning of the quantityn 2 w p (r 1 θ) is the excess number of population 2 around population 1.

Physical properties of photometric sources around spectroscopically identified sources
To obtain physical properties for galaxies, distance or redshift information is needed. Unfortunately, wide deep photometric surveys do not have measured redshift for most of galaxies. Photometric redshift z p is often used to approximate the distance of galaxies, but the errors of z p are very difficult to estimate especially for faint galaxies. In our method we do not use any information of z p . Instead, since those photometric sources correlated with the spectroscopically identified galaxies must share the same redshift, we can use the spectroscopic redshift for the physical properties for photometric sources around the spectroscopic object. Asn 2 w p (r 1 θ) is just the number excess of neighbors around population 1, we can extend our method to estimating the distributions and properties of Photometric objects Around Cosmic webs (PAC) delineated by a spectroscopic survey.
Assuming we only have a spectroscopic sample (population 1) and a large photometric sample, and we want to calculaten 2 w p (r 1 θ) between those sources in population 1 at a single redshift z s and sources with physical property X included in the photometric sample. To select population 2, we assume the whole photometric sample is at redshift z s and calculate their physical properties using SED fitting. As one may note, the calculation is correct only for sources around population 1. The physical properties calculated are incorrect for foreground and background sources. However, since foreground and background sources along the line-of-sight (LOS) to population 1 are distributed statistically in the same way as those along a random LOS, the foreground and background will be canceled out when we calculate w 12 (ϑ) (see Wang et al. (2011, Sec 3.3.1) for a more rigorous derivation). Therefore, population 2 is just selected as the sources with physical property X in the photometric sample even if we assume the whole sample is at redshift z s in the calculation. With this method, we can study the distribution of satellites and neighbors with specified physical properties around spectroscopically identified galaxies.

Spectroscopic sample with a redshift distribution
Till now, we have assumed that all the sources in population 1 have the same redshift. However, in reality, spectroscopic samples all have a redshift distribution. If the redshift range is relatively narrow and the evolution of the universe can be neglected,n 2 w p (r p ) will not vary much while the change of r 1 and θ = r p /r 1 may not be negligible from one spectroscopic source to another. Thus,n 2 w p (r p ) is a better statistic quantity thann 2 w p (r p )r 2 1 . We re-derive the Equation 3: Here D 1 D 2 w and D 1 R w are the weighted cross pair counts between population 1 and population 2, and between population 1 and the random sample. Given a set of sources in population 1 with positions {r 1,i } i=1,N1 , in order to measure the PCCF at the projected separation r p , we can measure the number dN 2,i of population 2 neighbors within the angular separation θ ±1/2dθ around object i in population 1. Then we measure the weighted count: Similarly we can get the weighted count for random sample: In an alternative way, we can also change the angular separation used to measure the number counts and make r p = r 1,i θ i the same for each galaxy in population 1. Then, we sum the counts around each galaxies without weighting. These two methods are equivalent.
To minimize the effects of complex survey geometries, we adopt an estimator in analogy to the Landy-Szalay estimator for two-point auto-correlation (Landy & Szalay 1993): where R 1 and R 2 are the random points for spectroscopic and photometric samples respectively.
The above method has accounted for the variance of θ with redshift for sources in population 1, while r 1 in S 2 w 12,weight /r 2 1 still varies with redshift. Therefore, we divide population 1 into narrower redshift bins and reduced the error from the change of r 1 .
For population 1 with a redshift distribution, we can divided them into m redshift bins. m can be as large as possible only if there are sufficient galaxies in each bin. If the mean redshifts for these bins are {z s,j } j=1,m , we calculate the physical properties for the whole photometric catalog and select a population 2 for each z s,j . Then, we calculateS 2,j w 12,weight,j (ϑ)/r 2 1,j for each redshift bin, and the meann 2 w p (r 1 θ) of population 1 can be obtained by averaging over these redshift bins:

APPLICATIONS TO CMASS AND HSC SAMPLES
In this section, we apply PAC to CMASS spectroscopic sample in the Baryon Oscillation Spectroscopic Survey (BOSS; Ahn et al. (2012); Bolton et al. (2012)) and Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP; Aihara et al. (2019)) PDR2 wide field photometric sample.

Observational data
We use the HSC-SSP PDR2 wide field photometric catalog (Aihara et al. 2019) as the photometric sample. To obtain more accurate physical properties, we choose sources in the footprints observed with all five bands (grizy) to ensure that there are enough bands for SED. Sources around bright objects are masked using {grizy} mask pdr2 bright objectcenter flag provided by the HSC collaboration (Coupon et al. 2018). And we use the {grizy} extendedness value flag to exclude stars in the sample. Finally, there are around 2 × 10 8 galaxies in our photometric sample. We can also construct a random point catalog (100/arcmin 2 ) with the same selection criteria from the HSC database for ACCF analysis. The effective area calculated from the random point number is 501 deg 2 .
To ensure the distribution of foreground and background galaxies is the same for all LOS directions, the easiest way is to make sure that the sample is complete for galaxies with specified physical properties at the required redshift. Since the survey depth is not uniform across the HSC survey, for a low stellar mass limit, some patches remain complete while others may not. Therefore, we use HSC-SSP PDR2 deep field catalog and DEmP photo-z (Nishizawa et al. 2020) to study the completeness of galaxies for different stellar mass. We select galaxies with photo-z between 0.5 and 0.7 and calculate the physical properties for these galaxies with five bands grizy using the SED code CIGALE (Boquien et al. 2019). The stellar population synthesis models of Bruzual & Charlot (2003) are used to compute the physical properties of galaxies. In these calculations, the Chabrier (2003) Initial Mass Function is adopted. We assume a delayed star formation history φ(t) ≈ t exp −t/τ where τ is taken from 10 7 to 1.258 × 10 10 yr with an equal logarithmic interval ∆lgτ = 0.1. Three metallicities, Z/Z = 0.4, 1, and 2.5, are considered, where Z is the metallicity of the Sun. We use the extinction law of Calzetti et al. (2000) with dust reddening in the range 0 < E(B − V ) < 0.5. As shown in Figure 1, stellar mass shows a clear correlation with z band magnitude at 0.5 < z p < 0.7, so magnitude limit can be used to derive a complete sample for specified stellar mass. Particularly, we define the z-band completeness limit C 95 (M * ) that 95% of the galaxies are brighter than C 95 (M * ) in the z-band for given stellar mass M * . (orange line in Figure 1). Therefore, for stellar mass bin [m l , m h ] of interest, we only use the data in the survey footprints Figure 1. The stellar mass-z-band magnitude relation for HSC deep field galaxies with photo-z between 0.5 and 0.7. The blue dots are galaxies in the sample and the orange line shows the 95% completeness limit C95(M * ) of the z-band magnitude.
deeper than C 95 (m l ), where the depth is defined as the z band limiting magnitude of 10σ detection for a point source. We would like to stress again that the photo-z is used only here to find the completeness limit, and it is not used in the following clustering analysis.
We use the CMASS sample in BOSS that consists of massive galaxies with i < 19.9 mag as spectroscopic sample (population 1). We first select galaxies in the redshift range of 0.5 < z s < 0.7, and then cross match it with the HSC photometric sample we constructed above and The DESI Legacy Imaging Surveys DR9 catalog (Dey et al. 2019). After that, we get the magnitudes in seven bands grizyW 1W 2 for each CMASS galaxy in the footprint of HSC. We calculate the physical properties for these galaxies using the same SED templates but with seven bands grizyW 1W 2. As noticed by previous studies (Maraston et al. 2013;Leauthaud et al. 2016;Guo et al. 2018), the CMASS sample is complete to stellar mass M * ≈ 10 11.3 M . Therefore, we adopt a stellar mass cut at 10 11.3 M in this study. Moreover, we only consider central galaxies in the spectroscopic sample, so we select galaxies which do not have more massive neighbors within the projected distance of 1M pc h −1 .
Finally, there are 8028 massive (> 10 11.3 M ) central galaxies left in the spectroscopic sample.

PAC for different mass bins
We first divide the spectroscopic sample into three mass bins [10 11.3 , 10 11.5 , 10 11.7 , 10 11.9 ] M . In each mass bin, galaxies are divided into four redshift bins [0.5, 0.55, 0.6, 0.65, 0.7]. Then, we perform SED for the whole photometric sample at redshifts 0.525, 0.575, 0.625 and 0.675 respectively. After that, the photometric sample is also divided into four mass bins [10 9.0 , 10 9.5 , 10 10.0 , 10 10.5 , 10 11.0 ] M at each redshift. Using Equation 8,n 2 w p (r p ) for twelve mass bins (three spectroscopic × four photometric) in the redshift range 0.5 < z s < 0.7 can be obtained. During the calculation, only footprints deeper than C 95 (m l ) are used for each mass bin.
To evaluate the statistical error, we use the Jackknife resampling technique (Efron 1982). Mean value and error of the mean value ofn 2 w p (r p ) for each mass bin can be calculated as: where N sub is the number of jackknifed realizations and n 2,k w p,k (r p ) is the excess of projected density of the kth realization. In this work, we adopt N sub = 50. The results are shown in Figure 2 with colored dots. Each panel shows the results for the same mass bin of the photometric sample. Sincen 2 is the same in the same panel, the difference ofn 2 w p (r p ) reflects the difference of w p (r p ) between different spectroscopic mass bins. Though the footprint become smaller for the lowest mass bin 10 9.0 M , the clustering signal is still quite good for all the three spectroscopic subsamples, so we can study the properties such as SMF and SHMR down to the low mass end. To test the robustness of PAC to the redshift bin used, we compare the results when different number of redshift bins are used to divide the spectroscopic subsample. As we can see from Figure 3, for a narrow redshift range as in our study (0.5 < z < 0.7), the measurement is nearly independent of the number of redshift bins used, indicating that our algorithm is robust.

ABUNDANCE MATCHING WITH N-BODY SIMULATION
With the density and clustering information of galaxies, we can study the galaxy-halo connection using Nbody simulation. HOD and AM are the two most commonly used methods for populating galaxies to dark matter halos. We follow Wang & Jing (2010) and use AM to study the galaxy-halo relation based onn 2 w p (r p ) for different mass bins obtained in the observation. As we will show, the stellar-halo mass relation (SHMR), the stellar mass function (SMF) and the conditional stellar mass function (CSMF) for satellites can be inferred from the AM results for galaxies in a wide range of stellar mass (10 9.0 − 10 12.0 M ).

CosmicGrowth Simulation
We use the CosmicGrowth Simulation (Jing 2019) for our studies. CosmicGrowth Simulation is a grid of high resolution N-body simulations run in different cosmologies using an adaptive parallel P 3 M code (Jing & Suto 2002). We use the ΛCDM simulation with cosmological parameters Ω m = 0.268, Ω Λ = 0.732 and σ 8 = 0.831. The box size is 600 M pc h −1 with 3072 3 dark matter particles and softening length η = 0.01 M pc h −1 . Groups are identified with the Friends-of-Friends algorithm with a linking length 0.2 times the mean particle separation. The halos are then processed with HBT+ (Han et al. 2012(Han et al. , 2018 to obtain subhalos and their evolution histories. We use the catalog of Snapshot 83 at redshift about 0.57 to compare with the observation. We also use the fitting formula in Jiang et al. (2008) to evaluate the merger timescale for subhalos with less than 20 particles (including orphans), which may be unresolved, and abandon those that have already merged into central subhalos. In Figure 4, we compared our subhalo mass function in halos with M 200 = 10 14.0 M h −1 to the universal subhalo mass function of Han et al. (2016) who used high resolution zoom-in simulations (mass resolution ∼ 10 3 M h −1 for the highest one) and carefully corrected for the resolution effect. In this comparison, the halos are defined as objects of the radius R 200 within which the average density equals 200 times the critical density of the universe, and the subhalo mass is defined as its peak halo mass m peak in its history, to be consistent with Han et al. (2016). Our subhalo mass function is in good agreement with Han et al. (2016) down to 10 10.5 M h −1 , which is good enough for this study. The SHMR can be described by a formula of double power-law form: where M acc is defined as the viral mass M vir of the halo at the time when the galaxy was last the central dominant object. We use the fitting formula in Bryan & Norman (1998) to find M vir . The scatter in log(M * ) at a given M acc is described with a Gaussian function of the width σ. We use the same set of parameters for centrals and satellites (unified model) as in many studies (Wang & Jing 2010;Behroozi et al. 2019).
Once the parameters {M 0 , α, β, k, σ} are fixed, galaxies can be assigned to each dark matter halo. To com-paren 2 w p (r p ) with observation, we define χ 2 as: where N p is the total number of points over which n 2 w p (r p ) is compared. We only consider the radius range of 0.1 < r p < 10 M pc h −1 , in order to avoid the deblending problem of the HSC catalog at the smaller r p (Wang et al. 2021a), and large errors at r p > 10M pc h −1 . In order to perform a maximum likelihood analysis, we use the Markov chain Monte Carlo (MCMC) sampler emcee (Foreman-Mackey et al. 2013).
Posterior PDFs of the parameters of the SHMR model from MCMC are shown in Figure 5 and Table 1  the correspondingn 2 w p (r p ) and errors for each mass bin is shown by solid lines with shadows in Figure 2. The fitting is overall good for all mass bins in our samples and all the parameters are constrained well.

SHMR and SMF
After we derive the parameters, we can calculate the SHMR and SMF at redshift 0.5 ∼ 0.7 for stellar mass range 10 9.0 ∼ 10 12.0 M covered by our observational samples.
In the left panel of Figure 6, we compare the SHMR from our results with other works. Blue line with shadow shows the SHMR from our work. Green line shows the unified SHMR model from Wang & Jing (2010) by fitting the SMF and GC at z s ∼ 0.8 from VVDS (Pozzetti et al. 2007;Meneux et al. 2008). Moster et al. (2013) uses a redshift-dependent parametrization of SHMR and constrains the parameters by fitting the SMFs from SDSS (Li & White 2009), Spitzer (Pérez-González et al. 2008 and Wide Field Camera 3 (WFC3) (Santini et al. 2012) varying from z s ∼ 4 to the present. We show their results at z s = 0.6 with red line in the figure. The SHMR from our and other works are in good agreement with each other.
SMF from our AM model and other observational measurements are shown in the right panel of Figure  6. Blue dashed line with shadow shows the results from AM of our work. Orange dots show our own measurement from the CMASS spectroscopic sample at the high mass end. Green triangles and red diamonds show the measurements from VIPERS (0.50 < z s < 0.60; Davidzon et al. (2013)) and the PRism MUlti-object Survey (PRIMUS; 0.50 < z s < 0.65; Moustakas et al. (2013)) with spectroscopic redshift. And purple squares are the results from McLeod et al. (2021) (0.25 < z p < 0.75), in which they measure the SMF using various deep photometric surveys (UKIDSS Ultra Deep Survey (UDS), COSMOS and CFHTLS-D1) with photo-z. At the high mass end, the results from AM and various observational measurements are consistent. While at the low mass end, the discrepancy between observations is still significant (nearly by a factor of two between Davidzon et al. (2013) and Moustakas et al. (2013)), which has been exhaustively discussed in literature (Li et al. 2006;Pozzetti et al. 2007;Davidzon et al. 2013). Since it is hard to detect the faint objects and the survey areas of the high redshift deep surveys are usually small, completeness, selection effects and cosmic variance should be carefully considered. Different weighting methods to compensate for the incompleteness and regions with different densities that the surveys observed can produce huge difference in SMF. Therefore, the SMF of the faint end at higher redshift is still very hard to measure, for which PAC is a very promising method that combines the advantages of both the photometric and spectroscopic surveys. Interestingly, the SMF from our work using PAC and AM is well consistent with a very recent measurement of McLeod et al. (2021) down to 10 9.0 M .

CSMF of satellites
The CSMF of satellites around a central galaxy of given stellar mass can be derived from the PAC measurement and the AM modeling. We define satellite to be the galaxies within R vir of the dark matter halos of the central galaxies. To get the number of satellites, the most straight forward way is to sum up the excess surface densityn 2 w p (r p ) weighted by area within R vir . However, two effects should be corrected in our study. The surface density includes all excess of galaxies along the line-of-sight direction rather than within R vir (see Jiang et al. (2012)) and the measurement within 0.1R vir is unreliable for the HSC PDR2 photometric catalog due to the deblending problem (Wang et al. 2021a).
Therefore, we use the results from AM to compensate these two effects. After populating the galaxies to halos, we measure the average number of satellites within the projected radius range 0.1R vir < r p < R vir and within virial radius R vir for each central and satellite mass bins. Then, we calculate the average number of satellites within 0.1R vir < r p < R vir in observation and infer the number within R vir using the ratio calculated from simulation. We show the CSMF from 10 9.0 M to 10 11.6 M for central galaxies with different mass bins in Figure 7. Colored dots show the results from the observation and the solid lines are the results from AM. The results are consistence with each other at all mass bins for observation and simulation.

CONCLUSION
In this paper, we provide a method for estimating the projected density distributionn 2 w p (r p ) from w 12 (ϑ) and extending this method to measure the distributions and properties of Photometric objects Around Cosmic web (PAC) traced by spectroscopic surveys. Basically, by assuming the whole photometric sample at the same redshift as spectroscopic sources, we can calculate the physical properties of the photometric sample. And through cross-correlation, foreground and background galaxies with wrong properties are canceled out and the true distribution of photometric sources with specified physical properties around the spectroscopic sources can be obtained.
We apply PAC to massive (> 10 11.3 M ) central galaxies in BOSS CMASS sample (0.5 < z < 0.7) and HSC-SSP PDR2 wide field photometric sample. We calculatē n 2 w p (r p ) for several stellar mass bins (three for CMASS × four for HSC) from 10 9.0 M to 10 12.0 M and the measurement is overall good at 0.1 M pc h −1 < r p < 10 M pc h −1 for all the mass bins. Then, we use abundance matching to modeln 2 w p (r p ) in N-body simulation with MCMC sampling. We use the same set of parameters for central and satellite galaxies to model the observation. All the parameters are constrained well and the fitting to observation is overall good for all mass bins. Our AM model can accurately reproduce the SMF comparing to observations. The SHMR from our results are also in good agreement with previous works. Using PAC and AM, we also calculate the CSMF of satellites for centrals with different mass.
We expect that PAC will have many applications with ongoing and upcoming photometric and spectroscopic surveys. However, since the wide spectroscopic surveys at higher redshifts only target specific populations of galaxies like ELGs and QSOs, the galaxy-halo connection of these populations, which is also one of the key challenges for galaxy formation and cosmological studies, should be well established to make full use of PAC. Recently, there has been some progress made in ELG HOD and AM modeling with spectroscopic or narrow band data (Guo et al. 2019;Okumura et al. 2021, H, Gao et al. 2021. PAC can provide more information for studying the galaxy-halo connection of ELGs. As in this work, the SHMR of normal galaxies can be obtained using PAC with a stellar mass limited spectroscopic sample such as Large Red Galaxies (LRGs). If the redshift range of ELGs is overlapped with that of LRGs, by applying PAC to ELGs and the same photometric sample, we can obtain galaxy bias of ELGs with respect to the underlying matter distribution. With the better understanding of the galaxy-halo connection for ELGs, we can extend PAC to a higher redshift where LRGs cannot be reached. It is even worth trying to simultaneously study the connections of the ELGs and the normal galaxy population to dark matter halos, since the projected density distributionn 2 w p (r p ) is expected to be precisely measured with next generation galaxy surveys.
With PAC and future surveys, we can study the galaxy-halo connection for galaxies with different phys-ical properties other than mass, such as SFR, color and morphology (Xu & Jing 2021). We can also push the understanding of SHMR, SMF and other properties to higher redshift and fainter luminosity end. With the properties and distribution of satellite, we can also study the galaxy evolution such as galaxy merger rate, merger timescale and environment quenching. Moreover, since PAC has very strong signal at the small scale, we can also quantify the fiber collision effect in spectroscopic samples. We may also apply the method to photo-z samples to quantify photo-z errors. We will explore these applications in our future studies.

ACKNOWLEDGMENTS
The work is supported by NSFC (12133006, 11890691, 11621303) and by 111 project No. B20019. This work made use of the Gravity Supercomputer at the Department of Astronomy, Shanghai Jiao Tong University.
The Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University.
This publication has made use of data products from the Sloan Digital Sky Survey (SDSS). Funding for SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the U.S. Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society, and the Higher Education Funding Council for England.
The BASS is a key project of the Telescope Access Program (TAP), which has been funded by the National Astronomical Observatories of China, the Chinese Academy of Sciences (the Strategic Priority Research Program "The Emergence of Cosmological Structures" Grant # XDB09000000), and the Special Fund for Astronomy from the Ministry of Finance. The BASS is also supported by the External Cooperation Program of Chinese Academy of Sciences (Grant # 114A11KYSB20160057), and Chinese National Natural Science Foundation (Grant # 11433005).
The Legacy Survey team makes use of data products from the Near-Earth Object Wide-field Infrared Sur-vey Explorer (NEOWISE), which is a project of the Jet Propulsion Laboratory/California Institute of Technology. NEOWISE is funded by the National Aeronautics and Space Administration.
The Legacy Surveys imaging of the DESI footprint is supported by the Director, Office of Science, Office of High Energy Physics of the U.S. Department of Energy under Contract No. DE-AC02-05CH1123, by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract; and by the U.S. National Science Foundation, Division of Astronomical Sciences under Contract No. AST-0950945 to NOAO.