We report evidence of galaxy assembly biasthe correlation between galaxy properties and biased secondary halo properties at fixed halo mass M SUBHSUBin the stellartohalo mass relation for red central galaxies from the Sloan Digital Sky Survey. In the M SUBHSUB  10SUP11.5SUP10SUP13.5SUP h SUP1SUP M SUBSUB range central galaxy stellar mass M SUBSUB is correlated with the number density of galaxies within 10 h SUP1SUP Mpc  SUB10SUB a common proxy for halo formation time. This galaxy assembly bias signal is also present when M SUBHSUB M SUBSUB and  SUB10SUB are substituted with group luminosity galaxy luminosity and metrics of the largescale density field. To associate differences in  SUB10SUB with variations in halo formation time we fitted a model that accounts for 1 errors in the M SUBHSUB measured by the J. L. Tinker group catalog and 2 the level of correlation between halo formation time and M SUBSUB at fixed M SUBHSUB. Fitting of this model yields that 1 errors in M SUBHSUB are 0.15 dex and 2 halo formation time and M SUBSUB are strongly correlated Spearmans rank correlation coefficient 0.85. At fixed M SUBHSUB variations of 0.4 dex in M SUBSUB are associated with 13 Gyr variations in halo formation time and galaxy formation time from stellar population fitting. These results are indicative that halo properties other than M SUBHSUB can impact central galaxy assembly.
['Galaxies', 'Extragalactic astronomy', 'Early-type galaxies', 'Large-scale structure of the universe', 'Observational cosmology', 'Dark matter', 'Galaxy dark matter halos', '573', '506', '429', '902', '1146', '353', '1880', 'Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Cosmology and Nongalactic Astrophysics']
Galaxy Assembly Bias in the Stellartohalo Mass Relation for Red Central Galaxies from SDSS
{'Galaxy assembly bias in the stellar-to-halo mass relation for red central galaxies from SDSS': "Grecco A. Oyarz'un, 1 Jeremy L. Tinker, 2 Kevin Bundy, 3 Enia Xhakaj, 4 and J. Stuart B. Wyithe 5, 6  \n1 The William H. Miller III Department of Physics & Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA 2 Center for Cosmology and Particle Physics, Department of Physics, New York University, New York, USA, 10003 3 University of California Observatories - Lick Observatory, University of California, Santa Cruz, CA 95064, USA 4 High-Energy Physics Division, Argonne National Laboratory, 9700 S Cass Ave, Lemont, IL 60439, USA  \n- 5 Research School of School of Astronomy and Astrophysics, Australian National University, Mt Stromlo, Act 2611, Australia 6 ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Australia  \n(Received May 20, 2024; Revised July 17, 2024; Accepted August 8, 2024)", 'ABSTRACT': "We report evidence of galaxy assembly bias - the correlation between galaxy properties and biased secondary halo properties at fixed halo mass (M H ) - in the stellar-to-halo mass relation (SHMR) for red central galaxies from the Sloan Digital Sky Survey. In the M H = 10 11 . 5 -10 13 . 5 h -1 M ⊙ range, central galaxy stellar mass (M ∗ ) is correlated with the number density of galaxies within 10 h -1 Mpc ( δ 10 ), a common proxy for halo formation time. This galaxy assembly bias signal is also present when M H , M ∗ , and δ 10 are substituted with group luminosity, galaxy luminosity, and metrics of the largescale density field. To associate differences in δ 10 with variations in halo formation time, we fitted a model that accounts for (1) errors in the M H measured by the Tinker (2021, 2022) group catalog and (2) the level of correlation between halo formation time and M ∗ at fixed M H . Fitting of this model yields that (1) errors in M H are ∼ 0 . 15 dex and (2) halo formation time and M ∗ are strongly correlated (Spearman's rank correlation coefficient ∼ 0 . 85). At fixed M H , variations of ∼ 0 . 4 dex in M ∗ are associated with ∼ 1 -3 Gyr variations in halo formation time and in galaxy formation time (from stellar population fitting; Oyarz'un et al. 2022). These results are indicative that halo properties other than M H can impact central galaxy assembly.", '1. INTRODUCTION': "In the current cosmological picture, the growth of dark-matter halos and the formation of galaxies are closely connected (e.g. White & Rees 1978; Blumenthal et al. 1984; Mo & White 1996; Wechsler & Tinker 2018; Vogelsberger et al. 2020). As halo mass (M H ) increases, so does the gravitational pull of the halo, facilitating the accretion of gas, fueling star-formation, and promoting the build-up of stellar mass (M ∗ ). The resulting correlation between M H and M ∗ is known as the stellar-to-halo mass relation of central galaxies (SHMR; e.g. Moster et al. 2013).  \nWhile M H emerges as the primary halo property driving galaxy growth, simulations of galaxy formation indicate that halo properties other than M H also play a role. One of these properties is thought to be halo formation time (e.g. Gao et al. 2005; Tojeiro et al. 2017). Varia-  \ngoyarzu1@jhu.edu  \ntions in halo assembly history lead to variations in halo gravitational potential with time, ultimately impacting the M ∗ assembly history of the galaxies that they host (e.g. Wechsler et al. 2002, 2006; Gao et al. 2005; Li et al. 2008). Similarly, halos of different concentrations feature different radial dark-matter density profiles, affecting the assembly histories of halos and galaxies alike, according to simulations (Ludlow et al. 2014; Matthee et al. 2017; Xu & Zheng 2020). These correlations between galaxy properties and biased secondary halo properties at fixed M H that are apparent in theoretical work are known as galaxy assembly bias (e.g. Croton et al. 2007; Mao et al. 2017; Zehavi et al. 2018; Hadzhiyska et al. 2021, 2023; Montero-Dorta & Rodriguez 2024).  \nThe major role played by dark-matter halos in simulations of galaxy formation has encouraged observational studies to search for signatures of galaxy-halo co-evolution. For instance, Obuljen et al. (2020) found the power spectrum of galaxies in the Baryon Oscillation Spectroscopic Survey of SDSS-III (BOSS; Dawson  \net al. 2013) to vary with galaxy central stellar velocity dispersion. These variations - of the order of 5 σ -appear at spatial scales that extend beyond the virial radii of dark-matter halos, indicating that galaxy properties correspond with the distribution of halos in the large scale structure. In a separate study, Wang et al. (2022) found that the spatial distribution of galaxies in the Sloan Digital Sky Survey (SDSS; Alam et al. 2015) was best reproduced with halo occupation models (e.g. Tinker et al. 2008) that account for galaxy assembly bias in halo concentration. The same conclusion was reached by Yuan et al. (2021) in BOSS and by Pearl et al. (2023) in the Dark Energy Spectroscopic Instrument (DESI) One-Percent Survey (DESI Collaboration et al. 2016). Taken together, these recent studies suggest that galaxy assembly bias is apparent in the spatial distribution of galaxies in the nearby Universe.  \nDespite the mounting observational evidence in support of galaxy assembly bias, the connection between secondary halo properties and galaxy properties (e.g. stellar mass or star-formation rate; Kauffmann 2015) remains unconstrained due to challenges affecting observational studies of the galaxy-halo connection. For instance, signatures of one-halo conformity, i.e., the similarity in the observable properties of galaxies within a halo (Weinmann et al. 2006), can be interpreted as evidence of galaxy assembly bias. However, because these signatures can also be explained by the co-evolution between central and satellite galaxies (e.g. Pasquali et al. 2010; Peng et al. 2012; Wetzel et al. 2013; Cortese et al. 2019; Gallazzi et al. 2020; Oyarz'un et al. 2023) and by systematic errors (e.g. Calderon et al. 2018), using one-halo conformity as a probe for galaxy assembly bias can be difficult. While other work (e.g. Kauffmann et al. 2013) have attempted to mitigate these degeneracies by searching for signatures of galaxy assembly bias in two-halo conformity instead (the similarity between the physical properties of galaxies in neighboring halos; e.g. Lacerna et al. 2022), biases in the isolation criterion emerge as a major concern (Sin et al. 2017; Tinker et al. 2018).  \nInstead of studying the impact of galaxy assembly bias on galaxy formation through galactic conformity, other work have turned to the clustering of galaxies. Studies have reported variations in the clustering signal with galaxy environment (Lehmann et al. 2017; Zentner et al. 2019) and/or galaxy color (Montero-Dorta et al. 2017). However, the reliability of these results has been questioned in light of systematic errors affecting the photometric group catalogs used for estimating halo masses (Campbell et al. 2015; Lin et al. 2016).  \nWith systematic errors hampering observational studies of the galaxy-halo connection, recent efforts have focused on mitigating some of these shortcomings. For instance, the group catalog by Tinker (2021, 2022) is calibrated with observations of color-dependent galaxy clustering and estimates of the total satellite luminosity, which allows for a more accurate reproduction of the central galaxy quenched fraction. In addition, Tinker (2021) used deep photometry from the DESI Legacy Imaging Survey (DLIS; Dey et al. 2019), allowing for precise accounting of galaxy M ∗ and total group luminosity (L group ).  \nFor instance, in Oyarz'un et al. (2022) we exploited the Tinker (2021) catalog to estimate the M H for over 2000 central galaxies in the SDSS-IV MaNGA survey (Mapping Nearby Galaxies at Apache Point Observatory; Bundy et al. 2015) and identified variations in the stellar populations of central galaxies at fixed M H . Specifically, central galaxies with high M ∗ for their M H feature old and Mg-enhanced stellar components. These results, along with equivalent trends recovered in spatially unresolved spectroscopy from the SDSS (ScholzD'ıaz et al. 2022, 2023, 2024), are an indication that variations in the assembly histories of central galaxies cannot be fully ascribed to variations in M H , thus hinting that halo properties other than M H may play a role in central galaxy formation.  \nIn this paper, we test whether these variations in the observable properties of central galaxies at fixed M H can be attributed to galaxy assembly bias. To this end, we quantify how halo formation time varies within the SHMR for ∼ 150,000 red central galaxies from the SDSS main Galaxy Survey (MGS; York et al. 2000; Gunn et al. 2006; Alam et al. 2015) by utilizing the number density of galaxies within a radius of R = 10 h -1 Mpc (Alpaslan &Tinker 2020). The association between halo formation time and the number density of galaxies is motivated by how higher density regions of the large scale structure form at earlier cosmic epochs (e.g. Hahn et al. 2007; Lee et al. 2017).  \nThis manuscript is organized as follows. We present our dataset in Section 2, methodology in Section 3, and results in Section 4. We discuss the implications of our work in Section 5 and summarize in Section 6. This paper makes use of the spectroscopic M ∗ measured by Chen et al. (2012) and the M H estimated by Tinker (2021). The M ∗ adopt a Kroupa (2001) IMF and h = 0 . 7. All magnitudes are reported in the AB system (Oke & Gunn 1983).", '2.1. Galaxy sample': "To study the SHMR for central galaxies in the nearby Universe, we turned to the largest sample of nearby galaxies publicly available, the SDSS main Galaxy Survey sample (MGS; York et al. 2000; Gunn et al. 2006; Alam et al. 2015). In particular, we used the dr72bright34 subsample of the MGS, which was defined in the construction of the NYU Value-Added Galaxy Catalog (NYU-VAC; Blanton et al. 2005). Spectroscopic and photometric data outputs are available as part of the eighth data release (DR8) of the SDSS (Brinchmann et al. 2004). To limit biases in the survey selection function, we restricted our sample of central galaxies to the redshift range 0 . 04 < z < 0 . 16, yielding 434,811 galaxies. To identify central galaxies in the MGS, we turned to the halo-based galaxy group catalog for SDSS by Tinker (2021). We selected all galaxies with satellite probabilities P sat < 0 . 1, resulting in 299,409 central galaxies. Further details on the group finder algorithm that we used to identify central galaxies are presented in Section 2.2.  \nBecause group catalogs can imprint systematic errors into the M H distributions of red and blue galaxies, separating the two populations is recommended (Lin et al. 2016). To separate red and blue galaxies in the MGS, we followed the Gaussian-mixture model employed by Tinker (2021). This model operates in the space composed of the 4000 ˚ A break index (D n 4000) and galaxy luminosity (L gal ). The outputs of the model were fitted, yielding a threshold of  \nD crit n 4000(L gal ) = 1 . 42 + 0 . 35 2 [ 1 + erf ( log L gal -9 . 9 0 . 8 )] (1)  \nto separate the two populations. Here, erf(x) is the error function. We refer the reader to Tinker (2021) for more details on how the threshold D crit n 4000(L gal ) was fitted to SDSS data.  \nThis characterization of D crit n 4000(L gal ) was used to identify and remove all blue central galaxies. The resulting sample contains 142,292 red central galaxies, which constitutes the working sample of this study. A sample of red central galaxies has two benefits over a sample of blue central galaxies. First, the majority of centrals are red at the high halo mass end (e.g. Yang et al. 2007), which is where the halo masses are the most reliable (e.g. Tinker 2022). Second, our prior characterization of the assembly histories of centrals within the SHMR focused on a passive sample (Oyarz'un et al. 2022), which better corresponds with a red sample.", 'SHMR for\\_red centrals from SDSS': 'Figure 1. Stellar-to-halo mass relation for red central galaxies from SDSS. The M ∗ were measured in spectra from the SDSS by Chen et al. (2012) and the M H were estimated with the group-finding algorithm by Tinker (2022) on deep photometry from the DLIS (Tinker 2021). We used the 33rd and 66th percentiles in M ∗ -to-M H ratio to divide centrals into high M ∗ (magenta) and low M ∗ (green) subsamples.  \n<!-- image -->', '2.2. Halo masses': 'For the identification of central galaxies and to obtain estimates of their M H , we utilized the halo-based galaxy group finder by Tinker (2022). In this algorithm, the probability of a galaxy being a central depends on both galaxy type and luminosity, enabling accurate reproduction of the central galaxy quenched fraction. After central galaxies and their corresponding groups are determined, halo masses are assigned through abundance matching with the Bolshoi-Planck simulation (Klypin et al. 2016). At this point, color-dependent galaxy clustering and total satellite luminosity are utilized as observational constraints until the best fitting model is found.  \nFor our analysis, we worked with the publicly available implementation of the Tinker (2022) algorithm on SDSS 1 (Tinker 2021). To construct the catalog, Tinker (2021) used deep photometry from the DESI Legacy Imaging Survey (DLIS; Dey et al. 2019), allowing for more accurate accounting of galaxy M ∗ and group L ∗ than in other group catalogs. The M H for our sample of red central galaxies is shown in Figure 1.  \n<!-- image -->  \n<!-- image -->  \nFigure 2. Left: The top panel shows the stellar-to-halo mass relation for red central galaxies from SDSS (same as Figure 1). The bottom panel shows the number density of galaxies within a radius R = 10 h -1 Mpc as a function of halo mass. High M ∗ centrals have higher galaxy number densities than low M ∗ centrals. Middle: Same as the left panel but with galaxy r-band luminosity instead of stellar mass and total group r-band luminosity instead of halo mass. Right: Same as the left panel, albeit with galaxy number density substituted with a metric of the density field measured with the MCPCM Slime Mold reconstruction of the Cosmic Web (Burchett et al. 2020). Apparent in all three panels is how galaxy stellar mass (or luminosity) varies with proxies for the density of the large scale structure at fixed halo mass (or at fixed group luminosity).  \n<!-- image -->', '2.3. Stellar masses': "The M ∗ were retrieved from the SDSS DR8 outputs. They were measured by Chen et al. (2012) as part of a PCA-based approach to constrain the star-formation histories and stellar populations for galaxies from SDSS DR7 and BOSS (Aihara et al. 2011). The small scatter at fixed halo mass shown by these stellar masses ( ∼ 0 . 18 dex) make them ideal for this work (Tinker et al. 2017). To measure their M ∗ , Chen et al. (2012) modeled the optical wavelength range between 3700 ˚ A and 5500 ˚ A of every spectrum with a linear combination of seven eigenspectra. The model space was constructed by generating a spectral library spanning a wide variety of star-formation histories. In this work, we retrieved the M ∗ that used the Maraston & Stromback (2011) stellar population synthesis models and Kroupa (2001) IMF for the construction of the library.  \nTo determine whether some of the scatter in the red SHMRis driven by secondary halo properties (e.g. Zent-  \nner et al. 2014; Xu & Zheng 2020), we defined two subsamples of red central galaxies at fixed M H : high-M ∗ and low-M ∗ centrals. These samples were computed as the 33rd and 66th percentiles in M ∗ -to-M H ratio, in similar fashion to Oyarz'un et al. (2022). In equation form,  \nHigh-M ∗ sample: M ∗ M H (M H ) > [ M ∗ M H (M H ) ] . 66 Low-M ∗ sample: M ∗ M H (M H ) < [ M ∗ M H (M H ) ] . 33 . (2)  \nThe two subsamples are plotted in Figure 1. High-M ∗ centrals are plotted in magenta and low-M ∗ centrals are plotted in green.", '3. METHODOLOGY': "3.1. Galaxy number densities and halo formation times  \nIn ΛCDM, higher density regions of the large scale structure collapse at earlier cosmic epochs (Lemson &  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 3. Model for interpreting variations in the number density of galaxies with M ∗ at fixed M H . The free parameters of the model are the random error in the halo mass estimates from the group catalog ( σ log M H ) and the Spearman's rank correlation coefficient ( r s ) between t halo and M ∗ at fixed M H (i.e., the galaxy assembly bias signal). Left: The SHMR for red central galaxies from SDSS for different values of σ log M H (top). Differences in ∆ log(1 + δ 10 ) between high-M ∗ and low-M ∗ centrals for the different values of σ log M H (bottom). Right: The correlation between t halo and M ∗ for different values of r s in the M H = 10 12 . 25 -10 12 . 5 h -1 M ⊙ bin (top). Variations in ∆ log(1 + δ 10 ) between high-M ∗ and low-M ∗ centrals for the different values of r s (bottom). We note how σ log M H and r s are only partially degenerate, allowing us to simultaneously quantify their impact on ∆log(1 + δ 10 ).  \n<!-- image -->  \nKauffmann 1999). We can take advantage of this property of ΛCDM to use the normalized number density of galaxies within a large radius R as our observational proxy for halo formation time.  \nIn principle, any radius larger than the sizes of halos (i.e., R ≳ 1 h -1 Mpc) should suffice. However, scales R ≲ 5 h -1 Mpc are not ideal because splashback subhalos can magnify the halo assembly bias signal (e.g. Sunayama et al. 2016). Moreover, our ability to decouple the effects of M H and halo formation time on central galaxy M ∗ improves as R increases. For these reasons, R ≳ 10 h -1 Mpc are preferred (Mansfield & Kravtsov 2020). In this work, we opted for R = 10 h -1 Mpc to keep the computation time manageable.  \nThe number density of galaxies within R = 10 h -1 Mpc (hereafter δ 10 ) was computed for every central as  \nδ 10 = ρ z ⟨ ρ z ⟩ -1 , (3)  \nwhere ρ z is the number density of galaxies around a top-hat sphere of radius R = 10 h -1 Mpc and ⟨ ρ z ⟩ is the average ρ z at the redshift of the central. These two  \nquantities were estimated in simulations that randomized the spatial distribution of galaxies, and they account for incompleteness and for edge effects. Further details on how δ 10 was computed can be found in Alpaslan & Tinker (2020). The average value of δ 10 as a function of M H for our red central galaxies is shown in Figure 2.  \nTo connect galaxy number density and halo formation time, we turned to the MultiDark Planck 2 Simulation (MDPL2; Prada et al. 2012). We first measured M H , δ 10 , and t halo for every host halo in the simulation, where t halo is the lookback time at which the host halo assembled half of its total mass at z = 0 (e.g. Tojeiro et al. 2017). Although all halo formation times reported in this manuscript were obtained according to the aforementioned definition, the main conclusions of this paper also stand if we instead use t 11 . 5 halo , which corresponds to the lookback time at which a critical halo mass of M H = 10 11 . 5 h -1 M ⊙ was reached (more in Section 4).  \nFor the measurement of the number density of galaxies in MDPL2, only subhalos with V peak > 100 kms -1 were counted, where V peak is the maximum peak circu-  \nlar rotation velocity across all redshifts. This choice can be interpreted as counting only subhalos with M H ≳ 10 11 h -1 M ⊙ (Kravtsov et al. 2004; Zehavi et al. 2019). Then, we ensured that the dynamic range of δ 10 in the data and in the simulation were comparable. To this end, we implemented an abundance matching approach, i.e., the halos were rank ordered in δ 10 at fixed M H and their δ 10 values were matched.", '3.2. Galaxy assembly bias model': "Here, we describe the analytical model we used to associate differences in δ 10 at fixed M H with variations in halo formation time. Our model has two parameters: the error in the M H reported by the group catalog and the strength of the correlation between halo formation time and stellar mass at fixed M H .  \nImprecise halo mass estimates can broaden the observed scatter of the SHMR, potentially biasing our high-M ∗ and low-M ∗ subsamples (see Figure 1). To account for this effect, we added random perturbations to the measured halo masses prior to the sample assignment step (Section 2.3). These random perturbations were drawn from a normal distribution characterized by a scatter with magnitude ' σ log M H '. The impact that different values of σ log M H have on the galaxy assembly bias signal is shown in the left panels of Figure 3. It is apparent in this figure that errors in the group catalog broaden the scatter of the SHMR and induce variations in δ 10 toward M H > 10 12 . 5 h -1 M ⊙ .  \nThe second free parameter in our model is the strength of the galaxy assembly bias signal at a given M H , which can be quantitatively expressed via the Spearman's rank correlation coefficient ( r s ). In the case of halo formation time and galaxy stellar mass, r s can be written as:  \nr s (M h ) = 1 -6 ∑ [ r (M ∗ (M h )) -r (t halo (M h ))] 2 n (M h )( n 2 (M h ) -1) . (4)  \nHere, r ( X ) symbolizes the ordered rank of the galaxies in variable X , which in our case takes the values X =M ∗ and X = t halo at the corresponding M H bin. The total number of galaxies in the M H bin is represented by ' n (M H )'. If r s = 1, then M ∗ and t halo are perfectly correlated (i.e., the perfect galaxy assembly bias case). In the case of no galaxy assembly bias, r s = 0. Values in the range r s = (0 , 1) correspond to different strengths of the galaxy assembly bias signal, as illustrated in the right panels of Figure 3.  \nIn the context of our work, we had at our disposal a collection of host halos with (t halo , δ 10 ) measured in MDPL2 and a collection of galaxies with (M ∗ , δ 10 ) measured in SDSS. The correspondence between M ∗ and t halo at fixed M H is determined by the value of r s ,  \nwhich in turn dictates how host halos and galaxies are mapped onto each other. The mapping between host halos and galaxies was performed in narrow M H bins, with a unique value for r s at all M H (i.e., we assume r s not to vary with M H ). We note that blue central galaxies were included in this process to ensure appropriate matching between the populations of halos and galaxies.", '3.3. Model fitting': "We used our two-dimensional ( σ log M H , r s ) model to fit the observed variations in log(1 + δ 10 ) at fixed M H . These variations - hereafter ∆ log(1+ δ 10 ) - were computed relative to the median value for log(1 + δ 10 ) in a given M i H bin.  \nHigh-M ∗ sample: ∆log(1 + δ 10 (high-M ∗ , M H )) = log(1 + δ 10 (high-M ∗ , M i H )) -Median[log(1 + δ 10 (M ∗ , M i ))] .  \nLow-M ∗ sample: ∆log(1 + δ 10 (low-M ∗ , M H )) = log(1 + δ 10 (low-M ∗ , M i H )) -Median[log(1 + δ 10 (M ∗ , M i ))]  \ni H i H .  \nThe uncertainties on  \n∆log(1 + δ 10 (high-M ∗ , M i H )) and ∆log(1 + δ 10 (low-M ∗ , M i H ))  \n(5)  \nwere propagated from the uncertainties on  \nδ 10 (high-M ∗ , M i H ) and δ 10 (low-M ∗ , M i H ) , (6)  \nwhich were computed by bootstrapping the galaxies selected for the computation of δ 10 . The errors on the average M H of every bin were not fitted for, and the different M H bins were assumed independent.  \nThis set of ∆ log(1 + δ 10 ) measurements as a function of M H compose the dataset that was fitted with the galaxy assembly bias model. Explicitly, the quantities that we fitted were  \n∆log(1 + δ 10 (high-M ∗ , M i H )) and (7)  \n∆log(1 +  \nδ  \n10  \n(low-M  \n∗  \n,  \nM  \ni  \nH  \n))  \n,  \nas they appear in the bottom panels of Figure 3. We highlight that ∆ log(1+ δ 10 ) is only explicitly dependent on three quantities: M H , M ∗ , and δ 10 . In observational space, these quantities were measured in SDSS. In model space, the M H correspond to the M H in MDPL2, the M ∗ were assigned according to the value of r s , and the t halo were measured in MDPL2 (see Sections 3.1 and  \nFigure 4. Posterior distribution (1 σ , 2 σ , 3 σ contours) of the model that fits for variations in galaxy number density within the SHMR for red central galaxies. The x-axis shows the random error in the halo mass estimates from the group catalog. The y-axis shows the Spearman's rank correlation coefficient between t halo and M ∗ at fixed M H (i.e., the galaxy assembly bias signal). The no-correlation case corresponds to r s = 0 and the perfect correlation case corresponds to r s = 1. The absence of assembly bias ( r s = 0) is inconsistent with the data.  \n<!-- image -->  \n3.2). Although the δ 10 in model space were measured in MDPL2, they were abundance matched to the values in SDSS (end of Section 3.1).  \nIn the computation of the posterior, we adopted uniform, bounded priors for both parameters, i.e.,  \np ( σ log M H ) = TopHat [0 , 1] p ( r s ) = TopHat [0 , 1] . (8)  \nThe Markov chain Monte Carlo sampler emcee (Foreman-Mackey et al. 2013) was used to constrain the posterior. We implemented a setup of 10,000 walkers, 2400-step chains, and 200 burn-in steps.", '4. RESULTS': "Plotted in the top left panel of Figure 2 is the SHMR for red central galaxies from SDSS separated into two bins of M ∗ at fixed M H (Section 3.1). The bottom left panel shows that centrals with high-M ∗ have systematically higher galaxy number densities within R = 10 h -1 Mpc than centrals with low-M ∗ , especially in the M H range between 10 11 . 5 -10 13 . 5 h -1 M ⊙ . This result is evidence that the M ∗ of galaxies vary with the density of the large scale structure at fixed M H .  \nWe can confirm that this result holds true after controlling for all three quantities involved in the measure-  \nment. Upon substitution of M H with group r-band luminosity and M ∗ with central galaxy r-band luminosity, differences in δ 10 remain significant (middle panel of Figure 2). We also controlled for δ 10 by utilizing reconstructions of the Cosmic Web with the Monte Carlo Physarum Machine (MCPM) Slime Mold algorithm (Burchett et al. 2020). In contrast to the number density of galaxies, which is sensitive to the abundance of luminous matter in the large scale structure, the MCPM Slime Mold catalog (Wilde et al. 2023) provides estimates for the density field in its totality, including dark-matter (Hasan et al. 2023a). This was achieved through the adaptation of an algorithm designed to characterize the growth of 'slime mold', enabling reconstruction of the tri-dimensional structure of the Cosmic Web. We used the publicly available 2 measurements in replacement of the number density of galaxies, achieving a consistent - and perhaps even more significant - signal in the right panel of Figure 2. Taken together, these results are evidence that the baryonic properties of central galaxies correlate with the density of the large scale structure after controlling for the mass at the halo scale.  \nWhile differences in δ 10 are suggestive of differences in halo formation time (e.g. Lee et al. 2017; Alpaslan & Tinker 2020), modeling is required in order to account for other possible sources of variation in δ 10 . The framework introduced in Section 3.2 serves this purpose by simultaneously modeling the effects of M H errors and galaxy assembly bias on δ 10 variations within the red SHMR. We fitted this model to the observed differences in δ 10 (left panel of Figure 2) and obtained the posterior distribution shown in Figure 4. Inspection of this figure reveals that random errors in the halo mass estimates and a galaxy assembly bias signal are both needed to explain the differences in δ 10 , with σ log M H ∼ 0 . 15 dex and r s ∼ 0 . 85 best reproducing the observations. Errors in the M H of this magnitude are within expectations from the group finder algorithm ( σ log M H ≲ 0 . 2; Tinker 2022) and the marginal posterior probability for r s rejects the absence of galaxy assembly bias ( r s = 0) with confidence > 5 σ .  \nThe model solutions are plotted along with the data in Figure 5. As revealed by the middle panel, the differences in galaxy number density at fixed M H are well reproduced by the converged model. Plotted in the right panel are the the halo formation times of low- and highM ∗ red centrals that have been inferred by our model.  \n<!-- image -->  \n<!-- image -->  \nFigure 5. Model solution after convergence. Left: The SHMR for red central galaxies from SDSS highlighting the two subsamples (same as Figure 1). High-M ∗ centrals are plotted in magenta ( > 66th percentile in the SHMR) and low-M ∗ centrals are plotted in green ( < 33rd percentile in the SHMR). Middle: Difference in galaxy number density within R = 10 h -1 Mpc between high- and low-M ∗ red centrals as a function of M H (datapoints). The shaded regions show the 1 σ difference in galaxy number density from a random sampling of the converged posterior. Right: 1 σ intervals of the t halo of high-M ∗ and low-M ∗ red centrals as constrained by the model. Also shown in black is the median t halo of host halos in the MDPL2 simulation as a function of M H . Note that the average between the t halo of high-M ∗ and low-M ∗ red centrals differs from the median t halo in MDPL2 because the median t halo was computed for all centrals (i.e., both red and blue).  \n<!-- image -->  \nMore massive dark-matter halos form at later times, and halos that are old for their M H host galaxies with higher M ∗ . We can confirm that this result also holds true for different definitions of halo formation time: at fixed M H , halos with older t 11 . 5 halo (as defined in Section 3.1) host higher M ∗ galaxies.", '5. DISCUSSION': "Variations in the baryonic properties of galaxies at fixed M H have been observed for a few years (e.g. Oyarz'un et al. 2022; Scholz-D'ıaz et al. 2022, 2023, 2024). However, diagnosing the source of these variations has been difficult due to the stochastic nature of galaxy formation and because of the systematic errors affecting observational halo diagnostics (e.g. Calderon et al. 2018; Tinker et al. 2018). In this context, we have determined that some of these variations at fixed M H , in particular differences in stellar mass and/or galaxy luminosity, are correlated with how luminous matter is arranged in large ( R = 10 h -1 Mpc) scales. In this section, we discuss the biases and/or the underlying physical processes that could be driving this connection.", '5.1. Systematic errors in the group catalog': "It is important to address how systematic errors in our approach could fabricate a signal in the galaxy number densities. Errors in the halo mass estimates from the group catalog are particularly relevant, which is the reason for their inclusion in the model. As we show in Figure 4, random errors in the group catalog alone can-  \nnot explain the differences in galaxy number densities at fixed M H . That said, should the errors in the group finder not be random (i.e., stellar- or halo-mass dependent), then any signal in the galaxy number densities could be fabricated. For instance, weak-lensing analyses have revealed differences between the underlying M H distributions of red and blue central galaxies in the Yang et al. (2007) group catalog, indicating that the scatter of the SHMR correlates with the physical properties of galaxies in an unaccounted for manner (Lin et al. 2016).  \nOur approach accounts for these difficulties in multiple ways. First, our detection is not reliant on differences between red and blue centrals, since the latter were removed from our sample. Second, the impact of the color-dependent spatial distribution of galaxies on halo mass estimation is accounted for, as the Tinker (2022) algorithm is designed to fit for the clustering of red and blue galaxies separately. It is also reassuring that the SDSS implementation of the algorithm can reproduce the weak-lensing derived SHMR (Mandelbaum et al. 2016), especially for red centrals in our M H regime (Tinker 2021). Finally, the fact that we also recover the signal when using total group luminosity instead of M H (Figure 2) is an indication that the signal is not fabricated by how the halo masses are assigned by the group finder algorithm.  \nSystematic errors in group catalogs can also stem from central and satellite galaxy mis-classification. Lin et al. (2016) showed that satellite contamination can affect the correlation functions and the dark matter profiles of  \nFigure 6. Variations in halo (this work) and galaxy (Oyarz'un et al. 2022) formation times within the SHMR for red central galaxies. High- and low-M ∗ centrals were defined according to the 33rd and 66th percentiles in M ∗ -to-M H ratio in the red SHMR (magenta and green, respectively). The differences for both high- and low-M ∗ centrals are relative to the median formation time at a given M H . The differences in halo formation time correspond to the 1 σ spread of a random posterior sampling (shaded regions). The differences in galaxy formation time correspond to the average massweighted stellar ages derived from stellar population fitting of MaNGA spectra (datapoints). Differences in the formation times are ∼ 1 -3 Gyr for both halos and galaxies.  \n<!-- image -->  \ncentral galaxy samples. Yet, as also pointed out by Lin et al. (2016), satellites that are mis-classified as centrals reflect the large-scale bias of their higher M H hosts. As the normalized densities converge at high M H (Figure 3), central and satellite galaxy mis-classification would weaken the assembly bias signal; a conclusion that we can confirm applies to our dataset after conducting simulations of galaxy miss-classifications. Moreover, adopting a central selection criterion of P sat < 0 . 01 instead of P sat < 0 . 1 does not impact our results.  \nAny systematic errors arising from our M ∗ are unlikely. At fixed M H , the measurement error of our stellar masses is ∆M ∗ < 0 . 18 dex (Tinker et al. 2017), much smaller than the ∼ 0 . 4 dex separating our two M ∗ subsamples (left panel of Figure 1). Moreover, the signal is also recovered after substituting M ∗ with L ∗ (Figure 2). The normalized number densities are similarly robust, with the flux limit of the DLIS reaching 22 . 8 mag at 5 σ (Dey et al. 2019), which is 5 -6 magnitudes deeper than the m ∗ of the galaxy luminosity function at z = 0 . 1 (Blanton et al. 2005). We also confirm that the signal is recovered with independent metrics for the density of the  \nlarge scale structure (i.e., the Slime Mold reconstruction of the Cosmic Web; Figure 2).  \nStill, quantifying the magnitude and significance of galaxy assembly bias with other probes for M H remains a priority. Weak lensing can be used to search for systematic differences between the total mass surface density profiles of high- and low-M ∗ centrals (e.g. Lin et al. 2016). In addition, host M H can also be estimated through the tight correlation between central M ∗ , M H , and M max ∗ (where M max ∗ is the the stellar mass of the central out to 100 kpc; Huang et al. 2020) or via the stellar-to-total dynamical mass relation for central galaxies (Scholz-D'ıaz et al. 2024). Also key will be alternative methods for estimating halo formation times in observations (e.g. machine learning algorithms; Lyu et al. 2024).", '5.2. Analytic model for galaxy assembly bias': "The magnitude of the difference in halo formation time between the two M ∗ -selected subsamples that has been inferred by our method is ∼ 1 -3 Gyr (Figure 6). Differences of this magnitude are consistent with the variations in stellar age that we found in Oyarz'un et al. (2022, also plotted in Figure 6) after similarly dissecting the red SHMR of nearby galaxies. In an approach that is analogue to that of Section 3.1, we also separated central galaxies according to the 33rd and 66th percentiles in M ∗ -to-M H ratio, recovering variations in the timescales of galaxy formation at fixed M H that are not drastically different from the variations in the timescales of halo formation. In this subsection, we will show that variations in the stellar ages and stellar masses of galaxies with the large-scale environment (at fixed M H ) can be interpreted as a natural outcome of linear structureformation theory and the regulation of star-formation through energy-conserving feedback.  \nWe will first discuss the difference in stellar age with galaxy environment in large scales (Figure 6). In this context, R = 10 h -1 Mpc is a natural scale to consider because that is the non-linear mass scale in the local Universe. The redshift at which a particular region of the Universe reaches a virialized or collapsed mass fraction (i.e., M H ) depends on the large scale overdensity, with overdense regions having accelerated structure formation. According to extended Press-Schechter theory, the variance in redshift at which a region of the Universe within a radius R will reach a particular collapsed fraction is  \nσ z = σ ( R ) 1 . 69 (1 + z ) , (9)  \nwhere σ ( R ) is the variance in the density field (Barkana & Loeb 2004; Wyithe & Loeb 2004). For R = 10 cMpc  \nat z = 0, the variance in redshift is σ z ≈ 0 . 6. These variations can be linked to differences in the formation of structure for an overdense Universe, i.e.,  \ndt = ∣ ∣ ∣ ∣ dt dz ∣ ∣ ∣ ∣ dz = σ z (1 + z ) H ( z ) ≈ 1 Gyr at z = 0 . (10)  \nVariations of this magnitude are consistent with the observed differences in age between the stellar populations of galaxies in large scale overdensities and underdensities at fixed M H (Figure 6).  \nWith this result in hand, we can predict the differences in M ∗ between old and young populations at fixed M H . Stellar populations with ages of 10 Gyr formed before quasars were dominant in regulating stellar growth (e.g. Croton et al. 2006a,b). Under the ansatz that stellar growth is set by the balance between the binding energy of the star-forming gas and the mechanical energy from supernovae (the base assumption in almost all semi-analytic models and numerical simulations), the ratio between stellar and halo mass corresponds to  \nM ∗ M H ∝ M 2 / 3 H (1 + z ) , (11)  \nas in Dekel & Woo (2003) and Wyithe & Loeb (2003). This implies that the range in M ∗ at fixed M H owing to variations in the formation time is  \n∆M ∗ = (1 + z + σ z ) (1 + z -σ z ) . (12)  \nFor the derived value of σ z ≈ 0 . 6 for R = 10 cMpc, ∆M ∗ ≈ 3 for stars forming at z = 0 and ∆M ∗ ≈ 2 for stars forming at z = 1. These values match the range of stellar mass that separates the high-M ∗ and low-M ∗ subsamples (∆logM ∗ ≈ 0 . 4; Figure 1).", '5.3. Galaxy assembly bias in simulations': "The results shown in Figure 5 and the analytical arguments presented in Section 5.2 can be interpreted as an underlying correlation between M ∗ and halo formation time. In this scenario, the properties of galaxies are dependent not only on the masses of dark matter halos, but also on secondary halo properties like formation time or assembly history. A fixed M H , older halos form deeper gravitational potentials that facilitate the assembly of higher M ∗ central galaxies (e.g. Matthee et al. 2017).  \nOur model associates variations in δ 10 across M H ∼ 10 11 -10 14 h -1 M ⊙ with a strong correlation between halo formation time and stellar mass ( r s ∼ 0 . 85) across this M H range. In broad agreement with this result, the correlation between halo age and M ∗ is present across the M H ∼ 10 11 -10 13 . 5 h -1 M ⊙ range in the IllustrisTNG simulation (Zehavi et al. 2018; Xu & Zheng  \n2020). Interestingly, the galaxy assembly bias in IllustrisTNG is also present when using the number density of particles within R ∼ 5 h -1 Mpc instead of t halo (Zehavi et al. 2018), which is consistent with our analysis in SDSS (albeit with the considerations regarding the choice of R outlined in Section 3.1).  \nWe note that the magnitude of the δ 10 differences peak around M H ∼ 10 12 M ⊙ . However, because we adopt a constant value for r s across all M H , we may not interpret any dependence of the δ 10 variations on M H as a signature of a M H -dependent galaxy assembly bias signal. Some theoretical work - e.g. IllustrisTNG (Vogelsberger et al. 2014; Nelson et al. 2015) -predict the galaxy assembly bias signal to decrease in magnitude toward M H > 10 13 . 5 h -1 M ⊙ (Zehavi et al. 2018; Xu & Zheng 2020). On the other hand, galaxy assembly bias remains strong in the UniverseMachine (Behroozi et al. 2019) for M H > 10 14 h -1 M ⊙ (Bradshaw et al. 2020). The EAGLE simulation (Crain et al. 2015; Schaye et al. 2015) is at the other end of the spectrum, predicting that halo formation time has little impact above M H ≳ 10 11 . 5 h -1 M ⊙ (Matthee et al. 2017).  \nSome of the discrepancy between these simulations may originate in how halo and galaxy mergers are handled, especially given that mergers become increasingly important as M H increases (Huang et al. 2013a,b, 2018, 2020; Oyarz'un et al. 2019; Angeloudi et al. 2024). In the case of our work, this is apparent in the little sensitivity of δ 10 variations to the correlation between t halo and M ∗ toward M H > 10 13 . 5 h -1 M ⊙ (Figure 4). This justifies the choice of fitting for a single, M H -independent value for r s , and it ultimately stems from how t halo and δ 10 are correlated at these M H in MDPL2.  \nFor a more quantitative comparison with simulations, we can turn to r s . We are particularly interested in comparing our recovered value of r s ∼ 0 . 85 - which is indicative of a strong correlation between t halo and M ∗ at fixed M H -with expectations from semi-empirical models. For this purpose, we analyzed data for the SHMR from the UniverseMachine. We calculated t halo following the same definition outlined in Section 3.1 and computed r s between t halo and M ∗ at fixed M H . We found r s to peak at M H ∼ 10 11 -10 12 h -1 M ⊙ with values of r s ∼ 0 . 5 -0 . 6. We conclude that semi-empirical models like the UniverseMachine predict a correlation between t halo and M ∗ that is only slightly weaker than what we inferred from SDSS.  \nSo far, we have discussed how variations in halo formation time can affect the M ∗ (this work) and stellar ages (Oyarz'un et al. 2022) of central galaxies (Figure 6). We should highlight that the stellar population variations that we identified in Oyarz'un et al. (2022) also  \nextend to chemical content of galaxies. Besides being older, centrals with high M ∗ -to-M H ratios also feature low [Fe/H] and high [Mg/Fe]. The natural interpretation for these trends is that centrals with high-M ∗ for their M H not only form their stars early, but that they also assemble them in rapid timescales, resulting in high [ α /Fe] (e.g. Faber & Jackson 1976; Faber et al. 1985; Thomas et al. 2005). Variations in the star-formation timescale may also be connected with variations in secondary halo properties, in particular with halo formation time and halo concentration. At fixed M H , older and higher-concentration halos are thought to bind the stellar and gaseous contents of centrals more strongly (e.g. Wechsler et al. 2002; Hearin et al. 2016). These results - along with the work by Winkel et al. (2021); Donnan et al. (2022); Zhou et al. (2024) - provide tentative observational evidence for a correlation between secondary halo properties and the quenching epochs and the metal content in galaxies (as seen in simulations; e.g. Booth & Schaye 2010; Matthee et al. 2017; MonteroDorta et al. 2020, 2021; Hasan et al. 2023b).  \nOur studies of how δ 10 and the stellar populations of galaxies vary within the SHMR were in great part enabled by publicly available datasets of high statistical power (i.e., SDSS and MaNGA). SDSS is the largest publicly available spectroscopic survey of galaxies in the nearby Universe, providing the necessary statistics for dissecting the SHMR, the required spectra for accurately measuring stellar masses, and the deep ancillary imaging needed for precise estimation of the halo mass. At the same time, MaNGA is the largest spatially resolved spectroscopic survey of nearby galaxies. In terms of spectroscopic constraining power, MaNGA surpasses even SDSS (see Oyarz'un et al. 2022), allowing for the most precise measurements of stellar age and chemical abundances possible in the SHMR of nearby galaxies. This will soon change, with surveys like the Bright Galaxy Sample of the DESI Experiment (DESI Collaboration et al. 2016) spearheading the next generation of datasets that will enable detailed characterization of the galaxy-halo connection at low redshifts.", '6. SUMMARY': "We used the group catalog from Tinker (2021, 2022) to study the stellar-to-halo mass relation for red central galaxies from SDSS across the halo mass range M H = 10 11 -10 15 h -1 M ⊙ . We showed that the stellar masses of red central galaxies vary not only with the mass of the dark-matter halo, but also with the number density of galaxies in large scales ( R = 10 h -1 Mpc). Our preferred interpretation for this result is that the stellar masses of central galaxies depend on both halo  \nmass and halo formation time. This interpretation is motivated by our current cosmological picture, in which high density regions of the large scale structure assemble earlier in cosmic history.  \nTo quantify the relation between the large-scale number density of galaxies and their halo formation time, we exploited the MDPL2 dark-matter-only cosmological simulation (Prada et al. 2012). With the halo formation times in hand, we constructed a model that reproduces differences in the observed number density of galaxies via two parameters: (1) random errors in the halo masses estimated by the group catalog and (2) the correlation between halo formation time and stellar mass at fixed halo mass (i.e., the galaxy assembly bias signal). We quantified the model posterior with SDSS and found that errors in the halo masses from the group catalog are ∼ 0 . 15 dex (in agreement with the expected ≲ 0 . 2 dex errors from the group finder algorithm; Tinker 2022) and that galaxy assembly bias is necessary for explaining our results at > 5 σ significance.  \nOur best-fit model indicates that variations in the stellar masses of red central galaxies at fixed halo mass correlate with differences in halo formation time of ∼ 1 -3 Gyr. Differences of this magnitude are remarkably consistent with differences in the star-formation histories of galaxies at fixed halo mass, which we measured previously in Oyarz'un et al. (2022). Taken together, these results are indicative that old halos (at fixed halo mass) are associated with more massive, older, and more rapidly assembled central galaxies.  \nWe also considered how systematic errors in the group finder algorithm could fabricate the signal that we observe in the galaxy number densities. While systematic and random errors are likely present, they alone cannot reproduce the signal. This is because we still detect a similarly significant difference in the galaxy number densities after substituting M H for total group luminosity and M ∗ for galaxy luminosity.", '7. DATA AVAILABILITY': "Underlying data in this article will be shared upon reasonable request to the corresponding author.  \nWe thank the reviewer of this paper for their insightful comments and suggestions. G.O. acknowledges support by the National Science Foundation through grants AST-2107989 and AST-2107990. This work made use of Astropy 3 , a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration et al. 2013, 2018, 2022).  \nFunding for the Sloan Digital Sky Survey (SDSS) has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Aeronautics and Space Administration, the National Science Foundation, the U.S. Department of Energy, the Japanese Monbukagakusho, and the Max Planck Society. The SDSS Web site is http://www.sdss.org/. The SDSS is managed by the Astrophysical Research Consortium (ARC) for the Participating Institutions. The Participating Institutions are The University of Chicago, Fermilab, the Institute for Advanced Study, the Japan Participation Group, The Johns Hopkins University, Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, Univer-  \nsity of Pittsburgh, Princeton University, the United States Naval Observatory, and the University of Washington. The MultiDark Database used in this paper and the web application providing online access to it were constructed as part of the activities of the German Astrophysical Virtual Observatory as result of a collaboration between the Leibniz-Institute for Astrophysics Potsdam (AIP) and the Spanish MultiDark Consolider Project CSD2009-00064. The Bolshoi and MultiDark simulations were run on the NASA's Pleiades supercomputer at the NASA Ames Research Center. The MultiDark-Planck (MDPL) and the BigMD simulation suite have been performed in the Supermuc supercomputer at LRZ using time granted by PRACE."}