Stellar and weak lensing profiles of massive galaxies in the Hyper-Suprime Cam survey and in hydrodynamic simulations Free

ABSTRACT

We perform a consistent comparison of the mass and mass profiles of massive (M_⋆ > 10^11.4 M_⊙) central galaxies at z ∼ 0.4 from deep Hyper Suprime-Cam (HSC) observations and from the Illustris, TNG100, and Ponos simulations. Weak lensing measurements from HSC enable measurements at fixed halo mass and provide constraints on the strength and impact of feedback at different halo mass scales. We compare the stellar mass function (SMF) and the Stellar-to-Halo Mass Relation (SHMR) at various radii and show that the radius at which the comparison is performed is important. In general, Illustris and TNG100 display steeper values of α where |$M_{\star } \propto M_{\rm vir}^{\alpha }$|⁠. These differences are more pronounced for Illustris than for TNG100 and in the inner rather than outer regions of galaxies. Differences in the inner regions may suggest that TNG100 is too efficient at quenching in situ star formation at M_vir ≃ 10¹³ M_⊙ but not efficient enough at M_vir ≃ 10¹⁴ M_⊙. The outer stellar masses are in excellent agreement with our observations at M_vir ≃ 10¹³ M_⊙, but both Illustris and TNG100 display excess outer mass as M_vir ≃ 10¹⁴ M_⊙ (by ∼0.25 and ∼0.12 dex, respectively). We argue that reducing stellar growth at early times in |$M_\star \sim 10^{9-10} \, \mathrm{M}_{\odot }$| galaxies would help to prevent excess ex-situ growth at this mass scale. The Ponos simulations do not implement AGN feedback and display an excess mass of ∼0.5 dex at r < 30 kpc compared to HSC which is indicative of overcooling and excess star formation in the central regions. The comparison of the inner profiles of Ponos and HSC suggests that the physical scale over which the central AGN limits star formation is r ≲ 20 kpc. Joint comparisons between weak lensing and galaxy stellar profiles are a direct test of whether simulations build and deposit galaxy mass in the correct dark matter haloes and thereby provide powerful constraints on the physics of feedback and galaxy growth. Our galaxy and weak lensing profiles are publicly available to facilitate comparisons with other simulations.

1 INTRODUCTION

In the Λ cold dark matter (ΛCDM) cosmological paradigm, dark matter haloes (and their associated galaxies) grow hierarchically, building up mass over time through mergers and accretion. However, developing a complete model for galaxy evolution which fully explains the properties of observed galaxies has historically been challenging. For a number of years, massive galaxies in simulations have typically been brighter, bluer, and contained more stellar mass than what is actually observed (e.g. Borgani & Kravtsov 2011). This has often been referred to as the “overcooling problem” because it results from excess cooling of gas leading to extended star-formation and overly massive galaxies (Benson et al. 2003). An energetic feedback process is necessary to quench star formation in massive galaxies and to reproduce the high mass end of the galaxy mass function (Borgani & Kravtsov 2011). Enhanced stellar feedback has been proposed as a solution (Springel & Hernquist 2003a,b), but recent work has shown that it is simply not enough to account for the necessary energy injection (e.g. Teyssier et al. 2011; Martizzi, Teyssier & Moore 2012a). The most popular scenario for this energy source is feedback by supermassive black holes (BHs) at the centres of active galactic nuclei (AGN).

With the introduction of feedback (from both stars and AGN), many modern hydrodynamic simulations of galaxy formation reproduce, to first order, the observed galaxy stellar mass function (SMF; e.g. Genel et al. 2014; Schaye et al. 2015; Beckmann et al. 2017; Pillepich et al. 2018b). The subgrid models of hydrodynamical simulations where these feedback processes are numerically approximated are usually calibrated with some set of observables, including: the galaxy SMF, the present-day stellar-to-halo mass relation (SHMR), the star formation rate density (SFRD), the BH mass to galaxy or halo mass relation, the halo gas fraction, galaxy stellar sizes, and the mass–metallicity relation. However, as the demand for even higher fidelity simulations increases, higher order effects must be taken into account. In particular, when calibrating to the galaxy SMF, consistent definitions for the masses of galaxies are not always adopted between the observations and the simulations (e.g. see discussion in Pillepich et al. 2018b).

In simulations, stellar mass is the direct outcome of the underlying physical recipe and can be precisely known. None the less, there is still a variety of stellar mass definitions that are adopted, the most common of which are: the sum of all the stellar particles gravitationally bound to a galaxy as defined by some halo finder (e.g. Dubois et al. 2016; Remus, Dolag & Hoffmann 2017), the sum of all stellar particles within 3D spherical apertures at fixed radius (e.g. Schaye et al. 2015), or some combination of those two definitions (e.g. Pillepich et al. 2018b). Other works have also created mock data from the simulations and performed observationally motivated measurements in order to recover the mass in a way that is more consistent with observations (e.g. Price et al. 2017; Laigle et al. 2019).

In data, stellar mass estimates depend on the mass-to-light ratio estimates (M/L), modelling of galaxy light profiles to extract galaxy luminosity (L), sky subtraction, knowledge of redshifts, and can also depend galaxy morphology (e.g. Bernardi et al. 2013; Muzzin et al. 2013; D’Souza, Vegetti & Kauffmann 2015). In this paper, however, we focus specifically on very massive galaxies (⁠|$M_\star \gt 10^{11.4}\, \mathrm{M}_{\odot }$|⁠), which present different challenges compared to the general galaxy population. Supermassive galaxies have redshifts that are generally well known. They are red and form a fairly homogeneous population with little dust and relatively shallow (M/L) gradients. Instead, the dominant challenge for super massive galaxies is the estimate of their luminosity due to their extended diffuse components that extend to 100 kpc and beyond. These outer low surface brightness regions of a galaxy can lead to important corrections in stellar mass (Bernardi et al. 2013; D’Souza et al. 2015; Huang et al. 2018b). However, recently, deep and wide multiband imaging from the Subaru Strategic Program (SSP; Aihara et al. 2018a,b) using Hyper Suprime-Cam (HSC; Miyazaki et al. 2012; Miyazaki et al. 2018) has advanced our understanding of the light profiles of individual very massive galaxies out to 100 kpc and beyond (Huang et al. 2018a,b). This allows for a robust determination of the stellar mass within 100 kpc (⁠|$M_{\star }^{100}$|⁠) as well as the full shape of the galaxy light profile from 10 to 100 kpc.

The goal of this paper is to outline a framework for performing a more direct comparison between observations and hydrodynamic simulations that capitalizes on new generation surveys that are wide, deep, and that have weak lensing capabilities. First, with surveys such as HSC, weak lensing enables comparisons at fixed halo mass (e.g. Huang et al. 2020). Secondly, in addition to making use of weak lensing, we also advocate for the comparison of galaxy mass profiles. This avoids having to define “galaxy mass” as a single number and can account for the fact that observations have finite depth. Finally, the radial shape of the galaxy mass profile also contains information about the assembly history and feedback processes that have shaped massive galaxies.

This paper is organized as follows: in Section 2, we describe our data for both observations and simulations; we then discuss how we measure mass density profiles in observations (Section 3) and simulations (Section 4); we present results from comparing measurements in observed and simulated galaxies in Section 5 and discuss potential reasons for disagreement in Section 6; finally, we summarize and conclude in Section 7. We provide median surface mass density profiles and weak-lensing ΔΣ profiles for the HSC galaxies here: https://github.com/f-ardila/HSC_versus_hydro-paper.

The following cosmological parameters are assumed throughout: H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, and Ω_Λ = 0.7. We use the virial mass for dark matter halo mass (⁠|$M_{\mathrm{vir}} = 4\pi \Delta (z) \rho _{\mathrm{c}}(z) R_{\mathrm{vir}}^{3} / 3$|⁠) as defined in Bryan & Norman (1998). All length units use physical units (not comoving).

2 DATA AND SIMULATIONS

2.1 HSC data

In this paper, we take advantage of new high quality imaging from the WIDE layer of the HSC SSP, which is a simultaneously wide (>1000 deg²) and deep (r ∼ 26 mag) imaging survey (with the DEEP and ULTRADEEP layers reaching r ∼ 27 and r ∼ 28 and a coverage of 27 and 3.5 deg², respectively; Aihara et al. 2018a,b), meaning that we can observe a large sample of massive galaxies and detect their faint extended stellar envelopes. The imaging is deep enough that we are able to measure surface brightness profiles of individual massive galaxies to ∼28.5 mag arcsec⁻² in i-band (Huang et al. 2018b).

This work uses the internal s16A data release (equivalent to “DR1”) which covers |$\sim 140\ \rm {deg}^2$| in all five broad-bands (grizy) to full WIDE depth. The combination of a wide (1|${_{.}^{\circ}}$|5) field of view and exceptional imaging depth and quality (median i-band seeing FWHM ∼0|${_{.}^{\prime\prime}}$|6) makes this survey ideally suited to study the surface brightness profiles of galaxies out to large radii. HSC i-band images were used to make our surface brightness (and stellar mass density) profile measurements because of the superior seeing in this band as a result of strict requirements imposed by weak-lensing science.

We use the HSC massive galaxy sample from Huang et al. (2018b) and refer the reader to that paper for the full details regarding the construction of this sample; we give only a brief summary here. Our sample consists of ∼15000 galaxies with reliable spectroscopic redshifts in the range 0.3 ≤ z ≤ 0.5. Our selection begins with an initial magnitude cut of i_{HSC, cModel} ≤ 21.5 mag to select massive galaxies (log (M_⋆/M_⊙) ≥ 11.5) at z < 0.5 based on Leauthaud et al. (2016). We limit our sample to regions that have full depth coverage in i-band in the WIDE layer, objects without deblending errors, with well-defined centroids, and usable cModel magnitudes in all five HSC bands. We also exclude objects affected by pixel saturation, cosmic rays, or other optical artefacts.

We then include only objects with a reliable spectroscopic redshift and restrict redshifts to the range 0.3 ≤ z ≤ 0.5. We focus on this range for several reasons: (1) to resolve and reliably measure the stellar mass within the inner 10 kpc of galaxies; (2) to limit background noise and cosmological dimming when measuring stellar mass out to 100 kpc and beyond; (3) to ignore the redshift evolution of the stellar populations in these galaxies; (4) to ensure stellar mass completeness of our sample; (5) to reduce the known issue of oversubtraction of the background that occurs at lower redshifts. Huang et al. (2018b) have shown that the light profiles of individual galaxies in this sample can be mapped to 100 kpc and beyond (without using stacking techniques).

2.2 Illustris and TNG100 simulation data

In this work, we focus on two uniform-volume (i.e. not zoom-in) cosmological hydrodynamic simulations: the Illustris-1 simulation of the Illustris Project (referred to hereafter as “Illustris”; Vogelsberger et al. 2014a,b; Genel et al. 2014; Sijacki et al. 2015) and the TNG100 simulation of the IllustrisTNG Project (referred to hereafter as “TNG100”; Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018; Pillepich et al. 2018b; Springel et al. 2018). Table 1 gives the main parameters of both simulations. They are also briefly described below.

Illustris consists of three full physics simulation boxes of the same size (75 h⁻¹ Mpc), but varying resolutions, run to z = 0 with the moving-mesh arepo code. There are also additional runs with dark matter only and adiabatic scenarios. For this project, we use the highest resolution, full physics run (Illustris-1). The Illustris galaxy formation model includes gravitational interactions among the different resolution elements (dark matter, stars, gas, and BHs); hydrodynamical equations for the gas component; a treatment of radiative cooling and heating processes; a mechanism for converting gas into stars; stellar evolution and the resulting chemical enrichment of the interstellar, circumgalactic and intergalactic media; stellar feedback induced outflows of gas; and the formation, growth and energetic feedback of supermassive BHs in distinct low- and high-accretion rate states (Vogelsberger et al. 2013; Torrey et al. 2014). The fiducial model was chosen to broadly reproduce the z = 0 galaxy SMF, the evolving cosmic SFR density, and the z = 0 relation between galaxy mass and gas-phase metallicity (Torrey et al. 2014).

IllustrisTNG (“The Next Generation”) is the successor to Illustris and has a larger range of box sizes and resolutions and several enhancements. For this project, we use the highest resolution of the medium size TNG100 run, which has the same 75 h⁻¹ Mpc box size and similar resolution as Illustris-1. The IllustrisTNG runs include improved numerical methods, new physics (e.g. magnetohydrodynamics), and modified sub-grid physical models to address the key shortcomings of the Illustris model. Most importantly, the new dual mode (thermal and kinetic) AGN feedback model helps to better regulate the stellar content of massive galaxies while preserving realistic halo gas fractions; and the galactic winds feedback model has been updated to better reproduce the abundance or mass of intermediate- and low-mass galaxies (Weinberger et al. 2017; Pillepich et al. 2018a). The fiducial IllustrisTNG model was chosen by comparing the outcome of many TNG model variations to that of the original Illustris model (and not directly to observations) and by ultimately choosing an implementation and parameter set that simultaneously had the potential to alleviate, at least qualitatively, the largest number of targeted tensions between Illustris and observations. For this purpose, the comparison focused on the global SFR density as a function of time; the current galaxy SMF; the stellar mass–halo mass relation at z = 0 in addition to the BH mass versus stellar or halo mass relation; the gas fraction within the virial radii of haloes; and the stellar half-mass radii of galaxies.

Galaxy stellar masses between Illustris and IllustrisTNG were compared by adopting the same operational definition, namely the sum of gravitationally bound stellar particles contained within twice the stellar half-mass radius of each SUBFIND (sub)halo. This was also the operational definition used for the original “calibration” of the Illustris model, against SMFs measured by Baldry, Glazebrook & Driver (2008) and Pérez-González et al. (2008). It should be kept in mind, however, that in the development of the Illustris and IllustrisTNG frameworks, the best model parameter values have not been obtained via an actual fit to a selection of galaxy observables, but simply by requiring consistent trends and overall normalizations. The adopted approach was intentionally kept simple because an actual fine tuning of galaxy population models against observational data would have required complex steps to be applied to the output of hundreds of model variations (e.g. transforming raw simulated data into realistic mock observations and applying selection functions tailored to each targeted observational data set). Furthermore, the very high-mass end of the galaxy population could not be taken into account during the development phase, because of the computational costs of simulating hundreds of large volumes at fixed resolution. The high mass end is therefore a regime where the models are de facto predictive.

For this paper, we choose all central galaxies in Illustris and TNG100 with total stellar masses |$M_* \gt 10^{11.2} \, \mathrm{M}_{\odot }$| at z = 0.4 (snapshot 108 in Illustris and 72 in TNG). This results in a sample of 339 galaxies in Illustris and 235 galaxies in TNG100. We choose this redshift in order to be consistent with the median redshift of our observations. Here, a central galaxy is defined by the subfind halo finder as the most bound subhalo within a larger friends-of-friends (FOF) group. All other bound haloes within the FOF group (and their galaxies) will be denoted as subhalos and satellites, respectively. This stellar mass cut results in a |$\sim 85{{\ \rm per\ cent}}$| (⁠|$\sim 91{{\ \rm per\ cent}}$|⁠) completeness of central galaxies with M_vir > 10¹³ M_⊙ for TNG100 (Illustris).

2.3 Ponos simlation data

We also compare the stellar mass profiles of our HSC galaxies to two highly resolved galaxies from the Ponos zoom-in numerical simulations of dark matter substructures in massive ellipticals (Fiacconi et al. 2016, 2017) including hydrodynamics using the TreeSPH code gasoline2 (Wadsley, Stadel & Quinn 2004), which employs a modern implementation of the SPH equations based on the geometric density formulation of the pressure force, sub-grid turbulent diffusion of thermal energy and metals, and the Wedland C4 kernel. At z = 0.4, the two Ponos galaxies (PonosV and PonosSB) have identical dark matter halo masses (⁠|$M_{\mathrm{\rm vir}}=1.04 \times 10^{13} \, \mathrm{M}_{\odot }$|⁠) but different stellar feedback mechanisms resulting in slightly different stellar masses between the two (⁠|$M_{\star }^{\mathrm{100}}\sim 10^{12} \, \mathrm{M}_{\odot }$|⁠).

PonosV¹ (‘V’ for violent merger history compared to more quiescent haloes) adopts the blastwave feedback sub-grid model (Stinson et al. 2006), in which thermal feedback from explosions of SN Type II is achieved by locally shutting off cooling for a time-scale comparable to the duration of the Sedov–Taylor and snowplough phases alltogether. Blastwave feedback can lead to realistic stellar masses and structural properties of disc galaxies provided that the internal parameters are set according to resolution (Sokołowska et al. 2017). PonosSB employs the more recent Superbubble (SB) feedback implementation by Keller et al. (2014), in which thermal energy and ejecta are deposited into multiphase particles, with the warm/hot phase returning to the cold phase over a few Myr owing to thermal conduction. Superbubble feedback does not require to shut-off cooling, and has been shown to be more efficient at regulating stellar masses in low mass and disc galaxies (up to M_vir ∼ 10¹²) in a resolution-independent way, matching easily the stellar mass-to-halo mass relation in this mass range (Keller, Wadsley & Couchman 2015).

Internal parameters of the two sub-grid feedback models, and of the underlying star formation recipe, are set as in Mayer et al. (2016). Neither of the two galaxies include AGN feedback, and both include radiative cooling from metal lines. Comparisons with the mass profiles from these Ponos galaxies and Illustris/TNG100 will shed light both on the impact of the resolution of the simulations, as well as the impact of AGN feedback (but see also Appendices and discussions on the effects of both resolution and AGN feedback within the TNG model itself in Pillepich et al. 2018a,b)

2.4 Mass maps

One of the main goals of this paper is to perform a consistent comparison between observations and theory. To treat simulated galaxies in a similar way as observations, we first translate the 3D distribution of their stellar population particles into 2D stellar mass maps. We choose to use projected mass maps rather than mock multiband observations for each galaxy because we wish to limit the number of assumptions we need to make. Because, we are interested in studying mass profiles, mock observations would require additional assumptions to convert these to stellar mass. Even if we were to study luminosity profiles rather than mass, we would still need assumptions about stellar population models to perform the k-corrections to build average luminosity profiles of our observed sample. Therefore, we choose to put these assumptions only on the observational side and keep the simulation data as is.

For each galaxy, we create a map by projecting the positions of stellar population particles on to a 300² grid of pixels. We have experimented with smoothing the maps using a Gaussian kernel, but find that such smoothing can lead to slight errors in the integrated stellar mass in the map. Thus, we simply assign the total mass of a stellar particle to the map pixel in which its centre lies.

Our maps have a pixel size of 1.0 kpc pixel⁻¹ and a physical side length of 300 kpc. These allow us to more accurately trace galaxy profiles on the scales that we are interested, below r < 100 kpc. We tested different resolution scales and found that with maps of 1.0 kpc pixel⁻¹, the measured mass profile was stable beyond 6 kpc. Within 6 kpc, we are limited by seeing in observations anyway, so we do not compare profiles below this scale.

We use all stellar population particles within the FOF group (i.e. we impose the particular projection depth determined by the FOF) to perform various tests. We distinguish between particles bound to the galaxy, particles bound to satellite galaxies in the same group, and particles that are not bound to any galaxy in the group. In practice, the latter component is negligible, meaning that all the stars in our maps are bound to either the central galaxy or its satellites. We also find the satellite component does not affect the measured mass profile within 100 kpc (§4.2) so we only include mass maps with the central galaxies in our analyses.

3 MEASURING MASSES AND MASS PROFILES IN DATA

3.1 One-dimensional surface mass density profiles

We measure the 1D stellar mass density profiles of HSC galaxies from their HSC i-band images using the method of Huang et al. (2018b). We use the galaxy surface brightness profile function (galsbp) in the publicly available kungpao package.² First, we apply an empirical background correction and mask out neighbouring objects based on their brightness and proximity. Then, we draw concentric elliptical isophotes with a fixed ellipticity on the target, and at a given radius (along the semimajor axis), we measure the median intensity along each isophote after iterative 3σ-clipping. The ellipticity parameters are computed from the intensity-weighted moments of flux distribution of the object in the map using SE xtractor (Bertin & Arnouts 1996, through the sep python library). Hereafter, in all figures, R corresponds to a distance along the semimajor axis of the elliptical isophote.

We then convert these profiles into surface stellar mass profiles (μ_*) by assuming a radially constant mass-to-light ratio measured from SED fitting. We note that massive elliptical galaxies are known to have shallow negative colour gradients (e.g. Carollo, Danziger & Buson 1993; Davies, Sadler & Peletier 1993; La Barbera et al. 2012; D’Souza et al. 2015), which may in principle underestimate the stellar mass in the centre, while overestimating the stellar mass in the outskirts. However, because the gradients are shallow, and they are smooth out to a few times the effective radius (e.g. La Barbera et al. 2010; D’Souza et al. 2014), an average M/L is unlikely to bias our stellar mass measurements (for details see Appendix C of Huang et al. 2018b).

Integrating μ_* profiles provides us with the stellar mass (M_*) within elliptical apertures. We use the stellar mass within a 100 kpc aperture (hereafter noted |$M_{\star }^{\mathrm{100}}$|⁠) as a value close to the total stellar mass of a galaxy. We also use the stellar mass within a 10 kpc aperture (⁠|$M_{\star }^{\mathrm{10}}$|⁠) as the inner stellar mass of the galaxy, and the difference (⁠|$M_{\star }^{\mathrm{100}} - M_{\star }^{\mathrm{10}}$|⁠) as the outer stellar mass. Estimates for exactly how much mass is missed by cutting off at 100 kpc will be given in a follow-up paper (Ardila et al., in preparation). As shown in Huang et al. (2018b) and Li et al. in preparation, at large scales we can reliably measure μ_* profiles for individual galaxies out to more than 100 kpc without being limited by the background subtraction, while on small scales our profiles are resolved down to ∼6 kpc. As a reference, 1|${_{.}^{\prime\prime}}$|0 corresponds to 3 and 6 kpc at z = 0.2 and 0.5, respectively, and the mean i-band seeing has FWHM |$=0{_{.}^{\prime\prime}}58$|⁠.

Our HSC WIDE data can reach >29 mag arcsec^–2 in surface brightness radial profile measurements in the i-band. It is known that the hscPipe deblending process is not optimized for these extended objects and that it tends to oversubtract background light. To avoid these problems, we use SE xtractor (Bertin & Arnouts 1996) background subtraction and object detection (using the sep pythoN library) to generate the appropriate masks. Appendix B in Huang et al. (2018b) gives more details about the masking performed in our HSC sample. As is shown in that paper, different masking methods do not affect the measured surface brightness profiles of HSC galaxies within 100 kpc. Even past 100 kpc, the difference resulting from different masking methods is small. We further show that the effect of satellite contamination is negligible with our simulated galaxies in Section 4.2 (see also Fig. 1).

We are unable to measure reliable 1D profiles for |$\sim 11 {{\ \rm per\ cent}}$| of the massive galaxies in our sample due to a complex inner structure (e.g. on-going a major merger), or substantial contamination from a bright star or foreground galaxy, but excluding them does not bias our sample, nor does it significantly affect the measured SMF (Huang et al. 2018b).

3.2 Halo masses from weak lensing

For our massive galaxies, we measure their ΔΣ profiles using the weak lensing shape catalogue defined by Mandelbaum et al. (2018). We follow the strategy outlined by Singh et al. (2017) and use the dsigma code ,³ which is optimized to work with HSC-SSP data. See Speagle et al. (2019) and Huang et al. (2020) for details of the g–g lensing measurements.

With the help of these g–g lensing profiles and the deep stellar mass density profiles, we are able to build an empirical model that connects the stellar mass distribution of massive galaxies to their underlying halo masses (Huang et al. 2020). This model (referred to as ASAP) is based on the Small MultiDark Planck (SMDPL; Klypin et al. 2016) N-body simulation of halo merging history and a modified version of the UniverseMachine semi-empirical model (Behroozi et al. 2019). ASAP is constrained by weak lensing data and by HSC observations of massive galaxies with |$M_{\star }^{\mathrm{100}}\gt 10^{11.6} \, \mathrm{M}_{\odot }$|⁠. In this mass range, the satellite fraction is very low (<10 per cent). ASAP provides the link between the halo masses (M_vir) of massive central galaxies, the stellar mass within 10 kpc (⁠|$M_{\star }^{\mathrm{10}}$|⁠), and the “total” stellar mass (⁠|$M_{\star }^{\mathrm{100}}$|⁠).⁴ The intrinsic scatter in this relation is 0.15–0.20 dex.

3.3 Comparisons at fixed halo mass

Given the ASAP model that best describes the HSC data (Huang et al. 2020), we can also calculate the full halo mass distribution, |$P(M_{\rm vir}|\overline{M}_{\rm vir})$|⁠, corresponding to any given cut in |$\overline{M}_{\rm vir}$|⁠. In practice, we measure |$P(M_{\rm vir}|\overline{M}_{\rm vir})$| from the full distribution of halo masses in the SMDPL simulation which is populated with galaxies drawn from the best-fit ASAP model. Additionally, for a given |$\overline{M}_{\rm vir}$| cut in HSC, we can select a distribution of haloes in Illustris/TNG that approximately matches |$P(M_{\rm vir}|\overline{M}_{\rm vir})$|⁠. Note that this matching relies on modelling the HSC data with the ASAP model. However, the ASAP model is more sophisticated than most models used in past work (e.g. Halo Occupation Models) and was specifically designed to interpret lensing data for supermassive galaxies. In the text, we will clarify when |$\overline{M}_{\rm vir}$| and |$P(M_{\rm vir}|\overline{M}_{\rm vir})$| are used.

4 MEASURING MASSES AND MASS PROFILES IN SIMULATIONS

4.1 Mass profiles in simulations

To measure the mass profiles of simulated galaxies, we adopt a nearly identical procedure as described in Section 3.1 for observed galaxies. The only difference is that we start with 2D projected stellar mass maps of each galaxy and measure the surface mass density profiles directly from those (i.e. we do not need to consider the mass-to-light ratio in the simulations). As with the HSC data, we draw concentric elliptical isophotes with a fixed ellipticity on each galaxy, and measure the mean intensity at a given radius along the semimajor axis of each isophote, after applying iterative 3σ clipping. The ellipticity and position angle of the isophotes are set to the average values from fitting the 2D shape of each galaxy. The techniques we applied to HSC data to subtract contributions from backgrounds and foregrounds are not required by our simulation analysis since these contaminants are not present in the stellar particle data extracted from the simulation output based on halo finders, as used here. This allows us to test the impact of masking satellites in the next section (Section 4.2).

4.2 Testing the effect of masking satellites

In order to reliably measure the surface brightness profiles of central galaxies, it is necessary to mask out nearby satellites, foreground, and background sources. This can be challenging at the depth of HSC, especially for our sample of massive ellipticals with very extended profiles and for which nearby satellites are common. As mentioned in Section 3.1 (details can be found in Huang et al. 2018b), different masking choices did not affect the measured surface brightness profile of HSC galaxies within 100 kpc.

Our simulated galaxies also allow us test the effects of masking satellites in the surface brightness profile measurements of HSC data. Given that we can easily include or exclude particles belonging to satellites, we can measure the effect satellites have on the surface mass density profile. Fig. 1 shows the mean, 1σ, and 2σ fractional difference in the mass density profiles between mass maps with and without satellites, when the same 3σ-clipping procedure to remove satellite contamination is applied. Even at the 2σ level, when including satellites (i.e. without any masking applied at all) the difference in the mass density profile will be less than |$2{{\ \rm per\ cent}}$| within 150 kpc. Since this test only compares perfect masking to non-masking, any realistic satellite masking performed on our data can be expected to perform better than the non-masking situation. This demonstrates that masking is not so important when using our profile measurement technique within 100 kpc.

While it has been shown that the stellar mass in satellites is not negligible compared to the total stellar mass of a cluster (e.g. Pillepich et al. 2018b), our results show that including satellites will not heavily affect the measurement of the mass density profile of a galaxy. This can be understood through the technique used to make the measurement. As described in Section 3.1, along each isophote, we measure the mean intensity after iterative 3σ clipping. This means that any satellites along an isophote that would heavily bias our measurement due to their high mass intensity would be excluded by the 3σclipping. These results show that the details of our masking strategy do not impact our HSC results, and also that our results are not affected by satellites below our detection threshold.

It is worth noting that this test only includes physically associated satellite galaxies. In reality, surface brightness profiles may also be contaminated by nearby foreground and background objects. Using mock galaxies on a more recent HSC data release, Li et al. (in preparation) uses a more realistic scenario to further demonstrate that our masking strategy can secure the surface brightness profile out to 100 kpc.

5 RESULTS

With a consistent methodology for measuring surface mass density profiles for massive central galaxies in both HSC observations and hydrodynamic simulations, we can now begin to compare and quantify similarities and differences between the data and the simulations. There are several ways to make these comparisons. We begin by investigating differences in the stellar mass function (SMF; Section 5.1.1) and the stellar-to-halo mass relation (SHMR; Section 5.1.2). We also compare the stellar mass density profiles of galaxies in HSC and simulations matched by their “total” stellar mass (⁠|$M_{\star }^{\mathrm{100}}$|⁠) in Section 5.2, and matched by their halo mass (⁠|$M_{\mathrm{\rm vir}}$|⁠) in Section 5.3. Finally, to demonstrate that we are indeed comparing at fixed halo mass, we compare weak-lensing profiles in bins of halo mass in Section 5.4.

5.1 Stellar mass functions and stellar-to-halo mass relation

In Fig. 2, we present both the z ∼ 0.4 SMF (left-hand column) and the stellar-to-halo mass relation (right-hand column) for our HSC sample compared to Illustris and TNG100, split by stellar mass definitions. The top two panels use |$M_{\star }^{\mathrm{r}}$| at |$r= $| 10 and 100 kpc for stellar mass. The bottom panel uses the stellar mass between 10 and 100 kpc (representing the outer mass of the galaxy).

5.1.1 SMF

We find that Illustris galaxies are overly massive at all radial scales and mass bins, particularly at higher masses where there can be an offset of ≳0.3 dex. This discrepancy (see also Genel et al. 2014) is likely related to the AGN feedback model, which in the Illustris implementation is now known to not be efficient enough to sufficiently quench star formation (Weinberger et al. 2017; Donnari et al. 2019). In TNG100, where the AGN feedback model was improved, we see that the SMF more closely matches that of HSC. While the amplitude of the TNG100 SMF compared to Illustris is reduced and results in a better match to observations, the shape of the SMF is nearly the same, much shallower than the observed SMF. In TNG100, this means too few lower mass galaxies and too many higher mass galaxies at all aperture masses (the HSC and TNG100 SMFs cross at |$M_{\star }^{\mathrm{10}}\sim 10^{11.5}$|⁠, |$M_{\star }^{\mathrm{100}}\sim 10^{11.7}$|⁠, and |$M_{\star }^{\mathrm{100}}-M_{\star }^{\mathrm{10}}\sim 10^{11.2}$|⁠). The bottom panel, which uses the outer galaxy mass (⁠|$M_{\star }^{\mathrm{100}}-M_{\star }^{\mathrm{10}}$|⁠) SMF, shows that both the Illustris and TNG100 galaxies are slightly offset from HSC across all masses, indicating that there is more outer mass in Illustris/TNG100 galaxies than in HSC galaxies.

Fig. 2 also highlights the fact that the level of agreement between simulations and observations depends on the adopted mass definition (inner mass, total mass, or outer mass). It expands upon the analysis of the galaxy SMFs of Illustris and IllustrisTNG presented in Pillepich et al. 2018b, where 3D based stellar mass definitions (within 10, 30, 100 kpc and twice the stellar half-mass radius) are contrasted at z = 0. Here, we further this work by contrasting with high-quality HSC data and by adopting 2D masses which are more directly comparable to data such as from HSC.

The volume of our HSC sample between 0.3 < z < 0.5 is approximately 420³ Mpc³, compared to ∼110³ Mpc³ for both Illustris and TNG100. For this reason, the errors on the HSC SMF are smaller than for Illustris and TNG100.

5.1.2 SHMR

Fig. 2 compares the SHMR between the observed and simulated galaxies. Rather than use HSC mass measurements directly, we use the values of the ASAP model, applied to the SMDPL simulation. This corresponds to the best-fitting model to the HSC data (HSC SMF, |$M_{\star }^{10}$|⁠, |$M_{\star }^{100}$|⁠, and weak lensing). This means that we use the full distribution of halo masses from SMDPL, |$P(M_{\rm vir}|\overline{M}_{\rm vir})$|⁠, and stellar masses are derived from those halo masses using ASAP. By using this best-fitting model, we can display not only the mean SHMR but also the scatter in the SHMR that best fits the HSC data. We have checked that using a top hat bin in M_vir versus using the full distribution |$P(M_{\rm vir}|\overline{M}_{\rm vir})$| does not impact these results.

Fig. 2 shows that Illustris galaxies in haloes with log (M_vir) ≳ 13.5 contain too much stellar mass for all aperture masses. For Illustris galaxies in haloes of log (M_vir) ∼ 13.5 (log M_vir ∼ 14), there is a constant offset of ∼0.1 dex (∼0.3 dex), corresponding to a |$\sim 30{{\ \rm per\ cent}}$| (⁠|$\sim 100{{\ \rm per\ cent}}$|⁠) difference in stellar mass. Galaxies in less massive Illustris haloes (log M_vir ∼ 13) tend to match SMDPL+ASAP well, and are offset by less than 0.04 dex (⁠|$\lt 10{{\ \rm per\ cent}}$|⁠).

TNG100 galaxies show a slightly different result, but similarly consistent across aperture mass definitions. In general, the lowest mass haloes host galaxies with less stellar mass than SMDPL+ASAP, the highest mass haloes have more stellar mass than SMDPL+ASAP, and galaxies in the middle halo mass bin have stellar masses consistent with SMDPL+ASAP (offset by less than 0.05 dex). For galaxies in haloes of log (M_vir) ∼ 13 there is a constant deficit of ∼0.1 dex in stellar mass, corresponding to a |$\sim 25{{\ \rm per\ cent}}$| difference compared to SMDPL+ASAP. In the most massive haloes (log M_vir ∼ 14), there is an excess of ∼0.15−0.30 dex which corresponds to a |$\sim 40\!-\!90{{\ \rm per\ cent}}$| difference from SMDPL+ASAP, this disagreement being larger in the inner regions of galaxies. In fact, TNG100 returns an outer stellar mass (⁠|$M_{\star }^{100}-M_{\star }^{10}$|⁠) that is consistent with the SMDPL+ASAP  results to better than 0.15 dex across the halo mass range studied here. The fact that Illustris and TNG100 show significant offsets in some halo mass bins in either the inner (⁠|$M_{\star }^{10}$|⁠) or outer (⁠|$M_{\star }^{100}-M_{\star }^{10}$|⁠) stellar mass, or both, suggests that the offsets may be associated to both in situ and ex-situ stellar mass processes. We will discuss these implications in Section 6.

If α is the slope of the high mass end of the SHMR with |$M_{\star }\propto M_{\rm vir}^{\alpha }$| then TNG100 appears to have a steeper value of α compared to SMDPL+ASAP, with differences most notable in the inner regions of galaxies. This could suggest that while the overall strength of AGN feedback better matches observations in TNG100, the halo mass dependence may still need adjusting. Accretion of low mass satellites also play a role at these mass scales (see discussion in Section 6.1).

We have also checked that the scatter in stellar mass at fixed halo mass (Fig. 3) is comparable for TNG100 and SMDPL+ASAP for most of the halo mass range which we investigate in this work (13.0 < log M_vir < 14.25). This suggests that the dominant driver in the difference we see in the SHMR (Fig. 2) is the change in slope rather than the scatter.

Both Ponos galaxies have inner stellar masses that are higher than any galaxy in Illustris, TNG100, or SMDPL+ASAP  at similar halo mass (a difference of >0.25 dex from the median SMDPL+ASAP  value). The outer mass, on the other hand, appears to be comparable to that of the other samples (a difference of ≲ 0.1 dex from the median SMDPL+ASAP  value). The significance of this will be discussed further in Section 6.1.

5.2 Density profile comparison at fixed aperture stellar mass

In this section, we compare the galaxy samples in bins of stellar mass (⁠|$M_{\star }^{\mathrm{100}}$|⁠). Fig. 4 shows the stellar mass density profiles of HSC galaxies and of simulated galaxies (Illustris in red and TNG100 in blue) in three bins of |$M_{\star }^{\mathrm{100}}$|⁠. The lighter thinner dotted lines are the profiles of individual simulated galaxies, while the thicker lines are the median profiles. Because of the large number of HSC galaxies in our sample, we only show the median profile (in black) and the 1σ width in grey. When matching by stellar mass, the overall amplitudes will match by design (since by construction integrated mass within 100 kpc is identical). This comparison therefore serves mainly the purpose of understanding whether or not the profiles have similar slopes. We also note that the completeness of our HSC sample is only |$\sim 65{{\ \rm per\ cent}}$| in the lowest mass bin (⁠|$10^{11.4}\lt M_{\star }^{\mathrm{100}}\lt 10^{11.6}$|⁠), but close to |$100{{\ \rm per\ cent}}$| at higher masses.

For both Illustris and TNG, in every mass bin, the median profiles of the simulated galaxies lie at or within the 1σ width of the distribution of HSC galaxy profiles, suggesting that the slopes of the profiles are in relatively good agreement. Our finding supports recent work on the density profiles of early-type galaxies in TNG100, which show that the power-law density slope of these galaxies is largely in agreement with different observations in different bins of stellar mass (Wang et al. 2020).

The offset of the median Illustris stellar mass from the median HSC mass in bins of stellar mass, is small (∼0.01 dex) in all mass bins and at all mass aperture radii. In TNG100, the mass offsets are of similar magnitude as Illustris in the inner regions of the galaxies (r < 10 kpc), but become larger in the outer regions (r > 10 kpc). Overall, Illustris and TNG100 tend to have slightly shallower profile slopes than HSC, with differences increasing towards the high mass end.

5.3 Density profile comparison at fixed halo mass

5.3.1 Illustris and TNG100

A more informative way to compare galaxies in HSC with those from simulations is to match them by dark matter halo mass. This is a key and unique novelty in this paper that is enabled by the weak lensing measurements from HSC. Fig. 5 shows the median stellar mass density profiles of HSC compared to Illustris (red lines in top panels) and TNG100 (blue lines in bottom panels) in bins of halo mass (same format as Fig. 4). We select three equally spaced bins in the mass range for which we expect to be complete in our HSC sample (⁠|$M_{\mathrm{\rm vir}} \gt 10^{13} \, \mathrm{M}_{\odot }$|⁠). For HSC, we use |$\overline{M}_{\rm vir}$| as described in Section 3.2. For Illustris/TNG100 galaxies, we use M_vir (labelled Group_M_TopHat200) as measured directly by the halo finder from the simulation output.

Fig. 5 shows that in general, both Illustris and TNG100 galaxy profiles agree relatively well with HSC galaxies. For both Illustris and TNG100, at larger halo masses there are larger offsets from HSC. In haloes of |$M_{\star } \sim 10^{14} \, \mathrm{M}_{\odot }$|⁠, the stellar mass profiles of TNG100 (Illustris) are offset from HSC by about 0.3 (∼0.5) dex in the very outskirts (∼100 kpc). In the same high mass bin, a similar significant offset in stellar mass is manifest also in the innermost regions (<6 kpc) for both Illustris and TNG100, while TNG100 exhibits stellar mass profiles in agreement with observed galaxies to better than 0.1 dex across the 10–80 kpc range. In the lower halo mass bins (⁠|$M_{\star } \sim 10^{13} \, \mathrm{M}_{\odot }$|⁠), Illustris galaxies have an excess of up to ∼0.1 dex in both the inner and outer regions. On the other hand, lower halo mass TNG100 galaxies display a notable deficit (∼0.1 dex) in stellar mass in the inner regions (around 6–10 kpc galactocentric distances). Overall, these differences translate into the result that Illustris/TNG100 display steeper SHMRs than HSC as is also seen in Fig. 2, but with a remarkable agreement in the stellar mass distributions between HSC and TNG100 galaxies in |$M_{\star } \sim 10^{13.5} \, \mathrm{M}_{\odot }$| haloes.

5.3.2 Ponos simulation

Our two Ponos galaxies are different hydrodynamic realizations of the same mass halo (⁠|$M_{\mathrm{\rm vir}} = 1.04 \times 10^{13} \, \mathrm{M}_{\odot }$|⁠). In this paper, they serve as an interesting counterpart to Illustris/TNG100 because they do not include AGN feedback and so help to inform our understanding of the impact of AGN feedback on galaxy mass profiles.

We compare Ponos galaxies with HSC at fixed halo mass. We select 75 HSC galaxies from our full sample such that |$1.03 \times 10^{13} \, \mathrm{M}_{\odot }\lt \overline{M}_{\mathrm{\rm vir}}\lt 1.05 \times 10^{13} \, \mathrm{M}_{\odot }$|⁠, where |$\overline{M}_{\mathrm{\rm vir}}$| is calculated using the mean M_vir−M_⋆ relation as given by the ASAP model presented in Huang et al. (2020; see Section 3.2). It should be noted that the median mass of our HSC comparison subsample (⁠|$\log M_{*}^{\mathrm{100}}/\mathrm{M}_{\odot }=11.35$|⁠) is below our completeness limit and so the comparison in this regard is not as accurate as in previous sections.

Fig. 6 shows a comparison of the stellar mass surface density profiles for HSC, Ponos, Illustris, and TNG100. Even though we only have two Ponos galaxies, it is clear from Figs 2 and 6 that Ponos galaxies are offset compared to HSC/Illustris/TNG100. The very inner regions of Ponos galaxies display an extreme excess of mass compared to HSC, likely resulting from the lack of AGN feedback in these simulations. This impacts total stellar mass and results in a ∼0.3−0.5 dex excess in |$M_{*}^{\mathrm{100}}$| compared to HSC. Ponos galaxies suffer from the traditional overcooling problem and this is primarily manifested in the inner region of their profiles. Other analyses of simulated galaxies without AGN feedback show a similar excess of stellar mass for a given halo mass (e.g. Martizzi et al. 2012a, 2014; Vogelsberger et al. 2013; Pillepich et al. 2018a). The stellar mass excess seen in Ponos is on scales with r ≲ 20 kpc but scales larger than r ≳ 20 kpc are in good agreement with HSC/Illustris/TNG100.

5.4 Comparison of weak-lensing profiles

The main results in this paper stem from our ability to perform comparisons at fixed halo mass. In order to demonstrate the validity of our halo mass matching procedure, we also show direct comparisons of the weak lensing signals in HSC and TNG100.

We bin galaxies in HSC according to |$\overline{M}_{\rm vir}$| and measure the stacked g–g lensing observable, ΔΣ. The results are displayed in Fig. 7. To demonstrate the validity of the ASAP model, we also show the weak lensing signal as measured for the best-fitting ASAP applied to the SMDPL simulation (see Section 3.2)

For TNG100, ΔΣ was measured with halotools (v0.7; Hearin et al. 2017) using particle data. To estimate the sample variance error on ΔΣ for TNG which is smaller in volume compared to HSC, we use samples drawn from the MultiDark Planck 2 simulation (MDPL2; Klypin et al. 2016). We divide the full 1 Gpc h^–1 volume of MDPL2 into several 75 Mpc h^–1 volumes corresponding to the box size of TNG100 and measure the standard deviation of the weak-lensing signal across the various smaller volumes in each halo mass bin. The red error bars shown in Fig. 7 corresponds to this standard deviation at each radial bin across the entire MDPL2 volume. The HSC weak lensing profile and the one from TNG100 are in good agreement given the sample variance errors associated with TNG100. This suggests that our method for assigning a halo mass to HSC galaxies is effective. More precise comparisons at fixed halo mass will necessitate larger volume hydro simulations.

It is worth noting that in TNG100, we only consider central galaxies, but in observation, we could have |$\sim 10{{\ \rm per\ cent}}$| satellite galaxy contamination in our sample, which could affect the inner |$\Delta \Sigma $| profile. However, we have checked that the impact of satellites on ΔΣ is negligible compared to the TNG errors on ΔΣ (Huang et al. in preparation).

6 DISCUSSION

We have shown that the full stellar mass surface density profiles of massive galaxies are more informative and general than summarizing galaxies by a single number for “stellar mass”. We have compared new data from HSC with galaxy formation models from three different suites of hydrodynamic simulations and find that while there is impressive agreement with some simulations, there are still important differences to point out. We now discuss some explanations for our results and how they may present a new way to constrain feedback models.

6.1 The dual impact of AGN and stellar feedback on the stellar mass profiles of massive galaxies

One way to explain the mass offsets of galaxies from Illustris/TNG100 is by appealing to the feedback models in the simulations. Feedback regulates the rate of star formation and the distribution of stellar mass within a galaxy. For massive galaxies (⁠|$M_{\rm halo} \gtrsim 10^{12} \, \mathrm{M}_{\odot }$|⁠), AGN feedback is thought to be the dominant mode (e.g. Croton et al. 2006), while in less massive galaxies (⁠|$M_{\rm halo} \lesssim 10^{12} \, \mathrm{M}_{\odot }$|⁠) stellar feedback is most important (e.g. Hopkins, Quataert & Murray 2012; Li et al. 2018; Wang et al. 2020).

Present day massive elliptical galaxies (like the ones studied in this paper) grow both via in situ star formation as well as from ex-situ growth via mergers with lower mass galaxies. The in situ and ex-situ components will dominate the light profiles of massive galaxies in different radial ranges (Fig. 8); in situ mass is expected to dominate the centre of the galaxy and the ex-situ mass becomes dominant at larger radii, though this transition radius moves inwards at higher masses (Rodriguez-Gomez et al. 2016; Pillepich et al. 2018b).

AGN feedback limits the growth from in situ star formation and thus impacts the inner regions of massive galaxies (Martizzi et al. 2014; Peirani et al. 2017). Stellar feedback is dominant in regulating the mass of the satellites that are accreted to the central galaxy and that deposit stellar mass in the outskirts, thus affecting the outer ex-situ component of the central galaxy’s stellar mass profile. By z = 0, in massive galaxies like the ones in our sample (⁠|$M_{\star } \gt 10^{11.2} \, \mathrm{M}_{\odot }$|⁠), more than 40 per cent (and up to 80 per cent in the most massive systems) of the stellar mass comes from mergers (i.e. from their ex-situ component). This has been shown to be the case in Illustris and IllustrisTNG. (e.g. Rodriguez-Gomez et al. 2015, 2016; Pillepich et al. 2018b; Tacchella et al. 2019). Furthermore, in Illustris, independent of galaxy mass, about half of that ex-situ mass will come from major mergers (1:4 mass ratio), |$\sim 40{{\ \rm per\ cent}}$| from minor and very minor mergers, and the remaining |$\sim 10{{\ \rm per\ cent}}$| is stripped from surviving galaxies (Rodriguez-Gomez et al. 2016).

In the original Illustris, comparisons with cosmic star formation density measurements compiled from various multiwavelength surveys (Behroozi, Wechsler & Conroy 2013) suggest that the AGN feedback recipe does not quench star-formation sufficiently at low redshift (Vogelsberger et al. 2014a). Particularly, in massive haloes, the radio-mode feedback, which outputs kinetic energy to mimic the effects of relativistic radio jets in X-ray cavities (“bubble” mode), is not effective at halting star formation. For this reason, in IllustrisTNG the implementation of low and high accretion feedback modes was modified, and the bubble mode was replaced with a kinetic mode that acts at the very centre of galaxies. Indeed, Fig. 2 shows that TNG100 SMF is improved compared to Illustris (also see Pillepich et al. 2018a,b). At the lower mass end (⁠|$M_{\star } \lesssim 10^{10} \, \mathrm{M}_{\odot }$|⁠), this improvement is mostly due to the changes in the implementation of stellar-driven winds, while at higher masses, it is due to the modifications to the AGN feedback, both of which suppress the SMF at their respective mass scales. Our new measurements create a new opportunity to yield qualitative new insight into the physical nature of feedback processes. In the following, we discuss how our findings can be interpreted in terms of feedback models. We focus on the results from TNG100 and propose two possible solutions to the relatively mild (dis)agreements between TNG100 and HSC galaxies:

• In centrals: AGN feedback in central galaxies needs to be adjusted to allow for a higher in situ fraction in lower mass haloes (M_vir ∼ 10¹³ M_⊙) while preventing an excess in higher mass haloes (M_vir ∼ 10¹⁴ M_⊙).

• In satellites that will eventually merge with centrals in dark matter haloes with M_vir ∼ 10¹⁴ M_⊙: feedback from stars and/or AGN in lower mass galaxies (M_* ∼ 10⁹ − 10¹⁰ M_⊙) needs to be modified to suppress their mass growth at early times. This would reduce both the inner regions and the stellar haloes for central galaxies of dark matter haloes with M_vir = 10¹⁴ M_⊙, bringing TNG100 into closer agreement with HSC.

6.1.1 In situ component

Previous work has demonstrated the significant impact AGN feedback can have on the density profiles of massive elliptical galaxies through its in situ component (Fan et al. 2008; Peirani, Kay & Silk 2008; Duffy et al. 2010; Martizzi et al. 2012a, b; Peirani et al. 2017; Wang et al. 2019). Once the galaxy is massive enough and the AGN turns on, radiation and winds from the accretion process at the galactic nucleus will be emitted outwards, injecting energy into the galaxy and preventing the cooling of gas to allow for star formation (e.g. Ishibashi & Fabian 2012). When the AGN feedback is absent or not efficient enough at quenching star formation at these masses, there is an overcooling in the central region of these galaxies resulting in increased star formation and excess mass in the inner regions of the galaxy (Martizzi et al. 2014; Peirani et al. 2017). Fig. 8 shows a schematic of how this might affect the stellar mass profiles of massive galaxies. AGN feedback that is not efficient at quenching star formation will result in overcooling and a prominent in situ component. More efficient AGN feedback would limit star formation and would result in a lower in situ component.

The inner density profiles of galaxies from HSC can be used as observational constraints on the strength of AGN feedback in central galaxies. In doing so, however, it is important to keep in mind that towards the highest galaxy masses (M_* ≳ 10^11.5 M_⊙), the ex-situ stellar mass becomes the dominant component at progressively smaller galactocentric distances (see next section). The effect of weak to no-AGN feedback is revealed in our comparison with the Ponos simulations (Section 5.3.2). As in Martizzi et al. (2014), the lack of AGN feedback resulted in galaxies that are too massive, particularly in their interior region (r ≲ 20 kpc), compared to observations. On larger scales, though, there is good agreement with HSC/Illustris/TNG100. This is interesting as it suggests that the inner profiles of massive galaxies contain information about the physical scale of the impact of AGN feedback on in situ star formation. The energy injected by the AGN can provide a mechanism for regulating star formation, either by preventing the cooling of gas or by expelling gas from the central regions of galaxies. The details of how and when AGN feedback occurs and its impact on galaxy evolution remain uncertain from both an observational and theoretical point of view. Our analysis allows us to estimate the galaxy scale over which the deposition of energy and momentum via outflows and radiation from the AGN affects the formation of stars in the host galaxy. The comparison shown in Fig. 6 would suggest that AGN feedback is limiting star formation in massive galaxies at roughly r ≲ 20 kpc.

Turning now to TNG100, Fig. 2 displays the |$M_{*}^{10}-M_{\rm vir}$| relation and reveals that TNG100 has a steeper relation than HSC with |$M_{*}^{10}$| being lower than HSC at M_vir = 10¹³ M_⊙ and |$M_{*}^{10}$| being higher than HSC at M_vir = 10¹⁴ M_⊙. This could be an indication that TNG100 is too efficient at quenching star formation at M_vir = 10¹³ M_⊙ but not efficient enough at M_vir = 10¹⁴ M_⊙. We will study this possibility further in future work by using information about the origin of the stellar particles in question (i.e. whether they were formed in situ or ex-situ; Chowdhury et al. in preparation).

6.1.2 Ex-situ component

Massive galaxies experience significant growth via merging. This growth can significantly impact the stellar mass profiles of massive galaxies given that stars that were formed in situ are generally closest to the centre of the galaxy, followed by stars accreted in major mergers, minor mergers, very minor mergers, and finally stars that were stripped from surviving galaxies (e.g. Rodriguez-Gomez et al. 2015, 2016; Pillepich et al. 2018b; Tacchella et al. 2019). Therefore, we can appeal to processes that add ex-situ stars to a galaxy to help explain what we see in the outer regions of a galaxy. Overly massive satellites merging with the central may explain the cases in which Illustris and TNG100 galaxies have more extended outskirts compared to HSC. We typically observe this in the most massive halo mass bin.

The most massive galaxies in our sample live in dark matter haloes of masses log (M_vir/M_⊙) ≳ 14. In TNG100 (and TNG300), 90 per cent of the total ex-situ mass of central galaxies (including the ICL component) in haloes of these masses is built from progenitors of masses log (M_⋆/M_⊙) ≳ 9−10 (Pillepich et al. 2018b). Most of that ex-situ growth happens through mergers at z ≲ 1 (Genel et al. 2018). While TNG100 has much better agreement with observations than Illustris, there is still a noted excess of lower mass (log M_⋆/M_⊙ ≲ 10) galaxies at z ∼ 1 compared to observations (e.g. fig. 14 in Pillepich et al. 2018b). It is, therefore, possible that the excess outer mass we see in the highest halo mass bin in Illustris/TNG100 galaxies compared to HSC is related to an excess ex-situ component built up from overly massive satellite galaxies.

6.2 Artificial stripping and disruption

Another possibility is that there are non-physical effects in simulations contributing to the excess buildup of mass in the outskirts of galaxies. For example, van den Bosch (2017); van den Bosch et al. (2018); van den Bosch & Ogiya (2018) have shown that at the numerical resolutions typical of large-volume cosmological simulations, DM subhaloes may undergo overdisruption and overmerging with their central hosts because of numerical effects. Admitting that those studies focus on idealized orbits and in the absence of gas and stars, they nevertheless inspire a critical approach for the interpretation of our findings. As the complete disruption of satellites is not of relevance for this analysis, we instead discuss two hypothetical scenarios: the first is when stellar mass is removed from satellites which should in reality remain gravitationally bound to satellites; the second is when stellar mass is stripped in the right amount from satellites but deposited in the wrong location (too early or too late) along their orbits.

In the first scenario, satellite galaxies are being artificially overstripped and this stripped mass adds to the stellar halo (whereas in reality this mass should still be counted to be part of the satellite). van den Bosch & Ogiya (2018) show, albeit in idealized configurations, that, as subhaloes are stripped of their DM, the first 99 per cent of their mass loss is accurately captured also at the numerical resolutions relevant for this paper. This may indicate that the largest majority of the stripped material (in this case stellar mass) is stripped in a quantitatively consistent manner in simulations like Illustris, IllustrisTNG, or Ponos as in extremely high-resolution controlled experiments, i.e that the lion’s share of the stellar mass in the outskirts of massive galaxies should be numerically converged when stripping is concerned (see also Lovell et al. in preparation).

In the second scenario, mass is being stripped in the right amounts, but deposited in the wrong location (e.g. earlier along orbits). Presumably, this would preferentially deposit mass in the outskirts of centrals (at tens or even hundreds of kpc from the centre), leading to an excess of mass at large radii and a deficit of mass at small radii, but the total mass would remain unchanged. In this paper, we find that in our highest mass bin (M_* ∼ 10¹⁴ M_⊙), galaxies in TNG100 are overly massive in total stellar mass (i.e. within 100 kpc) compared to HSC (e.g. Fig. 2). Hence, we believe that this is unlikely to explain the differences between TNG100 and HSC.

6.3 Uncertainties in mass measurements

It is also important to take into account the assumptions made in measuring masses from our observed profiles. To convert from surface brightness to mass density, we assume a constant M/L ratio for the entire galaxy. Huang et al. (2018b) discuss why this assumption is justified given that our sample of massive ellipticals is dominated by old stellar populations and known to have only shallow colour gradients (Appendix C in Huang et al. 2018b). None the less, there is a known shallow but negative M/L gradient driven by the metallicity gradient in these galaxies (e.g. La Barbera et al. 2012; D’Souza et al. 2015, Chowdhury et al. in preparation). If we instead assume a negative M/L gradient, that would result in decreased masses in the outer regions of our observed HSC galaxies. This would typically result in an even larger discrepancy with the simulated galaxies.

We also consider the limitations of our 1D profile measurement approach. If instead we considered a 2D functional form to approximate the distributions of light (and mass), we would be able to use multiple components to take into account variations in the isophotal shape. We could also integrate it to infinity to derive a true “total” mass. We could also take seeing into account. However, we would be limited by the fact that we do not know the perfect functional form (or combination of components) to describe the outskirts of massive galaxies. Multiple-Sersic components can be promising, but internal degeneracy and the lack of clear physical meaning are the downside. These will be explored in a future work (Ardila et al. in preparation).

7 SUMMARY AND CONCLUSIONS

In this paper, we perform a consistent comparison of the masses and mass profiles of massive (⁠|$M_\star \gt 10^{11.4}\, \mathrm{M}_\odot$|⁠) central galaxies at z ∼ 0.4 from deep HSC observations and from the Illustris, TNG100, and Ponos simulations. HSC weak lensing measurements further enable comparisons at fixed halo mass using a tight relation with |$M_{\star }^{\mathrm{10}}$|⁠, and |$M_{\star }^{\mathrm{100}}$| from the ASAP empirical model. We measure the stellar mass density profiles, the galaxy SMFs, and the stellar-to-halo mass relations using different operational definitions of stellar masses for massive galaxies in these simulations and compare with measurements from the sample of Huang et al. (2018b). Stellar profiles are measured by drawing concentric elliptical isophotes with a fixed ellipticity on the target galaxy, at given radii along the semimajor axis. In the simulations, we use stellar mass maps and perform measurements in a way that closely mimics the HSC methodology. Our main results are summarized below:

• The direct comparison of galaxy light profiles avoids having to define “galaxy mass” as a single number and can account for surface brightness limits in observations. The full mass profile also contains more information than mass at any fixed radius. This is shown in Fig. 2 which compares the galaxy SMF and the SHMR for various aperture definitions of stellar mass. It is clear that different simulations match HSC data better at different radii and in different mass bins. Understanding the radial ranges over which simulations and observations agree/disagree provides important clues about the assembly history and feedback processes that have shaped massive galaxies.

• We also compare the weak lensing signals of galaxies from HSC with Illustris/TNG100 (Fig. 7). This comparison, together with the ASAP model developed in Huang et al. (2020), enables us to perform comparisons at fixed halo mass. The ability to control for halo mass is a key component of this paper as it enables us to draw more informative conclusions about the strength and impact of feedback at different halo mass scales.

• Both Illustris and TNG100 display overall good agreement with the mass profiles between a few and 100 kpc of massive central galaxies in HSC, albeit with some interesting differences (Fig. 5). Generally, central galaxies in both Illustris and TN100 exhibit steeper stellar-to-halo mass relations, both within 10 and 100 kpc apertures (Fig. 2). This disagreement against the HSC observational inferences are more severe in the innermost regions (r < 10 kpc), at higher masses (⁠|$M_\star \sim 10^{14}\, \mathrm{M}_\odot$|⁠), and in Illustris than in TNG100. In TNG100, the stellar mass distributions of galaxies in |$M_\star \sim 10^{13.5}\, \mathrm{M}_\odot$| haloes are in excellent agreement with observations, while those in |$M_\star \sim 10^{13}\, \mathrm{M}_\odot$| haloes have a deficit of stellar mass (by 0.05–0.15 dex, depending on radius) and those in haloes of |$M_\star \sim 10^{14}\, \mathrm{M}_\odot$| have an excess (by 0.15–0.30 dex). In fact, TNG100 returns an outer stellar mass (⁠|$M_{*}^{100}-M_{*}^{10}$|⁠) that is consistent with the SMDPL+ASAP results to better than 0.12 dex across the halo mass range studied here.

• We interpret our results by assuming the two-phase formation scenario of massive galaxies in which the inner regions of |$M_\star \sim 10^{11}\, \mathrm{M}_\odot$| galaxies are dominated by in situ mass and the outer regions are dominated by ex-situ mass (Fig. 8). Given that both the inner and outer stellar masses show offsets in some halo mass bins in both Illustris and TNG100, the culprit for the (dis)agreement may be found in both in situ and ex-situ stellar mass processes. The steeper |$M_{*}^{10}-M_{\rm vir}$| relation in TNG100 than in HSC (TNG100 |$M_{*}^{10}$| is lower than HSC at |$M_{\rm vir}\sim 10^{13}\, \mathrm{M}_\odot$| and higher at |$M_{\rm vir}\sim 10^{14}\, \mathrm{M}_\odot$|⁠; Fig. 2) may be signaling that quenching in the TNG100 model is too efficient at |$M_{\rm vir}\sim 10^{13}\, \mathrm{M}_\odot$| but not efficient enough for the stars that end up assembling in M_vir ∼ 10¹⁴ haloes. On the other hand, central galaxies in |$M_{\rm vir}\sim 10^{14}\, \mathrm{M}_\odot$| haloes display too much outer mass in TNG100 compared to HSC (albeit by only ∼0.12 dex in |$M_{*}^{100}-M_{*}^{10}$|⁠). We hypothesize that feedback from stars and/or AGN needs to be adjusted to suppress the mass growth at early times of satellites that build up the stellar halo (⁠|$M_\star \sim 10^{9-10} \, \mathrm{M}_{\odot }$|⁠). Resolution effects could also be at play and warrant further investigation.

• Galaxies in the zoom-in Ponos simulation, which does not implement AGN feedback, display a substantial excess mass of ∼0.5 dex at r < 30 kpc compared to HSC galaxies of similar halo mass (though a comparable slope and amplitude for the rest of the profile). This is indicative of overcooling and excess star formation in the central regions due to the lack of regulation from the central BH in the Ponos simulations Fig. 6). The comparison between Ponos and HSC suggests that the physical scale over which the central AGN limits star formation is r ≲ 20 kpc.

• We have also performed a number of tests to validate our results. In particular, we showed the choice of method used to mask neighbouring satellite galaxies does not impact our measurements. This is due to the iterative 3σclipping that we apply along each isophote prior to measuring the mean intensity. Even when a satellite mask is not applied, the impact of satellites represents less than a 2 per cent effect for 95 per cent of our sample at 150 kpc.

• We provide median surface stellar mass density profiles and weak-lensing ΔΣ profiles for the HSC galaxies (https://github.com/f-ardila/HSC_versus_hydro-paper) to facilitate comparisons with other simulations.

In the era of large surveys that are wide, deep, and that have lensing capability (HSC, LSST, Euclid, WFIRST), joint comparisons between weak lensing and galaxy stellar profiles offer a powerful method for comparing simulations with observations. In particular, this combination provides a direct test of whether simulations build and deposit galaxy mass in the correct dark matter haloes and therefore constrains the physics of feedback and galaxy growth.

ACKNOWLEDGEMENTS

FA is supported by the National Science Foundation Graduate Research Fellowship Program under grant no. DGE-1842400. This research was supported in part by the National Science Foundation under grant no. NSF PHY-1748958. This material is based upon work supported by the National Science Foundation under grant no. 1714610. This research was also supported in part by National Science Foundation under grant no. NSF PHY11-25915. AL acknowledges support from the David and Lucille Packard Foundation, and from the Alfred P. Sloan Foundation. We acknowledge use of the lux supercomputer at UC Santa Cruz, funded by NSF MRI grant AST 1828315. PM acknowledges a NASA contract supporting the WFIRST-EXPO Science Investigation Team (15-WFIRST15-0004), administered by the GSFC. Work done at Argonne National Laboratory was supported under the DOE contract DE-AC02-06CH11357.

Based [in part] on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by the Subaru Telescope and Astronomy Data Center at National Astronomical Observatory of Japan.

The Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University.

The Pan-STARRS1 Surveys (PS1) have been made possible through contributions of the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under grant no. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation under grant no. AST-1238877, the University of Maryland, and Eotvos Lorand University (ELTE).

This paper makes use of software developed for the Large Synoptic Survey Telescope. We thank the LSST Project for making their code available as free software at http://dm.lsst.org.

The following software and programming languages made this research possible: python (v2.7); astropy (v2.0; Astropy Collaboration et al. 2013, 2018), a community-developed core python package for Astronomy; Source Extraction and Photometry in python (SEP v1.0; Barbary, Boone & Deil 2015); halotools (v0.7; Hearin et al. 2017); COsmology, haLO and large-Scale StrUcture toolS (COLOSSUS; Diemer 2018).

FA would like to express deep gratitude to his grandfather, Andrés Rosselli Quijano, for all the knowledge and passion he passed on, as well as the inspiration to pursue this subject area.

DATA AVAILABILITY

The observational HSC data underlying this article are publicly available at https://github.com/f-ardila/HSC_versus_hydro-paper. Illustris/TNG simulation data are publicly available at https://www.tng-project.org/data/.

Footnotes

REFERENCES