Cluster Cosmology Without Cluster Finding

We propose that observations of super-massive galaxies contain cosmological constraining power similar to conventional cluster cosmology, and we provide promising indications that the associated systematic errors are comparably easier to control. We consider a fiducial spectroscopic and stellar mass complete sample of galaxies drawn from the Dark Energy Spectroscopic Survey (DESI) and forecast how constraints on Omega_m-sigma_8 from this sample will compare with those from number counts of clusters based on richness. At fixed number density, we find that massive galaxies offer similar constraints to galaxy clusters. However, a mass-complete galaxy sample from DESI has the potential to probe lower halo masses than standard optical cluster samples (which are typically limited to richness above 20 and halo mass above 10^13.5); additionally, it is straightforward to cleanly measure projected galaxy clustering for such a DESI sample, which we show can substantially improve the constraining power on Omega_m. We also compare the constraining power of stellar mass-limited samples to those from larger but mass-incomplete samples (e.g., the DESI Bright Galaxy Survey, BGS, Sample); relative to a lower number density stellar mass-limited samples, we find that a BGS-like sample improves statistical constraints by 60% for Omega_m and 40% for sigma_8, but this uses small scale information which will be harder to model for BGS. Our initial assessment of the systematics associated with supermassive galaxy cosmology yields promising results. The proposed samples have a 10% satellite fraction, but we show that cosmological constraints may be robust to the impact of satellites. These findings motivate future work to realize the potential of super-massive galaxies to probe lower halo masses than richness-based clusters and to avoid persistent systematics associated with optical cluster finding.


INTRODUCTION
The current ΛCDM paradigm successfully explains numerous cosmological phenomena, including the formation and growth of largescale structure and the cosmic expansion of the Universe (e.g., Weinberg et al. 2013;Lahav & Liddle 2014).However, the physical origin of cosmic acceleration is still poorly understood.Cosmic acceleration is studied through probes such as Type Ia Supernovae (SNe), baryon acoustic oscillations (BAO), weak gravitational lensing, and galaxy clusters (e.g.Albrecht et al. 2006;Dawson et al. 2018;Abdalla et al. 2022, and references therein).Type Ia Supernovae, and ★ E-mail: exhakaj@ucsc.edubaryon acoustic oscillations probe the expansion history, while weak gravitational lensing and galaxy clusters probe the growth of dark matter structure.In turn, the expansion history of the Universe and the growth of dark matter structure shed light on cosmic expansion.
Galaxy clusters are formed from the highest peaks of the matter density field.Their abundance and spatial distribution contain a wealth of information on the growth of dark matter structure (Schuecker et al. 2003;Schmidt et al. 2004;Bonamente et al. 2006;Albrecht et al. 2006;Allen et al. 2008;Ettori et al. 2009;Vikhlinin et al. 2009;Henry et al. 2009;Mantz et al. 2010;Rozo et al. 2010;Allen et al. 2011;Weinberg et al. 2013;Mantz et al. 2015;Planck Collaboration et al. 2016;Bocquet et al. 2019;Costanzi et al. 2021;To et al. 2021;Chiu et al. 2022;Salvati et al. 2022;Lesci et al. 2023;Park et al. 2023).Recently, optical cluster counts and weak gravitational lensing have offered competitive constraints on  8 and Ω m (e.g.Rozo et al. 2010;Abbott et al. 2020;To et al. 2021;Costanzi et al. 2021;Giocoli et al. 2021;Lesci et al. 2022;Park et al. 2023).Optical surveys such as the Dark Energy Survey (DES, The Dark Energy Survey Collaboration 2005) and Hyper Suprime-Cam (HSC, Aihara et al. 2018) have detected larger cluster samples than surveys in other wavelengths providing better statistics for cosmological analyses.Three important systematics associated with optical cluster finding are selection effects, projection effects (the failure to separate halos along the line of sight into distinct objects), and cluster miscentering.These effects can bias cosmology if not modeled correctly (e.g.Erickson et al. 2011;Zhang et al. 2019;Costanzi et al. 2019;Sunayama et al. 2020;Wu et al. 2021).A great deal of effort has been devoted to building mocks that study projection effects (DeRose et al. 2019;Korytov et al. 2019;Sunayama et al. 2020;DeRose et al. 2022;Wu et al. 2022;Wechsler et al. 2022;Park et al. 2023).However, mocks with accurate red galaxy populations have proven challenging to build.In addition to galaxy color, the spatial distribution of cluster satellites is also challenging to model (DeRose et al. 2019;Korytov et al. 2019;DeRose et al. 2022;Wechsler et al. 2022;To et al. 2023).
Another popular method to constrain the growth of structure is the combination of galaxy-galaxy lensing and projected clustering.Galaxy-galaxy lensing is proportional to Ω m  82 at large radial scales, where  is the galaxy bias (e.g., Yoo et al. 2006), Ω m is the mean matter density of the Universe, and  8 is the amplitude of the power spectrum.Projected clustering, on the other hand, scales as  2  8 2 .When fit jointly, the bias term is constrained, allowing for a measurement of Ω m  8 (Baldauf et al. 2010;More et al. 2013;Mandelbaum et al. 2013;Cacciato et al. 2013;More et al. 2015;Leauthaud et al. 2017;Lange et al. 2019;Singh et al. 2020).Cosmological constraints from lensing and clustering analyses have often been carried out using red galaxy samples.For example, the Baryon Oscillation Spectroscopic Survey (BOSS, Dawson et al. 2013) provided a large number of massive galaxies with spectroscopic redshifts that have been used for lensing plus clustering studies (Miyatake et al. 2015;Leauthaud et al. 2017;Singh et al. 2020;Sunayama et al. 2020;Lange et al. 2022;Tröster et al. 2022;Leauthaud et al. 2022;Amon et al. 2022;Lange et al. 2023).The DES survey instead used photometric samples of red galaxies (Abbott et al. 2019(Abbott et al. , 2022)).However, these galaxy samples have typically been selected with color cuts and are rarely complete, except at the highest galaxy masses (Leauthaud et al. 2016).As such, unlike studies of galaxy clusters where abundance is a key constraint, traditional lensing plus clustering analyses do not use the observed number density for cosmological constraints 1 .
Here we introduce the idea of using mass-complete samples of super-massive galaxies to constrain cosmology by applying a methodology that unifies aspects of cluster cosmology and traditional lensing plus clustering analyses.Table 1 summarizes the key observational probes typically used for the two classic approaches.As highlighted in Table 1, cluster cosmology has typically not used information from the clustering of clusters (but see Wu et al. 2021;Park et al. 2021;To et al. 2021) whereas lensing plus clustering typically does not use number density as a constraint.As explained in the previous paragraph, this is because galaxies samples from surveys 1 For example, the  Γ parameter in Lange et al. (2022) allows the amplitude of the central halo occupation float as a free parameter which removes any constraining power from the observed number density.such as BOSS are highly incomplete.However, the Dark Energy Spectroscopic Survey (DESI, DESI Collaboration et al. 2016) will change this picture (DESI Collaboration et al. 2016).Indeed, DESI will be deep enough to detect large samples of massive galaxies which are complete above  * = 11.5 M ⊙ at intermediate redshifts (e.g.,  < 0.6).We propose that mass complete samples from DESI can be used in a similar fashion to traditional cluster cosmology but with two added advantages: 1)it may be possible to bypass some of the complications associated with cluster finding (see below), and 2) it will be possible to take advantage of spectroscopic galaxy clustering measurements from DESI.
New results from Huang et al. (2022, hereafter H+22) show that super-massive galaxies can be used to trace dark matter halos with scatter comparable to state-of-the-art red sequence methods.This work suggests that the stellar mass measured within  = [50, 100] kpc (the outer envelope of the stellar mass distribution -hereafter the outskirt mass) yields a halo mass tracer with scatter similar to red-sequence-based cluster finders such as redMaPPer (e.g., Rykoff et al. 2014;Rozo & Rykoff 2014;Rozo et al. 2015a,b;Rykoff et al. 2016) and CAMIRA (e.g., Oguri 2014;Oguri et al. 2018).The outskirt mass has already been measured in multiple surveys (Li et al. 2021) and will be well measured with future surveys like the Rubin Observatory Legacy Survey of Space and Time (LSST, Ivezic et al. 2008).
Given the results of H+22, we explore here the idea that massive galaxies from DESI can be used to trace halos, thus bypassing the cluster-finding process.This method will not rely on prior knowledge about red galaxies in clusters and will not suffer from cluster selection effects such as those described in Abbott et al. (2020).However, we expect this method to have other systematics that must be studied further.In particular, about 10 − 20% of massive galaxies will be satellite galaxies.For example, Mandelbaum et al. (2006) predict  sat = 0.16 ± 0.09 for early type galaxies of mass  * = 10 11.3 M ⊙ .Reddick et al. (2013) find  sat = 0.13 ± 0.05 for galaxies of mass  * = 10 11.6 M ⊙ .Finally, Saito et al. (2016) find  sat = 0.11 for CMASS galaxies at  = 0.55 and of mass  * > 10 11.6 M ⊙ .In order to use our proposed methodology, it will be important to model the impact of satellite galaxies.However, satellite modeling may be more straightforward than understanding and making mocks of red galaxy populations.This paper aims to study how super-massive galaxies from the DESI survey could be used for cosmological constraints and how this methodology compares with other traditional analyses (see Table 1).Our main probes are galaxy-galaxy lensing (ΔΣ), projected galaxy clustering ( p ), and galaxy number density ( ngal ).We adopt three narrow mass bins in the outer stellar mass range 2 10 10.8 − 10 12.1 M ⊙ and a redshift bin of  ∈ [0.3, 0.6].We assume the full 1000 deg 2 area of HSC.We show how our constraints compare with traditional methods such as cluster cosmology and lensing plus clustering analyses.It is important to note that because this is a new method, systematics (e.g., satellite modeling) must be quantified and understood.This will be the goal of future work.Here we focus first on the relative statistical constraining power of different techniques and analysis choices.
This paper is structured as follows.In Section 2, we introduce our model.In Sections 3 and 3.2, we introduce our model fitting routine and the derivation of the covariance matrix.We present our results in Section 4. We discuss our results and possible future direc-tions in Section 5. Finally, we summarize and conclude in Section 6.We adopt the cosmology of the Covariance AbacusSummit Boxes (Maksimova et al. 2021), namely, a flat, ΛCDM cosmology with Ω m = 0.307, Ω b = 0.0482,  8 = 0.829, ℎ = 0.678,  s = 0.9611, corresponding to the best-fit Planck cosmology (Planck Collaboration et al. 2020).

THEORETICAL FRAMEWORK
Our goal is to study the cosmological constraining power of complete samples of super-massive galaxies.Our model is based on the halo model in which all dark matter is clumped into spherically over-dense regions known as halos.The following sections describe how the data vectors are modeled.We first introduce our mass to observable relations.We then describe how we model the data vector.

Philosophy
Traditionally, the correlation between halo mass and stellar mass has been modeled through the stellar-halo mass relation (SHMR, e.g., Behroozi et al. 2010;Leauthaud et al. 2012;Tinker et al. 2017;Kravtsov et al. 2018; also see Wechsler & Tinker 2018 for a recent review).This relation has generally been calibrated using the measured stellar masses of central galaxies.The total stellar mass of galaxies can be divided into inner and outer components.Recently, Huang et al. (2020Huang et al. ( , 2022) ) showed that the inner light of massive galaxies has a very poor correlation with halo mass (scatter of halo mass at fixed stellar mass,  log 10 M h |log 10 M * , is about 0.8 dex), unlike the outskirt mass of massive galaxies ( log 10 M h |log 10 M * ∼ 0.4 dex) 3 .The latter is comparable to state-of-the-art red sequence cluster finders and may outperform red sequence methods below  = 10 (Figure 8, H+22).
For this reason, instead of using total galaxy mass, we adopt the outskirt mass.More specifically, this is the stellar mass measured at  = 50 − 100kpc (hereafter  * , [50,100] ).Although an SHMR using  * , [50,100] has not yet been calibrated, this paper assumes an unbroken power-law relation with log-normal scatter.This is justified because we are focusing on just the very high mass range (log 10  * , [50,100] ∈ [10.8, 12.1]) well above the pivot mass (Leauthaud et al. 2012) where the SHMR can be approximated as a simple log-normal relation.Furthermore, we assume that our massto-observable relation models only central galaxies and neglects satellites.The impact of satellites will be discussed in Section 4.4.
Our log-normal model (N ) consists of three free parameters: the slope of the power law (), the scatter in stellar mass or richness ( log 10  * or  log 10  ), and the amplitude of the stellar mass function ( 0 ).Given a halo proxy O, the mass to observable relation4 is: (1) H+22 do not provide the calibrated mass to observable relation for  * , [50,100] .Here, we begin by deriving an SHMR for  * , [50,100] and  that matches the bins in  * , [50,100] and richness from H+22 along with their respective underlying halo distributions.In the following subsections, we describe how this goal is achieved.Note that we use the total scatter in the mass to observable ratio.This includes the intrinsic scatter of the relation convolved with the statistical uncertainties of stellar mass and richness measurements.

2.1.2
The halo mass to  * , [50,100] relation From measurements based on the halo mass predictions of the ASAP model, we derive the slope of the halo mass to  * , [50,100] relation to be  [50,100] = 0.7 (Huang et al. 2018(Huang et al. , 2020(Huang et al. , 2022)).We tune the remaining free parameters in Eq. ( 1) to reproduce the halo mass distribution of H+22 The left panel of Figure 1 displays this relation.

The halo mass to richness relation
Previous studies have modeled the mass-richness relation with a simple log-linear function (e.g.McClintock et al. 2019) with a slope of   = 0.9.To calibrate the other free parameters ( log 10  and  0 ), we follow the same methodology as above.Assuming a pivot halo mass of  pivot = 3 × 10 14 M ⊙ , the mass to richness relation (Figure 1, right panel) is: log 10  = N 0.9 log 10 ( ℎ / pivot ) + 1.15, 0.46 . (3) Note that   is larger than  logM * , [50,100] (Equation 2).This is the result of the degeneracy between  and  log 10 ( O ) (H+22).Different combinations of  and  log 10  * can yield similar distributions in halo mass.Because the slopes of the mass to observable relations are different, a relative comparison of  log 10  * , [50,100] and   does not provide any meaningful information.For this project, we calibrate both mass to observable relations to replicate the same halo mass distributions as in H+22.We illustrate this further in Section 3.3 and Figure 3.

The Stellar Mass Function
The stellar mass function (Φ (O |  ℎ ), hereafter SMF) describes the number of galaxies (clusters) 5 with halo proxy, O, in the range Table 1.Typical observables and analysis choices for cluster cosmology and lensing plus clustering.Cluster cosmology uses ΔΣ and cluster abundance (which offers the same information as ngal ) but typically does not use the clustering of clusters (except in more recent studies such as Mana et al. 2013;Park et al. 2021;To et al. 2021).Cluster samples are also studied in bins of cluster richness (for example, ), defined as the number of red galaxies associated with a cluster within a certain radius.The traditional lensing plus clustering technique uses ΔΣ and  p (but not ngal ) and often uses one large sample without any binning in galaxy mass.Here we propose to unify the two techniques using complete samples of super-massive galaxies from the DESI survey.In this approach, all observables are used as constraints, and the data can be divided into mass bins.O ± O/2, at constant halo mass.To measure the SMF, we start by computing the conditional SMF:

ΔΣ
The average count of galaxies (clusters) within the bin hosted by a halo of mass  ℎ can be computed through the halo occupation function: Finally, the SMF, the abundance of galaxies in a stellar mass bin, can be computed from the halo mass function and the halo occupation function: where as a halo mass function, d d ℎ , we use Tinker et al. (2008).

Power Spectra
In the halo model, the cross-galaxy-matter power spectrum can be written as the sum of the one and two-halo terms defined as follows: Here, ũm (, ) is the Fourier transform of the radial matter density profile of a halo mass  ℎ .ũm (, ), is modeled assuming a Navarro, Frenk & White profile (NFW, Navarro et al. 1996) with the Diemer & Joyce (2019) mass-concentration relation (, ).
In the two halo term,  lin (, ) is the linear matter power spectrum, while  O () is the mean linear bias of galaxies (clusters) computed as: where  h () is the halo bias relation in Tinker et al. (2010).The galaxy-galaxy power spectrum can be modeled as follows:

Observable Data Vectors
We derive our three data vectors (ΔΣ,   and ngal ) 6 from the power spectra and stellar mass function introduced previously.ΔΣ and   are obtained from the galaxy-matter and galaxy-galaxy two-point correlation functions respectively.These are just the inverse Fourier transform of the galaxy-matter and galaxy-galaxy power spectra: Here,  and  represent the galaxy and matter components of the power spectrum.
To compute ΔΣ, we begin with the surface density, Σ.The surface density is the integral of the galaxy-matter correlation function along the line of sight: ΔΣ is the excess surface density along the line of sight.Therefore it can be computed via: where Σ(< |) is the mean surface density spanning from  ′ = 0 to  ′ = : To compute   , we integrate the galaxy-galaxy correlation function along the line of sight (Davis & Peebles 1983;van den Bosch et al. 2013): Finally, the galaxy (cluster) number density, ngal is equal to the SMF measured between O  1 and O  2 : (15)

Model Implementation
Our theoretical framework is implemented in CosmoLike (Krause & Eifler 2017).CosmoLike is a code package that aims to constrain 6 Derivations in this section are adapted from van den Bosch et al. ( 2013).
cosmology through fast likelihood analysis of joint cosmological probes.In our case, these probes are ΔΣ,  p , n, and are derived assuming the halo model and a fiducial cosmology as described in Sec. 2. To constrain cosmological parameters, we follow a Bayesian approach when fitting a fiducial data vector given a joint covariance.
We describe the joint covariance in Section 3.2 and our fiducial sample in Section 3.3.

Covariance matrix
To calculate a covariance matrix, we use the AbacusSummit covariance simulation suite (Garrison et al. 2018).These are a suite of dark matter-only simulations with box size 500 Mpc/ℎ and particle mass  p = 2•10 9 M⊙/ℎ.The set of covariance boxes includes 2000 different realizations, 15% of which are unavailable at the redshift of interest.Thus, for this project, we use 1700 boxes at  = 0.2.We generate galaxy and cluster mock catalogs using the CompaSo halo catalogs available for this suite (Hadzhiyska et al. 2021).First, we select all host halos.Then, we assign a central galaxy or cluster to host halos using a mass-to-observable relation, including scatter (Section 2.1).We create bins in stellar mass and richness, corresponding to our fiducial samples (Table 2).Then, we compute our data vectors from all 1700 mocks.We build galaxy-galaxy lensing profiles through the excess surface density of dark matter particles in cylinders surrounding host halos.This is implemented in the mean delta sigma function of Halotools (Hearin et al. 2017).We compute  p by integrating the redshift space two-point correlation function up to a  max = 100Mpc/ℎ implemented in the Halotools function w p.
We build a covariance matrix using the mock data vectors from all 1700 mock realizations and show the results in Figure 2.There are two additional steps to obtain the covariance matrix for our analysis.First, for simplicity, we assume that ΔΣ and  p are drawn from the same survey area of 1000 deg 2 (although, in practice, for DESI, these areas will be different).We rescale the covariance matrix to the HSC area.The scaling factor is computed following: where  1 is the comoving volume to  = 0.3,  2 is the comoving volume to  = 0.6,  HSC is the HSC survey area (1000 deg 2 ), and  sky is the full sky area,  sky = 41253 deg 2 .Second, we include lensing shape noise in the ΔΣ data vector.We compute the shape noise analytically as in Singh et al. (2017).Finally, we add the diagonal terms of the shape noise to the simulation covariance.

Fiducial Galaxy and Cluster Samples
We define several samples in a single redshift bin  ∈ [0.3, 0.6].For most of this paper, we assume samples drawn from the intersection of HSC wide and DESI and assume a survey area of 1000 deg 2 .The primary goal of this paper is a relative comparison of the statistical constraining power of different techniques, and so the exact area assumed is not of great importance.The use of the fiducial samples is twofold: We define the fiducial samples in our analytical model used to perform the cosmological analysis; we also define them in the mock realizations from the Abacus suite.The latter is necessary in order to build the covariance.We consider three sets of samples.The first set corresponds to the super-massive galaxy and cluster samples used in H+22, hereafter S  * [50,100] ,bins .This set consists of 3 samples binned by outskirt mass and  (see Table 2 for the bins' properties).We illustrate the underlying halo mass distribution of S  * [50,100] ,bins in Figure 3.These distributions are built with a mock realization of the AbacusSummit (specifically Phase 3000 at  = 0.2; see Section 3.2 for the description of our mocks).The  * , [50,100] and  bins in Figure 3 have similar underlying mass distributions.This is by design and is motivated by the work of H+22.As shown in H+22, galaxy outskirt mass is an excellent proxy of halo mass.Hence, it is possible to construct samples binned by outskirt mass with similar underlying halo mass distributions to ones binned by .Based on early analysis of DESI data, we expect that DESI will obtain complete samples of galaxies in the mass range probed by these bins.
S  * [50,100] ,bins extends to a lower halo mass than the fiducial sample of the DES cluster cosmology analysis (Abbott et al. 2020).The DES cluster sample ranges from  = 20 to  = 200.This corresponds to our highest richness bin.The lowest limit on our cluster sample is  = 9.22.We believe pushing to lower halo mass samples may be possible with DESI, as the DES sample was limited by cluster-finding systematics.In particular, these affect the low richness regime the most and can be notoriously difficult to accurately model and marginalize over (DeRose et al. 2019;Korytov et al. 2019;Sunayama et al. 2020;DeRose et al. 2022;Wu et al. 2022;Park et al. 2023).In contrast, pushing to lower halo masses with the outskirt mass tracer may be more straightforward, although future work is needed to determine if this is indeed possible.Pushing to lower halo masses (higher number densities) will improve the signal-to-noise of our data vector (especially  p ), leading to tighter cosmological constraints.
Second, we consider a galaxy sample that combines all galaxies of S  * [50,100] ,bins into a single bin.We name this sample S  * [50,100] ,full .This will be a cumulative threshold sample that will include all super-massive galaxies with stellar mass above  * = 11.5M⊙ .Because DESI will be complete above  * = 11.5M⊙ , S  * [50,100] ,full will also be a mass complete sample.This sample aims to study the impact of binning on our cosmological constraints.
Finally, we also assume a sample of galaxies with the same number density as the DESI Bright Galaxy Survey (BGS) sample at  ∼ 0.3, hereafter S BGS .BGS is a magnitude-limited sample ( < 19.5) at low redshift with ngal = 1.70 × 10 −3 ℎ 3 /Mpc 3 (Hahn et al. 2022).Therefore S BGS has a much higher number density than S  * [50,100] ,full ( ngal = 2.36 × 10 −5 ℎ 3 /Mpc 3 ).However, the full BGS sample will not be mass-complete.We do not expect to measure outskirt masses for the full BGS sample for two reasons.[1] To measure the outskirt mass, we use a fixed physical kpc boundary.This means we cannot apply it to the high number density regime since 50 kpc is too large for low-mass galaxies.[2] Extending our sample to low-mass galaxies would mean that we need to introduce a break in our SHMR model, and we are not currently doing this.The purpose of our BGS-like sample is solely to compare the constraining power of our massive galaxies sample (S  * [50,100] ,bins ) with a larger but mass incomplete sample (S BGS ).Thus, in the context of this work, the BGS-like sample is a comparison sample that aims to understand the impact of mass completeness on cosmological constraints.All three samples and their properties are displayed in Table 2.
Figure 4 illustrates the different samples.The S  * [50,100] ,bins is designed to study how samples of galaxies binned by outskirt mass compared with samples binned by .The S  * [50,100] ,full sample is designed to study the impact of binning.Finally, S BGS allows us to study the trade-off between samples that are large but incomplete and for which number density cannot be used in the cosmological parameter fitting versus samples that have overall lower number densities but which are complete, and for which number densities can be used for cosmological parameter fitting (similar to cluster cosmology analyses).Throughout this paper, the S  * [50,100] ,bins is our fiducial sample, and we study how cosmological constraints from this sample compare with constraints from the other samples.

Radial range of data vectors
We use CosmoLike to predict the projected observables, ΔΣ and  p , in 10 logarithmically distributed radial bins ( p ) ranging from 0.4 to 40 Mpc/ℎ.The number of radial bins is determined such that the Hartlap factor is minimal (∼ 3%) (Hartlap et al. 2007).Since we use data vectors of different lengths throughout the paper, the Hartlap factor will change.However, we expect this change to be at the percent level due to the large number of realizations we are using (1700 different phases of the Abacus Suite).Both the upper and lower radial limits are determined by the AbacusSummit Simulations (Maksimova et al. 2021), which we used to compute the Gaussian and non-Gaussian components of the covariance (Section 3.2).The lower limit is restricted by the resolution of the simulations, while the upper limit is set by the box size (500 Mpc/ℎ).

Model parameters
Model parameters are divided into two categories.The first set consists of cosmological parameters.We vary Ω m and  8 as they are the main parameters influencing the amplitude or the shape of the halo mass function.Other cosmological parameters are fixed to the values in Table 3. Figure 5 shows how the halo mass function changes when we vary Ω m and  8 .Ω m affects the overall amplitude, whereas  8 affects the massive end of the halo mass function.At the halo mass range of interest, we should, in principle, be able to constrain both Ω m and  8 .
The second category of model parameters describes the massto-observable relation.These are the slope, the y-intercept, and the scatter in the mass to observable relation.In Figure 6, we study how Table 2. Fiducial richness and outer stellar mass bins for low, medium, and high halo masses (top row), a cumulative threshold sample (mid rows), and a DESI BGS-like sample (bottom row) along with the respective number densities.Note that the stellar mass values shown here are measured in outskirt mass rather than the total stellar mass.A galaxy with  * , [50,100] ≈ 10 10.85 M ⊙ resides in a halo with mass  ℎ ≈ 10 13.5 M ⊙ /ℎ.[50,100] and , respectively.The underlying halo mass distribution for massive galaxies and clusters is similar across all three bins.
the SMF changes as we vary each parameter.All three parameters affect the shape and amplitude of SMF at the stellar mass range of interest.Therefore, we vary all these parameters in our fits.

Model Fitting
We follow a Bayesian approach to fit our fiducial data vectors.We sample posteriors from the parameter space following a Markov Chain Monte Carlo (MCMC) approach with emcee (Foreman-Mackey et al. 2013).We follow the same likelihood analysis as in Krause & Eifler (2017), in which the likelihood of the cosmological and nuisance parameters is parametrized jointly as a multivariate Gaussian .
(17) Here, D is the multi-probe data vector computed at the fiducial cosmological (p  ) and nuisance (p  ) parameters chosen for this paper.C is the covariance matrix measured with the Abacus simulations (Figure 2).Finally, the model M (p c , p n ) is a function of the free parameters.
We use broad, uninformative top hat priors for all free parameters (Table 3).We assess convergence by computing the mean auto-correlation time across all dimensions for each MCMC step.We assume that convergence is achieved when the change of the mean auto-correlation time between steps is smaller than 1% after at least 100 correlation times (Sokal 1996) 7 .

RESULTS
This paper introduces massive galaxies as competitive halo tracers for future cosmological analyses.Our data vector consists of ΔΣ,  p , and n.We assume 1000 deg 2 of HSC lensing data with spectroscopic redshifts from DESI and a set of fiducial parameters given in Table 3.This section presents our main results.First, we study the impact of binning on cosmological constraints.We then ,full , is one threshold-limited sample that combines all three bins from S  * [50,100] ,bins .The third sample is S BGS (pink).This sample has a higher number density and extends to lower halo masses but is incomplete in galaxy mass.
Table 3. Model parameters and priors used in the joint likelihood analysis introduced in Section 3.6.The top row describes the cosmological parameters, while the two bottom rows describe the free parameters of the mass to observable relations (Section 2.1).

Parameter
Fiducial Prior study the impact of ngal and  p on  8 and Ω m .We then compare how the massive galaxy methodology compares with two other traditional methods: 1) cluster cosmology and 2) joint lensing plus clustering analyses.Finally, we investigate one potential systematic associated with this methodology: the impact of satellites.

Impact of binning
We start by analyzing the impact of binning on our constraints.We consider two scenarios: one in which the super-massive galaxies are divided into three narrow bins by outskirt mass (S  * [50,100] ,bins ) and another in which all galaxies reside in a single cumulative bin (S  * [50,100] ,full ; see Table 2 for details on each sample).Both scenarios assume the same survey area (HSC 1000 deg 2 ), fiducial parameters (Table 3), and data vectors (ΔΣ,  p , and, ngal ).
Figure 7 shows how constraints on Ω m and  8 compare for S  * [50,100] ,bins and S  * [50,100] ,full .We find that binning by stellar mass does not improve constraints on Ω m but yields stronger constraints for  8 .The S  * [50,100] ,bins sample yields  8 = 0.86 +0.08 −0.06 while S  * [50,100] ,full yields  8 = 0.85 +0.12 −0.11 .Binning by galaxy mass, thus, yields a 34% improvement over a threshold sample.We can understand this result by looking at the impact of  8 and Ω m on the HMF in Figure 3. Ω m and  8 have different impacts on the HMF: Ω m mainly influences the overall amplitude, while  8 induces mass-dependent changes.At  h ∈ [10 13.2 , 10 14 ]M ⊙ ℎ −1 , the slope of HMF .This translates into tighter constraints in  8 .
Our result qualitatively agrees with Wu et al. (2021), who demonstrate that when the bins are correlated with a mass to observable relation, binning does indeed improve constraints in  8 .One important assumption here, however, is that the mass-to-observable relation follows a simple log-linear relation with constant scatter and satellite fraction throughout the three bins.While it needs to be studied in greater detail, this is a reasonable assumption for our work as the massive galaxy sample that we adopt covers a very narrow mass range (e.g.McClintock et al. 2019) and well above the pivot mass scales at which the SHMR bends (e.g., Leauthaud et al. 2012).

Impact of ngal
Next, we study the impact of number density on cosmological constraints.Utilizing the sample S  * [50,100] ,bins , we consider two scenarios: one in which we include ngal in our data vector and another in which we exclude it.The results are displayed in Figure 8 (purple contours versus teal contours).Including number density improves  8 constraints by 33% and Ω m by 23%.It is clear that ngal affects both cosmological parameters.This has a simple explanation.By definition, ngal drives the amplitude of the SMF and, thus, the amplitude of the HMF.As both parameters impact the amplitude of the HMF, ngal helps to constrain both.

Impact of 𝑤 p
Here we analyze the impact of  p on cosmological constraints.We adopt S  * [50,100] ,bins as our baseline sample and consider two cases: one in which we include  p as part of the data vector along with ΔΣ and ngal and another in which we exclude it.Figure 9 displays the results.Clearly,  p substantially impacts Ω m , improving the constraints by as much as 84%.This is because the combination of lensing and clustering is historically known to be a powerful cosmology tool, and using the two in combination, helps to break the degeneracy with galaxy bias (Yoo et al. 2006;Baldauf et 3).All three nuisance parameters affect the amplitude and the shape of SMF.
Traditionally, cluster cosmology with optically selected clusters only uses lensing and number densities as their main observables (Abbott et al. 2020;Costanzi et al. 2021).When samples are limited to high-mass halos, the measurement of the clustering of clusters can be noisy and, therefore, may not lead to improvements in cosmological parameter constraints.Only a handful of analyses have explored the clustering of clusters as part of the data vector (Mana et al. 2013;Baxter et al. 2016;Park et al. 2021;To et al. 2021), but in these analyses, additional data are required to constrain cosmological parameters.For example, To et al. (2021) carried out a 4x2 analysis, including the autocorrelation of galaxies and clusters, the cross-correlation of galaxies and clusters, cluster lensing, and cluster abundances.When comparing their constraints with the traditional cluster cosmology results from Costanzi et al. (2021), they report a 38% improvement in constraints of Ω m and mild improvements in constraints on  8 .Qualitatively this agrees with our results, although this improvement is less significant quantitatively.This is mainly due to the different halo mass cuts employed.Indeed, To et al. (2021) use a cluster sample with a richness cut of  ∈ [20,235].Instead, our cuts are defined by our estimates for the galaxy mass completeness limit for DESI, which corresponds to  ∈ [9.22, 150].This points to a significant advantage in using super massive galaxies over traditional clusters.Namely, supermassive galaxies may enable us to push down in halo mass into the group regime, where the number densities are higher and  p offers larger constraining power.

Massive Galaxies vs. Clusters
Here we compare constraints from clusters versus super-massive galaxies.We assume the same data vector (ΔΣ,  p , and n) and adopt three bins of richness and outskirt mass (Figure 3).We use the covariance predicted for a 1000 deg 2 survey (Section 3.2). Figure 10 compares traditional cluster cosmology and massive galaxy cosmology.While constraints in Ω m are similar, massive galaxies yield tighter constraints in  8 by 36%.3), and the same binning system (Figure 3).This assumes a galaxy sample that selects halos with a similar number density as a  > 9.22 selection.Including  p greatly improves constraints on Ω m , with an improvement of 84%.This figure presents one of the key potential advantages of massive galaxies over clusters: the ability to push to lower halo masses where  p is better-measured.
The fact that both methods offer similar statistical constraining power can be explained by Figure 3.Although the mass-toobservable relations for massive galaxies and clusters are different, they present similar values for the scatter in halo mass at fixed observable (Figure 8, H+22) and can therefore be used to select similar bins in halo mass.Because they trace the halo mass function in a similar fashion, they also derive comparable constraints for  8 and Ω m .
We investigate why massive galaxies yield tighter constraints in  8 .The only difference between the cluster and massive galaxies analyses in Figure 10 is the mass to observable relation.To test which parameters of the mass to observable relation shift the contours, we vary each of them in their respective prior range.We recompute the cosmological constraints for each variation.Our tests conclude that the shift in the contours observed in Figure 3 is related to differences in the power-law index of the mass-observable relation.
We conclude that massive galaxies can be used as alternatives to clusters to perform similar cosmological analyses.However, massive galaxies may offer advantages over traditional cluster cosmology.These arguments will be presented further in Section 5.1.

The Trade-off Between Complete versus Incomplete Samples
In lensing and clustering cosmological analyses with spectroscopic samples, number density is seldom used as a constraint (e.g.van Uitert et al. 2018;Joudaki et al. 2018;Abbott et al. 2018;McClintock et al. 2019;Singh et al. 2020;Heymans et al. 2021;Miyatake et al. 2021;Abbott et al. 2022;Amon et al. 2022;Lange et al. 2023).This is because galaxy samples from surveys such as BOSS are incomplete, except in specific mass and redshift ranges (Leauthaud et al. 2016).In contrast to BOSS, however, DESI will be complete .Constraints for Ω m and  8 with massive galaxies (teal) and clusters (pink) as halo tracers.Both analyses assume the same survey area (HSC 1000 deg 2 ) and fiducial data vector (Table 3).Massive galaxies yield competitive constraints while potentially avoiding systematics such as projection effects and miscentering.
above log 10  * = 11.5 M ⊙ over a wide range in redshift which will allow number density to be used as a constraint in a similar fashion to galaxy clusters.Figure 8 shows that including ngal as a third observable largely improves cosmological constraints when using high-mass galaxy samples with number densities similar to those from cluster studies.However, Figure 8 was limited to high-mass galaxies with low number densities.In contrast, lensing and clustering studies using BOSS have often used the full BOSS samples (Leauthaud et al. 2017;Singh et al. 2020;Sunayama et al. 2020;Lange et al. 2022;Tröster et al. 2022) which are an order of magnitude larger in count than complete samples.Lensing and clustering studies with photometric red galaxies typically have higher number densities (Abbott et al. 2022).This raises the question of whether or not it is more advantageous to use low-number density samples that are complete or higher number density samples which will afford better signal to noise on ΔΣ and  p but which sacrifice ngal .We now study this question using as our baseline S BGS , a single-bin DESI-like BGS sample consisting of many more galaxies than our S  * [50,100] ,full and S  * [50,100] ,bins samples.Concretely, S BGS has a number density of ngal = 1.70 × 10 −3 ℎ 3 /Mpc 3 , while S  * [50,100] ,full has a number density of ngal = 2.36 • 10 −5 ℎ 3 /Mpc 3 .Figure 11 compares constraints from S BGS using ΔΣ and  p to constraints from S  * [50,100] ,bins using ΔΣ,  p , and ngal .We find that adding ngal does not compensate for the high signal-to-noise of S BGS .Specifically, the S BGS sample yields Ω m constraints that are 60% tighter than constraints from S  * [50,100] ,bins , and  8 constraints that are 40% tighter than constraints from S  * [50,100] ,bins .
While incomplete samples at face value offer stronger statistical constraining power, it is essential to remember that Figure 11 only compares the statistical constraining power of the two approaches.Systematic errors will play an important role in determining the relative merits of both approaches.For DESI, we do not expect to measure  * , [50,100] for the full BGS sample.Indeed, we would have to adopt a different method of measuring outskirt mass for low-mass galaxies, and it is not clear if outskirt mass would still be a good proxy for low-mass galaxies.The model we have used here is also not optimized for the full BGS sample.For this, we would have to model the galaxy halo connection with a break in the SHMR as we extend the sample to lower masses.In the context of this work, S BGS is a comparison sample that aims to replicate samples that apply the traditional lensing+clustering cosmology (e.g.van Uitert et al. 2018;Joudaki et al. 2018;Abbott et al. 2018;McClintock et al. 2019;Singh et al. 2020;Heymans et al. 2021;Miyatake et al. 2021;Abbott et al. 2022;Amon et al. 2022;Lange et al. 2023), which would give us a chance to understand the statistical importance of completeness on cosmological constraints.In this regard, smaller but mass-complete samples may offer significant advantages, particularly concerning modeling systematics (galaxy halo connection, baryonic effects, assembly bias, satellites, etc.).This will be discussed further in Section 5.2.

Impact of satellite galaxies on parameter constraints
A key difference between traditional cluster cosmology and supermassive galaxy cosmology is that samples binned by galaxy mass may be contaminated with satellites at the ∼10% level (Reid et al. 2016;Saito et al. 2016, etc.).There are multiple ways in which this can be avoided.For instance, one could improve the selection technique of massive galaxy samples to minimize satellite contamination.Alternatively, one could include satellites in the modeling and marginalize over the satellite contribution.However, to begin with, and as a worst-case scenario, here we consider the bias introduced if satellites are ignored altogether.We use a mock galaxy catalog created from a snapshot of MDPL2.This mock catalog was created using a subhalo abundance matching (SHAM, Vale & Ostriker 2004) type approach that matches the HSC stellar mass function derived in Huang et al. 2022. DeMartino et al., in prep will describe this mock catalog in more de-tail.We use this mock galaxy catalog simply to extract semi-realistic ΔΣ and  p profiles for satellite galaxies.We then vary the satellite fraction and evaluate the impact on our cosmological fitting.
Previous work has shown that the satellite fraction among galaxies in this mass range ranges between 8-20%.For example, Reid et al. (2016) find a satellite fraction of 10%, while Saito et al. (2016) find a satellite fraction of 11%.We roughly expect galaxies in the S M * ,bin samples to have a 10% satellite fraction (DeMartino et al., in prep).Here, we choose to study the effect of satellites in our analysis with an upper limit of 20% on the satellite fraction.
Figure 12 displays the effect of satellites on the shape of our key data vectors.We increase the satellite fraction in our selection from 0 to 20%.As the satellite fraction increases, ΔΣ changes mainly in the transition of the one to two halo term.On the other hand,  p is primarily affected in the one halo term.
We create mock data vectors with varying levels of satellite contamination (Figure 12) and analyze these mock data vectors with our fiducial CosmoLike model, which does not include a model for satellites.We constrain the model using the same covariance as in the massive galaxies analysis (Figure 2).The results are shown in Figure 13.Interestingly, even at the 20% level, the bias due to satellites is almost negligible compared to our 1 cosmological constraints.This may seem surprising given Figure 12.We, therefore, analyze why the cosmology parameters are unaffected.
We study how the galaxy-halo parameters are impacted when ignoring satellites.Figure 14 displays the constraints on the galaxyhalo parameters when  sat = 0, 0.1 and 0.2.Clearly, satellite contamination directly impacts the slope  and the −intercept of the mass to observable relation.We conclude that ignoring satellites mainly leads to biases in the galaxy-halo parameters but leaves the cosmology parameters relatively unchanged.In conclusion, while proper satellite modeling is needed to calibrate the mass to observable relation correctly, it is promising to find that satellites do not seem to be a major systematic for super-massive galaxy cosmology.

Massive Galaxies versus Clusters
H+22 showed that the outskirt mass of massive galaxies provides a mass proxy that traces halos with comparable scatter to red sequence richness.Given this information, we have demonstrated that the lensing, clustering, and number densities of massive galaxies constrain the growth of structure with similar statistical power to traditional cluster-finding techniques.This is of interest for several reasons.First, massive galaxies could enable us to push to a lower halo mass range where richness becomes uncertain.Second, by pushing down to lower mass, we find that clustering can be used to tighten constraints.Third, we can avoid cluster-finding systematics that are hard to model.Given the statistical precision of current data sets, systematic issues with traditional cluster-finding techniques have become an increasing concern.For example, Zhang et al. (2019) show that cluster richness estimates tend to be biased low due to miscentering.Furthermore, they indicate that this richness bias may affect cosmology and that future surveys should explicitly take this into account in their cosmological analyses (also see Costanzi et al. 2019;Sunayama et al. 2020).
Tackling systematics due to selection and projection effects is of primary importance for optical cluster-finding techniques.Massive galaxies present an alternative approach that could potentially bypass some of the more difficult systematics.For example, by selecting halos by their outskirt mass, we avoid the issue of projection effects which have proven difficult to model.This is not to say that this alternative approach will be free of systematics.Additional work will be needed to characterize and understand issues associated with this technique.Here, we discuss three effects that will need to be considered.
The first is satellite contamination and its impact on cosmology.Indeed, we expect that super-massive galaxy samples will have a ∼10% satellite fraction (e.g., Li et al. 2014;Reid et al. 2016;Saito et al. 2016;Leauthaud et al. 2017;Kumar et al. 2022).However, Figure 13 is promising and suggests that even a 20% satellite contamination has a less than 1 effect, even with a model that completely ignores satellites.This is because the bias induced by ignoring satellites gets absorbed into the galaxy halo parameters leaving cosmological parameters intact.In addition, in practice, one would attempt to model satellites and not ignore them as we have here.It is plausible that modeling the effects of a 10% satellite fraction could be more straightforward than modeling red galaxies.This is because the emergence of the red sequence depends on several complex galaxy formation processes that are yet to be well characterized.The second potential systematic is the impact of assembly bias.Indeed, it is possible that galaxies selected by outskirt mass might preferentially select older halos and thus be subject to assembly bias (Bradshaw et al. 2020).This effect is being studied using the TNG simulations and will be presented in upcoming work (Xu et al. in prep).Finally, we expect baryonic effects (e.g., Schneider et al. 2019Schneider et al. , 2020a,b;,b;Huang et al. 2021;Shao et al. 2022;Chen et al. 2023) to impact cosmology, especially since we are using both small and large scales in this forecast.However, this systematic is common to both cluster cosmology and super-massive galaxy cosmology, so it needs to be tackled regardless.

Massive Galaxies versus Lensing+Clustering
BAO surveys such as BOSS have collected spectra for millions of massive galaxies (Dawson et al. 2013).The samples of interest are generally incomplete, except at the highest galaxy masses (Leauthaud et al. 2016).This is the main reason why number densities cannot be used as part of the data vector but are instead left to float as a free parameter.The overall number density of a sample is often introduced to scale the normalization of the overall central halo occupation distribution (e.g.Lange et al. 2022;Yuan et al. 2022).However, DESI will be complete in certain mass and redshift range (DESI Collaboration et al. 2016).This allows for restricting analyses to regions of mass and redshift where DESI will be complete.This work analyzes some of the trade-offs in these choices.In Section 4.3, we showed the trade-off between using the full DESI BGS sample compared to smaller  * limited samples.We found that a larger galaxy sample has better constraining power despite not including number density (Figure 11).
However, while at face value, the cosmological constraints from the overall sample outperform the complete samples, our comparisons only consider the statistical constraining power of these different methodologies.When considering systematics, massive  * limited samples may prove advantageous for several reasons: (i) The galaxy halo modeling of  * complete super massive galaxy samples is more straightforward than incomplete samples where the effects of color cuts are poorly understood (e.g., Saito et al. 2016).Therefore,  * complete super massive galaxy samples offer reduced systematic errors associated with the galaxy halo modeling.This may allow smaller radial scales to be used for these samples.
(ii) Along similar lines, the impact of assembly bias should be reduced for  * complete super massive galaxy samples compared to color-selected incomplete samples.
(iii) Baryons are expected to modify the galaxy-dark matter con- . impact of satellites on the mass to observable relation parameters.
Higher satellite fraction translates into higher  and lower −intercept.This strong correlation demonstrates that proper satellite modeling is needed to correctly calibrate the mass to observable relation for massive galaxies.
nection, the dark matter distribution, and the gas distribution, on scales below a few Mpc.Stellar complete super massive galaxy samples offer a more straightforward way of modeling baryonic ef-fects and their dependence on halo mass (e.g., Schneider et al. 2019Schneider et al. , 2020a,b) ,b) than color-selected incomplete samples.Furthermore, the fact that the proposed galaxy samples have a simple selection function means that comparison with hydrodynamic simulations will be more straightforward.
While previous work either uses large scales only or ignores baryonic effects and/or satellite contamination, we expect that masscomplete samples from DESI will be better at simultaneously constraining cosmology, the impact on baryons on the galaxy-dark matter cross-correlation, as well as the halo mass dependence of these effects.Therefore, it could be that the decrease of 30% in constraining power on Ω m (Figure 11) is a worthwhile trade-off for gaining a more simple and accurate galaxy halo modeling as well as constraints on baryonic effects.
Recently, Dvornik et al. (2022) presented constraints on  8 and Ω m using the same observational probes as considered here, namely: ΔΣ,  p , and galaxy number densities using the KiDS survey and covering an area of 1006 deg 2 .Dvornik et al. (2022) use a sample spanning 0.1 <  < 1.3.Their analysis leads to  8 = 0.781 +0.033 −0.029 and Ω m = 0.290 +0.021 −0.017 .Here we discuss some differences between this work and our proposed approach.First, Dvornik et al. (2022) adopt photometric samples, whereas we are assuming spectroscopic samples such as those that will be delivered by DESI, which offer high signal-to-noise measurements of  p .Furthermore, DESI will avoid all systematics associated with photo-'s (for the lens samples).Second, they use total stellar mass as a proxy which has a larger scatter compared to  * , [50,100] (Huang et al. 2022(Huang et al. , 2021)).Third, they push to lower mass than advocated for in this paper.In Figure 13, we showed that the effect of satellites in the higher mass range of interest is negligible, but this may no longer be the case for low-mass samples with higher satellite fractions.In summary, while Dvornik et al. (2022) and this paper are similar in spirit, the proposed implementation details are different.
Finally, there is a possibility of using traditional lensing and clustering cosmology jointly with super-massive cosmology.This would utilize the sensitivity that super-massive galaxies have on  8 and Ω m in the higher mass range (Figure 5) and allow the use of ngal (Figure 8), while also exploiting the high statistical power of the DESI BGS sample (Figure 11).Of course, an accurate understanding of the systematics associated with both samples is needed to reach this point.We will discuss this in more detail in future work.

SUMMARY AND CONCLUSIONS
In this work, we introduced the idea that super-massive galaxies can be used to trace the growth of structure.We studied how cosmological constraints using massive galaxies as halo tracers compare with those from cluster cosmology and traditional lensing plus projected clustering analyses.We used CosmoLike to model our data vector by applying the halo model described in Section 2. We build a covariance matrix for the observables ΔΣ,  p , and n assuming 1000 deg 2 of HSC lensing data with spectroscopic redshifts from DESI.
Our key findings include the following: • We studied the impact of binning and compared how cosmological constraints from three narrow bins in outskirt mass compare to constraints from one cumulative  * , [50,100] bin.We found that binning the data in narrow outskirt mass sub-samples improves constraints on  8 by 34% (Figure 7).This is because  8 mainly impacts the high mass slope of the HMF.Because binning helps to constrain this slope, this translates into tighter constraints on  8 .
• We studied the impact of number density and compared constraints with and without, including ngal in the data vector.Our results show that including ngal as a third observational probe improves constraints on  8 by 33% and Ω m by 23% (Figure 8).
• We compared constraints from a stellar mass limited and mass complete sample to those from a larger but mass incomplete sample (e.g., the DESI BGS).Constraints for a BGS-like sample are tighter in Ω m by 60% and  8 by 40% (Figure 11).However, we note that these forecasts only consider the statistical constraining power of these different methodologies.We present arguments to suggest that stellar mass limited and mass complete samples may offer distinct advantages when considering the inclusion of systematic effects.
• We study the impact of including projected clustering in our data vector.We find that  p strongly impacts our constraining power on Ω m , representing an 84% improvement on Ω m (Figure 9).An analysis using ΔΣ and ngal but omitting clustering is similar to the approach used in cluster cosmology.While  p is not often included as an observational probe in cluster cosmology due to the low number density of clusters, it will be of great utility for super-massive galaxy cosmology where pushing down to lower halo masses may be more straightforward.
• We compare the cosmological constraining power of clusters and massive galaxies as halo tracers.We calibrate a mass to observable relation for each tracer (Figure 1).Assuming the stellar mass and richness bins in H+22, we obtain similar cosmological constraints from both clusters and massive galaxies (Figure 10).This results from similar underlying halo mass distributions of our two samples (Figure 3).
• One of the main caveats of working with massive galaxies as halo tracers will be satellite contamination.We study this effect and find that for the survey parameters assumed, this is less than a 1 effect on cosmological parameters (Figure 13).Instead, we find that satellite contamination is absorbed into the parameters of the mass to observable relation (Figure 14).
In this paper, we have shown that massive galaxies present an excellent avenue for performing precision cosmology.Massive galaxies offer competitive constraints to traditional cluster cosmology and will allow us to bypass systematics associated with cluster-finding systematics, such as miscentering and projection effects which can be hard to model and will bias cosmological constraints.This is especially important in the current era of the  8 tension, which could point to new physics or unaccounted systematics.Super-massive galaxies from DESI will be complete down to lower halo masses than  = 20 cluster samples.By pushing to lower halo masses,  p adds strong constraints.There are some caveats to working with massive galaxies.Satellite contamination is inevitable.However, we have shown that cosmological constraints are robust to the impact of satellites.Assembly bias can introduce another potential systematic.Finally, we must consider baryonic effects, which is a systematic for many low redshift probes of the growth of structure.In conclusion, while there is further work to be carried out to turn super-massive galaxies into a full-fledged cosmological probe, this paper has demonstrated that this approach holds tremendous promise because it can push to lower halo masses while simultaneously avoiding systematics associated with cluster finding.

Figure 1 .
Figure1.The mass to observable relations for massive galaxies (left) and clusters (right).The dashed lines show the mean relations while the darker and lighter colors show the 1 and 2 spread of the relation.The scatter is a combination of intrinsic and measurement error from the observations of H+22.The left-hand panel shows the relation for the galaxy outskirt mass.

Figure 2 .
Figure 2. The assumed covariance matrix for the lensing, projected clustering, and number density measurements across three different  * bins for the S  * [50,100] ,bins samples as measured from the Abacus simulations.

Figure 3 .
Figure3.The underlying halo mass distribution for the three fiducial bins in  * ,[50,100]  (filled) and  (step) from Phase 3000 of AbacusSummit simulation mocks at  = 0.2.Dotted and dashed lines show the mean halo mass of each bin for  * ,[50,100]  and , respectively.The underlying halo mass distribution for massive galaxies and clusters is similar across all three bins.

Figure 4 .
Figure 4. Illustration of the stellar mass function for our three sets of galaxies.The S  * [50,100],bins samples (blue, orange, green) are binned by galaxy mass in a range where DESI will be complete in stellar mass.In this mass range, we can use number densities for constraints, similar to how cluster abundances are used in cluster cosmology.The second sample, S  * [50,100] ,full , is one threshold-limited sample that combines all three bins from S  * [50,100] ,bins .The third sample is S BGS (pink).This sample has a higher number density and extends to lower halo masses but is incomplete in galaxy mass.
8 is steeper than the slope of HMF Ω m .Binned samples will better constrain the slope of HMF  8

Figure 5 .βσFigure 6 .
Figure 5. Sensitivity of the HMF to cosmological parameters.Each panel shows how the HMF changes when we vary Ω m (left) and  8 (right).Ω m changes the overall amplitude of HMF, while  8 affects only the higher mass end.The blue, orange, and green vertical lines indicate the mean values of the stellar mass bins from S  * [50,100] ,bins .The HMF is sensitive to both  8 and Ω m in the halo mass range of interest.

Figure 9 .
Figure9.Constraints for Ω m and  8 with massive galaxies when including (teal) and omitting (gray)  p .Both analyses assume the same survey area (HSC 1000 deg 2 ), fiducial parameters (Table3), and the same binning system (Figure3).This assumes a galaxy sample that selects halos with a similar number density as a  > 9.22 selection.Including  p greatly improves constraints on Ω m , with an improvement of 84%.This figure presents one of the key potential advantages of massive galaxies over clusters: the ability to push to lower halo masses where  p is better-measured.
Figure10.Constraints for Ω m and  8 with massive galaxies (teal) and clusters (pink) as halo tracers.Both analyses assume the same survey area (HSC 1000 deg 2 ) and fiducial data vector (Table3).Massive galaxies yield competitive constraints while potentially avoiding systematics such as projection effects and miscentering.

Figure 12 .Figure 13 .
Figure12.Changes in the shape of ΔΣ (left) and  p (right) due to satellite contamination.Higher satellite fraction impact ΔΣ in the transition between the one and two halo terms, while  p in the one halo term.This impact becomes more dominant as  sat increases both for ΔΣ and  p .Predicted HSC error bars are overlaid in grey.
Constraints for Ω m and  8 with massive galaxies when considering a single cumulative bin (navy) and three narrow bins (teal).Both analyses assume the same survey area (HSC 1000 deg 2 ) and the same fiducial data vector (Table3).The maroon dot and lines here and throughout the paper show the fiducial values of Ω m and  8 .Three narrow bins yield 34% tighter constraints on  8 .Constraints for Ω m and  8 with massive galaxies when including (teal) or omitting (orange) ngal .Both analyses assume the same sample (S  * [50,100] ,bins ), survey area (HSC 1000 deg 2 ), and fiducial parameters (Table3).Including ngal improves constraints of  8 by 33% and Ω m by 23%.
Constraints for Ω m  with massive galaxies when S  * [50,100] ,bins (teal) and S BGS (red).The S  * [50,100] ,bins analysis includes all three ΔΣ,  p , and ngal as data vectors, while the S BGS includes only ΔΣ and  p .Both analyses assume survey area (HSC 1000 deg 2 ) and fiducial parameters (Table3).The S BGS sample yields tighter constraints for  8 and Ω m .Specifically, it improves  8 constraints by 40% and Ω m constraints by 60%.