Cosmology with galaxy-galaxy lensing on non-perturbative scales: Emulation method and application to BOSS LOWZ

TL;DR

This paper introduces a nonlinear emulation framework combining HOD models and N-body simulations to accurately predict galaxy clustering and lensing signals, applied to BOSS LOWZ data to measure cosmological parameters.

Contribution

The paper presents a novel emulation method that achieves high accuracy in modeling galaxy signals on non-perturbative scales, enabling improved cosmological parameter estimation.

Findings

Achieved <2% interpolation accuracy across parameter space.

Recovered true S_8 value from mock data.

Measured S_8=0.847±0.037 from BOSS LOWZ data, in tension with Planck results.

Abstract

We describe our nonlinear emulation (i.e., interpolation) framework that combines the halo occupation distribution (HOD) galaxy bias model with $N$-body simulations of nonlinear structure formation, designed to accurately predict the projected clustering and galaxy-galaxy lensing signals from luminous red galaxies (LRGs) in the redshift range $0.16 < z < 0.36$ on comoving scales $0.6 < r_p < 30$ \hMpc. The interpolation accuracy is $\lesssim 1-2$ per cent across the entire physically plausible range of parameters for all scales considered. We correctly recover the true value of the cosmological parameter $S_8 = ({\sigma_8}/{0.8228}) ({\Omega_{\text{m}}}/{0.3107})^{0.6}$ from mock measurements produced via subhalo abundance matching (SHAM)-based lightcones designed to approximately match the properties of the SDSS LOWZ galaxy sample. Applying our model to Baryon Oscillation Spectroscopic Survey (BOSS) Data Release 14 (DR14) LOWZ galaxy clustering and galaxy-shear cross-correlation measurements made with Sloan Digital Sky Survey (SDSS) Data Release 8 (DR8) imaging, we perform a prototype cosmological analysis marginalizing over $w$CDM cosmological parameters and galaxy HOD parameters. We obtain a 4.4 per cent measurement of $S_8 = 0.847 \pm 0.037$, in $3.5\sigma$ tension with the Planck cosmological results of $1.00 \pm 0.02$. We discuss the possibility of underestimated systematic uncertainties or astrophysical effects that could explain this discrepancy.

Tables

Table 1. Table 1: Emulator parameters. Each non-cosmological parameter sample is assigned to a cosmological parameter sample, with 10 non-cosmological samples assigned to a given cosmology. The emulator training set therefore has a total sample size N = 400 𝑁 400 N=400 .

Parameter	Sampling range	Fiducial value	Units
$n_{\text{gal}}$	[$2.5$, $3.5$] $\times {10}^{-4}$	$3\times {10}^{-4}$	${\left(h^{-1}\text{Mpc}\right)}^{-3}$
${\sigma }_{\log M}$	[0.01, 0.8]	0.4	dimensionless
$M_0/M_1$	[0.0, 0.4]	0.1	dimensionless
$M_1/M_{\text{min}}$	[7.5, 20]	8.5	dimensionless
$\alpha$	[0.6, 1.5]	1.0	dimensionless
$A_{\text{conc}}$	[0.5, 3.0]	1.0	dimensionless
$R_{\text{rescale}}$	[0.5, 2.0]	1.0	dimensionless
${\sigma }_8$	[0.65, 1.0]	0.830	dimensionless
${\Omega }_{\text{CDM}}{h}^2$	[0.1045, 0.1322]	0.1199	dimensionless
${\Omega }_b{h}^2$	[0.0209, 0.0235]	0.02222	dimensionless
$H_0$	[61.567, 74.793]	67.26	km s-1 Mpc-1
$n_s$	[0.9300, 0.9898]	0.9652	dimensionless
$w_0$	[-1.370, -0.655]	-1.0	dimensionless
$N_{\text{eff}}$	—	3.046	dimensionless
$z_{\text{eff}}$	—	0.300	dimensionless

Table 2. Table 2: Priors on parameters for our cosmological analyses. Units are given in Table 1 .

Parameter	Prior range
$n_{\text{gal}}$	[$2.5$, $3.5$] $\times {10}^{-4}$
${\sigma }_{\log M}$	[0.01, 0.8]
$M_0/M_1$	[0.0, 0.4]
$M_1/M_{\text{min}}$	[7.5, 20]
$\alpha$	[0.5, 1.5]
$A_{\text{conc}}$	[0.5, 3.0]
$R_{\text{rescale}}$	[0.5, 2.0]
${\sigma }_8$	[0.65, 1.0]
${\Omega }_m$	[0.26, 0.35]
${\Omega }_b$	[0.0394, 0.0602]
$H_0$	[61.57, 74.80]
$n_s$	[0.93, 0.9898]
$w_0$	[-1.37035, -0.6548324]
$A_{\text{lensing}}$	$𝒩\left(\mu =1.0,\sigma =0.06\right)$

Table 3. Table 3: Posterior means ( ± 68 plus-or-minus 68 \pm 68 per cent credible intervals) of the parameters resulting from our cosmological analyses on both lightcones (indicated with non-boldface column names) and data (indicated with boldface column names). Units are given in Table 1 . Upper or lower limits indicate that the 68 per cent credible interval abuts the prior minimum or maximum value of a parameter. For A conc subscript 𝐴 conc A_{\text{conc}} in the ‘baryons’ column, the 68 per cent posterior interval is identical to the prior, so no values are given in the table for this entry. For entries marked ‘N/A’, the analysis did not include A lensing subscript 𝐴 lensing A_{\text{lensing}} as a free parameter. Due to the asymmetry of the marginalized posteriors, the posterior means are not necessarily identical to the posterior maxima shown in Figures 8 and 9 .

Parameter	Fiducial	Lower $f_{\text{sat}}$	Higher $f_{\text{sat}}$	15% Incompl.	‘Baryons’	LOWZ	$\left[>\mathrm{𝟐}{𝒉}^{-\mathrm{𝟏}}\text{Mpc}\right]$
$n_g\times {10}^4$	${2.90}_{-0.29}^{+0.19}$	$3.00\pm 0.23$	${3.13}_{-0.14}^{+0.32}$	${3.04}_{-0.22}^{+0.34}$	$2.98\pm 0.23$	${3.08}_{-0.16}^{+0.38}$	${2.97}_{-0.31}^{+0.25}$
$M_0/M_1$		${0.20}_{-0.13}^{+0.11}$		${0.178}_{-0.16}^{+0.073}$			$0.20\pm 0.10$
$M_1/M_{\text{min}}$	${12.8}_{-3.9}^{+2.0}$	${11.9}_{-2.8}^{+1.3}$	${12.0}_{-3.2}^{+1.6}$	${11.08}_{-3.3}^{+0.73}$	${13.0}_{-3.9}^{+2.1}$		${13.5}_{-3.8}^{+2.9}$
${\sigma }_{\log M}$	${0.165}_{-0.15}^{+0.049}$		${0.174}_{-0.16}^{+0.053}$	${0.207}_{-0.14}^{+0.099}$		${0.291}_{-0.082}^{+0.11}$	${0.186}_{-0.13}^{+0.081}$
$\alpha$	${0.90}_{-0.27}^{+0.19}$		${0.82}_{-0.25}^{+0.14}$	${0.85}_{-0.30}^{+0.13}$	${0.89}_{-0.27}^{+0.17}$		${1.04}_{-0.24}^{+0.34}$
$A_{\text{conc}}$	$>1.47$	$1.74\pm 0.65$	${1.71}_{-0.83}^{+0.69}$	${1.79}_{-0.68}^{+0.82}$	$---$	${1.78}_{-0.69}^{+0.82}$	${1.70}_{-0.78}^{+0.65}$
$R_{\text{rescale}}$	$0.91\pm 0.14$	${0.82}_{-0.12}^{+0.13}$	$1.03\pm 0.16$	${0.87}_{-0.15}^{+0.17}$	${0.92}_{-0.13}^{+0.16}$	${1.17}_{-0.095}^{+0.17}$	$>1.18$
$S_8$	${1.019}_{-0.055}^{+0.040}$	$1.017\pm 0.032$	${0.985}_{-0.043}^{+0.036}$	$0.981\pm 0.032$	${1.006}_{-0.049}^{+0.038}$	$0.847\pm 0.037$	${0.853}_{-0.057}^{+0.034}$
${\sigma }_8$	${0.818}_{-0.054}^{+0.036}$	$0.826\pm 0.033$	${0.803}_{-0.045}^{+0.036}$	${0.798}_{-0.035}^{+0.030}$	${0.808}_{-0.052}^{+0.036}$	${0.692}_{-0.037}^{+0.014}$
${\Omega }_m$	${0.324}_{-0.0089}^{+0.023}$	${0.318}_{-0.016}^{+0.020}$	$0.316\pm 0.014$	${0.317}_{-0.016}^{+0.020}$	$>0.317$	${0.315}_{-0.017}^{+0.020}$	${0.313}_{-0.016}^{+0.019}$
$H_0$	${67.1}_{-3.1}^{+2.2}$	${67.5}_{-2.7}^{+2.4}$	${70.3}_{-0.96}^{+4.4}$	$68.1\pm 2.7$	${66.9}_{-3.4}^{+2.1}$	${67.7}_{-4.3}^{+3.3}$	${67.1}_{-4.4}^{+2.3}$
$n_s$	${0.961}_{-0.016}^{+0.021}$	$0.960\pm 0.016$	${0.959}_{-0.019}^{+0.016}$	${0.957}_{-0.021}^{+0.013}$	${0.9559}_{-0.024}^{+0.0096}$	$>0.961$	${0.962}_{-0.013}^{+0.021}$
${\Omega }_b$	${0.0505}_{-0.0039}^{+0.0049}$	$0.0493\pm 0.0042$	$>0.0535$	$0.0490\pm 0.0044$	${0.0498}_{-0.0042}^{+0.0051}$	${0.0485}_{-0.0064}^{+0.0043}$	${0.0473}_{-0.0076}^{+0.0023}$
$w_0$	$-{1.02}_{-0.18}^{+0.21}$	$-{0.98}_{-0.16}^{+0.22}$	$-{0.97}_{-0.15}^{+0.24}$	$-{1.00}_{-0.18}^{+0.21}$	$-1.03\pm 0.18$	$-1.02\pm 0.16$	$-{1.01}_{-0.17}^{+0.19}$
$A_{\text{lensing}}$	$1.020\pm 0.049$	N/A	N/A	N/A	$1.039\pm 0.048$	${0.970}_{-0.043}^{+0.048}$	$0.967\pm 0.048$

Equations

$${⟨z⟩}_{\text{clustering}}=\frac{\int dz{p}^2(z)(d{V}_c/dz{)}^{-1}z}{\int dz{p}^2(z)(d{V}_c/dz{)}^{-1}},$$

$${⟨z⟩}_{\text{lensing}}=\frac{\int d{z}_l{p}_l(z_l)w_l(z_l)z_l}{\int d{z}_l{p}_l(z_l)w_l(z_l)},$$

$$w_l(z_l)=D_L^{-2}(z_l)(1+z_l{)}^{-2}{\int }_{z_l}^{\infty }d{z}_s{p}_s(z_s){\Sigma }_c^{-2}(z_l,z_s),$$

$$\text{ratio}{w}_p(r_p;\mathbf{p})=\frac{w_{p,\text{sim}}(r_p;\mathbf{p})}{w_{p,\text{analytic}}(r_p;\mathbf{p})},$$

$$w_p(r_p)={\int }_{r_p}^{{\Pi }_{\text{max}}}{\xi }_{\text{gg}}\left(\sqrt{r_p^2+{\Pi }^2}\right)d\Pi ,$$

$$w_{\text{p,rsd}}(r_p)=\frac{1}{2{\pi }^2}{\int }_0^{\infty }d{k}_z{\int }_0^{\infty }d{k}_{\perp }{P}_{\text{gg}}(k){\left(1+\beta {\mu }^2\right)}^2$$

$$\times \cos (k_z{r}_{\pi })J_0(k_{\perp }{r}_p),$$

$$w_{\text{p,norsd}}(r_p)=\frac{1}{2{\pi }^2}{\int }_0^{\infty }d{k}_z{\int }_0^{\infty }d{k}_{\perp }{P}_{\text{gg}}(k)$$

$$\times \cos (k_z{r}_{\pi })J_0(k_{\perp }{r}_p),$$

$$\text{RSD corr.}(r_p)=\frac{w_{\text{p,rsd}}(r_p)}{w_{\text{p,norsd}}(r_p)}$$

$$\Delta \Sigma (r_p)=\bar{\rho }\left[\frac{4}{r_p^2}{\int }_{r_{p,\text{min}}}^{r_p}r{w}_{\text{gm}}(r)dr-2w_{\text{gm}}(r_p)\right],$$

$$w_{\text{gm}}(r_p)={\int }_0^{{\Pi }_{\text{max}}}{\xi }_{\text{gm}}\left(\sqrt{r_p^2+{\Pi }^2}\right)d\Pi .$$

$$||\text{LOSOE}|{|}^2=\sum _i^N{\left({\hat{y}}_{(\text{sim}(i)),i}-y_i\right)}^2,$$

$${\text{corr}}_{ij}=C_{ij}/\sqrt{C_{ii}{C}_{jj}}.$$

$${\lambda }_i\gtrsim \sqrt{\frac{2}{N_{\text{subvol}}}},$$

$$\text{Cov}($$

$$\times \left[\left(P_{\text{gg}}(k)+\frac{1}{n_g}\right)\left(P_{\gamma \gamma }^{2D}(k)+\frac{{\sigma }_{\gamma }^2}{{\Sigma }_s}\right)+{\left(P_{\text{g}\gamma }^{2D}(k)\right)}^2\right]$$

$$P_{\gamma \gamma }^{2D}={\int }_0^{{\chi }_s}d\chi {\left(\frac{\bar{\rho }}{{\Sigma }_c(\chi ,{\chi }_s)}\right)}^2{P}_{\text{mm}}\left(k\frac{{\chi }_l}{\chi }\right),$$

$$(P_{\text{g}\gamma }^{2D}{)}^2\approx {\Pi }_{\text{lens}}{\left(\frac{\bar{\rho }}{{\Sigma }_c({\chi }_l,{\chi }_s)}\right)}^2{P}_{\text{gm}}^2(k).$$

$${\Pi }_{\text{lens}}\equiv {\int }_0^{{\chi }_s}d\chi {\left(\frac{{\Sigma }_c({\chi }_l,{\chi }_s)}{{\Sigma }_c(\chi ,{\chi }_s)}\right)}^2,$$

$${\tilde{C}}^{-1}=(C+{\sigma }^{2}v{v}^T{)}^{-1}=C^{-1}-\frac{{\sigma }^2{C}^{-1}v{v}^T{C}^{-1}}{1+{\sigma }^2{v}^T{C}^{-1}v},$$

$${\tilde{C}}^{-1}=\underset{{\sigma }^2\to \infty }{\lim }(C+{\sigma }^{2}v{v}^T{)}^{-1}=C^{-1}-\frac{C^{-1}v{v}^T{C}^{-1}}{v^T{C}^{-1}v},$$

$$\ln \mathcal{L}=-\frac{1}{2}\Delta y^T{\tilde{C}}^{-1}\Delta y-\frac{1}{2}\ln \det \tilde{C}-\frac{1}{2}{N}_{\text{obs}}\ln 2\pi .$$

$$\Delta \Sigma (r_p)=\Delta \tilde{\Sigma }(r_p)+{\left(\frac{1}{r_p}\right)}^2\frac{M_0}{\pi },$$

$$\Delta \Sigma (r_p)$$

$$=\bar{\rho }\left[\frac{4}{r_p^2}{\int }_{r_{\text{p,min}}}^{r_p}r{w}_{\text{gm}}(r)dr-2w_{\text{gm}}(r_p)\right]$$

$$+\bar{\rho }\left[\frac{4}{r_p^2}{\int }_0^{r_{\text{p,min}}}r{w}_{\text{gm}}(r)dr\right],$$

$$=\Delta \tilde{\Sigma }(r_p)+\bar{\rho }\left[\frac{4}{r_p^2}{\int }_0^{r_{\text{p,min}}}r{w}_{\text{gm}}(r)dr\right],$$

$$=\Delta \tilde{\Sigma }(r_p)+{\left(\frac{r_{\text{p,min}}}{r_p}\right)}^2\bar{\Sigma }(r_{\text{p,min}}),$$

$$\mathcal{Z}=\int \mathcal{L}(\mathbf{p})p(\mathbf{p})d\mathbf{p},$$

$${\xi }_{\text{gg}}(r_p,\Pi )=\frac{DD-2DR+RR}{RR},$$

$$w_p(r_p)=2{\int }_0^{{\Pi }_{\text{max}}}{\xi }_{\text{gg}}(r_p,\Pi )d\Pi$$

$${\xi }_{\text{gm}}(r_p,\Pi )=\frac{D_1{D}_2-D_{1}R-D_{2}R+RR}{RR}$$

$$w_{\text{gm}}(r_p)={\int }_0^{{\Pi }_{\text{max}}}{\xi }_{\text{gm}}(r_p,\Pi )d\Pi ,$$

Full text

††pubyear: 2019††pagerange: Cosmology with galaxy-galaxy lensing on non-perturbative scales: Emulation method and application to BOSS LOWZ–D

1 Introduction

Weak gravitational lensing has emerged as the most powerful probe of matter clustering in the low redshift universe, critical to testing whether cosmic acceleration is caused by a cosmological constant, by an alternative form of dark energy, or by a breakdown of General Relativity on cosmological scales. Cosmic shear measures the clustering of foreground dark matter from the correlations induced in the shape of background source galaxies. Galaxy-galaxy lensing (GGL) uses a background shear map to measure the clustering of matter around a foreground galaxy population, which can be combined with the foreground galaxy clustering itself to infer the underlying dark matter clustering. In applications to state-of-the-art weak lensing data sets, the two approaches have comparable statistical uncertainties, with systematics that are partly in common and partly distinct (e.g. Hildebrandt et al. 2017; DES Collaboration et al. 2017). Building on our previous work (Wibking et al., 2019), this paper presents a numerical approach to modeling galaxy clustering and GGL into the deeply non-linear regime and applies it to weak lensing and galaxy clustering measurements (Singh et al., 2018) from the LOWZ galaxy sample of the Baryon Oscillation Spectroscopic Survey (BOSS; Eisenstein et al. 2011; Dawson et al. 2013), including tests on the LOWZ mock catalogs used by Singh et al. (2018).

On asymptotically large scales, where linear theory and scale-independent galaxy bias should be exact, one can think of galaxy clustering + GGL as measuring $\xi_{\text{gg}}=b_{g}^{2}\xi_{\text{mm}}$ and $\xi_{\text{gm}}=b_{g}\xi_{\text{mm}}$, allowing cancellation of the unknown $b_{g}$ and inference of $\xi_{\text{mm}}$. The accuracy of this method can be improved by using higher order perturbative models of galaxy bias (reviewed by Desjacques et al. 2018). However, these models break down on comoving scales below $\sim 10\,h^{-1}\,\text{Mpc}$, so non-perturbative models are needed to exploit clustering and GGL data on the $\sim$ Mpc scales where they are most precise. Demands on the accuracy of model predictions will become more stringent with the completion of current generation weak lensing surveys such as the Kilo-Degree Survey (KiDS; Hildebrandt et al. 2017), the Dark Energy Survey (DES; DES Collaboration et al. 2017), and the Hyper-Suprime Camera Strategic Survey Program (HSC SSP; Aihara et al. 2018), and with the advent of future surveys from the Euclid mission (Laureijs et al., 2011), the Large Synoptic Survey Telescope (LSST; The LSST Dark Energy Science Collaboration et al. 2018), and the Wide Field Infrared Survey Telescope (WFIRST; Doré et al. 2019). The dilemma of scales is already illustrated by existing analyses. Most studies of clustering + GGL on large scales infer an amplitude of matter clustering that is lower than predicted by a $\Lambda$CDM cosmological model (cold dark matter with a cosmological constant) normalized to Planck cosmic microwave background data (e.g. Mandelbaum et al. 2013; Hildebrandt et al. 2017; DES Collaboration et al. 2017; but see More et al. 2015 for a compatible result). However, the significance of the discrepancy with any one data set is limited because the statistical errors on these scales are large. Leauthaud et al. (2017) find a much stronger discrepancy on Mpc scales between measured GGL for BOSS CMASS galaxies and predictions of Planck-normalized mock galaxy catalogs tuned to CMASS galaxy clustering, but they are hesitant to draw strong conclusions because of theoretical uncertainties in the clustering models. Recently, Singh et al. (2018) used information down to $1\,h^{-1}\,\text{Mpc}$ scales with a nonparametric model of scale-dependent galaxy bias and inferred a lower amplitude of matter clustering than Planck at $>3\sigma$ significance. They likewise caution that uncertainties about modeling the galaxy population prohibit strong conclusions about cosmological physics.

In this work, we adopt the halo occupation distribution (HOD) model of galaxy bias (Jing et al., 1998; Peacock & Smith, 2000; Seljak, 2000; Scoccimarro et al., 2001; Berlind & Weinberg, 2002), which is commonly used as a model of galaxy clustering on $\sim$Mpc to sub-Mpc scales (e.g., Zehavi et al. 2005, 2011; Coupon et al. 2012; Sinha et al. 2018) and also in combination with GGL (e.g. Zu & Mandelbaum 2015). Several previous papers have advanced the idea of modeling non-linear galaxy clustering and GGL with HODs, in effect allowing the HOD to provide nuisance parameters that one can marginalize over when deriving cosmological constraints (e.g. Yoo et al. 2006; Cacciato et al. 2009, 2013; Yoo & Seljak 2012). Cacciato et al. (2013) apply this approach to Sloan Digital Sky Survey (SDSS) data and find $\sigma_{8}$ and $\Omega_{m}$ values in good agreement with the WMAP7 results (Komatsu et al., 2011) but low compared to recent values from Planck (Planck Collaboration et al., 2016). More et al. (2015) likewise apply this approach to CMASS galaxy clustering and Canada-France-Hawaii Telescope Lensing Survey (CFHTLS) shear catalogs and find $\sigma_{8}$ and $\Omega_{m}$ values in good agreement with both final WMAP and early Planck results. These papers have relied on analytic formulations of the HOD/halo model, which are accurate at the $\sim 5$ per cent level relative to numerical predictions from cosmological $N$-body simulations. It seems unlikely that a first-principles analytic approach can achieve the per cent-level accuracy demanded by current data sets, in part because of uncertainties in the effects of halo exclusion and scale-dependent halo bias (see e.g., van den Bosch et al. 2013).

By directly interpolating results from populated $N$-body simulations, on the other hand, we can compute the predictions of the halo model for the projected galaxy clustering $w_{p}$ and galaxy-galaxy lensing $\Delta\Sigma$ with sub-percent accuracy into the deeply non-linear regime. The main disadvantage of this approach is that it requires a large library of $N$-body simulations to sample the cosmological parameter space and many repopulations and pair counting computations of each of these simulations to sample the HOD parameter space. The range of galaxy samples that one can model is limited by the minimum mass of a well-resolved halo for a simulation of a given resolution ($\sim 300-500$ particles per halo; e.g. Trenti et al. 2010; Reed et al. 2013). In this paper, we use the AbacusCosmos suite of simulations (Garrison et al., 2017) sampling the parameters of $w$CDM cosmology, extending our previous work (Wibking et al., 2019) that used a grid of $(\sigma_{8},\Omega_{m})$ values within $\Lambda$CDM. A similar effort, using different simulations and a somewhat different approach to computing their predictions, has been undertaken by Zhai et al. (2019).

Interpolation across the outputs of simulations has become popularly known within the cosmology community as emulation (Heitmann et al., 2009; Kwan et al., 2015). As with many previous efforts, we use Gaussian processes, which are commonly adopted for relatively small multidimensional datasets that are expensive to obtain and where it is desired to propagate quantitative uncertainty estimates from the input training data to the resulting predictions.111Spline interpolation, for instance, becomes increasingly difficult in high dimensions, except when implemented as a Gaussian process kernel, in which case it generally is not a good choice of kernel function. The underlying idea of using Gaussian process regression methods to smoothly predict the output of computer simulations is at least three decades old (Sacks et al., 1989), and the general theory of using Gaussian processes for regression is even older (O’Hagan & Kingman, 1978), with its original one-dimensional time-series formulation due to Kolmogorov (1941) and Wiener (1949). To obtain sufficient accuracy for our interpolation of projected galaxy clustering $w_{p}$ across our parameter space, we emulate the ratio of the numerically computed $w_{p}$ to an approximate analytic calculation using the same cosmological and HOD parameters. For $\Delta\Sigma$, we find that direct emulation is sufficiently accurate.

Before applying our emulator to the BOSS LOWZ data, we test its ability to recover the correct cosmological parameters from the light-cone mock catalogs created by Nuza et al. (2013) and Rodríguez-Torres et al. (2016), as used in Singh et al. (2018). These mock catalogs use subhalo abundance matching (SHAM) with parameters designed to reproduce the clustering and number density properties of the LOWZ samples from $z=0.16-0.36$. Crucially, they are not created with an HOD prescription, so these tests provide at least some evidence that our parameterization is flexible enough to represent a range of possible scenarios for galaxy formation physics. We compute samples from the posterior conditioned on these mock data and show that we can recover the parameter $S_{8}\propto\sigma_{8}\Omega_{m}^{0.6}$ with bias $\leq 0.5\sigma$, even when using scales down to $\sim 0.6\,h^{-1}\,\text{Mpc}$, which is the minimum usable scale because of fiber collision systematics in the clustering measurements.

In section 2, we describe the simulations and methodology used to construct our emulator and the technical details necessary to obtain our target accuracy. Section 3 describes the construction of the covariance matrix for the BOSS LOWZ analysis. Section 4 presents the mock catalogs tests and section 5 the results from application to the Singh et al. (2018) GGL and clustering measurements. We conclude in section 6 with a discussion of the limitations of our method, proposed extensions and modifications for future work, and prospects for application to DES, HSC, WFIRST, and LSST.

We use comoving densities and distances throughout and assume flat $\Lambda$CDM for computations of these quantities, unless specified otherwise.

2 Emulator Construction

Figure 1 presents an overview of our emulation and inference framework. The critical inputs are the AbacusCosmos simulations, which we combine with HOD populations to train the Gaussian process emulator. We use 20 Abacus simulations of the fiducial cosmology given in Table 1 to compute part of the covariance matrix for the BOSS LOWZ data and also to correct the mean model predictions for cosmic variance of a single simulation volume. Applying the emulator to the data with these covariances yields our cosmological parameter constraints marginalized over HOD parameters and nuisance parameters describing possible weak lensing systematics.

2.1 Parameter space

We initially chose a halo occupation distribution (HOD) parameter space designed to encompass the posterior parameter ranges for the LOWZ sample found by Parejko et al. (2013) when fitting a simulation-based HOD model to the projected galaxy correlation function $w_{p}$ at fixed cosmology ($\log M_{\text{min}}=13.25\pm 0.26$, $\sigma_{\log M}=0.43\pm 0.25$, $\log M_{0}=13.27^{+0.49}_{-0.76}$, $\log M_{1}=14.18\pm 0.39$, $\alpha=0.94\pm 0.49$), recast in the parameterization used by Wibking et al. (2019) (hereafter Paper I). However, we found that this parameter space was not sufficient to fit the mock galaxy catalogs we used for tests in section 4, and so we extended the parameter space to include two additional parameters, $A_{\text{conc}}$ and $R_{\text{rescale}}$, which account for differences between galaxy and halo profiles. We also extended the minimum of the ranges of the dimensionless HOD parameters $M_{0}/M_{1}$ and $\sigma_{\log M}$ to extend to zero and increased the maximum range of the satellite galaxy mass parameter $M_{1}/M_{\text{min}}$ to 20. We determine the characteristic mass scale of halos $M_{\text{min}}$ from the desired number density of galaxies $n_{\text{gal}}$ by performing an integral over the mass function tabulated from a given simulation we wish to populate with galaxies This means that the halo occupation $\langle N_{\text{gal}}|M_{\text{halo}}\rangle$ is not fully specified by the HOD parameters alone, and must be emulated as a function of the full set of parameters in Table 1.

The $A_{\text{conc}}$ parameter is a multiplicative correction parameter to the concentrations used for the galaxy number density profile within halos, where these concentrations are taken from the $v_{\text{max}}$-based concentrations (Klypin et al., 2011) determined for each dark matter halo within the simulations by the Rockstar halo finding code (Behroozi et al., 2013). As shown in Paper I, this parameter is essentially equivalent on the scales considered in this work to varying the power-law slope of the NFW profile, but it has the advantage that its convolution with itself can be written in closed form in configuration space (Zheng & Weinberg 2007; see also Appendix A). We replace the power-law slope variation parameter $\Delta\gamma$ with $A_{\text{conc}}$ in this paper.

The parameter $R_{\text{rescale}}$ is a multiplicative correction factor to the halo radii used for the galaxy number density profiles, where the halo radii are the virial radii likewise determined for each halo by Rockstar. While we only have halo catalogs for spherical overdensity halos determined by the virial overdensity criterion of Bryan & Norman (1998), and thus cannot marginalize over halo definition, marginalizing over the $R_{\text{rescale}}$ parameter allows us to take into account the main effect of varying halo definitions, namely the dependence of halo radii on halo definition. Since the extent and spatial distribution of ‘satellite’ galaxies is uncertain (and empirically, for luminous galaxy samples, is not the same as that of the dark matter; e.g., Watson et al. 2010; Piscionere et al. 2015), we populate halos with satellite galaxies according to an NFW profile (Navarro et al., 1997) with concentration and radius corrected by $A_{\text{conc}}$ and $R_{\text{rescale}}$ from the concentration and radius determined by Rockstar for a given halo.

We use a cosmological parameter space that approximately encompasses the union of the WMAP9 (Hinshaw et al., 2013) and Planck (Planck Collaboration et al., 2016) $w$CDM posteriors, as given by the design of the AbacusCosmos simulations (Garrison et al., 2017). The ranges of these parameters are shown in Table 1 (these ranges are not uniformly sampled, but are rather sampled along the principal components of the combined WMAP and Planck posteriors). We use the $(720\,h^{-1}\,\text{Mpc})^{3}$ set of simulation boxes run with identical phases of their initial conditions, with particle mass $\sim 1\times 10^{10}\,M_{\odot}h^{-1}$ at the fiducial (Planck) cosmology. Since LOWZ galaxies typically live in halos of mass $\sim 10^{13}\,M_{\odot}h^{-1}$, these halos are well resolved with $\sim 10^{3}$ particles per halo (although it is important to resolve halos down to the mass scale $\log M_{\text{halo}}\sim(\log M_{\text{min}}-\sigma_{\log M})$ in order to avoid biases in clustering predictions; see e.g. Sinha et al. 2018).

We use simulation outputs at redshift $z=0.3$, which is close to the effective redshifts of the LOWZ sample. There are two distinct effective redshifts, one for clustering and one for lensing, due to the differing line-of-sight weight functions for each signal. We compute the effective clustering redshift as

[TABLE]

where $p(z)$ is the weighted redshift distribution of the number of galaxies $dN_{g}/dz$ (Mandelbaum et al., 2011). For the LOWZ sample we use in this work, $\langle z\rangle_{\text{clustering}}\approx 0.27$. We compute the effective lensing redshift as

[TABLE]

with

[TABLE]

where $D_{L}(z)$ is the luminosity distance at redshift $z$, $p_{s}(z)$ is the source redshift distribution (weighted by the inverse variance of the shape measurements), $\Sigma_{c}$ is the lensing critical density, and the factor $(1+z_{l})^{-2}$ accounts for the fact that in comoving coordinates lower redshift lens galaxies have larger effective apertures (Nakajima et al., 2012) This window function is equivalent to that given by Singh et al. (2019), assuming that differences between photometric and true source galaxy redshifts are negligible. For the spectroscopic sample and lensing catalog used here, we find $\langle z\rangle_{\text{lensing}}\approx 0.24$.

Since these redshifts are both close to the AbacusCosmos simulation outputs at $z=0.3$ (Garrison et al., 2017), we adopt $z_{\text{eff}}=0.3$ as the effective redshift for all of our emulator predictions and neglect the difference between the clustering and lensing effective redshift. The adequacy of this effective redshift approximation for clustering is tested by fitting to mock catalogs (section 4), while the effective redshift approximation for lensing is only approximately tested by our mock catalog fits, due to the approximation we use to compute the lensing signal from the mock catalogs. However, the tests of Singh et al. (2019), conducted by applying the lensing weights (eq. 3) to particle-galaxy pairs in simulations, indicate that computing lensing predictions at a single effective redshift is accurate to better than $1-1.5$ per cent in this redshift range.

2.2 Sampling strategy

We use a modified stratified sampling design for the training dataset, generating sampling designs separately for the cosmological and non-cosmological parameters. The sampling design for the cosmological parameters, described by Garrison et al. (2017), is a min-max Latin hypercube (a variant of the classical Latin hypercube with additional space-filling properties described by Heitmann et al. 2009) with sampling dimensions given by the principal components of the combined WMAP and Planck posteriors and $N=40$ sub-cells per dimension (therefore yielding 40 samples). For the non-cosmological parameters, we use a Latin hypercube sampling design (McKay et al., 1979) with dimensions given by the nominal parameters in Table 1 and $N=400$ cells per dimension. Once we have generated the discrete Latin hypercube, we sample uniformly within each occupied cell of the hypercube to obtain the parameter samples, yielding 400 total samples.

We then have 400 non-cosmological parameter samples to distribute among 40 simulations. We carry this out without duplicating any non-cosmological parameter samples. For each cosmological simulation (corresponding to a realization of one of the cosmological parameter samples), we assign 10 non-cosmological parameter samples to the given simulation box from the set of unassigned non-cosmological samples. There is therefore no special relationship assumed by the sample design between the cosmological subspace of parameters and the non-cosmological subspace of parameters. We experimented with various other sampling strategies for the subspace of non-cosmological parameters, including Latin hypercubes with additional symmetry or volume-filling properties as well as uniform random samples, and found none that offered an improvement in emulator accuracy. We note that a hybrid combination of clustered and space-filling sampling designs shows promise in other application domains (Zhu & Zhang, 2006; Zimmerman, 2006), and we may consider such designs in future work.

2.3 Emulated quantities

2.3.1 Projected galaxy correlation function

We emulate the ratio of the projected galaxy-galaxy clustering $w_{p}$ relative to its analytic halo model prediction:

[TABLE]

where $\mathbf{p}$ is the vector of parameters described in Table 1. We compute $w_{p,\text{sim}}(r_{p})$ via a projection integral over the real-space $\xi_{gg}$ computed from the simulations:

[TABLE]

where we choose $\Pi_{\text{max}}=100\,h^{-1}\,\text{Mpc}$ and use a piecewise linear integration scheme to minimize finite bin-size effects (see Appendix D). This quantity is averaged over 20 stochastic realizations of the halo occupation in order to reduce noise; for a given set of parameters $\mathbf{p}$, the mean and variance of these realizations is used as the input datapoint to the Gaussian process model (see Appendix C). The analytic $w_{p}(r_{p})$ is computed by the same projection integral over the analytic $\xi_{\text{gg}}$.

We correct for residual RSD effects on $w_{p}$ (which are $\approx 15$ per cent at the largest scales considered in this work) by computing the ratio between $w_{p}$ computed with and without the Kaiser (1987) model for redshift-space distortions on the galaxy power spectrum (van den Bosch et al., 2013):

[TABLE]

[TABLE]

where $k^{2}=k_{\perp}^{2}+k_{z}^{2}$, $\mu=k_{z}/k$, and $\beta=f(z)/b$, with $f(z)$ as the growth factor at redshift $z$ and $b$ as the large scale galaxy bias (measured from $\sqrt{P_{\text{gg}}(k)/P_{\text{mm}}(k)}$ averaged over $k\lesssim 0.1$). The galaxy power spectrum $P_{\text{gg}}$ is computed via an integral over the real-space correlation function $\xi_{\text{gg}}$ measured in the simulations at the fiducial parameter values given in Table 1. We then multiply the emulated value of $w_{p}(r_{p})$ by the ratio

[TABLE]

to obtain a scale-dependent correction factor for residual RSD effects. While this is not strictly applicable except for the fiducial parameter values, we find that this ratio is relatively insensitive to parameter changes and adopt it as a fixed correction factor for $w_{p}(r_{p})$ at all parameter values.222Computed with traditional multidimensional cubature, the integrals involved are relatively expensive to compute. In future work, we recommend applying the FFTLOG method (Hamilton, 2000) to compute this integral (for a more complicated application, see McEwen et al. 2016). We test the accuracy of this approximation by fitting mock catalogs with observables computed in redshift space (section 4).

2.3.2 Galaxy-galaxy lensing

For $\Delta\Sigma$, we directly emulate the observable, since it is substantially smoother (and therefore the training data have higher signal-to-noise). $\Delta\Sigma$ is computed by integrating over $w_{\text{gm}}$, which in turn is computed by projecting $\xi_{\text{gm}}$ from the simulations:

[TABLE]

where

[TABLE]

In practice both the lower limit $r_{p,\text{min}}$ and the upper limit $\Pi_{\text{max}}$ of these integrals are taken to be finite values such that the integrals converge within an acceptable precision (in our case, we use a minimum integration limit of $0.01\,h^{-1}\,\text{Mpc}$ and a maximum integration limit of $100\,h^{-1}\,\text{Mpc}$). We likewise average this quantity over 20 realizations of the halo occupation in order to reduce noise in the training data, and use the mean and variance of these realizations as the emulator input datapoint $(y_{\text{obs,i}},\,\sigma_{\text{obs,i}}^{2})$ for a given set of parameters $\mathbf{p}_{i}$.

Finally, we correct $\Delta\Sigma$ as measured in the simulations to the amplitude and radial binning an observer making the assumption of an $\Omega_{m}=0.3$ cosmology would measure. As noted in Paper I, this correction for projected distances is negligible, but the assumed-cosmology correction to the critical lensing density $\Sigma_{c}$ changes the degeneracy direction between $\Omega_{m}$ and $\sigma_{8}$ by $\Omega_{m}^{-0.1}$ (due to geometric effects on the line-of-sight lensing distance).

2.3.3 Ensemble mean correction

Since our varying-cosmology simulations are computed with matched phases of their initial conditions (see Garrison et al. 2017), we first compute the fiducial matched-phases signal from the fiducial parameters (given in Table 1) with the same phases as those used to construct the emulation training data. We then compute the ensemble mean fiducial signal by averaging over 20 additional simulations that provide realizations of the fiducial cosmology with varying phases. We then divide the ensemble mean signal by the fiducial matched-phases signal to obtain a multiplicative correction for the observables due the sample variance of our simulations.

This procedure reduces the sample variance of our training set data, the accuracy of which depends on the assumption that the derivatives of the simulated correlation functions with respect to cosmology and HOD parameters are independent of the phases of the initial conditions. Our lightcone mock tests indicate that this correction slightly reduces bias of the recovered cosmological parameters, but that our inferences of cosmological parameters are relatively unbiased even without this correction. For future datasets, a more sophisticated approach may be needed, such as the “fixed-paired” approach to initial conditions advocated by Angulo & Pontzen (2016).

2.3.4 Analytic halo model for $w_{p}$

The analytic halo model for $w_{p}$ consists of a 1-halo term given by the sum of the central-satellite and satellite-satellite terms and a simple two-halo term given by the square of the large-scale bias multiplied by the linear matter correlation function computed with the ‘no-wiggles’ fitting formula for the linear matter power spectrum of Eisenstein & Hu (1998). For speed, we compute the one-halo terms in configuration space, using the closed-form expression for the self-convolution of the NFW profile given by Sheth et al. (2001) and Zheng & Weinberg (2007). A complete description of the equations used in our analytic halo model is given in Appendix A.

We emulate the ratios of the observables using Gaussian process regression (Rasmussen & Williams, 2006) with a squared-exponential kernel. For each radial bin of each of the observables, we determine the hyperparameters of the kernel by maximizing the leave-one-out cross-validation pseudo-likelihood (Appendix C). This pseudo-likelihood directly minimizes the prediction error, which may be more robust to the (arbitrary) choice of kernel function than maximizing the marginal likelihood of the Gaussian process model itself (Rasmussen & Williams, 2006). We thus obtain separate Gaussian process emulators for each radial bin of each emulated quantity ($w_{p,\text{sim}}/w_{p,\text{analytic}}$ and $\Delta\Sigma$).

2.4 Interpolation Accuracy

2.4.1 Cross-validation error

We measure our emulator accuracy by using the leave-one-simulation-out cross validation error (computed separately for each radial bin of each emulated quantity). This is defined as

[TABLE]

where $y_{i}$ is the $i$-th training datapoint value ($i=1,\dots,400$) and $\hat{y}_{(\text{sim}(i))}$ is the prediction of the emulator for point $i$ trained on all datapoints except for those derived from the same simulation as point $i$. This ‘leaves out’ the training data from a given simulation when predicting the observables at the same cosmological parameters as those of the given simulation, thereby providing a conservative estimate of emulator accuracy. Note that we fix the hyperparameters of the emulator (see Appendix C) when computing eq. 11.

We show emulator accuracy in Figure 5. We find that emulating in logarithmic variables does not improve the emulator accuracy, so we use the linear (i.e. untransformed) variables, in contrast to Paper I. We find that accuracy substantially increases when directly emulating the (projected) observables $w_{p}$ and $\Delta\Sigma$ instead of the three-dimensional correlation functions. In order to further improve accuracy, we emulate in the ratio of $w_{p}$ with respect to its analytic halo model prediction as described in section 2.3.1. We hypothesize that this improves emulator accuracy for $w_{p}$ at a fixed number of training points because the variance of the training data is reduced due to approximate prior knowledge of how the observable should respond to a parameter change.

We also compute the covariance matrix $C$ for the leave-one-simulation-out prediction errors of our emulators across different $r_{p}$ bins. While one may naively expect that the errors in predicting the values of the correlation function would not be correlated between bins because we train separate emulators for each bin, this is not the case. In Figure 6 we show the correlation matrix of the prediction errors for $w_{p}$, defined as

[TABLE]

We find that nearest-neighbor bins have highly correlated prediction errors, with a sharp drop-off in correlations for more distant bins, with the exception of bins $\gtrsim 3$ $h^{-1}\,\text{Mpc}$, which are all moderately correlated with each other. We emphasize that although the underlying reason for the correlated prediction errors is due to correlated uncertainties in the correlation function itself, the resulting pattern of correlated errors of the emulators is not identical to the form of statistical correlations in $w_{p}$ itself.

2.4.2 Emulator efficiency

Since we have virtually eliminated the noise contributions from sample variance and from the stochasticity of the HOD (via the averaging over multiple populations described in section 2.3.1), we can construct an emulator of a given accuracy using a much smaller training set than other emulators (e.g. Zhai et al. 2019, who used 2000 samples to train their GP-based emulator). This directly translates into two orders of magnitude (a factor of $[2000/400]^{3}=125$) smaller computational requirements for the hyperparameter optimization, since we gain efficiency by a factor $N^{3}$ in the kernel matrix inversion (eq. 98).

3 Covariance matrix

We use a combination of simulation-based and analytic covariance matrices for our mock cosmological analyses in section 4. For $w_{p}$, we use a simulation-based covariance to capture the non-Gaussian part of the covariance (excess variance compared to the variance of a Gaussian field with the nonlinear galaxy power spectrum), which we find is important on 1-halo scales. For $\Delta\Sigma$, we use an analytic covariance matrix to capture the Gaussian part of the covariance (the variance purely due to a Gaussian field with the nonlinear galaxy power spectrum), including shape noise and large-scale-structure noise. From simulations, we find that the non-Gaussian part of the $\Delta\Sigma$ covariance is negligible, and we therefore neglect it in our analysis.333The non-Gaussian $\Delta\Sigma$ covariance would not be neglible for cluster mass halos with low shape noise weak lensing data (Wu et al., 2019). Since the $\Delta\Sigma$ covariance matrix is dominated by shape noise, we find that the cross-covariance between $w_{p}$ and $\Delta\Sigma$ is negligibly small and so we do not include it in our analysis. For completeness, we derive the Gaussian part of the cross-covariance between $w_{p}$ and $\Delta\Sigma$ in Appendix B.

3.1 Simulation-based covariance for $w_{p}$

Due to the presence of a large non-Gaussian component on 1-halo scales, we use bootstrap-estimated simulation covariances to estimate the $w_{p}$ covariance matrix, populating simulation boxes with galaxies according to the fiducial HOD parameters shown in Table 1. We use 20 volumes of the Planck cosmology with varying phases from the ($1100$ $h^{-1}\,\text{Mpc}$)3 boxes of Garrison et al. (2017) and divide each volume into $5\times 5$ subvolumes which tesselate the $x$-$y$ plane in projection, for a total of $20\times 5\times 5=500$ subvolumes. We compute $w_{p}$ in projection within each subvolume, and compute 500 bootstrap resamples by choosing 500 subvolumes with replacement and averaging $w_{p}$ for each resample. Our estimate for the covariance of $w_{p}$ is then the covariance among the bootstrap resamples rescaled to the effective volume $V_{\text{eff}}$ of the LOWZ galaxy sample.

Finally, we compute the eigenvalues and eigenvectors of this covariance matrix, in order to check for noisy modes that may bias the inverse covariance matrix. A common noise threshold is

[TABLE]

where $\lambda_{i}$ is a given eigenvalue of the covariance matrix and $N_{\text{subvol}}=500$ is the number of semi-independent subvolumes used to estimate the covariance (e.g., Gaztañaga & Scoccimarro 2005; Sinha et al. 2018; Zhai et al. 2019). While it is the most conservative choice to eliminate any potentially noisy modes when computing the inverse covariance matrix by setting the eigenvalues below the noise threshold (eq. 13) to zero (then taking the reciprocal of the non-zero eigenvalues to yield the inverse covariance matrix), we did not do so in our analysis. We show the eigenvalues in Figure 7 in order to illustrate that no modes are particularly noisy and therefore any noise bias in the inverse covariance matrix should be minimal.

3.2 Analytic covariance for $\Delta\Sigma$

We use a similar form for the Gaussian covariance for $\Delta\Sigma$ as in Paper I, but we use a slightly more accurate treatment of the large-scale-structure noise contribution (equivalent to eq. A36 of Singh et al. 2017b). As we derive in Appendix B, the full (Limber-approximate) form of this covariance is

[TABLE]

where

[TABLE]

[TABLE]

Here $\Pi_{\text{lens}}$ is the effective line-of-sight depth of the squared lensing weight function, defined as

[TABLE]

and $\Sigma_{c}$ is the lensing critical surface density, $\Sigma_{s}$ is the projected surface density of source galaxies (in units $[h^{-1}\,\text{Mpc}]^{-2}$), $\sigma_{\gamma}$ is the shape noise per galaxy, $P_{\text{gg}}$ is the 3D galaxy power spectrum, $P_{\text{mm}}$ is the 3D matter power spectrum, $P_{\text{gm}}$ is the 3D galaxy-matter cross spectrum, $\chi_{l}$ is the effective distance of the lens galaxy population, $\chi_{s}$ is the effective distance of the source galaxy population, and $\bar{J}_{2}$ denotes the bin-averaged Bessel function of order 2. The galaxy power spectra are computed from integrals over the real-space galaxy correlation functions computed at the fiducial parameters (Table 1). We note that it is necessary to average the Bessel functions over the bin widths before computing the integral (eq. 3.2) in order to avoid divergence (see Appendix B for details).

3.3 Comparison to survey jackknife

We find that our simulation-based covariance for $w_{p}$ agrees in magnitude almost perfectly with a jackknife covariance produced from the LOWZ mocks used in section 4.

The $\Delta\Sigma$ covariance cannot be directly compared to the mocks, because raytraced lensing was not done for the mocks and so any estimates from the mocks do not include shape noise or any line-of-sight large scale structure noise. Instead, we compare the amplitude of the jackknife covariance from the LOWZ data to that of our analytic covariance, finding that the analytic covariance (at our fiducial parameter choices) has error bars that are $10-15$ per cent larger than that of the data jackknife, which we deem to be acceptable agreement given the difference of methodology.

3.4 Point-mass marginalization

Since our simulations strictly only provide a prediction for $\Delta\Sigma$ up to a constant in enclosed projected mass (due to unresolved substructure at very small scales), it is necessary to include a point mass term in our cosmological analyses. This term can also absorb baryonic physics effects (such as dissipation and feedback) that are not represented in our $N$-body calculations. We use the Sherman-Morrison matrix identity in the form

[TABLE]

where $\mathbf{v}$ is a column vector with values $[r_{p,0}^{-2}\,,r_{p,1}^{-2}\,,...\,,r_{p,N}^{-2}]$, and take the limit

[TABLE]

in order to modify our inverse covariance matrix in a way that analytically marginalizes the contribution from a point mass (e.g., MacCrann et al. 2019). The modified inverse covariance matrix $\tilde{C}^{-1}$ is then used in our likelihood function as such:

[TABLE]

This procedure is equivalent to explicitly marginalizing a point mass term of the form:

[TABLE]

with an infinitely wide Gaussian prior on the amplitude of $M_{0}$. We favor this procedure over explicitly adding a point mass parameter, since no additional computational expense is added to the posterior sampling compared to neglecting the point mass term.

Adopting a point mass is equivalent to marginalizing over the uncertainty in the un-modeled inner region of the galaxy-mass cross correlation:

[TABLE]

where we can identify the enclosed mean projected mass $\pi(r_{\text{p,min}})^{2}\,\bar{\Sigma}(r_{\text{p,min}})$ with the point mass $M_{0}$. A related (but not equivalent) approach is that of Baldauf et al. (2010), who construct an estimator that effectively subtracts the point mass term from the observations at the cost of increased noise. This approach is used by Singh et al. (2018). Our approach does not increase the noise of the lensing signal itself but instead requires marginalization over the amplitude of the point mass.

Since our procedure makes the updated covariance matrix formally singular (i.e. $\det\tilde{C}=0$), we likewise update the $\frac{1}{2}\ln\det\tilde{C}$ term of eq. 20 according to eq. 10 of Bridle et al. (2002) in order to correctly normalize the likelihood and therefore obtain a correct value of the Bayesian evidence integral

[TABLE]

where $\mathbf{p}$ is the vector of parameters, $p$ is the prior (defined in Table 2), $\mathcal{L}$ is the likelihood (implicitly dependent on the observed data; eq. 20), and $\mathcal{Z}$ is the evidence, or normalization constant for the posterior distribution.

4 Cosmological analysis on mocks

We test the ability of our emulator-based model to recover the correct cosmological parameters from a mock galaxy catalog produced from an independent simulation (Klypin et al. 2016; run with a different $N$-body code, GADGET-2; Springel 2005) with an independent method for populating halos with galaxies (subhalo abundance matching; e.g. Kravtsov et al. 2004; Vale & Ostriker 2004). This is the same galaxy sample lightcone as used in Singh et al. (2018), produced with subhalo abundance matching (SHAM) tuned to match the number density and clustering of LOWZ galaxies in three disjoint redshift bins in the range $0.16<z<0.36$ (originally described in Nuza et al. 2013; Rodríguez-Torres et al. 2016). The lightcone was generated by populating subhalos found with Rockstar (Behroozi et al., 2013), using its default halo definition corresponding to the virial mass given by the fitting function of Bryan & Norman (1998). Because we have only a single lightcone volume, we can only test for biases of our parameter constraints at the level of the statistical uncertainty of the LOWZ data set. However, we do construct variants of the lightcone catalog to test for specific possible systematic effects as discussed below.

The projected clustering $w_{p}$ is computed by converting the mock catalog galaxy coordinates into redshift space comoving distances and counting galaxy pairs to tabulate the statistic $\xi_{\text{gg}}(r_{p},\Pi)$ in transverse bins $r_{p}$ and line-of-sight bins $\Pi$ with the Landy-Szalay estimator (Landy & Szalay, 1993),

[TABLE]

then integrating along the line of sight:

[TABLE]

with $\Pi_{\text{max}}=100$ $h^{-1}\,\text{Mpc}$.

The galaxy-galaxy lensing signal $\Delta\Sigma$ is computed from the mock by computing the galaxy-matter correlation function

[TABLE]

where $D_{1}$ corresponds to galaxies and $D_{2}$ corresponds to matter particles and integrating along the line-of-sight to obtain

[TABLE]

where $\Pi_{\text{max}}=100\,h^{-1}\,\text{Mpc}$, then integrating once more (using eq. 9) to obtain $\Delta\Sigma$, using $r_{p,\text{min}}=0.1\,h^{-1}\,\text{Mpc}$. We add an additional term to $\Delta\Sigma$ based on a power-law extrapolation of $w_{\text{gm}}(r_{p})$ in order to account for the projected mass at scales $<r_{p,\text{min}}$:

[TABLE]

where we assume $w_{\text{gm}}(r_{p})\propto 1/r_{p}$.

This calculation does not precisely respect the true redshift weighting of the galaxy-shear signal over the range of the LOWZ sample, which has a redshift-dependent lensing kernel. A fully correct calculation requires raytraced simulations with resolved halo substructure and an explicit source galaxy population, which are not available to us. However, the difference between the true $\Delta\Sigma$ and the approximate signal used here has been calculated (using the lensing weights in eq. 3 applied to simulation particle-galaxy pairs) to be $\lesssim 1-1.5$ per cent at any projected scale used here (Singh et al., 2019).

4.1 Test of Emulator

For the fiducial mock, we adopt the emulator model described in section 2 and the covariance matrix described in section 3 and compute the Bayesian evidence of the model via nested sampling (Skilling, 2004) with the MultiNest implementation (Feroz & Hobson, 2008). This method integrates Eq. 23 by sampling from the prior and then rejecting samples that lie outside of successively-smaller nested likelihood contours. The MultiNest implemention of nested sampling approximates each likelihood contour with a set of bounding ellipsoids which attempt to fully enclose the given likelihood contour, thereby gaining efficiency compared to rejection sampling from the entire prior volume. The sampling stops when the code estimates that the log evidence lower bound computed via nested sampling is within 0.05 of the log evidence upper bound estimated via the size of the remaining likelihood contour. The most significant failure mode of this method is to fail to sample one of the modes of a multi-modal posterior in the first iteration of the algorithm, which we attempt to guard against by checking that the posterior and evidence are stable against increasing the initial number of sampling points (so-called ‘live points’) by a factor of 2.444We use at least 400 initial samples for all posteriors and evidence computations used here. We do not use importance nested sampling but instead use the more conservative classical MultiNest algorithm as implemented in MultiNest version 3.11. As a byproduct of the Monte Carlo evaluation of the evidence integral, we also obtain the posterior distribution of the parameters (Skilling, 2004). This method has compared favorably with Markov Chain Monte Carlo posterior estimation in the DES Collaboration et al. (2017) cosmological analysis.

As a conservative test of the accuracy of our emulator, we do not add any additional uncertainties to our analysis covariance matrix, with the understanding that any emulator inaccuracies (due to interpolation error or the finite volume of our simulations, or due to systematic errors in our methodology) could cause biases in parameter recovery. In this case, a lack of identifiable bias will indicate that any emulator inaccuracies are subdominant to the statistical precision of the measurements, as we will show is the case in the following sections.

We adopt the prior ranges shown in Table 2.

We show the projections of our parameter posteriors in Figure 8 and tabulate the posterior means and 68 per cent credible intervals in the leftmost column of Table 3. The parameter of greatest interest is the combined cosmological parameter $S_{8}$, defined as

[TABLE]

which is chosen such that it is approximately the best-constrained combination of $\Omega_{m}$ and $\sigma_{8}$ (as determined in Paper I), and is normalized to unity for the true values of $\sigma_{8}$ and $\Omega_{m}$ used in the lightcone simulation, so as to be consistent with the definition used by Singh et al. (2018).555We note that this definition is not consistent with that used by DES Collaboration et al. (2017). For the fit to the fiducial mock, we recover this parameter to within 1.7 per cent of its true value, which corresponds to approximately half of its posterior uncertainty. Because we have only a single mock realization of the LOWZ volume, we cannot test our recovery method at a precision higher than the statistical error of the LOWZ data.

Since our model is nonlinear, we cannot use the classical $\chi^{2}$ goodness-of-fit test to assess whether our model is consistent with the mock data.666This is because the number of degrees of freedom of a nonlinear model is not well-defined. See Andrae et al. (2010) for a pedagogical explanation. Instead, we use the Bayesian posterior predictive discrepancy with a $\chi^{2}$-like test statistic (Guttman, 1967; Gelman et al., 1996) to quantify the tension of the posterior with the mock data. The test statistic $\Delta\chi^{2}$ is computed for the data by finding the minimum $\chi^{2}$ between the data vector and all of the posterior samples. In the limit of a large number of posterior samples, this is equivalent to the classical $\chi^{2}$ statistic computed from the data. However, the reference distribution differs from the classical test, and we compute the reference (i.e. expected) distribution of $\Delta\chi^{2}$ via a Monte Carlo procedure. For each posterior sample, we add noise by sampling with our fiducial covariance matrix, then compute the minimum $\chi^{2}$ between the noisy sample and all other posterior samples.777We choose a ‘minimum discrepancy’ statistic rather than an ‘average discrepancy’ statistic (Gelman et al., 1996) because the latter is more sensitive to the normalization of the covariance matrix, which we do not attempt (and do not need, for the purpose of parameter inference) to compute to the same accuracy as our model predictions. By comparing the $\Delta\chi^{2}$ from the data with this reference distribution, we obtain a convenient measure of how consistent the data are with being drawn from the posterior distribution. We find that the posterior test statistic $\Delta\chi^{2}$ for the mock data is smaller than that of $68.3$ per cent of the posterior samples. Interpreting this result as a frequentist statistical test, we conclude that there is no significant evidence against the null hypothesis that the data are drawn from the posterior. We emphasize that this test cannot give the probability that our model is correct (for which Bayesian model comparison is required), but rather only suggests that our model (considered without reference to other models) is an adequate description of the mock data given the uncertainties.

4.2 Varying satellite fractions and incompleteness

We vary the satellite fraction in the mock, as described by Singh et al. (2018), and recompute the posteriors with MultiNest. For both a 15 per cent lower satellite fraction and a 15 per cent higher satellite fraction, we recover the true value of $S_{8}$ within $\sim 0.5\sigma$.

We test the effect of a 15 per cent incompleteness fraction (applied uniformly to both central and satellite galaxies) by increasing the number density and then stochastically removing 15 per cent of all galaxies such that the observed number density of the mock catalog is consistent with the observed number density of LOWZ galaxies. For this mock, we recover the true value of $S_{8}$ within $\sim 0.5\sigma$ (see Table 3).

4.3 Parameterized lensing systematics

For the previous tests, we assumed that neither our emulator nor our mock data contributed any additional systematic errors, and so the resulting biases would be conservative. However, we computed an additional mock recovery test including a multiplicative lensing calibration parameter to account for the combined effects of multiplicative shear calibration bias and photometric redshift calibration bias (which becomes multiplicative in $\Delta\Sigma$ through its multiplicative effect on the critical lensing surface density $\Sigma_{c}$). The prior on this calibration parameter was set equal to $\mathcal{N}(\mu=1,\sigma=0.06)$, consistent with the calibration prior used on the analysis with LOWZ data (see following section). We did not change the mock data themselves, so the true value of this parameter for all of our lightcone test scenarios is $A_{\text{lensing}}=1.0$. In this scenario, we find the recovered $S_{8}$ to be within 2 per cent of the true value (within $\sim 0.3\sigma$), and we find the recovered calibration parameter $A_{\text{lensing}}$ to be $\sim 3$ per cent lower than the true value (within $\sim 0.5\sigma$), as shown in Figure 8. Adding the $A_{\text{lensing}}$ parameter increases the $S_{8}$ uncertainty from 0.029 to $0.048$, implying that 40 per cent of the error budget is due to systematic uncertainties in the lensing signal. We show the posterior mean values of our parameters for all of the tests of this section in Table 3.

5 Cosmological analysis on data

We determined the design of the emulator, computed the covariance matrices, and performed all tests on the mock catalogs discussed in the previous sections prior to computing any posteriors from the data. We formally blinded our analysis by multiplying the observed $\Delta\Sigma$ (on all scales) by a constant drawn from $\mathcal{N}(\mu=1.0,\sigma=0.06)$ and only known to one of us before unblinding.888We note that in the context of Bayesian model averaging, which we advocate in section 6, blinding is superfluous, inasmuch as confirmation bias is caused by allowing for post-hoc model selection. In practice, however, there may be some benefit. The width of this distribution was chosen to be moderately larger than the expected statistical precision of our posterior inference on $S_{8}$. From tests on our fiducial mock, we find that this constant is strongly degenerate with $\sigma_{8}$, so it meaningfully blinded our posterior inference of the combined parameter $S_{8}$, which is the only cosmological parameter with strong constraints in this analysis.

We use the measurements of LOWZ clustering and galaxy-galaxy lensing ($w_{p}$ and $\Delta\Sigma$) from Singh et al. (2018). To summarize the data sources: for the galaxy sample, we use the Baryon Oscillation Spectroscopic Survey (BOSS) Data Release 12 (DR12) LOWZ sample (Alam et al., 2015; Reid et al., 2016), which has an effective area of 8,337 deg2 (including NGC and SGC regions, rejecting the small part of the NGC region that had incorrect targeting for the LOWZ sample, and weighting the resulting sky area by observational completeness). We select a subsample of this catalog to include only galaxies in the redshift range $0.16<z<0.36$ and further reject galaxies within regions that do not pass the photometric quality cuts for shape measurements used by Singh et al. (2017a) (originally defined by Reyes et al. 2012), removing an additional 8 per cent of galaxies. The weights applied to each galaxy are the combined systematic, fiber collision, and redshift failure weights developed by the BOSS large-scale structure working group (Ross et al., 2014). We only use scales for which the fiber collision weights provide a correction that is accurate at the $\lesssim 1$ per cent level, above approximately twice the projected scale of fiber collisions at the maximum redshift of the sample (section 6.2).

For the galaxy shear catalog, we use a catalog derived from the Sloan Digital Sky Survey (SDSS) Data Release 8 (DR8) imaging (Aihara et al., 2011) with the photometric calibration method of Padmanabhan et al. (2008). Galaxy shapes were measured with the re-Gaussianization method (Hirata & Seljak, 2003), with the shear response calibrated via image simulations that included the effects of nearby neighbors (Mandelbaum et al., 2012; Mandelbaum & Hyper Suprime-Cam (HSC) Collaboration, 2017). Photometric redshifts of this catalog were measured by Nakajima et al. (2012) and Reyes et al. (2012) with the ZEBRA code (Feldmann et al., 2006) and tested with clustering redshifts by Singh et al. (2019). For a more complete description, we refer the reader to Singh et al. (2019). We allow for lensing systematic uncertainties by including the $A_{\text{lensing}}$ parameter with a 6 per cent prior.

In Figure 9, we show the posterior parameter contours inferred from our fiducial cosmological analysis for the parameters $S_{8}$, $\sigma_{8}$, $\Omega_{m}$, and $A_{\text{lensing}}$. However, we marginalize over all parameters shown in Table 2, including all $w$CDM cosmological parameters (except $N_{\text{eff}}$). We see that $S_{8}$ is degenerate with $A_{\text{lensing}}$, as expected from our tests on mocks, and that this degeneracy substantially increases the uncertainty of our final measurement on $S_{8}$. The parameter $S_{8}$, while almost perfectly decorrelated with $\Omega_{m}$ in our tests on mocks (Figure 8), is not exactly the best-constrained combination of $\sigma_{8}$ and $\Omega_{m}$ for the data, since $\Omega_{m}$ is not fully decorrelated with it. Our fiducial measurement of $S_{8}$ is $0.847\pm 0.037$, a 4.4 per cent measurement. With the results of Planck Collaboration et al. (2018) giving $S_{8}=1.00\pm 0.02$ for their fiducial flat $\Lambda$CDM model, this represents a $3.5\sigma$ tension with Planck data. We show the posterior predictive distribution for $w_{p}$ in Figure 10 and for $\Delta\Sigma$ in Figure 11. While there is very good agreement between the data and the model overall, there may be a hint of tension on scales $\gtrsim 5-10\,h^{-1}\,\text{Mpc}$. In Figure 11, we show a model computed without a point mass term in addition to our fiducial model with a point mass term, illustrating that the data favor a negative value of the point mass term and that this term is important to obtain a model that is in good agreement with the data on small scales.

We additionally compute the posterior resulting from only using scales in the 2-halo regime ($\gtrsim 2\,h^{-1}\,\text{Mpc}$). The marginalized parameter constraints for this analysis are shown in the rightmost column of Table 3. We find that the posterior mean of $S_{8}$ is consistent in this analysis with the posterior mean obtained from our fiducial analysis that includes two-point information from the 1-halo regime, but with a larger, 6 per cent uncertainty, $S_{8}=0.85\pm 0.05$ (quoting the standard deviation instead of the 68 per cent credible interval as used in Table 3). Incorporating the uncertainty in the Planck analysis, this is a $2.6\sigma$ tension between our results and the fiducial Planck cosmological analysis. While there are small ($\sim 1\sigma$) parameter shifts for some of the HOD parameters compared to our fiducial analysis, the agreement of the cosmological parameter values between these analyses suggests that any possible tension between small and large scales does not affect the cosmological parameters of interest.

In Paper I, we predicted that including 1-halo information down to $\sim 0.5\,h^{-1}\,\text{Mpc}$, compared with using only 2-halo information, could improve constraints on $S_{8}$ by a factor of $\sim 2$. This is not achieved in our LOWZ analysis because much of the total uncertainty in both cases comes from the weak lensing systematics parameterized by $A_{\text{lensing}}$. To explicitly verify this point, we analyzed the data on the fiducial set of scales ($0.6\lesssim r_{p}\lesssim 30$) with the parameter $A_{\text{lensing}}$ fixed to 1 and without including the point mass term, yielding $S_{8}=0.803\pm 0.023$, a 2.8 per cent measurement, and additionally with $A_{\text{lensing}}$ fixed to 1 but with the point mass term, yielding $S_{8}=0.831\pm 0.029$, a 3.5 per cent measurement (in contrast to 4.4 per cent precision when both including a point mass term and marginalizing over $A_{\text{lensing}}$). We do not consider these modified analyses as part of our science results, except to illustrate that if the $\Delta\Sigma$ signal could be fully modeled without additional nuisance parameters, then precision in good agreement with our predictions from Paper I can be obtained (where we predicted 2 per cent uncertainty on $S_{8}$ in the case where the only cosmological parameters marginalized over were $\Omega_{m}$ and $\sigma_{8}$). Current weak lensing surveys (e.g. DES Collaboration et al. 2017) aim to reduce photometric redshift zero-point and shear calibration uncertainties to the percent or sub-percent level so that they can take full advantage of their smaller statistical errors.

Figure 12 shows the posterior predictive distribution of the mean halo occupation. The shape of this distribution is fairly well constrained, but the occupation in rare high mass halos is uncertain and many of the individual HOD parameters are poorly constrained (Table 3). These uncertainties may appear surprising when compared to the per cent-level constraints on some halo occupation parameters reported previously in the literature (e.g. Zehavi et al. 2011; Sinha et al. 2018), but it was in fact predicted by our forecasts in Paper I. Some of the additional uncertainty is undoubtedly due to marginalizing over cosmological parameters in addition to HOD parameters, which has been seldom done in previous analyses (and never, to our knowledge, in the context of making halo model predictions by populating $N$-body simulations). The uncertainties in individual parameters also depends on the choice of parameterization because of degeneracies, so a representation like Figure 12 is more informative in terms of both the mean and uncertainties of the inferred halo occupation. A parameterization which more tightly follows the posterior predictive distribution of halo occupation as a function of mass (Figure 12) may be worth exploring in future work.

Finally, we compute the posterior discrepancy distribution for our fiducial analysis on data just as we did previously for the analyses on mock lightcones (Figure 13), following the methodology advocated by Gelman et al. (1996). The $\text{min}\,\Delta\chi^{2}$ discrepancy computed for the data lies near the peak of the expected distribution, computed assuming that the data is drawn from the posterior samples. As before, we compute the minimum $\Delta\chi^{2}$ between the data and the posterior samples to obtain $\text{min}\,\Delta\chi^{2}$ of the data. For the expected distribution of $\text{min}\,\Delta\chi^{2}$, we choose posterior samples at random and add noise according to our fiducial covariance matrix, and then compute the minimum $\Delta\chi^{2}$ between this synthetic observation and the set of posterior samples. Interpreted as a frequentist statistical test, the value of $\text{min}\,\Delta\chi^{2}$ of our data is such that we do not reject the null hypothesis that our data is drawn from the posterior distribution. While this test does not imply anything about the correctness of our model, it does suggest that our model is an adequate description of the data given its statistical precision, a desirable and nontrivial property of any forward modeling framework.

6 Discussion

6.1 Comparison with previous work

Our posterior constraints on $S_{8}$ are in good agreement with those of Singh et al. (2018), who use the same data but using a minimum scale of $1\,h^{-1}\,\text{Mpc}$ and with a parameterized model of the cross-correlation coefficient between galaxies and matter inspired by that obtained in simulations instead of a full forward model of the galaxy clustering and galaxy-galaxy lensing, as used in this work. While they conduct multiple analyses with various estimators and scale cuts, their statistical precision is comparable to ours, as they obtain (for their tightest-constrained value) $S_{8}=0.823\pm 0.035$, a $4.2$ per cent uncertainty. They do not include systematic uncertainties directly in deriving posterior parameter constraints but rather quote an additional systematic uncertainty of 6 per cent in addition to their quoted statistical uncertainty. As a result, it is difficult to directly compare whether our analysis results in more precise constraints, but we note that when excluding systematic uncertainties from our analysis (i.e. removing the point mass term and fixing $A_{\text{lensing}}=1$) we obtain a precision of 2.8 per cent on our posterior value of $S_{8}$, suggesting that including sub-Mpc scales and explicitly modeling both the clustering and lensing signal may provide additional cosmological constraining power in our analysis.

Our results may be expected to be very similar, as our halo occupation models predict similar values of the scale-dependent galaxy-matter cross-correlation coefficient assumed in the parameterization of Singh et al. (2018), given that the galaxy-galaxy lensing signal can be exactly decomposed into a product of the projected galaxy-matter cross-correlation and the projected galaxy bias (Baldauf et al., 2010). Galaxy assembly bias within the subhalo abundance matching framework does not significantly alter (at the percent level) the scale-dependent galaxy-matter cross-correlation coefficient (McEwen & Weinberg, 2016) and thus both approaches are relatively robust to the possible presence of galaxy assembly bias.

Our results also exhibit essentially no tension with the combined galaxy clustering, galaxy-galaxy lensing and cosmic shear analysis of DES Collaboration et al. (2017), as readily seen in Figure 9. This analysis includes all $w$CDM parameters and additionally marginalizes over the effective number of neutrino species $N_{\text{eff}}$ as well of the sum of neutrino masses $\sum m_{\nu}$, whereas we fix $N_{\text{eff}}$ to the standard model value (3.046) and assume massless neutrinos. We also compare with an alternative analysis performed by DES Collaboration et al. (2017) (section VIID) wherein they assume $\Lambda$CDM and fix the sum of neutrino masses $\sum m_{\nu}$ to the minimal value allowed by neutrino oscillation experiments (0.06 eV). We have explicitly verified from their posterior samples that these priors leave the $S_{8}$ (as defined in this work) constraint essentially unchanged compared to their fiducial $w$CDM priors, which marginalize over neutrino mass, and our results are likewise fully compatible with their alternative analysis with fixed neutrino mass.

Our results are also qualitatively similar to those of Leauthaud et al. (2017), who in addition to conducting a qualitative investigation of the galaxy-galaxy lensing signal in Planck-normalized galaxy mocks of the BOSS CMASS sample (with effective redshift $z_{\text{eff}}\approx 0.57$), perform a fitting-function-based (van den Bosch et al., 2013) HOD analysis (without including possible central galaxy miscentering or incompleteness) to obtain a constraint on $S_{8}$ that is lower than that of the fiducial analysis of Planck data at $2-3\,\sigma$ significance (see their Figure 9 and accompanying text). Our HOD model is fully self-consistent and does not require any additional nuisance parameters to describe the ‘halo exclusion’ effect, unlike their analysis, which may contribute to the comparatively greater statistical significance of the tension of our value of $S_{8}$ with that of Planck. Our result also uses the BOSS LOWZ galaxy sample instead of the BOSS CMASS galaxy sample, which may be more robust to galaxy selection effects altering the form of the halo occupation, since the stellar mass completeness of LOWZ is greater over a wider redshift range than that of CMASS (Leauthaud et al., 2016), although a definitive statement would require convincing modeling of the complicated redshift- and color-dependent selection effects.

In Figure 14, we show the posterior predictive distribution for $\Delta\Sigma$ for three separate analysis cases. The first (dashed black line) is the posterior mean derived from our fiducial analysis, with a free point mass parameter and a lensing amplitude parameter $A_{\text{lensing}}$ with a 6 per cent prior uncertainty derived from our systematic error estimates. The second (dashed-dotted blue line) is the posterior mean derived from assuming that the cosmological parameters are fixed to the Planck best-fit $\Lambda$CDM values, but still allowing the lensing amplitude to vary subject to the same prior as before. The third (dotted green line) is the posterior mean derived from assuming that the Planck cosmological parameter values are correct and fixing the lensing amplitude parameter $A_{\text{lensing}}=1$. The first two analyses are in qualitative agreement with the data (red line), but the third analysis is not. This comparison illustrates that the most significant (in terms of fractional amplitude) discrepancy between the data and the model is at large scales when the cosmological parameters are forced to agree with the Planck results. We hypothesize that the contrary conclusion drawn by Lange et al. (2019) may be due to their lack of point mass term in modeling the lensing signal.

In Figure 15, we present an alternative method of characterizing the tension between our results and the fiducial Planck cosmology. We use the discrepancy distance measure described in section 4 to compute the expected distribution of $\Delta\chi^{2}$ assuming that the data is drawn from the posterior of a Planck-normalized model. The $\Delta\chi^{2}$ computed for the data lies at the 97th percentile of the expected distribution. Interpreted as a $p$-value, this corresponds to $p=0.024$ or a $\sim 2\sigma$ tension when interpreted as a single-tailed statistical test. This measure of tension between the data and the Planck cosmology may be more conservative than the standard comparison of the marginalized posterior of $S_{8}$ but it may be more robust to the (somewhat arbitrary) choice of priors on individual parameters and related effects due to the prior volume scaling with the number of nuisance parameters used in the analysis.

6.2 Observational systematics

Since we restrict our analysis to scales $\lesssim 30\,h^{-1}\,\text{Mpc}$, the dominant systematic uncertainty for $w_{p}$ is fiber collisions, which we avoid the need to explicitly model by adopting the standard fiber collision weights given in the LOWZ catalogs (Ross et al., 2014) and only using projected scales greater than twice the fiber collision scale at the maximum redshift of the LOWZ sample ($\gtrsim 0.6\,h^{-1}\,\text{Mpc}$). This is shown to be sufficient for $\lesssim 1$ per cent accuracy of the projected correlation function in tests against mocks with synthetic fiber collisions by Guo et al. (2012) and Yang et al. (2019). We neglect relativistic effects on clustering measurements, which are likely negligible for current datasets.

As discussed in Singh et al. (2018), the dominant source of uncertainty for lensing measurements used here is photo-z calibration ($\sim 5$ per cent), followed by shear calibration ($\sim 2$ per cent), intrinsic alignments (consistent with zero at $\sim 1-2$ per cent precision, according to the analysis of Blazek et al. 2012 on the same dataset used here), and magnification bias (estimated by Singh et al. 2019 to be $\lesssim 1$ per cent). We have not attempted to estimate reduced shear corrections or any additional higher-order lensing effects. We note that fourth-order and higher lensing effects have been shown to be subdominant to shape noise for cosmic shear measurements even at LSST precision (Petri et al., 2017). We have used the multiplicative calibration parameter $A_{\text{lensing}}$ on the amplitude of $\Delta\Sigma$ in order to account for the combined effects of photo-z calibration bias, shear calibration bias, intrinsic alignments, and magnification bias, with a Gaussian prior centered at unity with a 6 per cent standard deviation.

6.3 Baryonic and neutrino effects

Singh et al. (2018) conducted extensive parameter recovery tests by fitting to hydrodynamic simulations with strong AGN feedback and found that baryonic effects only affected $S_{8}$ by at most $\sim 2-3$ per cent, with their fits using a minimum scale of $1\,h^{-1}\,\text{Mpc}$. However, when including shape noise in their analysis covariance matrix, they find essentially no bias in recovered $S_{8}$. We conducted identical tests, but down to our minimum scale of $0.6\,h^{-1}\,\text{Mpc}$, by modifying the ‘observed’ $\Delta\Sigma$ from the fiducial mock lightcone by the ratio of $\Delta\Sigma$ (measured with $\Pi_{\text{min}}=0.01\,h^{-1}\,\text{Mpc}$; eq. 9) with and without baryonic effects as measured in the original Illustris simulation by Singh et al. (2018). We find that when the point mass is allowed to be negative, the effect of baryons on $S_{8}$ is consistent with zero ($\lesssim 0.6$ per cent difference from the true value of $S_{8}$ in our lightcones). Red contours in Figure 8 show the parameter constraints recovered from the mock catalog with this treatment of baryonic effects.

For future surveys with increased statistical precision, the best way to test for baryonic effects is likely to be cross-correlations with probes of hot gas, such as thermal Sunyaev-Zeldovich (tSZ) maps from Planck and tSZ and kinetic Sunyaev-Zeldovich (kSZ) maps from Stage IV cosmic microwave background experiments (CMB-S4; Abazajian et al. 2016). With additional model flexibility, combining this cross-correlation with the cosmic shear correlation functions should enable self-calibration of baryonic effects on the matter distribution, or otherwise provide a convincing null test of their absence.

Clustering effects of massive neutrinos are not included in our emulator, as it has only recently become feasible to routinely simulate the effects of massive neutrinos on nonlinear matter clustering thanks to algorithmic advances (Bird et al., 2018; Banerjee et al., 2018). The tests of Leauthaud et al. (2017) at slightly higher redshift ($z\sim 0.57$) suggest that massive neutrinos may affect the matter clustering at the $5-10$ per cent level on the scales used in this work, with an effect size depending on the (currently unknown) value of the sum of the neutrino masses. The best current limits are $\lesssim 0.1$ eV, assuming $\Lambda$CDM cosmology, but these relax to $\lesssim 0.3$ eV within $w$CDM cosmology (Alam et al., 2017; Choudhury & Choubey, 2018). The effect of neutrino mass is strongly, but not exactly, degenerate with $S_{8}$ (Ichiki & Takada, 2012). Future emulators would be well advised to include such effects, especially given the strong discrepancy we infer between the low-redshift amplitude of matter clustering and that of Planck CMB measurements in the absence of massive neutrino effects.

6.4 Halo model uncertainties

Central galaxies may not be at the centers of their halos (e.g., Ho et al. 2009), and thus the lensing signal may be mismatched with the clustering signal for this reason. A model for miscentering and galaxy sample incompleteness was included in a fitting-function-based HOD cosmological analysis of CMASS galaxy clustering and galaxy-galaxy lensing by More et al. (2015), with posterior constraints indicating a $\sim 2\sigma$ detection of a miscentered population of central galaxies. Since the LOWZ sample also consists of luminous red galaxies, it may be worthwhile to parameterize the miscentering effects in LOWZ galaxies in future work (but see Lange et al. 2019 for a counterargument regarding the plausibility of such effects).

In Paper I we introduced a parameter $Q_{\text{env}}$ to represent the possible effects of galaxy assembly bias (e.g. Zentner et al. 2016). We found that $Q_{\text{env}}$ was not needed to obtain unbiased cosmological results from the SHAM lightcone (consistent with the result of McEwen & Weinberg 2016), and we therefore did not include it in our analysis. The real universe may exhibit galaxy assembly bias in a form unlike that of SHAM, so it may be desirable to include an explicit parameterization of assembly bias in future cosmological analyses. Several alternate forms of assembly bias were explored by Yuan et al. (2019), who found that none of the candidate forms of assembly bias could explain the lensing amplitude of BOSS CMASS galaxies in a $\Lambda$CDM cosmology with Planck-compatible parameters. Taken together, these results disfavor assembly bias as a solution to the observed lensing tensions.

6.5 Possible empirical tests of robustness

Confidence in small scale measurements of cosmological parameters as described here will require either Bayesian model averaging over a wide range of phenomenological models of galaxy bias (i.e., for a set of models which are considered equally likely a priori, weighting the posteriors for the parameters in common between all models under consideration by the evidence integral of each model; see Marshall et al. 2003; Liddle et al. 2006; Parkinson & Liddle 2010; Vardanyan et al. 2011 for applications to cosmological parameters) or will require the development of null tests on data, not just simulations.999As justification, we invoke Cromwell’s dictum (as named by Bayesian statisticians): “I beseech you, in the bowels of Christ, consider it possible that you are mistaken” (quoted in Rasmussen & Williams 2006). While it is not possible to compute all models, we hope that our study will encourage model comparison and model averaging to be attempted with a finite number of models on LOWZ galaxy-galaxy lensing.101010However, we note that any such comparison is subjective, since the evidence of a given model depends sensitively on the prior adopted on the parameters of that model, much more so than the posterior distribution itself (e.g., Gelman et al. 2017). There exist non-Bayesian model averaging methods which are more robust to prior choices, e.g. using hyperparameter weights (Lahav et al., 2000; Trotta, 2008).

When one is limited to evaluating a single model, or when it is unknown whether the correct model lies within the set of models under consideration, a reasonable test of model robustness is to fit multiple galaxy samples, showing one can recover the same cosmological parameters and also fit the cross-correlation functions between the galaxy samples without significant discrepancies between model and data. This would likely require luminosity- and color-dependent modeling of the galaxy population. (In principle, one would also need to include subsamples split on all of the galaxy properties used for sample selection.) This is an ambitious program of research, which will no doubt require the work of many people.111111For usage of this phrase in an entirely different context, see Teller (1955).

7 Summary and Outlook

We have extended the methodology of Wibking et al. (2019) (Paper I) to enable emulation of galaxy clustering and GGL for $w$CDM cosmologies, sampling the cosmological parameter space allowed by WMAP and Planck CMB measurements and HOD parameters characteristic of the BOSS LOWZ sample of massive galaxies at $z=0.16-0.36$. We use the Garrison et al. (2017) suite of 40 $w$CDM $N$-body simulations with matched Fourier phases, supplemented by 20 simulations of a fiducial $\Lambda$CDM cosmology with varying phases which we use for covariance matrix calculations and for a sample variance correction to the mean model predictions. We construct a Gaussian process emulator for the GGL observable $\Delta\Sigma(r_{p})$ and for the ratio of the projected galaxy correlation function $w_{p}(r_{p})$ to the prediction of an analytic halo model calculation. This ratio changes much more slowly with parameters than $w_{p}(r_{p})$ itself, enabling us to achieve percent-level emulator accuracy over a wider range of parameters. We train the emulator by maximizing the leave-one-out pseudo-likelihood (Appendix C).

We test our approach on the light-cone mock catalogue of the BOSS LOWZ sample constructed by Singh et al. (2018). This catalogue uses subhalo abundance matching rather than an HOD prescription to populate the simulated dark matter distribution with galaxies, so it tests the ability of an HOD-based emulator method to derive unbiased results from a catalogue that incorporates a different model of nonlinear galaxy bias as well as evolution over the redshift range of the LOWZ sample. When fitting the mock catalogue we adjust the covariance matrix in a way that effectively marginalizes over a central point mass (MacCrann et al., 2019), i.e., an additive $1/r_{p}^{2}$ contribution to $\Delta\Sigma(r_{p})$ with arbitrary normalization. This parameter accounts for finite resolution in our N-body simulations, and it can also account for baryonic physics effects and for the possibility that central galaxies do not lie at the potential minima of their parent halos. Our HOD parameterization also includes a novel parameter $R_{\text{rescale}}$ that multiplies all halo virial radii by a constant factor; this parameter reduces sensitivity to the halo definition.

With fiducial parameter choices in the mock catalogue, we recover the true input value of the cosmological parameter combination $S_{8}\equiv(\sigma_{8}/0.8228)(\Omega_{m}/0.3107)^{0.6}$ to within 1.7 per cent, about half of the estimated $1\sigma$ statistical uncertainty. We constructed alternative mock catalogues with higher or lower satellite fractions, central galaxy incompleteness, or baryonic effects on $\Delta\Sigma(r_{p})$, and we again found unbiased recovery of $S_{8}$. Because we have only a single mock LOWZ volume, we cannot average over many realizations to test for biases that are small compared to the statistical error.

After completing all mock catalogue tests and finalizing modeling choices, we applied our method to the Singh et al. (2018) measurements of LOWZ $w_{p}$ and $\Delta\Sigma$. We find $S_{8}=0.847\pm 0.037$, a $3.5\sigma$ tension with the Planck-normalized value of $1.00\pm 0.02$ for a $\Lambda$CDM cosmology. A significant fraction of the 4.4 per cent error budget comes from systematic uncertainty in the weak lensing measurements, modeled by a multiplicative parameter $A_{\text{lensing}}$ with a 6 per cent Gaussian prior to represent estimated uncertainties in photometric redshift zero points and shear calibration. With $A_{\text{lensing}}$ fixed to one, our posterior weighted estimate would be $S_{8}=0.831\pm 0.029$, a $3.5$ per cent statistical error. The point mass marginalization also contributes signficantly to the error budget. If we fixed both $A_{\text{lensing}}$ and the point mass parameter, our statistical error on $S_{8}$ would improve to 2.8 per cent. If we instead force $S_{8}=1$ and adopt our fiducial prior on $A_{\text{lensing}}$, we can fit the clustering and GGL data with $A_{\text{lensing}}\approx 0.85$, implying that reproducing the Planck $S_{8}$ normalization with our modeling choices requires a lensing systematic error that is about 2.5 times higher than estimated by Singh et al. (2018).

Our results are in excellent agreement with those of Singh et al. (2018), who model the same data with a different method, adopting a parameterized galaxy-matter cross-correlation coefficient constrained by their mock catalogues. They are in qualitative agreement with the findings of Leauthaud et al. (2017), who analyzed the higher redshift BOSS CMASS sample using deeper weak lensing data over much smaller area. Our normalization of $S_{8}$ is consistent with (but lower than) that derived from the DES Y1 $3\times 2$pt analysis (DES Collaboration et al., 2017). Our measurement joins others (e.g., Leauthaud et al. 2017; Singh et al. 2018; Mandelbaum et al. 2013; Hildebrandt et al. 2017) that suggest an amplitude of low redshift matter clustering significantly below that of a Planck-normalized $\Lambda$CDM model.

In future work, several aspects of our methodology could be improved:

• •

the calculation of the redshift-space distortion correction to $w_{p}$ could be made more efficient and properly cosmology-dependent via the use of the FFTLOG algorithm (Hamilton, 2000) to compute the integrals required in the Kaiser (1987) model;

• •

the sample variance of the emulator could be further reduced via $N$-body simulations initialized with ‘fixed and paired’ phases, as advocated by Angulo & Pontzen (2016);

• •

for a fixed number of samples, the emulator accuracy might be improved with a sampling design that incorporates both clustered and volume-filling subsamples of the parameter space (Zhu & Zhang, 2006; Zimmerman, 2006);

• •

future datasets may be precise enough as to require explicit redshift-dependent modeling, which may be straightforwardly incorporated via an additional label dimension of the training data;

• •

halo definition may be marginalized over in a fully self-consistent way by introducing the spherical overdensity threshold as an emulator parameter and sampling from this parameter space jointly with the cosmological parameters;

• •

the robustness of the cosmological results may be empirically tested by analyzing disjoint subsamples split by galaxy properties used in target selection, such as luminosity and color;

• •

the effect of massive neutrinos on nonlinear matter clustering may be significant at the $\sim 5$ per cent level in $\Delta\Sigma$ (Leauthaud et al., 2017) and should be included in $N$-body simulations with the benefit of recently-developed numerical techniques (Bird et al., 2018; Banerjee et al., 2018); and

• •

although fourth-order and higher effects are likely negligible for any planned weak lensing experiment (Petri et al. 2017), weak lensing effects arising at third order in the gravitational potential such as reduced shear and magnification bias corrections (e.g. Krause & Hirata 2010) may need to be explicitly modeled in future analyses via emulation of the three-point function.

In the near future, expanded data sets from DES, KiDS, and HSC will allow this lensing amplitude tension to be tested at higher precision using both cosmic shear and galaxy clustering + GGL. The principal challenge will be controlling weak lensing systematics at the necessary level, though further tests of the modeling methods will also be needed as the statistical precision improves. GGL and cosmic shear have different quantitative response to biases in photometric redshifts, shear calibration, and galaxy intrinsic alignments, so they provide useful cross-checks as well as increased statistical precision in combination. GGL + clustering analyses of galaxy samples with distinct clustering properties and at multiple redshifts allow further tests for the robustness of cosmological parameter constraints. If future observations indeed show the need for physics beyond that of GR+$\Lambda$CDM, then GGL+clustering analysis will be a powerful tool for measuring its redshift, scale, and environment dependence.

Acknowledgements

We thank Niall MacCrann and the DES Collaboration for providing the posterior samples for the DES ‘fixed-neutrino’ analysis used in section 6.1.

BDW thanks Chris Hirata for helpful discussions of covariance matrices and many other topics, Joe McEwen and Xiao Fang for explaining the FFTLOG algorithm and its application to one-loop standard perturbation theory, Martin White for suggesting empirical tests of the robustness of small-scale measurements of cosmological parameters, and Andreas Berlind for first suggesting to him the possibility of using small-scale galaxy clustering to measure cosmological parameters.

BDW is supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE-1343012. ANS is supported by the Department of Energy Computational Science Graduate Fellowship Program of the Office of Science and National Nuclear Security Administration in the Department of Energy under contract DE-FG02-97ER25308. BDW, ANS, and DHW are supported in part by NSF grant AST-1516997. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. DJE is supported by U.S. Department of Energy grant DE-SC0013718 and as a Simons Foundation Investigator. LG is supported by the Simons Foundation.

Simulations were analyzed in part on computational resources of the Ohio Supercomputer Center (Center, 1987), with resources supported in part by the Center for Cosmology and AstroParticle Physics at the Ohio State University.

Software: matplotlib (Hunter, 2007), GNU Scientific Library (Galassi et al., 2009), Corrfunc (Sinha & Garrison, 2017). BDW especially thanks the authors of the Corrfunc pair-counting code, without which this project (and his Ph.D. thesis) would not have been feasible.

This research has made use of NASA’s Astrophysics Data System.

Appendix A Analytic halo model

Although we use a numerical halo model for our predictions, we need an analytic halo model in order to reduce the dynamic range of the emulated quantity and increase our emulation accuracy. For this purpose, we construct a simple analytic halo model that does not attempt to model halo exclusion or residual RSD effects. Although we only use the analytic halo model in order to compute the ratio $w_{\text{p,sim}}/w_{\text{p,analytic}}$ in this work, we include the equations for $\Delta\Sigma$ for completeness. For notational clarity, we drop the ‘analytic’ subscripts used elsewhere in this work in order to indicate the use of the halo model described in this Appendix.

We write the configuration-space correlation functions as a the sum of a 1-halo term and a 2-halo term

[TABLE]

where the 1-halo term $\xi_{\text{1h}}$ arises from galaxy pairs (or galaxy-matter ‘pairs’) within a given halo and the 2-halo term $\xi_{\text{2h}}$ arises from galaxy pairs (or galaxy-matter ‘pairs’) between two distinct halos.

A.1 1-halo term $\xi_{\text{1h}}$

In writing the 1-halo term, we follow the real-space formulation of Zheng & Weinberg (2007). This allows for faster computation of the 1-halo term than its Fourier-space formulations (e.g., van den Bosch et al. 2013).

A.1.1 Galaxy autocorrelation

The 1-halo galaxy autocorrelation is given by the normalized differential pair counts of galaxies within halos (Berlind & Weinberg, 2002), assuming Poisson satellite counts:

[TABLE]

where

[TABLE]

and

[TABLE]

For an NFW profile, the differential (w.r.t dimensionless radius $r/R_{\text{vir}}$) pair count functions $I^{\prime}$ and $F^{\prime}$ are (Sheth et al., 2001):

[TABLE]

and

[TABLE]

where $s=xc$ and $a=1/c$, and the galaxy number density is

[TABLE]

Note that here $\langle N_{\text{sat}}\rangle$ is the fully marginalized satellite halo occupation (i.e., not conditioned on having a central in a given halo).

The virial radius is defined as

[TABLE]

where $\rho_{m}=\rho_{\text{crit}}\Omega_{m}$ ($\rho_{\text{crit}}$ is the critical density of the universe) and we choose $\Delta_{\text{SO}}=200$ for consistency with the halo mass-concentration relation described below (although this choice is one only of convenience).

We use the halo mass-concentration relation of Correa et al. (2015)

[TABLE]

where

[TABLE]

and the mass function is the fitting formula of Tinker et al. (2008):

[TABLE]

where

[TABLE]

and

[TABLE]

A.1.2 Galaxy-mass cross-correlation

For the 1-halo term, we have

[TABLE]

We write the convolution of an NFW profile with another NFW profile of differing concentration $F^{\prime}(x,c_{1},c_{2})$ as (Appendix A, Zheng & Weinberg 2007), with $s=2x$:

[TABLE]

where

[TABLE]

where $A_{\star}$ is approximately

[TABLE]

where

[TABLE]

A.2 2-halo term $\xi_{\text{2h}}$

For the linear matter power spectrum $P_{\text{lin}}(k)$, we use the ‘no-wiggles’ fitting formula of Eisenstein & Hu (1998), which is substantially cheaper to compute than using Boltzmann codes, which makes our posterior inferences faster.

We use the halo mass-bias relation of Tinker et al. (2010):

[TABLE]

where $\delta_{c}=1.686$. We obtain the overall large-scale bias of the galaxy sample as

[TABLE]

A.2.1 Galaxy autocorrelation

For the galaxy autocorrelation between distinct halos, we have

[TABLE]

A.2.2 Galaxy-matter cross-correlation

For the galaxy-matter cross-correlation between distinct halos, we have

[TABLE]

where we further assume that $r_{\text{gm}}=1$ on all 2-halo scales.

Appendix B Analytic covariance matrices

Rewriting equation A2 of Krause & Eifler (2017) for the three-dimensional power spectrum, we have:

[TABLE]

where $\delta_{ij}$ refers to the Kronecker delta function and $N^{i}_{A}$ is the appropriate noise term for probe $A$ in redshift bin $i$ (e.g., for the 3D galaxy density field, this is $1/n^{i}_{g}$, where $n^{i}_{g}$ is the number density of galaxies in redshift bin $i$).

Specializing to the case of a single redshift bin (i.e., $i=j=k=l$), we have

[TABLE]

where the sum of products of power spectra (in brackets) is equal to the four-point function $\langle ABCD\rangle$ when $A$, $B$, $C$, and $D$, are Gaussian fields (compare with Eq. 4 of Cooray & Hu 2001).

For projected two-point functions with pairwise line-of-sight weight functions $W_{AB}$ and $W_{CD}$, we have

[TABLE]

where $\Pi$ is the relative line-of-sight distance between the points of a pair used in computing the two point statistic of interest.

Ignoring the finite-size correlation function bin width and survey boundary effects, for scalar projected correlation functions $w_{AB}(r_{p})$ and $w_{CD}(r_{p})$, we have

[TABLE]

where $i,j$ here refer to the bin indices $r_{p,i}$ of the projected correlation functions. When applying this expression to compute a covariance matrix, we stress that it is necessary to average the integrand over the bin width before performing the outer integrals. Otherwise, integrals of this form diverge whenever $N_{A}$ or $N_{B}$ are nonzero, due to the identity (Jackson, 1975; Gradshteyn & Ryzhik, 2007):

[TABLE]

Physically, this divergence represents the variance tending toward infinity as the number of galaxy pairs contained within an infinitesimal radial bin goes to zero. The required bin-averaging of the integrand can be carried out by replacing each Bessel function $J_{\nu}(kr_{i})$ with the bin-averaged Bessel function $\bar{J}_{\nu}(kr_{i})$, where the average is taken over an annulus of inner radius $r_{i}$ and outer radius $r_{i+1}$.

B.1 Covariance of $w_{p,gg}$

For the projected two point function $w_{p,gg}(r_{p})$, the pairwise line of sight weight function is

[TABLE]

(compare with equation A21 of Singh et al. 2017b) and the integral over the weight function is

[TABLE]

Therefore eq. B becomes

[TABLE]

B.2 Covariance of $\gamma_{t}$

For galaxy-galaxy lensing, after modifying the Fourier transform to account for the fact that the $\gamma_{t}$ is a component of a spin-2 tensor, eq. B reduces to

[TABLE]

where

[TABLE]

[TABLE]

and $\Pi_{\text{lens}}$ is the effective line-of-sight depth of the squared lensing weight function, given by

[TABLE]

B.3 Covariance between $\gamma_{t}(r_{p})$ and $w_{gg}(r_{p})$

For the cross-probe covariance between the galaxy-galaxy lensing and the galaxy-galaxy projected correlation function, we have

[TABLE]

B.4 Projected radius of a survey

For various computations related to the covariance matrix (e.g., the effective source density $\Sigma_{s}$), it is necessary to compute the effective projected radius. For this purpose, we assume a survey at a single redshift $z$ with spherical cap geometry in a flat universe, where the effective survey radius $R_{s}$ is such that the area of a flat circle with radius $R_{s}$ is equivalent to that of a spherical cap survey. Then we have

[TABLE]

where $f_{\text{sky}}$ is the fraction of the sky covered by the survey and $\chi$ is the comoving radial distance to the effective survey redshift $z$. Then we have

[TABLE]

For a survey of 9736 sq. deg. at an effective redshift $z=0.3$, $R_{s}\approx 812.9\,h^{-1}\,\text{Mpc}$.

Appendix C Gaussian process implementation

The Gaussian process predictor for the expected value of scalar-valued process (i.e., function) $y(\mathbf{x_{\star}})$ is

[TABLE]

where the coefficients $\alpha_{i}$ are computed by

[TABLE]

where $K_{ij}=k(\mathbf{x}_{i},\mathbf{x}_{j})$, $\mathbf{y}$ is the vector of training data observations, $k$ is the kernel (i.e., covariance) function, and the indices $i$,$j$ run over the training data of size $N$. This is simply a Wiener filter applied to the training data (Rybicki & Press, 1992). As such, a key assumption of this method is that noise of the process that we seek to predict is Gaussian. In our application, the training data consists of the values of $w_{\text{p,sim}}/w_{\text{p,analytic}}$ (or $\Delta\Sigma$) in a given radial bin and the parameter values at which these were computed. We estimate the noise $\sigma_{ii}^{2}$ of each datapoint via the sample variance estimator of the 20 stochastic HOD realizations computed at each value of the parameters.

The surge in recent applications is due to the realization that such a smoothing filter can be applied to arbitrary machine learning (i.e., curve fitting) problems by marginalizing over the unknown function $k$ (O’Hagan & Kingman, 1978; Sacks et al., 1989; Rasmussen & Williams, 2006). Due to computational expense, and as is standard practice, one can find a single covariance function $k$ (over some well-defined class of functions) which maximizes the marginal likelihood of the Gaussian process.

However, instead of maximizing the marginal likelihood itself, we instead maximize a closely-related function, namely, the leave-one-out pseudo-likelihood

[TABLE]

where $N$ is the number of training points, $||\text{LOOE}_{i}||$ is the leave-one-out cross-validation error on training point $i$, $\sigma_{i}^{2}$ is the noise variance of on training point $i$, and $\lambda$ is a free parameter that rescales all of the variances. This function is obtained by taking the product over $i=1,2,\dots,N$ of the conditional likelihoods for training point $y_{i}$ given that all training points except $y_{i}$ are known, which is precisely the situation when performing leave-one-out cross-validation. This pseudo-likelihood has the useful property that it is independent of the form of the kernel function $k$, unlike the true marginal likelihood of the Gaussian process (which depends on the log determinant of $K$), which ensures that the maximum of this function attempts to minimize the prediction error regardless of the complexity of the kernel function (Rasmussen & Williams, 2006).

In order to efficiently compute the leave-one-out predictions $\hat{y}_{(i)}$, we use the identity for the inverse of a $2\times 2$ block matrix (e.g., Press et al. 1992), assuming without loss of generality that $i=N$. In this case, this identity reduces to the covariance matrix $K$ partitioned into sub-blocks, with $A$ an $(N-1)\times(N-1)$ matrix, $b$ an $(N-1)$ column vector, and $c$ a scalar:

[TABLE]

where

[TABLE]

With repeated application of the Sherman-Morrison-Woodbury formula, we then solve for $A^{-1}$ in terms of $\tilde{b}$, $b$ and $\tilde{A}$ to obtain

[TABLE]

Recognizing that $\hat{y}_{(i)}=b^{T}A^{-1}\mathbf{y}_{(i)}$, we obtain the leave-one-out prediction

[TABLE]

which can be rewritten in terms of the true value $y_{i}$, the matrix elements of $K^{-1}$, and the matrix elements of $K^{-1}\mathbf{y}$ as (Rasmussen & Williams 2006, Eq. 5.12)

[TABLE]

The leave-one-out error on point $i$ is therefore

[TABLE]

For our kernel function $k$, we adopt the squared exponential family of kernel functions

[TABLE]

where $\mathbf{d}=\mathbf{x}_{i}-\mathbf{x}_{j}$, $\Lambda$ is a diagonal matrix of hyperparameters over which to optimize, and $\sigma^{2}_{\text{signal}}$ and $\sigma^{2}_{\text{mean}}$ are hyperparameters over which to optimize. This kernel is a convenient choice of an infinitely differentiable and translation invariant function and is perhaps the most common kernel used in machine learning applications (Rasmussen & Williams, 2006).

To optimize the hyperparameters ($\Lambda$, $\sigma^{2}_{\text{signal}}$, $\sigma^{2}_{\text{mean}}$, $\lambda$), we use the BOBYQA derivative-free quadratic surface optimization method (Powell, 2009) included in the nlopt software package (Johnson, 2018). We find that the performance of the emulator is very sensitive to the initial guess of the hyperparameters, and that an initial guess for the hyperparameters which implies a high signal-to-noise of the training data is crucial in order to avoid being trapped in a sub-optimal local maximum of eq. 100. In practice, we accomplish this by de-dimensionalizing all inputs and outputs of the emulator, rescaling the range of each input parameter $x_{i}$ to $[0,1]$ and rescaling the training data outputs $y_{i}$ so that $y_{i}$ has mean zero and variance unity. Then we choose an initial guess for $\lambda\equiv 0.01$ and $\sigma_{\text{signal}}\equiv 1$, which implies that the training data have a signal-to-noise $S/N$ greater than implied by the input uncertainties by a factor $\lambda^{-1/2}=10$. For our training data, even allowing $\lambda$ to vary, this is sufficent to ensure that the optimization finds a relatively high signal-to-noise solution rather than running away toward a low signal-to-noise solution with low predictive power.

Appendix D Projection integrals with finite bin size

In order to take advantage of statistical isotropy to increase the signal-to-noise of our simulation training data, we directly measure the real-space (cross-)correlation functions, and then transform to projected quantities. A slight complication arises when computing these from simulations rather than from theory since we must tabulate these quantities as bin-averaged correlation functions when counting pairs of particles in simulations, but the projection integrals are defined in terms of the raw correlation functions.

We compute the projection integrals from the tabulated correlation functions as follows. Starting from the Abel integral transform

[TABLE]

we assume $\xi_{\text{xy}}$ is tabulated such that correlation function values are exact at midpoints of the tabulated bins $r_{i}$ (a more accurate approach would be to compute the mean pair-weighted separation within each bin). We then analytically integrate, using piecewise linear elements on the interval between the midpoints of adjacent tabulated bins [$r_{i,-}$, $r_{i,+}$], to obtain a second-order accurate (in adjacent bin separation $\Delta r$) sum for $w_{p}(r_{p})$:

[TABLE]

where

[TABLE]

and

[TABLE]