Can Assembly Bias Explain the Lensing Amplitude of the BOSS CMASS Sample in a Planck Cosmology?

TL;DR

This study tests if galaxy assembly bias can resolve the discrepancy between observed galaxy clustering and lensing signals in the BOSS CMASS sample within a Planck cosmology, finding it insufficient to explain the difference.

Contribution

The paper demonstrates that incorporating galaxy assembly bias into extended HOD models does not reconcile the observed clustering and lensing signals at Planck cosmology.

Findings

Standard HOD models over-predict lensing signal by 20-40%.

Assembly bias models do not resolve the clustering-lensing discrepancy.

A persistent 31-34% mismatch remains despite extended models.

Abstract

In this paper, we investigate whether galaxy assembly bias can reconcile the 20-40% disagreement between the observed galaxy projected clustering signal and the galaxy-galaxy lensing signal in the BOSS CMASS galaxy sample reported in Leauthaud et al. (2017). We use the suite of AbacusCosmos Lamda-CDM simulations at Planck best-fit cosmology and two flexible implementations of extended halo occupation distribution (HOD) models that incorporate galaxy assembly bias to build forward models and produce joint fits of the observed galaxy clustering signal and the galaxy-galaxy lensing signal. We find that our models using the standard HODs without any assembly bias generalizations continue to show a 20-40% over-prediction of the observed galaxy-galaxy lensing signal. We find that our implementations of galaxy assembly bias do not reconcile the two measurements at Planck best-fit cosmology. In fact, despite incorporating galaxy assembly bias, the satellite distribution parameter, and the satellite velocity bias parameter into our extended HOD model, our fits still strongly suggest a 31-34% discrepancy between the observed projected clustering and galaxy-galaxy lensing measurements. It remains to be seen whether a combination of other galaxy assembly bias models, alternative cosmological parameters, or baryonic effects can explain the amplitude difference between the two signals.

Figures

Tables

Table 1. Table 1: The standard HOD parameters and their baseline values as quoted from Alam et al. ( 2017 ) .

Parameter name	Baseline value
${\log }_{10}\left(M_{\text{cut}}/h^{-1}{M}_{\odot }\right)$	$13.248$
${\log }_{10}\left(M_1/h^{-1}{M}_{\odot }\right)$	14.179
$\sigma$	0.897
$\alpha$	1.151
$\kappa$	0.137

Table 2. Table 2: The parameter ranges for the 3 tiers of HOD pairs we have generated for training the generalized HOD emulator. Each bracket shows the minimum and maximum bounds in that parameter for that tier. Our mocks are generated from HODs uniformly sampled from within these bounds. As discussed in the main text, each tier increases the bounds in all the HOD parameters, and this tiered approach ensures a higher density of samples in the region closest to the baseline HOD. The log 10 ⁡ M cut subscript 10 subscript 𝑀 cut \log_{10}M_{\textrm{cut}} notation is a shorthand for log 10 ⁡ ( M cut / h − 1 ​ M ⊙ ) subscript 10 subscript 𝑀 cut superscript ℎ 1 subscript 𝑀 direct-product \log_{10}(M_{\textrm{cut}}/h^{-1}M_{\odot}) for brevity.

Parameter	Tier 1 (15 pairs)	Tier 2 (75 pairs)	Tier 3 (120 pairs)
${\log }_{10}{M}_{\text{cut}}$	[13.225, 13.269]	[13.151, 13.327]	[12.748, 13.748]
${\log }_{10}{M}_1$	[14.157, 14.200]	[14.082, 14.258]	[13.679, 14.679]
$\sigma$	[0.847, 0.947]	[0.697, 1.097]	[0.497, 1.297]
$\alpha$	[1.101, 1.201]	[0.951, 1.351]	[0.751, 1.551]
$s$	[-0.05, 0.05]	[-0.2, 0.2]	[-0.3, 0.3]
$A$	[-0.05, 0.05]	[-0.2, 0.2]	[-0.3, 0.3]

Table 3. Table 3: The parameter ranges for the 4 tiers of HOD pairs we have generated for training the decorated HOD emulator. Each bracket shows the minimum and maximum bounds in that parameter for that tier. Our mocks are generated from HODs uniformly sampled from within these bounds. The log 10 ⁡ M cut subscript 10 subscript 𝑀 cut \log_{10}M_{\textrm{cut}} notation is a shorthand for log 10 ⁡ ( M cut / h − 1 ​ M ⊙ ) subscript 10 subscript 𝑀 cut superscript ℎ 1 subscript 𝑀 direct-product \log_{10}(M_{\textrm{cut}}/h^{-1}M_{\odot}) for brevity.

Parameter	Tier 1 (15 pairs)	Tier 2 (45 pairs)	Tier 3 (120 pairs)	Tier 4 (120 pairs)
${\log }_{10}{M}_{\text{cut}}$	[13.225, 13.269]	[13.151, 13.327]	[13.048, 13.448]	[12.748, 13.748]
${\log }_{10}{M}_1$	[14.157, 14.200]	[14.082, 14.258]	[13.979, 14.379]	[13.679, 14.679]
$\sigma$	[0.847, 0.947]	[0.697, 1.097]	[0.697, 1.097]	[0.497, 1.297]
$\alpha$	[1.101, 1.201]	[0.951, 1.351]	[0.951, 1.351]	[0.751, 1.551]
$A_{\text{cent}}$	[-0.05, 0.05]	[-0.2, 0.2]	[-0.2, 0.2]	[-0.3, 0.3]
$A_{\text{sat}}$	[-0.05, 0.05]	[-0.2, 0.2]	[-0.2, 0.2]	[-0.3, 0.3]

Table 4. Table 4: The prior information for the generalized HOD model. We choose the priors to be Gaussians centered on the baseline values with broad non-informative width. The baseline values are taken from Alam et al. ( 2017 ) .

Parameter name	${\mu }_{\text{prior}}$	${\sigma }_{\text{prior}}$
${\log }_{10}\left(M_{\text{cut}}/h^{-1}{M}_{\odot }\right)$	13.2	0.3
${\log }_{10}\left(M_1/h^{-1}{M}_{\odot }\right)$	14.2	0.3
$\sigma$	0.897	0.1
$\alpha$	1.151	0.1
$s$	0	0.1
$A$	0	0.1

Table 5. Table 5: The prior specification for the decorated HOD model. We choose the priors to be Gaussians centered on the baseline values with broad non-informative width.

Parameter name	${\mu }_{\text{prior}}$	${\sigma }_{\text{prior}}$
${\log }_{10}\left(M_{\text{cut}}\left[h^{-1}{M}_{\odot }\right]\right)$	13.2	0.3 (truncated at $\pm 1$)
${\log }_{10}\left(M_1\left[h^{-1}{M}_{\odot }\right]\right)$	14.2	0.3 (truncated at $\pm 1$)
$\sigma$	0.897	0.1
$\alpha$	1.151	0.1
$A_{\text{cent}}$	0	0.1
$A_{\text{sat}}$	0	0.1

Table 6. Table 6: Summary of the evidence log ⁡ 𝒵 𝒵 \log\mathcal{Z} and their uncertainties of the different models we tested. The two columns correspond to the model without a scaling parameter and with a scaling parameter. The first three rows show the generalized HOD implementation in Section  5.1 , with the first row showing the standard 5-parameter model without any generalized parameter, and the next two showing the standard model plus the parameters in parentheses. The next two rows show the evidence of the decorated HOD implementation provided in Halotools .

	No scaling	scaling
generalized HOD (standard)	-69.46$\pm$0.07
generalized HOD ($s$, $A$)	-65.09$\pm$0.08	-16.00$\pm$0.08
generalized HOD ($s$, $A$, $s_v$)	-65.18$\pm$0.08	-16.20$\pm$0.08
decorated HOD (standard)	-74.91$\pm$0.07
decorated HOD ($A_{\text{cent}}$, $A_{\text{sat}}$)	-72.89$\pm$0.08	-17.92$\pm$0.08

Equations

$${\bar{n}}_{\mathrm{cent}}=\frac{1}{2}\mathrm{erfc}\left[\frac{\ln (M_{\mathrm{cut}}/M)}{\sqrt{2}\sigma }\right],$$

$${\bar{n}}_{\text{sat}}={\left[\frac{M-\kappa M_{\text{cut}}}{M_1}\right]}^{\alpha }{\bar{n}}_{\mathrm{cent}},$$

$$c=\frac{r_{\text{vir}}}{r_{s,\text{Klypin}}},$$

$$\log M_{\text{pseudo}}=\{\begin{array}{ll}\log M+A& \text{if}c>\bar{c}(M),\\ \log M-A& \text{if}c<\bar{c}(M),\end{array}$$

$${\bar{n}}_{\text{sat}}={\left[\frac{M-\kappa M_{\text{cut}}}{M_1}\right]}^{\alpha },$$

$$⟨n_{\text{gal}}|M,c>c_{\text{piv}}⟩=$$

$$⟨n_{\text{gal}}|M,c<c_{\text{piv}}⟩=$$

$$\delta n_{\text{cent}}=A_{\text{cent}}\text{min}\left[⟨n_{\text{cent}}|M⟩,1-⟨n_{\text{cent}}|M⟩\right]$$

$$\delta n_{\text{sat}}=A_{\text{sat}}⟨n_{\text{sat}}|M⟩$$

$$w_p(r_{\perp })=2{\int }_0^{{\pi }_{\text{max}}}\xi (r_{\perp },\pi )d\pi ,$$

$$\Delta \Sigma (r_{\perp })=\overline{\Sigma }(<r_{\perp })-\overline{\Sigma }(r_{\perp }),$$

$$-0.2<{\log }_{10}\left(\frac{w_p(r_{\perp }<1h^{-1}\text{Mpc})}{w_{p,0}(r_{\perp }<1h^{-1}\text{Mpc})}\right)<0.2,$$

$$-0.15<{\log }_{10}\left(\frac{w_p(r_{\perp }<1h^{-1}\text{Mpc})}{w_{p,0}(r_{\perp }<1h^{-1}\text{Mpc})}\right)<0.15,$$

$$\text{centering}{r}_{\perp }{w}_{p,i}({\mathbf{p}}_0+\delta \mathbf{p})=\text{@add@centering}$$

$$+\frac{1}{2}{\Sigma }_{j,k=0}^7\frac{{\partial }^2{r}_{\perp }{w}_{p,i}}{\partial p_j\partial p_k}\delta p_j\delta p_k,$$

$$\text{centering}{r}_{\perp }\Delta {\Sigma }_i({\mathbf{p}}_0+\delta \mathbf{p})=\text{@add@centering}$$

$$+\frac{1}{2}{\Sigma }_{j,k=0}^7\frac{{\partial }^2{r}_{\perp }\Delta {\Sigma }_i}{\partial p_j\partial p_k}\delta p_j\delta p_k,$$

$$\log L(D|\Theta ,M)=\log L(w_p|\Theta ,M)+\log L(\Delta \Sigma |\Theta ,M),$$

$$\log L(w_p|\Theta ,M)=-\frac{1}{2}\delta {\mathbf{w}}_p^T{\mathbf{C}}^{-1}\delta {\mathbf{w}}_p,$$

$$\delta {\mathbf{w}}_p=\delta {\mathbf{w}}_p^{\text{obs}}-\delta {\mathbf{w}}_p^{\text{emulated}}(\Theta ).$$

$$\mathcal{Z}=P(D|M)={\int }_{{\Omega }_{\Theta }}P(D|\Theta ,M)P(\Theta |M)d\Theta .$$

Full text

††pagerange: Can Assembly Bias Explain the Lensing Amplitude of the BOSS CMASS Sample in a Planck Cosmology?–LABEL:lastpage

1 Introduction

Weak gravitational lensing refers to the distortions in the images of distant galaxies by intervening mass along the line of sight. Because it directly measures the total mass distribution of the Universe, it has long been considered a powerful yet unique cosmological probe. Galaxy-galaxy lensing (hereafter “g-g lensing") refers to the cross-correlation between foreground lens galaxy positions and the lensing shear of background source galaxies (Tyson et al., 1984; Brainerd et al., 1996; dell’Antonio & Tyson, 1996; Prat et al., 2018). At small scales, it provides a measure of the radial distribution of total mass around galaxies, presenting an unique opportunity to directly probe the properties of dark matter halos.

Recent surveys such as the Sloan Digital Sky Survey (SDSS, York et al., 2000), Dark Energy Survey (DES, Dark Energy Survey Collaboration et al., 2016), the Canada-France-Hawaii Telescope Lensing Survey (CFHTlenS, Heymans et al., 2012; Erben et al., 2013), the Kilo-Degree Survey (KiDS, de Jong et al., 2013; Kuijken et al., 2015), and the Hyper Suprime Cam survey (HSC, Aihara et al., 2018) have generated thousands of square degrees of high signal-to-noise g-g lensing data (e.g. Mandelbaum et al., 2013; Velander et al., 2014; Abbott et al., 2018; van Uitert et al., 2018). Upcoming missions such as Euclid (Laureijs et al., 2011), the Wide Field Infrared Survey Telescope (WFIRST, Spergel et al., 2013), and the Large Synoptic Survey Telescope (LSST, LSST Science Collaboration et al., 2009) promise to bring in even higher precision data over a vast fraction of the sky.

In parallel to these efforts, surveys such as the Baryon Oscillation Spectroscopic Survey (BOSS, Eisenstein et al., 2011; Dawson et al., 2013), have collected optical spectra for more than one million massive galaxies at $z<1$. These spectra enabled accurate measurements of galaxy clustering in the form of the 2-point correlation function (2PCF), placing tight cosmology constraints. Upcoming experiments such as the Dark Energy Spectroscopic Instrument (DESI, Levi et al., 2013), the Prime Focus Spectrograph (PFS, Takada et al., 2014), and Euclid (Laureijs et al., 2011) will measure the redshifts of tens of millions of galaxies, yielding exquisite measurements of galaxy clustering and also providing excellent lens samples for g-g lensing studies.

While galaxy clustering and g-g lensing represent two independent yet complementary cosmology probes, Leauthaud et al. (2017) find discrepancies of 20-40 percent between their measurements of g-g lensing for CMASS galaxies and a model predicted from mock galaxy catalogs generated at Planck cosmology that match the CMASS projected correlation function (Reid et al., 2014; Saito et al., 2016). Lange et al. (2019) extended this result by finding a similar $\sim 25\%$ discrepancy between the projected clustering measurement and the g-g lensing measurement in the BOSS LOWZ sample. They also found that this discrepancy is independent of redshift ($0.1<z<0.7$ and stellar mass ($11<\log M_{\star}/M_{\odot}<12$) in the BOSS CMASS and LOWZ sample. This discrepancy is well above the statistical error of the lensing signal and calls for a detailed re-examination of the forward model used to predict the g-g lensing signal.

Leauthaud et al. (2017) found that lowering the cosmological parameter $S_{8}=\sigma_{8}\sqrt{\Omega_{m}/0.3}$ by 2-3 $\sigma$ from the Planck 2015 value can reconcile the difference. However, cosmological effects are entangled with other effects due to details of galaxy-halo connection, baryons, and massive neutrinos. Before one can draw an inference about cosmological models, one must control for uncertainties in the astrophysical modeling of these other effects. One source of modeling uncertainties is the assumed connection between the observed galaxies and their dark matter halos. Leauthaud et al. (2017) use a standard Halo Occupation Distribution (HOD) model of the Zheng & Weinberg (2007) form. The mis-match in the amplitude of the g-g lensing signal may point to the failures of such empirical models.

In particular, one important aspect that these models neglect is galaxy assembly bias: the fact that in addition to halo mass, galaxy occupation depends on other properties such as halo age, spin, and concentration (e.g. Gao et al., 2005; Wechsler et al., 2006; Gao & White, 2007; Zentner, 2007; Dalal et al., 2008; Lacerna & Padilla, 2011; Zentner et al., 2014). Numerous recent studies have attempted using clustering data and simulations to detect galaxy assembly bias and constrain its effects (e.g. More et al., 2016; Miyatake et al., 2016; Saito et al., 2016; Lehmann et al., 2017; Xu & Zheng, 2018; Zehavi et al., 2018; Contreras et al., 2019; Zentner et al., 2019). Galaxy assembly bias is especially relevant for the g-g lensing discrepancy because the clustering measurements tightly constrain the large-scale galaxy bias, whereas the lensing measurement is mostly sensitive to the dark matter halo mass profile. Thus, the g-g lensing discrepancy is fundamentally a discrepancy between halo mass and large-scale galaxy bias, the signature of assembly bias.

In this paper, we explore the possibility of reconciling the g-g lensing discrepancy with two different extended HOD models incorporating galaxy assembly bias plus other halo scale physics. Specifically, we apply the generalized HOD model (GRAND-HOD Yuan et al., 2018), which incorporates a novel implementation of galaxy assembly bias plus other generalizations to the Zheng & Weinberg (2007) HOD, and the decorated HOD model (Hearin et al., 2016), which incorporates both a central assembly bias and a satellite assembly bias. We build emulator models of the projected galaxy correlation function and the g-g lensing as a function of the extended HODs at Planck cosmology. We present joint fits to the observed projected galaxy correlation function and g-g lensing to evaluate how well the generalized HOD model can reconcile the two measurements.

This paper is organized as follows. In Section 2, we introduce the standard 5-parameter HOD, the generalized HOD, and the decorated HOD. In Section 3, we present the galaxy clustering and weak lensing observables that we fit our extended HOD models to, and in Section 4 we present our forward model for emulating the observables as a function of the extended HOD parameters. Then we use these emulator models to fit the observables and present the results in Section 5. We discuss the limitations and implications of our results in Section 6. Finally, in Section 7, we summarize the main conclusions of our analysis.

Throughout this paper, unless otherwise specified, we assume a Planck 2015 (Planck Collaboration et al., 2016) cosmology with $H_{0}=67.26$ km/s/Mpc, $\Omega_{\textrm{m}}=0.3141$, and $\sigma_{\textrm{8}}=0.83$. We use halo mass definition $M_{200\textrm{b}}$, which we simply quote as $M$ for the rest of this paper.

2 Theoretical Background

The standard 5-parameter HOD model (Zheng & Weinberg, 2007) is a popular empirical framework to populate dark matter halos with mock central and satellite galaxies as a function of halo mass. However, for cosmology, this model may also be a source of systematics. In this section, we briefly review the standard HOD formalism and discuss physically motivated extensions to the standard HOD.

2.1 The standard HOD

The standard HOD (Zheng & Weinberg, 2007) gives the mean number of central galaxies and satellite galaxies as a function of halo mass

[TABLE]

where halo mass $M$ is again defined as $M_{200\textrm{b}}$. The five parameters of this model are $M_{\textrm{cut}},M_{1},\sigma,\alpha,\kappa$. However, for the rest of the paper, we only consider the first four parameters of the model as our tests show that the predicted clustering and lensing observables depend very weakly on $\kappa$. The actual number of central galaxies in a halo follows the Bernoulli distribution. The actual number of satellites follows the Poisson distribution with the mean equal to $\bar{n}_{\textrm{sat}}$. The central assumes the location and velocity of the center-of-mass of the halo. The satellites are assigned to halo particles to track the dark matter distribution, with each particle of the halo having an equal probability of hosting a satellite galaxy.

2.2 Generalized HOD

To generate predictions for galaxy clustering, we populate dark matter halos with mock galaxies using the publicly available GRAND-HOD package111https://github.com/SandyYuan/GRAND-HOD. The routine introduces five new parameters to the standard HOD, including a novel implementation of galaxy assembly bias. Yuan et al. (2018) describes the generalizations in detail. We highlight the three parameters relevant for this study.

We introduce the satellite distribution parameter $s$, which deviates the satellite spatial distribution away from the halo profile. In our standard HOD implementation, satellites are placed on halo particles with equal probabilities. However, a positive $s$ leads the algorithm to favor particles further from halo center, effectively decreasing the concentration of satellites. Figure 2 of Yuan et al. (2018) shows how $s$ affects the predicted 2PCF. The range of $s$ is defined to be between $-1$ and 1.

Similarly, we introduce the satellite velocity bias parameter $s_{v}$, which biases the satellite velocity distribution away from that of the halo. In our implementation, the satellites always assume the velocities of the particles they are placed on. A positive $s_{v}$ simply favors particles with higher velocity relative to the halo center to host satellite galaxies. By making sure that each galaxy still tracks the velocity and position of a dark matter particle, this implementation guarantees that the satellite galaxies still obey Newtonian physics in the halo potential. While peculiar velocities do not directly affect the projected correlation function, our implementation of velocity bias does change the clustering because particle velocity relative to halo center is correlated to particle position in the halo, and the higher velocity subsample of particles do not evenly trace the density profile of the halo. Figure 4 of Yuan et al. (2018) shows how $s_{v}$ affects the predicted 2PCF. The range of $s_{v}$ is defined to be between $-1$ and 1.

We also introduce a galaxy assembly bias parameter $A$. We set the NFW (Navarro et al., 1997) halo concentration as the secondary dependence for the galaxy occupation besides halo mass. The NFW concentration is defined as

[TABLE]

where $r_{\textrm{vir}}$ is the virial radius of the halo and $r_{s,\textrm{Klypin}}$ is the Klypin scale radius (Klypin et al., 2011). In our implementation, we first rank all halos by halo mass and calculate the number of galaxies $n_{\textrm{gal}}$ for each halo. We save the list of $n_{\textrm{gal}}$. Then we re-rank the halos according to a “pseudo-mass" defined as

[TABLE]

where $A$ is the assembly bias parameter which governs the strength of assembly bias in our model, and $\bar{c}(M)$ is the median concentration within a mass bin at mass $M$. Note that we do not re-rank the $n_{\textrm{gal}}$ list, just the halos themselves. Finally, we assign the numbers in the $n_{\textrm{gal}}$ list to the re-ranked list of halos in order. Effectively, for a positive $A$, we are swapping galaxies in a more massive less concentrated halo to a less massive more concentrated halo. This swapping routine ensures that the total number of galaxies is preserved in the catalog when only the assembly bias parameter is varied. However, it does not preserve the expected number of galaxies for a given halo mass $\langle\bar{n}_{g}|M\rangle$, in contrast to the standard assembly bias interpretation. Figure 6 of Yuan et al. (2018) shows the effect of $A$ on the predicted 2PCF. The range of $A$ is technically between $-\infty$ and $\infty$, but we expect $A$ to be on the order of $10^{-1}$.

In this paper we use a generalized HOD with 7 parameters: $M_{\textrm{cut}}$, $M_{1}$, $\sigma$, $\alpha$, $s$, $s_{v}$, and $A$. Again, we have ignored parameter $\kappa$. Our goal is to emulate the projected 2PCF and the g-g lensing signal as a function of these HOD parameters to search for good fits to both measurements within this generalized model space.

Again we stress that for this study we use $M_{200b}$ as our halo mass, whereas the publicly available GRAND-HOD code uses the virial mass $M_{vir}$. We also use a slightly modified formula for the mean number of satellites $\bar{n}_{\textrm{sat}}$ compared to that available on GRAND-HOD. Namely the GRAND-HOD uses

[TABLE]

whereas for this study we modulate the number of satellites with the number of centrals $\bar{n}_{\textrm{cent}}$ (see Equation 1). We make these adjustments to be more consistent with the HOD implementation of Alam et al. (2017), which fitted the standard HOD model to the CMASS projected 2PCF, though our implementation still differs in key aspects such as the satellite profile.

2.3 Decorated HOD

The decorated HOD (Hearin et al., 2016) provides another way of incorporating assembly bias into the standard HOD. The decorated HOD is conveniently implemented in the Halotools code package (Hearin et al., 2017). We borrow a slightly modified implementation of the decorated HOD implementation from Wang et al. (2019), which also uses the NFW halo concentration (Equation 2) as the secondary property of the halo. Again, they assume that $P(n_{\textrm{cen}}|M,c)$ is a Bernoulli distribution and that $P(n_{\textrm{sat}}|M,c)$ is a Poisson distribution, but that these distributions have first moments of

[TABLE]

where the $n_{\textrm{gal}}$ notation applies to both the centrals and satellites. The pivotal value $c_{\textrm{piv}}$ is chosen to be the median concentration at a given halo mass. This implementation can be conceptualized as the top $50\%$ of halos in concentration in a mass bin takes some galaxies away from the bottom $50\%$. The magnitude of $\delta n_{\textrm{cent}}$ and $\delta n_{\textrm{sat}}$ are characterized by

[TABLE]

where $A_{\textrm{cent}}$ and $A_{\textrm{sat}}$ are the central assembly bias and satellite assembly bias parameters, respectively. The two assembly bias parameters both range between $-1$ and 1. This decorated HOD implementation also introduces modifications to the standard Halotools implementation, in that it preserves the total number of galaxies in the mock catalogs that differ only in their assembly bias parameter values. This is achieved by conditioning the decorated HOD on the total number of galaxies before populating each individual halos with galaxies. This modification shifts the number of centrals from a Bernoulli distribution but preserves the Poisson distribution in the number of satellites in each halo.

3 Simulations and Data

In this section, we introduce the galaxy projected clustering and weak lensing observables that we fit our different HOD models to. We also discuss the simulations used to predict the observables.

3.1 Galaxy clustering and g-g lensing data

The first observable we consider in this paper is the projected galaxy 2PCF, commonly referred to as $w_{p}$. It is defined as

[TABLE]

where $\xi(r_{\perp},\pi)$ is the anisotropic 2PCF, and $r_{\perp}$ and $\pi$ are transverse and line-of-sight (LOS) separations in comoving units. In this paper, we will be matching our theory $w_{p}$ to the observed $w_{p}$ by Saito et al. (2016), which was measured on the BOSS DR10 CMASS sample in the redshift range $0.43<z<0.7$. The associated covariance matrix was determined by Reid et al. (2014). The covariance matrix is computed from 5,000,000 realizations drawn from 200 bootstrap regions in the survey of roughly equal size and shape. Note that the $w_{p}$ measurement of Saito et al. (2016) assumes a different fiducial cosmology with $H_{0}=70$ km/s/Mpc and $\Omega_{\textrm{m}}=0.274$. Our models generate the predicted $w_{p}$ in Planck 2015 cosmology but we convert the prediction to the Saito et al. (2016) cosmology using a set of simple conversion formulas presented in More (2013). Our $w_{p}$ plots throughout this paper are shown assuming the Saito et al. (2016) fiducial cosmology.

The g-g lensing observable we use is the mean surface mass density contrast profile $\Delta\Sigma$, defined as

[TABLE]

where $\overline{\Sigma}(r_{\perp})$ is the azimuthally averaged and projected surface mass density at radius $r_{\perp}$ and $\overline{\Sigma}(<r_{\perp})$ is the mean projected surface mass density within radius $r_{\perp}$ (Miralda-Escude, 1991; Wilson et al., 2001; Leauthaud et al., 2017). The observed $\Delta\Sigma$ signal is presented in Leauthaud et al. (2017) for a flat $\Lambda$CDM cosmology with $H_{0}=100$ km/s/Mpc and $\Omega_{\textrm{m}}=0.31$, where the lens galaxy sample is the BOSS DR10 CMASS sample between $0.43<z<0.7$ and the background galaxy sample is a combination of two datasets: the Canada France Hawaii Telescope Lensing Survey (CFHTLenS, Heymans et al., 2012; Miller et al., 2013) and the Canada France Hawaii Telescope Stripe 82 Survey (CS82, Leauthaud et al., 2017, Erben et al. in prep). The covariance matrix, computed via bootstrap, is presented in Leauthaud et al. (2017).

3.2 Simulation and Mocks

For the purpose of this paper, we use a series of galaxy mocks generated from the AbacusCosmos N-body simulation suite, generated by the fast and high-precision Abacus N-body code (Garrison et al., 2018, 2016, Ferrer et al., in preparation; Metchnik $\&$ Pinto, in preparation). We use 20 boxes of comoving size 1100 $h^{-1}$Mpc with Planck 2015 cosmology (Planck Collaboration et al., 2016) at redshift $z=0.5$. These boxes are set to different initial phases to generate unique outputs. Each box contains 14403 dark matter particles of mass $4\times 10^{10}$ $h^{-1}M_{\odot}$. The force softening length is 0.06 $h^{-1}$ Mpc. Dark matter halos are found and characterized using the ROCKSTAR (Behroozi et al., 2013) halo finder.

To generate mock galaxies, we implement 421 generalized HODs: 1 baseline HOD (the values of the 5 standard parameters are taken from Alam et al., 2017, as shown in Table 1, with the generalized parameters set to 0) and 210 pairs of perturbed HODs. Each perturbed pair consists of two generalized HODs symmetrically perturbed around the baseline values. The perturbation is uniformly sampled from a high-dimensional box in the generalized HOD parameter space, hard bounded in each parameter direction. The 210 pairs of perturbed HOD are divided into three tiers, with increasing bounds in parameter space. Tier 1 consists of the first 15 HOD pairs and sample within the tightest bounds around the baseline HOD in parameter space; tier 2 consists of the next 75 pairs of HODs and sample from a larger box in parameter space; tier 3 consists of the final 120 pairs of HODs and sample from the largest box in parameter space. This tiered approach guarantees a higher density of samples in the region closest to the baseline value, effectively increasing the weight of the closest region in the fits. Table 2 lists the specific bounds in each parameter for each tier.

Each HOD is run over 20 simulation boxes and in each box repeated 4 times with 4 different random number generator seeds. We take the average over the 4 runs with different seeds and 20 boxes to reduce sample variance.

To generate mock projected correlation function $w_{p}$, we run Corrfunc (Sinha, 2016) with $\pi_{\textrm{max}}=77.6h^{-1}$ Mpc (this value in Planck cosmology matches $\pi_{\textrm{max}}=80h^{-1}$ Mpc used in Saito et al. (2016) in their fiducial cosmology) on the mock galaxies to obtain the $w_{p}$. We use 18 evenly spaced logarithmic bins in $r_{\perp}$ between $r_{\perp}=0.165h^{-1}$ Mpc and $r_{\perp}=29.3h^{-1}$ Mpc. These values are again chosen to match those of Saito et al. (2016). To compute the mock g-g lensing signal, we use the $\Delta\Sigma$ functionality provided in Halotools, with 10 logarithmically spaced bins between $r_{\perp}=0.157h^{-1}$ Mpc and $r_{\perp}=15h^{-1}$ Mpc, matched with those used Leauthaud et al. (2017). We do not go below $0.157h^{-1}$ Mpc due to limited force softening resolution in our simulations.

Before we move on to the emulator fits, we perform a rejection step for our samples of paired HODs. Specifically, having computed the $w_{p}$ of each of the 421 HODs, we reject the ones whose $w_{p}$ deviate too far from the baseline $w_{p}$. The reason for this rejection step is that there are directions in the HOD space that drastically change the $w_{p}$, way beyond the capabilities of our emulator model (described in the following section). Even within our HOD parameter sample, we find examples where a combination of high $M_{\mathrm{cut}}$ and low $M_{1}$ leads to a very high satellite fraction, resulting in a $w_{p}$ that is 10 times higher in amplitude than our baseline $w_{p}$. To constrain ourselves to regions of the HOD space with a moderately behaved $w_{p}$, we choose an arbitrary but generous cut of

[TABLE]

where $w_{p}(r_{\perp}<1h^{-1}\textrm{Mpc})$ is the $w_{p}$ integrated up to $r_{\perp}=1h^{-1}\textrm{Mpc}$, and $w_{p,0}$ denotes the baseline $w_{p}$. This cut roughly corresponds to a percentage range of $63\%-158\%$ of the baseline $w_{p}$, more than sufficient to cover the a few percent level uncertainties in the observed $w_{p}$. By applying this cut, we remove 155 HODs from our sample of mocks, with the remaining 266 HODs moving on to the emulator fits.

To generate mocks for the decorated HOD emulator, we follow the same procedure. We first generate 1 baseline HOD (using the HOD parameters taken from Alam et al. (2017) with the assembly bias parameters set to 0), and 240 pairs of perturbed HODs around the baseline. The perturbed HODs are divided into four tiers, with each HOD applied to all 20 simulation boxes with 4 random number seeds. We summarize the parameter ranges in Table 3.

Then we perform the same selection step to select only the mocks whose $w_{p}$ is within a reasonable range of the baseline $w_{p}$. We employ the following cut

[TABLE]

and we select our final training set of 251 HODs.

4 Methods

In this section, we present our methodology for constructing the $w_{p}$ and $\Delta\Sigma$ emulators from simulations. Then we discuss the use of nested sampling to explore the generalized HOD parameter space when fitting the observed $w_{p}$ and $\Delta\Sigma$.

4.1 The $w_{p}$ and $\Delta\Sigma$ emulator

We first construct an emulator that models the galaxy projected 2PCF $w_{p}$ as a function of the generalized HOD parameters ($\log_{\textrm{10}}M_{\textrm{cut}},\log_{\textrm{10}}M_{1},\sigma,\alpha,s,s_{v},A$). We model the weighted projected 2PCF $r_{\perp}w_{p}$ because it has more moderate behavior, resulting in a more balanced covariance matrix.

There is a variety of ways to construct the model such as neural nets and gaussian processes. Here we are using a model based on first and second derivatives. Specifically, we can write down our model analytically as

[TABLE]

where $i$ is the bin number and $\mathbf{p}=[p_{1},p_{2},...,p_{7}]$ are the 7 HOD parameters. Thus, for each bin, we have a total of 36 coefficients to fit for: the intercept $r_{\perp}w_{p,i}(\mathbf{p}_{0})$, the 7 first derivatives $\partial r_{\perp}w_{p,i}/\partial p_{j}$, and the 28 second derivatives $\partial^{2}r_{\perp}w_{p,i}/\partial p_{j}\partial p_{k}$.

Similarly, we can construct our g-g lensing model with the following formula

[TABLE]

where we have again modulated the lensing signal $\Delta\Sigma$ with $r_{\perp}$ to produce a more flat behavior as a function of scale and a more balanced covariance matrix.

We fit the models given in Equation 11 and Equation 12 to the mock observables presented in Section 3.2 using a standard least-squares routine. The two emulators are fitted individually and each $r_{\perp}$ bin is also fitted independently. For each emulator, each bin in $r_{\perp}$ has 36 unknowns and 266 data points in the fit.

We run a set of cross validation tests to examine the performance of our fit. We perform 10 cross validation runs. For each run, we train our model on a randomly selected $\sim 90\%$ of HODs plus the baseline HOD and then test on the remaining $\sim 10\%$. This way, we can compute the in-sample and out-sample errors of our fit. The in-sample errors for each validation fun are computed by testing our best-fit emulator for that run against the set of generalized HODs that the model is trained on, while the out-sample errors are computed by testing against the set of generalized HODs that the model is not trained on.

We show the average of the absolute values of the in-sample errors and out-sample errors of our best-fit emulators as a function of $r_{\perp}$ in Figure 1 and Figure 2, respectively. The average is calculated by averaging the fractional error across all generalized HODs in the test set, then averaged over all 10 cross-validation runs. The errors shown are in fraction of $w_{p}$ itself.

To reduce overfitting in the $w_{p}$ emulator, we calculate the mean and standard deviation of the best-ft of each of the 36 coefficients by combining the best fit values from the 10 cross-validation runs. We set to zero the 2 coefficients that have the lowest signal noise. These 2 coefficients correspond to the following second derivatives: $\partial s_{v}\partial A$ and $\partial\log_{10}M_{1}\partial s_{v}$. Dropping these 3 second derivatives gives us a slight reduction in out-sample error by $\approx 2\%$. The errors shown in Figure 1 are the errors after removing these 2 second derivatives from the emulator. We do not repeat this procedure for the g-g lensing emulator as it has low out-sample errors and does not appear to be overfit.

We repeat the same procedure for the decorated HOD to construct emulators of $w_{p}$ and $\Delta\Sigma$ as a function of 6 decorated HOD parameters ($\log_{\textrm{10}}M_{\textrm{cut}},\log_{\textrm{10}}M_{1},\sigma,\alpha,A_{\textrm{cent}},A_{\textrm{sat}}$). We run cross-validation tests and remove second derivatives $\partial\log_{\textrm{10}}M_{\textrm{cut}}\partial A_{\textrm{sat}}$ from the $w_{p}$ emulator due to their low signal-to-noise. We recover a relative out-sample error of $<2.1\%$ in the $w_{p}$ emulator and $<0.5\%$ in the $\Delta\Sigma$ emulator.

4.2 The likelihood functions

The joint log-likelihood function is computed as

[TABLE]

where the likelihood functions for $w_{p}$ and $\Delta\Sigma$ are assumed to be Gaussian. Thus, we have

[TABLE]

where $\mathbf{C}$ is the observed covariance matrix of $w_{p}$, and $\delta\mathbf{w}_{p}$ is the difference vector between the observed $w_{p}$ and the emulated $w_{p}$ given the HOD parameters,

[TABLE]

We construct the log-likelihood of $\Delta\Sigma$ in analogous manner.

4.3 Nested Sampling

We need a sampling algorithm to explore the generalized HOD posterior space, which we compute from the likelihood function and the priors. In this paper, we would also like to compare different generalized HOD models using Bayesian evidence. The Bayesian evidence is defined as

[TABLE]

where $M$ represents the model, $D$ represents the data, and $\Theta$ represents the model parameters. $P(D|\Theta,M)$ is the likelihood of the data given the parameters of our model, and $P(\Theta|M)$ is the prior for the parameters of our model. The evidence can simply be interpreted as the marginal likelihood of the data given the model, and serves as an important metric in Bayesian model comparisons.

The nested sampling technique, first developed by Skilling (2006) gives us an effective way to compute Bayesian evidence integrals. For this paper, we use the publicly available nested sampling code dynesty (Speagle & Barbary, 2018; Speagle, 2019). This code computes the Bayesian evidence while generating samples of the posterior parameter space. It replaces the multi-dimensional evidence integral over model parameters in Equation 16 with a 1D integral over the prior mass contained within nested isolikelihood contours. For a detailed description of the code, refer to Speagle (2019).

In our dynesty runs, we use a nested sampler with 1500 live points and a random walk sampler conditioned on the bounding distribution. The stopping criterion is set to $d\log\mathcal{Z}>0.01$.

5 Results

In this section, we use the $w_{p}$ and $\Delta\Sigma$ emulators to simultaneously fit the observed $w_{p}$ and $\Delta\Sigma$ signals and to explore the generalized/decorated HOD parameter space, using the likelihood function described in Section 4.2 and the nested sampling technique described in Section 4.3.

In the following subsections, we first showcase the marginalized parameter constraints for the generalized HOD model, then the marginalized parameter constraints for the decorated HOD model. We compare the two models and discuss whether either of them provides a good simultaneous fit of the observed projected clustering and g-g lensing signal.

5.1 Joint fits with generalized HOD

We first present the joint fit using the generalized HOD with the 6 parameters listed in the first column of Table 4. We impose broad Gaussian priors centered around their baseline values, as listed in the second and third column of Table 4. The baseline values of all generalized parameters are set to 0. The width of the Gaussian priors are chosen to be broad and non-informative, but constraining enough to mostly fit inside the emulators’ training range and reject unphysical values.

Figure 3 shows the 1D and 2D marginalized posterior constraints on the generalized HOD parameters. The black lines show the baseline values as reference. The blue contours show the 0.5, 1, 1.5, and 2 $\sigma$ uncertainties. The values displayed above the 1D marginals are posterior medians with the upper/lower bounds associated with the 0.025 and 0.975 quantiles. Compared to the baseline values in Table 1, the posterior shows a preference for a lower $M_{\textrm{cut}}$, $\sigma$, and $A$, while the other parameters remain consistent with the prior mean.

We point out that, even though the posterior marginals look Gaussian, the posterior medians shown on the 1D marginals are different from the maximum a posteriori (MAP) point due to the complex shape of the posterior surface. We illustrate this point by showing an example trace plot for parameter $s$ in the top panel of Figure 7, where we plot the posterior samples as a function of log prior volume, colored by their importance weight. The nested sampler marches from the left to the right and goes from being prior dominated to being likelihood dominated. The shape of the posterior marginal is determined by the most heavily weighted samples, colored in yellow, whereas the samples diverge into multiple local maxima towards the right of the plot. The MAP point corresponds to the long spike at $s\approx 0.15$, whereas the posterior median corresponds to the center of the yellow region, which is close to 0.

For completeness, the MAP point is given by $\log_{10}(M_{\textrm{cut}}/h^{-1}M_{\odot})=13.27$, $\log_{10}(M_{1}/h^{-1}M_{\odot})=14.34$, $\sigma=0.634$, $\alpha=1.144$, $s=-0.065$, and $A=-0.484$. Compared to the posterior median point shown on the 1D marginals, we see that the data actually prefer an even more extreme $A$. An assembly bias parameter of -0.5 is beyond our training range and also physically unlikely considering it would mean the shuffling of satellites across 1.0 dex in mass range, which is comparable to the entire halo mass range in the simulations. This fit suggests that the generalized HOD model failed to find a consensus between the 2PCF and g-g lensing signal within a reasonable parameter range.

The Bayesian evidence for this generalized HOD model is $\log\mathcal{Z}=-65.09\pm 0.08$. The Bayesian evidence for the standard HOD model without any generalizations is $\log\mathcal{Z}=-69.46\pm 0.07$. Thus, the generalized model is preferred over the standard model, with a Bayesian evidence roughly 4.4 e-fold higher than that of the standard model. Thus, the observed $w_{p}$ and $\Delta\Sigma$ somewhat favor the generalized HOD model despite the introduction of 2 new parameters, which increases the prior volume by order of a few. We summarize the evidence values for all the different models in Table 6 for comparison.

We show the maximum a posteriori (MAP) $w_{p}$ and $\Delta\Sigma$ in Figure 4. With nested sampling, the MAP is simply the last point in the chain. Again note that the MAP point is not the posterior median shown on the 1D marginal posterior plots. We then pass the MAP values back to the emulators to generate the emulated $w_{p}$ and $\Delta\Sigma$, which we plot in green. The red curves show the actual $w_{p}$ and $\Delta\Sigma$ computed from the mocks generated from the MAP HOD. We show the observed signals plotted in blue for comparison. It is clear that this generalized HOD model does not provide a good fit for either $w_{p}$ or $\Delta\Sigma$. Compare with Leauthaud et al. (2017), which shows that the HOD that fits $w_{p}$ predicts a g-g lensing signal $20-40\%$ higher than the observed g-g lensing signal. This fit shows that we do not find a point in the generalized HOD space that fits both observables, but instead the model is trying to find a comprise between the observed $w_{p}$ and g-g lensing, by bringing down the predicted g-g lensing while losing consistency with the observed $w_{p}$.

Comparing the emulated and mock observables in Figure 4, we see that our $\Delta\Sigma$ emulator successfully predicts the mock $\Delta\Sigma$ with remarkable accuracy. The $w_{p}$ emulator performs fairly well at $r_{\perp}>1h^{-1}$Mpc, but less so at smaller scales. This is due to the fact that the MAP point is somewhat outside of our training range due to an extreme best-fit assembly bias value.

We also repeat this joint fit with prior widths on all the non-mass parameters extended from 0.1 to 0.2, allowing us to explore regions of the HOD parameter space beyond our training. We recover the same qualitative result. We do not find a consensus fit between the two observables even though the best fit HOD now gives more extreme values, with $\sigma=0.369$ and $A=-0.848$. The best fit still significantly under-predicts $w_{p}$ while significantly over-predicting g-g lensing, similar to Figure 4.

We then test the 7-parameter generalized HOD model with the addition of satellite velocity bias $s_{v}$. The 7 parameters are $\log_{\textrm{10}}M_{\textrm{cut}},\log_{\textrm{10}}M_{1},\sigma,\alpha,s,s_{v}$, and $A$. We choose a broad Gaussian prior for $s_{v}$, with $\mu_{\textrm{prior}}=0$ and $\sigma_{\textrm{prior}}=0.1$. It is possible that the addition of $s_{v}$ can improve the model, considering that the generalized HOD implementation of satellite velocity bias is based on particle selection so it does affect the radial distribution of the satellites in the halo. Figure 4 of Yuan et al. (2018) shows that $s_{v}=\pm 0.2$ changes $w_{p}$ by approximately $1\%$.

We fit the 7-parameter generalized HOD and sample the posterior with dynesty. We find the same modes in the posterior space as the 6-parameter generalized HOD without $s_{v}$. The Bayesian evidence of this model is $\log\mathcal{Z}=-65.18\pm 0.08$, which is statistically the same as that of the 6-parameter model despite the addition of a new parameter. Thus, the 7-parameter model with $s_{v}$ is not preferred over the 6-parameter model without $s_{v}$. We omit the $s_{v}$ parameter in the following models we consider.

5.2 Generalized HOD with an amplitude scaling parameter

In this section, we introduce a floating amplitude parameter “scaling" to the lensing signal. This is a free parameter multiplied onto the emulated $\Delta\Sigma$ to allow its overall amplitude to shift up and down. We set a log-normal prior on “scaling” with mode equal to 1 and a scale of 0.2.

Figure 5 shows the 1D and 2D marginalized posterior constraints on the 6-parameter generalized HOD, but now with the “scaling" parameter. The posterior median is less than 1$\sigma$ deviation in most HOD parameters from the prior mean, except for $\log_{10}M_{1}$, which gives a $\sim 2\sigma$ deviation from prior mean. The MAP HOD is $\log_{10}(M_{\textrm{cut}}/h^{-1}M_{\odot})=13.36$, $\log_{10}(M_{1}/h^{-1}M_{\odot})=14.48$, $\sigma=0.590$, $\alpha=1.102$, $s=-0.221$, and $A=-0.098$, again somewhat different from the posterior medians shown on the 1D marginals. However, unlike the “unscaled" fit, all the best-fit parameters of this fit fall inside the training range of our emulator, with moderate assembly bias and satellite distribution parameter values. The only parameter whose best-fit value is significantly different from the prior mean is $M_{1}$, approximately 0.28 dec higher than the baseline value, which calls for fewer satellites given a halo mass but this value is not beyond reason. The “scaling" factor strongly favors a value of $\sim 0.67$ instead of 1. The Bayesian evidence of this model is $\log\mathcal{Z}=16.00\pm 0.08$, a 49.1 e-fold increase over the 6-parameter model with no floating amplitude. Thus, the inclusion of a floating amplitude parameter for the g-g lensing signal is strongly favored and we now find a reasonable consensus HOD between the two observables.

We show the trace plot for parameter $s$ in Figure 7, where the top and bottom panels show the trace plots for the fit without and with the scaling parameter, respectively. The nested sampler marches from the left to the right through the prior volume. The points denote the dead points over the course of the run, colored by their importance weights. For technical details, refer to the dynesty code release paper (Speagle, 2019). We see that for both runs, the sampling is initially prior dominated but becomes more and more likelihood dominated as we move to the right of the panel. For the “unscaled" case, the sample diverges and several local modes develop, suggesting disagreement in the dataset. However, for the “scaling" case, the sample remains compact throughout the prior volume, giving rise to one maximum likelihood mode, suggesting the discovery of a consistent fit.

Figure 6 shows the emulated and mock $w_{p}$ and $\Delta\Sigma$ of the MAP estimate of the 6-parameter + scaling model posterior. The emulated $\Delta\Sigma$ shown is after applying the best-fit scaling parameter, thus appearing in excellent agreement with the observed $\Delta\Sigma$. We see that both the $w_{p}$ and $\Delta\Sigma$ predictions are consistent with the observations, deviating by $<1\sigma$ in almost all bins. Unlike for the unscaled fit, the emulated observables and the mock observables are also in excellent agreement on all scales, showcasing the high accuracy of our emulators within its training range. The fact that the introduction of a flexible amplitude on $\Delta\Sigma$ produces a remarkably good fit with reasonable best-fit HOD parameters suggests that the inconsistency between the observed $w_{p}$ and $\Delta\Sigma$ measurements is well described by a scale independent amplitude shift. The best-fit value of 0.67 for the scaling parameter corresponds to a $33\%$ deficiency in bias, which is consistent with the $20-40\%$ disagreement between the observed and predicted $\Delta\Sigma$ found in Leauthaud et al. (2017). This suggests that assembly bias alone cannot reconcile the observed $w_{p}$ and $\Delta\Sigma$ measurements.

5.3 Joint fits with decorated HOD

The decorated HOD implementation provided in Halotools provides another framework to extend the standard 5-parameter HOD model. The decorated HOD presents several key differences compared to the generalized HOD. The decorated parameters include two separate assembly bias parameters for the central and satellite galaxies. A key difference in implementation is that it distributes satellite galaxies on NFW profiles given halo concentration whereas the generalized HOD distributes the satellites on halo particles, leading to differing observable predictions on the small scale. While we showed that the generalized HOD in its current form does not reconcile the observed projected clustering and g-g lensing signal, it is possible that the decorated HOD model provides a better joint fit of the two.

The construction of the decorated HOD emulator was discussed in Section 4.1. Then we use dynesty to perform a joint fit of the observed $w_{p}$ and g-g lensing following the same procedure as for the generalized HOD framework. The 6 HOD parameters in this joint fit are $\log_{\textrm{10}}M_{\textrm{cut}},\log_{\textrm{10}}M_{1},\sigma,\alpha,A_{\textrm{cent}}$, and $A_{\textrm{sat}}$. We again choose broad non-informative Gaussian priors for these parameters centered around their baseline values. However, we truncate the Gaussians at -1 and 1 for the two assembly bias parameters, as they are only defined between -1 and 1. We summarize these priors in Table 5.

The MAP HOD values are $\log_{10}(M_{\textrm{cut}}/h^{-1}M_{\odot})=13.01$, $\log_{10}(M_{1}/h^{-1}M_{\odot})=14.10$, $\sigma=0.562$, $\alpha=1.299$, $A_{\textrm{cent}}=-0.168$, and $A_{\textrm{sat}}=-0.387$. We see that the data prefer a very high satellite assembly bias value of almost $-0.4$, which is likely unphysical as it suggests the less concentrated half of halos have more than twice as many galaxies as the more concentrated half. Figure 8 shows the emulated $w_{p}$ and g-g lensing of the MAP HOD using the 6 parameter decorated HOD framework. We do not find a good fit, and we see behavior very similar to Figure 4 in that the joint fit attempts to bring the emulated lensing signal in better agreement with the observations by shifting the emulated $w_{p}$ out of agreement with the observations. The trace plot for this fit shows strong multi-modality, again suggesting that the model failed to find a consensus between the two observables within a reasonable parameter range. We do not show it for brevity.

The evidence of this model is given by $\log\mathcal{Z}=-72.89\pm 0.08$, compared to the evidence of the standard 5-parameter model without the 2 assembly bias parameters $\log\mathcal{Z}=-74.91\pm 0.07$. This suggests that the decorated HOD model with assembly bias is moderately favored over the standard 5-parameter HOD. The evidence is also 7.8 e-fold lower than that of the generalized HOD model, possibly due to the inclusion of satellite distribution parameter $s$ and differences in the implementation. Again we have compiled all the evidence values in Table 6 for comparison.

Following the same procedure as for the generalized HOD, we then add the scaling parameter to the decorated HOD model to allow the overall amplitude of the lensing signal to shift. We re-fit the observed $w_{p}$ and $\Delta\Sigma$ and find a set of best-fit values that seem more reasonable and give emulated signals that closely match the observables. The best-fit HOD parameters are $\log_{10}(M_{\textrm{cut}}/h^{-1}M_{\odot})=13.19$, $\log_{10}(M_{1}/h^{-1}M_{\odot})=14.28$, $\sigma=0.601$, $\alpha=1.370$, $A_{\textrm{cent}}=0.136$, and $A_{\textrm{sat}}=-0.079$. The best-fit value for the scaling parameter is $0.66^{+0.03}_{-0.03}$, which represents a $34\%$ inconsistency between the observed $w_{p}$ signal and the observed $\Delta\Sigma$ signal, consistent with the discrepancy reported in Leauthaud et al. (2017). We show the best-fit emulated $w_{p}$ and g-g lensing signal in Figure 8, where we see excellent agreement between the emulated signals and the observed signals. The integrated Bayesian bias for this model is given by $\log\mathcal{Z}=-17.92\pm 0.08$, which is a 55 e-fold increase compared to the decorated model without scaling. This again shows that the inclusion of a flexible lensing amplitude is strongly favored, and that the decorated HOD model by itself does not provide a good fit of the projected clustering and g-g lensing observables. We recover the same qualitative results when we increase the priors on the non-mass parameters ($\sigma$, $\alpha$, $A_{\textrm{cent}}$, and $A_{\textrm{sat}}$) from 0.1 to 0.2, though the joint fit in the unscaled case returns more extreme best-fit values in $\sigma$ and the assembly biases.

It is worth noting that Figure 8 shows great consistency between the emulated $w_{p}$ in green and the mock $w_{p}$ in red, except at smaller scales for the unscaled case, where we get approximately $2\%$ error. We achieve such good accuracy due to our large training range and restricting the joint fit to the training range. When we increase the HOD parameter priors and let the joint fit explore regions that the emulators are not trained in, the accuracy quickly collapses, especially in $w_{p}$ due to its highly non-linear dependence on HOD parameters. We see a minor manifestation of this in Figure 4, where the best-fit assembly bias value is moderately beyond the training range.

It is also interesting that while the two HOD frameworks produce the same qualitative behavior in both the scaled and unscaled case and noticeably yielding the same best-fit scaling parameter, the best-fit decorated HOD and the best-fit generalized HOD give rather different HOD parameter values. This is possibly due to differences in the base HOD implementation, the assembly bias implementation, and the inclusion of the satellite distribution parameter $s$ in the generalized HOD. This also showcases the effect of differences in HOD implementations and highlights the importance of comparing different HOD models in future cosmological inference work.

6 Discussion

In this paper, we have tested whether two different extended HOD models can fit the observed projected clustering signal $w_{p}$ and the g-g lensing signal $\Delta\Sigma$. We find that neither model provides a good joint fit, with Bayesian evidence more than 49 e-fold lower than that of the models with an extra scaling parameter. The best-fit scaling parameter is consistently showing a roughly $34\%$ discrepancy between the observed $w_{p}$ and $\Delta\Sigma$, in agreement with the conclusions of Leauthaud et al. (2017).

Table 6 summarizes the integrated evidence of the models we have tested in this paper. We see that in both the generalized HOD framework and the decorated HOD framework, the addition of generalized parameters — $s,A$ for the generalized HOD and $A_{\textrm{cent}},A_{\textrm{sat}}$ for the decorated HOD — is favored, with a 4.3 and 2.0 e-fold increase to the model evidence, respectively. The addition of satellite velocity bias parameter is dis-favored, however, with no significant change to the model evidence. The decorated HOD models also show significantly lower evidence than the generalized HOD models. This is due to different implementations of the standard HOD, different implementations of assembly bias, and the inclusion of a satellite distribution parameter $s$ in the generalized HOD that introduces additional flexibility.

The fact that the joint fit with generalized/decorated parameters give complex and multi-modal posteriors (see Figure 3) are most likely a result of the clustering and lensing data sets being incompatible, but it may also suggest that our HOD parameter space has moderate degeneracies with respect to the projected clustering and g-g lensing data sets. It is possible that if we fit to the anisotropic correlation function $\xi(r_{\perp},\pi)$ instead of the $w_{p}$, we would recover more constraining power on the parameter posteriors.

Our conclusions are limited in several ways. First of all, we only considered galaxy assembly bias implementation that uses halo concentration as the secondary dependency. In the more general sense of the term, galaxy assembly bias is not just limited to halo concentration as the only secondary dependence. It is possible that a galaxy assembly bias implementation whose secondary dependence is linked to the merger history, or local environment, or other halo properties, may reconcile the observed discrepancy. Lange et al. (2019) additionally tested the decorated HOD model with halo spin as the secondary dependence and found a decrease of at most $10\%$ to the predicted lensing signal on the small scale and almost no impact on the large scale, insufficient to reconcile the $20-40\%$ discrepancy.

We have encountered significant challenges when constructing our $w_{p}$ emulator over a wide range of HOD parameters. We find that $w_{p}$ has a highly non-linear dependence on the HOD parameters, making the second-order Taylor expansion model insufficient in modeling $w_{p}$ over the full range of HOD parameters. Thus, we implement the rejection sampling step to only emulate in regions of HOD space that are well-behaved. Even then, we still have difficulties in the accuracy of our decorated HOD emulator in certain regions of parameter space that are of potential interest. In this work, we have limited ourselves to exploring a large but not exhaustive range of HOD parameters. In the future, we intend to further expand our training set and adopt a neural network style emulator that can more easily incorporate these non-linearities. (Refer to Wibking et al. (2019) for a non-linear emulator implementation using Gaussian processes.) We also plan on changing to a different HOD parameter basis, possibly using parameters such as total galaxy number density and satellite fraction, to reduce degeneracies in the parameters.

Another limitation of our analysis is that we do not marginalize over Planck cosmological posteriors. We used a fixed cosmology in building our emulators and our fits. We believe the cosmological dependence of this problem is an interesting one, and in future work we will use the AbacusCosmos simulation boxes that sample the Planck posterior to extend our emulator to also emulating the cosmological parameters. We defer that discussion to a future paper.

Baryonic effects provide another possible explanation of the discrepancy between the projected clustering measurements and the g-g lensing measurements. We have used dark matter only N-body simulations for our model predictions. While the introduction of the satellite distribution parameter may have made our generalized HOD models flexible enough to marginalize over baryonic effects in the projected clustering predictions, it does not change the distribution of dark matter particles themselves and thus fails to account for any baryonic effects on the g-g lensing signal. Several studies have suggested that baryonic effect can impact the halo profile and influence subhalo properties (e.g. van Daalen et al., 2014; Velliscig et al., 2014; Chaves-Montero et al., 2016). Leauthaud et al. (2017) compared the g-g lensing signals from the full-physics Illustris simulations (Vogelsberger et al., 2014a, b; Genel et al., 2014; Sijacki et al., 2015; Nelson et al., 2015) and from the corresponding gravity-only simulations to estimate the effect of baryonic physics on CMASS-like samples, and found that baryonic effects can induce a noticeable increase in $\Delta\Sigma$, in the direction of reconciling the $20-40\%$ discrepancy. Refer to Figure 12 of Leauthaud et al. (2017), we see that baryonic physics can induce an increase of $\sim 20\%$ in $\Delta\Sigma$ at a radius of 0.2$h^{-1}$Mpc, and an essentially zero increase at larger radius of approximately $10h^{-1}$Mpc. We find in our analyses that the $33-34\%$ discrepancy in the lensing signal exists on all scales, so baryonic effects, at least as detected using the Illustris simulations, might not fully explain the lensing discrepancy. Moreover, Weinberger et al. (2017) found a smaller impact of the baryonic effects using an improved AGN feedback model with the Illustris simulations. Lange et al. (2019) used the newly-released and improved IllustrisTNG simulation (Pillepich et al., 2018; Springel et al., 2018; Nelson et al., 2019; Naiman et al., 2018; Marinacci et al., 2018; Nelson et al., 2018) and found an at most $10\%$ decrease to the lensing signal on the small scale due to baryonic effects, even in the $95\%$ posterior range. However, beyond just affecting the halo density profile, baryonic effects can also produce more complicated assembly bias beyond that using the halo concentration proxy.

7 Summary and Conclusions

In this paper, we address the question of whether galaxy assembly bias, defined as the secondary dependence of galaxy occupation on halo concentration, can explain the $20-40\%$ discrepancy between the galaxy projected clustering measurement $w_{p}$ and the g-g lensing measurement $\Delta\Sigma$ in a Planck cosmology. We apply two different extensions to the standard 5-parameter HOD of the Zheng & Weinberg (2007) form to fit the clustering and lensing observables. The first is a generalized HOD which adds the satellite distribution parameter and the assembly bias parameter. The second is the decorated HOD, which adds the central assembly bias parameter and the satellite assembly bias parameter. We find that neither model can reconcile the discrepancy between the two observables. However, allowing a $\sim 34\%$ increase to the measured lensing amplitude would yield a very good fit using either models. This result suggests that galaxy assembly bias, in the ways we have implemented it, does not reconcile the $20-40\%$ discrepancy between the galaxy projected clustering measurement and the g-g lensing measurement in Planck cosmology. While other implementations of galaxy assembly bias might explain the discrepancy, our findings suggest that galaxy assembly bias is in-plausible as the main explanation of the discrepancy. It is also possible that a slightly different cosmology or baryonic effects can explain the discrepancy, but these discussions are beyond the scope of this paper.

Acknowledgements

We thank our referee Andrew Hearin for helpful feedback on this paper. We also thank him and Kuan Wang for help with Halotools. We would also like to thank Josh Speagle and Lehman Garrison for their support with dynesty and AbacusCosmos throughout this paper. We are grateful to Shun Saito for allowing us to use their projected clustering measurement data.

DJE is supported by U.S. Department of Energy grant DE-SC0013718 and as a Simons Foundation Investigator. AL is supported by the U.D Department of Energy, Office of Science, Office of High Energy Physics under Award Number DE-SC0019301, by the Packard Foundation, and by the Sloan Foundation.

The ABACUS simulations used in this paper are available at https://lgarrison.github.io/AbacusCosmos.