Cosmic Shear Cosmology Beyond 2-Point Statistics: A Combined Peak Count and Correlation Function Analysis of DES-Y1

Mon. Not. R. Astron. Soc. 000, 1–26 (2020) Printed 7 December 2020 (MN LATEX style file v2.
2)
Cosmic Shear Cosmology Beyond 2-Point Statistics: A Combined Peak

Count and Correlation Function Analysis of DES-Y1
Joachim Harnois-Déraps1,2? , Nicolas Martinet3 , Tiago Castro4,5,6,7 , Klaus Dolag8,9 ,

Benjamin
1
Giblin2 , Catherine Heymans 2,10 , Hendrik Hildebrandt10 & Qianli Xia2
Astrophysics Research Institute, Liverpool John Moores University, 146 Brownlow Hill, Liverpool, L3 5RF, UK
2 ScottishUniversities Physics Alliance, Institute for Astronomy, University of Edinburgh, Blackford Hill, Scotland, UK
3 Aix-Marseille Univ, CNRS, CNES, LAM, Marseille, France
4 Dipartimento di Fisica, Sezione di Astronomia, Università di Trieste, Via Tiepolo 11, I-34143 Trieste, Italy
arXiv:2012.02777v1 [astro-ph.CO] 4 Dec 2020
5 INAF - Osservatorio Astronomico di Trieste, via Tiepolo 11, I-34131 Trieste, Ital
6 IFPU - Institute for Fundamental Physics of the Universe, via Beirut 2, 34151, Trieste, Italy
7 INFN - Sezione di Trieste, I-34100 Trieste, Italy
8 University Observatory Munich, Scheinerstr. 1, 81679 Munich, Germany
9 Max-Planck-Institut fur Astrophysik, Karl-Schwarzschild Strasse 1, 85748 Garching, Germany
10 Ruhr-University Bochum, Faculty of Physics and Astronomy, Astronomical Institute (AIRUB), German Centre for Cosmological Lensing, 44780 Bochum, Germany
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We constrain cosmological parameters from a joint cosmic shear analysis of peak-counts and the
two-point shear correlation function, as measured √ from the Dark Energy Survey (DES-Y1). We
find the structure growth parameter S 8 ≡ σ8 Ωm /0.3 = 0.766+0.033 −0.038 , which at 4.8% precision,
provides one of the tightest constraints on S 8 from the DES-Y1 weak lensing data. In our simulation-
based method we determine the expected DES-Y1 peak-count signal for a range of cosmologies
sampled in four wCDM parameters (Ωm , σ8 , h, w0 ). We also determine the joint peak-correlation
function covariance matrix with over 1000 realisations at our fiducial cosmology. With mock DES-
Y1 data we calibrate the impact of photometric redshift and shear calibration uncertainty on the
peak-count, marginalising over these uncertainties in our cosmological analysis. Using dedicated
training samples we show that our measurements are unaffected by mass resolution limits in the
simulation, and that our constraints are robust against uncertainty in the effect of baryon feedback.
Accurate modelling for the impact of intrinsic alignments on the tomographic peak-count remains
a challenge, currently limiting our exploitation of cross-correlated peak counts between high and
low redshift bins. We demonstrate that once calibrated, a fully tomographic joint peak-count and
correlation function analysis has the potential to reach a 3% precision on S 8 for DES-Y1. Our
methodology can be adopted to model any statistic that is sensitive to the non-Gaussian information
encoded in the shear field. In order to accelerate the development of these beyond-two-point cosmic
shear studies, our simulations are made available to the community, on request.
Key words: Gravitational lensing: weak – Methods: numerical – Cosmology: dark matter, dark
energy & large-scale structure of Universe
1 INTRODUCTION dedicated Stage-III weak lensing experiments, namely the Kilo De-
gree Survey1 , the Dark Energy Survey2 and the Hyper Suprime Cam-
Over the last decade, weak gravitational lensing has emerged as one
era Survey3 , were launched and aimed at constraining properties of
of the most promising techniques to investigate the properties of our
dark matter to within a few percent. These are now well advanced or
Universe on cosmic scales. Based on the analysis of small distortions
have recently completed their data acquisition, and the community is
between the shapes of millions of galaxies, weak lensing by large scale
preparing for the next generation of Stage IV experiments, notably the
structures, or cosmic shear, can directly probe the total projected mass
distribution between the observer and the source galaxies, as well as
place tight constraints on a number of other cosmological parameters
(for recent reviews of weak lensing as a cosmic probe, see Kilbinger
2015). Following the success of the Canada-France-Hawaii Telescope
Lensing Survey (Heymans et al. 2012; Erben et al. 2013), a series of
1 KiDS:kids.strw.leidenuniv.nl
2 DES:www.darkenergysurvey.org
? E-mail: jharno@roe.ac.uk 3 HSC:www.naoj.org/Projects/HSC/
© 2020 RAS
2 J. Harnois-Déraps, N. Martinet et al.
Rubin observatory4 , and the Euclid5 and Nancy Grace Roman6 space M18 hereafter) and of the DES Science Verification data (Kacprzak
telescopes. et al. 2016, K16 hereafter) calibrated their signal on a suite of full N-
The central approach adopted by these surveys for constraining body simulations spanning the [Ωm − σ8 ] plane described in Dietrich
cosmology is based on two-point statistics – mostly either in the form & Hartlap (2010). The accuracy of this suite has however been later
of correlation functions (e.g. Kilbinger et al. 2013; Troxel et al. 2018; shown to be only ∼10% (Giblin et al. 2018). Significant improvements
Hamana et al. 2020; Asgari et al. 2020a) or its Fourier equivalent, the on the simulation side are therefore critical for the new generation of
power spectrum, estimated using pseudo-C` (Hikage et al. 2019), band data analyses based on non-Gaussian statistics.
powers (Becker et al. 2016; van Uitert et al. 2018; Joachimi et al. 2020) This paper aims to address this issue: we present a cosmological
and quadratic estimators (Köhlinger et al. 2017). By definition, these re-analysis of the DES-Y1 cosmic shear data (Abbott et al. 2018b),
two-point functions can potentially capture all possible cosmological exploiting a novel simulation-based cosmology inference pipeline cal-
information contained in a linear, Gaussian density field, and are thus ibrated on state-of-the-art suites of N-body runs that are specifically
highly efficient at analysing large scale structure data. They have been designed to analyse current weak lensing data beyond two-point statis-
thoroughly studied in terms of signal modelling (Kilbinger et al. 2017), tics. In this work, the incarnation of our pipeline is tailored for the peak
measurement (Jarvis et al. 2004; Schneider et al. 2002; Alonso et al. count analysis of the DES-Y1 survey, however it is straightforward to
2019) and systematics (Mandelbaum 2018). extend it to alternative non-Gaussian probes. Our pipeline first cali-
With the improved accuracy and precision provided by current brates the cosmological dependence of arbitrary non-Gaussian mea-
and upcoming surveys, it becomes increasingly appealing to probe surements with the cosmo-SLICS (Harnois-Déraps et al. 2019), a seg-
small angular scales, where the signal is the strongest. In doing so, ment of the Scinet LIght-Cone Simulations suite that samples Ωm , σ8 ,
the measurements are intrinsically affected by the non-Gaussian nature w0 and h (the Hubble reduced parameter). We next estimate the co-
of the matter density field, and it is natural to seek analysis techniques variance from a suite of fully independent N-body runs extracted from
that can extract the additional cosmological information that two-point the main SLICS sample7 (Harnois-Déraps et al. 2018). We further use
functions fail to capture. A variety of alternative methods have been the cosmo-SLICS to generate systematics-infused control samples that
applied to lensing data with this in mind, including three-point func- we use to model the impact of photometric redshift and shear calibra-
tions (Fu et al. 2014), Minkowski functionals and lensing moments tion uncertainty. We study the impact of galaxy intrinsic alignment
(Petri et al. 2015), peak count statistics (Liu et al. 2015b,a; Kacprzak with dedicated mock data in which the ellipticities of central galaxies
et al. 2016; Martinet et al. 2018; Shan et al. 2018), density split statis- are aligned (or not) with the shape of their host dark matter haloes,
tics (Gruen et al. 2018), clipping of the shear field (Giblin et al. 2018) following the in-painting prescription of Heymans et al. (2006) (see
and convolutional neural networks (Fluri et al. 2019). Other promising also Joachimi et al. 2013b, for a more recent application). We finally
techniques are also being developed, notably the scattering transform use a suite of high-resolution simulations (SLICS-HR, presented in
(Cheng et al. 2020), persistent homology (Heydenreich et al. 2020a), Harnois-Déraps & van Waerbeke 2015) to investigate the impact of
lensing skew-spectrum (Munshi et al. 2020), lensing minimas (Coul- mass-resolution on the non-Gaussian statistics, and full hydrodynami-
ton et al. 2020) and moments of the lensing mass maps (van Waerbeke cal simulation light-cones from the Magneticum Pathfinder8 to assess
et al. 2013; Gatti et al. 2020). the effect of baryon feedback. All of the above are fully integrated with
While existing non-Gaussian data analyses revealed a constrain- the cosmoSIS cosmological inference pipeline (Zuntz et al. 2015) and
ing power comparable to that of the two-point functions, it is expected therefore interfaces naturally with the two-point statistics likelihood,
that the gain will drastically increase with the statistical precision of enabling joint analyses with the fiducial DES-Y1 cosmic shear cor-
P
the data. For example, constraints on the sum of neutrino mass ( mν ), relation function measurements presented in Troxel et al. (2018, T18
on the matter density (Ωm ) and on the amplitude of the primordial hereafter), with the 3 × 2pts analysis presented in DES Collaboration
power spectrum (As ), in a tomographic peak count analysis of LSST, et al. (2017), or any other analysis implemented within cosmoSIS.
are forecasted to improve by 40%, 39%, and 36% respectively, com- The current document is structured as follow: In Sec. 2 we present
pared to a power-spectrum analysis of the same data (Li et al. 2019). the data and the simulation suites on which our pipeline is built; Sec.
Upcoming measurements of the dark energy equation of state (w0 ) will 3 describes the theoretical background, the weak lensing observables
also benefit from these methods, with a forecasted factor of three im- and the analysis methods. A detailed treatment of our systematic un-
provement expected on the precision when combining two-point func- certainties is presented in Sec. 4, the results of our DES-Y1 data anal-
tions with aperture mass map statistics (Martinet et al. 2020). Similar ysis are discussed in Sec. 5, and we conclude afterwards. The Appen-
results are found in the context of a final Stage-III lensing experiment dices contain additional validation tests of our simulations and further
such as the (upcoming) DES-Y6 data release, where the combination details on our cosmological inference results.
of non-Gaussian statistics with the power spectrum√ method reduces the
error on the parameter combination S 8 ≡ σ8 Ωm /0.3 by about 25%
compared to a two-point function (Zürcher et al. 2020), where σ8 is 2 DATA AND SIMULATIONS
the normalisation amplitude of the linear matter power spectrum.
In the absence of accurate theoretical predictions for the signal, We present in this section the data and simulations included in our
the covariance, and the impact of systematics, non-Gaussian statistics analysis. We exploit multiple state-of-the-art simulation suites in order
must be carefully calibrated on numerical simulations specifically tai- to conduct our cosmological analysis, including a Cosmology training
lored to the data being analysed, which are generally expensive to run. set to model the response of our measurement to variations in cos-
Faster approximate methods exist (e.g. Izard et al. 2018), however they mology, as well as a Covariance training set and multiple Systematics
typically suffer from small scale inaccuracies exactly in the regime training sets. These DES-Y1-specific simulation products are created
where the lensing signal is the strongest, introducing significant biases from four suites of simulations, which we describe after introducing
in the inferred cosmological parameters. Previous peak count analy- the data.
ses of the third KiDS data release (KiDS-450) (Martinet et al. 2018, The total computing cost of the SLICS, cosmo-SLICS and
SLICS-HR are 12.3, 1.1 and 1.3 million cpu hours, respectively. They
4 LSST: www.lsst.org
5 Euclid: sci.esa.int/web/euclid 7 slics.roe.ac.uk
6 WFIRST: roman.gsfc.nasa.gov 8 www.magneticum.org
© 2020 RAS, MNRAS 000, 1–26

DES-Y1 Cosmology Beyond 2-Point Statistics 3
IM3SHAPE (Zuntz et al. 2013) that were both fully implemented

(see Zuntz et al. 2018, for details). While they provide consistent re-
sults, the former method has a larger acceptance rate of objects with
good shape measurements, and thereby results in measurements with
higher signal-to-noise. Following Troxel et al. (2018), we also adopt
the Metacalibration shear estimates in our analysis. This method pro-
vides a shear response measurement per galaxy, Rγ , a 2 × 2 matrix that
must be included to calibrate any measured statistics (we refer to Zuntz
et al. 2018, for more details on this calibration technique in the con-
text of shear two-point correlation functions). Additionally, the galaxy
selection itself can introduce a selection bias, which can be captured
by a second 2 × 2 matrix, labelled RS in T18, which we choose not to
include due to the small relative contribution. We compute from these
matrices the shear response correction, defined as S = Tr(Rγ )/2. As
explained in T18, the method imposes a prior on an overall multiplica-
Figure 1. Tiling strategy adopted to pave the full DES-Y1 data (black) with tive shear correction of m ± σm = 0.012 ± 0.013, which calibrates the
flat-sky 10 × 10 deg2 simulations (red squares). The squares overlap owing to
galaxy ellipticities as → (1 + m), with ≡ 1 + i2 .
the sky curvature, hence we separate the data at the mean declination in the
overlapping regions. In our pipeline, measurements are carried out in each tile
separately, then combined at the level of summary statistics.
Table 1. Survey properties. The effective number densities neff [in gal/arcmin2 ]
and shape noise σ listed here assume the definition of Chang et al. (2013).
The column ‘ZB range’ refers to the photometric selection that defines the four
DES-Y1 tomographic bins, while the mean redshift in each bin is listed under
hzDIR i.
tomo ZB range No. of objects neff σ hzDIR i

The galaxy sample is further divided into four tomographic red-
bin1 0.20 − 0.43 6,993,471 1.45 0.26 0.403 ± 0.008 shift bins based on the photometric redshift posterior estimated from
bin2 0.43 − 0.63 7,141,911 1.43 0.29 0.560 ± 0.014
griz flux measurements (Hoyle et al. 2018). The redshift distribution
bin3 0.63 − 0.90 7,514,933 1.47 0.26 0.773 ± 0.011
in these bins, ni (z), must then be estimated, and a number of methods
bin4 0.90 − 1.30 3,839,717 0.70 0.27 0.984 ± 0.009
are proposed to achieve this. The fiducial cosmic shear results pre-
sented in T18 are based on the Bayesian photometric redshift (bpz)
methodology described in Benı́tez (2000), which are consistent with
were produced on a system of IBM iDataPlex DX360M2 machines a n(z) estimated by resampling the COSMOS2015 field (Laigle et al.
equipped with one or two Intel Xeon E5540 quad cores, running at 2016) with objects of matched flux and size (Hoyle et al. 2018). How-
2.53GHz with 2GB of RAM per core. Every simulation was split into ever, the accuracy of these two methods has been questioned in Joudaki
64 mpi processes, each further parallelised with either four or eight et al. (2020, J20 hereafter), where it is argued that even though both
openmp threads. Modern compilers and cpus would likely bring the to- the bpz and COSMOS resampling estimates account for statistical un-
tal computing cost down if similar simulations had to be run again in certainty, residual systematics effects could significantly affect the in-
the future. ferred n(z) distributions.In particular, the COSMOS sample could be
populated with outliers and/or an overall bias that would affect the
calibration (e.g. figure 11 of Alarcon et al. 2020), and J20 proposes
2.1 DES-Y1 data
instead to calibrate with redshifts from matched spectroscopic cata-
In this paper we present cosmological constraints obtained from a re- logues10 . The direct reweighted estimation method (Lima et al. 2008,
analysis of the public lensing catalogues of the Year-1 data release9 of DIR hereafter) was selected for the fiducial cosmic shear analysis of
the Dark Energy Survey (Abbott et al. 2018b). These catalogues were the third KiDS data release (Hildebrandt et al. 2017, 2018), and for
obtained from the analysis of millions of galaxy images taken by the the DES-Y1 data re-analyses of J20 and Asgari et al. (2020b), where it
570 megapixel DECam (Flaugher et al. 2015) on the Blanco telescope is found that this calibration brings both DES-Y1 and KV450 results
at the Cerro Tololo Inter-American Observatory, observed in the grizY in excellent agreement, affecting the constraints on S 8 by only 0.8σ.
bands. The specific selection criteria of the DES-Y1 cosmic shear data It should be noted also that DIR has inherent systematic uncertainties
used in this paper exactly match those of the cosmic shear analysis that are hard to quantify. In particular, incomplete spectroscopy and
presented in Troxel et al. (2018): they consist of 26 million galaxies colour pre-selection (Gruen & Brimioulle 2017) can potentially bias
that pass the flags select, Metacal and the redMaGiC filters (Zuntz the DIR n(z). Despite these issues that can in principle be addressed
et al. 2018), thereafter covering a total unmasked area of 1321 deg2 , for by a pre-selection of sources via the self organising map technique
an object density of 5.07 gal arcmin−2 . The footprint of the DES-Y1 (Wright et al. 2020), we choose to adopt this DIR methodology for
data is presented in Fig. 1, which shows in black the galaxy positions simplicity and to be able to easily relate our findings to previous work.
from the selected sample. We use the same tomographic redshift distribution ni (z) and uncer-
The galaxy shears in the DES-Y1 data are estimated by two tainty about the mean redshift hziDIR i as in J20 here. In this method,
independent methods, Metacalibration (Sheldon & Huff 2017) and the uncertainties on the mean redshifts, σiz , are estimated from a boot-
strap resampling of the spectroscopic samples. The density, the mean
redshifts and the shape noise of the galaxies in individual tomographic
9 des.ncsa.illinois.edu/releases/dr1 bins are presented in Table 1.
© 2020 RAS, MNRAS 000, 1–26

Table 2. Summary of key properties from the four simulations suites used in our 2.3 Covariance training set
pipeline. Lbox is the box side [in h−1 Mpc], np is the number of particles evolved,
Our covariance matrix is estimated from the SLICS (Harnois-Déraps
Nsim is the number of N-body runs, NLC is the number of light-cones in the
full training set and Ncosmo is the number of cosmology samples. The bottom
et al. 2018, HD18 hereafter), a public simulation suite in which the
section summarises the range in cosmological parameters that is covered by the cosmology is fixed for every N-body run, but the random phases in the
cosmo-SLICS. initial conditions are varied, offering a unique opportunity to estimate
the uncertainty associated with sampling variance. The volume and
Sim. suite Lbox np Nsims NLC Ncosmo number of particles are the same as for the cosmo-SLICS, achieving
a particle mass of 2.88 h−1 M (see the properties summary in Table
cosmo-SLICS 505 15363 52 520 26 2). The light-cones are constructed in the same way as the cosmo-
SLICS 505 15363 124 124 1 SLICS, except that in this case the mass sheets are sampled only once
SLICS-HR 505 15363 5 50 1
per N-body run, generating 124 truly independent realisations. The
Magneticum 2 352 2 × 15833 1 10 1
accuracy of the SLICS has been quantified in Harnois-Déraps & van
Magneticum 2b 640 2 × 28803 1 10 1
Waerbeke (2015) by comparing their matter power spectrum to that of
parameter Ωm S8 h w0 the Cosmic Emulator (Heitmann et al. 2014), which match to within
sampling [0.1, 0.55] [0.6, 0.9] [0.6, 0.82] [-2.0, -0.5] 2% up to k = 2.0 hMpc−1 ; smaller scales progressively depart from the
emulator. The cosmo-SLICS have a similar resolution.
2.4 Systematics training set: mass resolution

2.2 Cosmology training set
Numerical simulations are inevitably limited by their intrinsic mass
The training set is constructed from the cosmo-SLICS (Harnois- and force resolution, and it is critical to understand how these af-
Déraps et al. 2019, HD19 hereafter), a suite of wCDM N-body sim- fect any measurements carried out on the simulated data. We employ
ulations specifically designed for weak lensing data analysis targeting for this purpose a series of ‘high-resolution’ runs, first introduced in
dark matter and dark energy. These simulations cover a wide range of Harnois-Déraps & van Waerbeke (2015) and labelled ‘SLICS-HR’
values in (Ωm , σ8 , h, w0 ). They sample the parameter volume at 25+1 therein. These consist of 5 independent N-body simulations similar
coordinates organised in a latin hypercube (25 wCDM plus one ΛCDM to the main SLICS suite, but in which the force accuracy of cubep3 m
point), and further include a sample variance suppression technique, has been increased significantly such as to resolve smaller structures,
achieving a sub-percent to a few percent accuracy depending on the even though the particle number is fixed. These have been shown to
scales involved. This is comparable to the accuracy of many widely- reproduce the Cosmic Emulator power spectrum to within 2% up to
used two-point statistics models based on non-linear power spectra k = 10.0 h−1 Mpc, indicating that even those small scales are correctly
from Halofit (Takahashi et al. 2012) or from HMcode (Mead et al. captured by the simulations. The SLICS-HR are post-processed with
2015, 2020). The full training range is detailed in Table 2, which also a strategy similar to that adopted for the cosmo-SLICS, re-sampling
influences our choice of priors when sampling the likelihood (see Sec. the projected mass sheets in order to generate 10 pseudo-independent
3.6). light-cones per run.
Each run evolved 15363 particles inside a 505h−1 Mpc co-moving
volume with the public cubep3 m N-body code (Harnois-Déraps et al.
2013), generating on-the-fly multiple two-dimensional projections of 2.5 Systematics training set: baryon feedback
the density field. These flat-sky mass planes were subsequently ar-
ranged into past light-cones of 10 degrees on the side, from which Another important systematic we investigate in this analysis is the im-
lensing maps were extracted at a number of redshift planes (see Sec. pact of strong baryonic physics that modifies the clustering property
2.6.1). This process was repeated multiple times after the mass planes of matter. As noted in multiple independent studies, AGN feedback
were randomly selected from a pool of six different projected sub- has a particularly important effect on the matter power spectrum but
volumes, then their origins were randomly shifted. In total, 50 pseudo- is challenging to calibrate. Simulations often struggle to reproduce the
independent light-cones per cosmology are available for the generation correct baryon fraction in haloes of different masses, and these dif-
of galaxy lensing catalogues (see HD19 for a complete description). In ferences in turn cause major discrepancies in the clustering properties
the end we include 10 light-cones per cosmology out of 50, after ver- (see Chisari et al. 2018, for example). In this paper, we examine one of
ifying that our results do not change when training on only five of these models and inspect which parts of our peak count measurements
them. Indeed, 1000 deg2 is enough to reach convergence on our statis- are affected by baryons.
tics, largely due to the sample suppression technique implemented in We used for this exercise a series of light-cones ray-traced from
HD19. a subset of the Magneticum Pathfinder hydrodynamical simulations,
Two of these models (cosmology-fid and -00, see HD19) are used which are designed to study the formation of cosmological structures
to infuse photometric redshift and shear calibration uncertainty, which in presence of baryonic physics and which were recently described in
we describe in Sec. 4.1 and 4.2, respectively. Castro et al. (2020). These are based on the smoothed particle hydro-
dynamics (SPH) code p-gadget3 (Springel 2005), in which a number
of baryonic processes are implemented, including radiative cooling,
star formation, supernovae, AGN and their associated feedback on the
matter density field. The Magneticum reproduce a number of key ob-
10
servations such as statistical properties of the large-scale, inter-galactic
Both the DIR and the COSMOS resampling methods have been shown to be
and inter-cluster medium, but also central dark matter fractions and
consistent with other n(z) estimation techniques such as the cross-correlation
between photometric and overlapping spectroscopic surveys (Morrison et al. stellar mass size relations (see Hirschmann et al. 2014; Teklu et al.
2017; Johnson et al. 2017; Hildebrandt et al. 2018; Gatti et al. 2020; Hoyle 2015; Castro et al. 2018, 2020, for more details). What is especially
et al. 2018). J20 also show that the DIR method is robust against the specific important in our case is that the total baryonic feedback on the matter
choice of spectroscopic calibration sample, provided that the combination is field is comparable to that of the BAHAMAS cosmological hydrody-
sufficiently wide and deep. namical simulations (McCarthy et al. 2017), in particular in terms of
© 2020 RAS, MNRAS 000, 1–26

the strength of the effect on the matter power spectrum. This derives must extend over tile boundaries. Both data and mock data are sepa-
from the similar baryon fractions produced by Magneticum and BA- rated in tiles at the catalogue level; these are then analysed individu-
HAMAS, which are in reasonable agreement with observations. This ally, and the data vectors are combined at the end12 .
validates the Magneticum as a good representation for the impact of Another subtle difference that needs to be taken into account is
baryon feedback, given the current uncertainty on the exact impact that the position coordinates (RA, Dec) and the galaxy ellipticities ()
(see Sec. 4.3 for further discussion). from the data are provided on the (Southern) curved sky, whereas all of
Among the various runs, we use a combination of the high- our simulations assume a (X, Y) cartesian coordinate system. Since the
resolution Run-2 (Hirschmann et al. 2014) and Run-2b (Ragagnin physics are independent of our choice of coordinate system, and since
et al. 2017), which both co-evolve dark matter particles of mass we analyse every tile individually, we apply a coordinate transforma-
6.9 × 108 h−1 M and gas particles with mass 1.4 × 108 h−1 M in co- tion to centre every tile onto the equator, where both coordinate frames
moving volumes of side 352 and 640 h−1 Mpc, respectively; the smaller coincide13 . The weak lensing statistics of a given tile are unaffected by
(larger) box is used at lower (higher) redshift, and the transition oc- this rotation, a fundamental fact that we verify with two-point correla-
curs at z = 0.31. The input cosmology is consistent with the SLICS tion functions in Sec. 3.1.
but slightly different, with Ωm = 0.272, h = 0.704, ns = 0.963 and As easily noticed from looking at Fig. 1, some of the galaxies fall
σ8 = 0.809. Both Run-2 and Run-2b also exist in pure gravity mode outside the tiles, which slightly affects the total number of galaxies in
(i.e. dark matter-only) with otherwise identical initial conditions, al- the sample. It is not ideal, but adding multiple simulated tiles for such
lowing us to isolate the impact of the baryonic sector on our observ- a small fraction (1.9%) of the data is arguably not worth the effort.
ables. The number of objects listed in Table 1 reflects this final selection and
amounts to a total of 25.5 million of galaxies.
2.6 Simulation post-processing

2.6.3 Mock galaxy shapes and redshifts
2.6.1 Light-cones
As mentioned above, the position and the intrinsic ellipticities of indi-
The simulation suites used in this paper all work under the flat sky
vidual galaxies in the simulated catalogues are taken from the obser-
approximation, which assumes that the maps are far enough from the
vations. Redshifts are assigned to every object in a given tomographic
observer so that cartesian axes can be used instead of angles and ra-
bin ‘i’ by sampling randomly the ni (z) described in Sec. 2. Therefore,
dial distances. At preselected redshifts z, the N-body/hydrodynamical
variations in survey depth are not included in our training sets. This
codes assign the particles onto a three-dimensional grid, select a sub-
induces a systematic difference with the data, but we expect that this
volume to be projected with pre-determined co-moving thickness, and
has a minor effect on our cosmological measurement, given the current
collapse the mass density along one of the axis. This procedure is
size of the statistical error. At this stage, every galaxy has position and
repeated with different projection directions and sub-volumes, creat-
a redshift, which are used to extract the lensing quantities (κ, γ) from
ing a collection of mass sheets at every redshift. These are next post-
the simulation light-cones.
processed to generate a series of past light-cone mass maps, δ2D (θ, z),
We finally include the intrinsic galaxy shapes and Metacal shear
each of 100 deg2 , which are then used to generate convergence κ(θ, zs )
response correction in the simulations by randomly rotating the ob-
and shear γ(θ, zs ) maps at multiple source redshift planes, zs , (see
served galaxy shapes such as to undo the cosmological correlations
HD18 and HD19 for full details), where γ = γ1 + iγ2 , the two com-
from the data, and we then combine the new ellipticity int with the
ponents of the spin-2 shear field. From these, mock lensing quantities
simulated lensing signal as:
(κ, γ) can be computed for any galaxy position provided its (RA, Dec)
coordinates and a redshift. int + g
= . (1)
1 + ∗int g
2.6.2 Assembling the simulated surveys Here g is the reduced shear, defined as g = γS /(1 + κ), and the bold-
font symbols g, γ, and int are again spin-2 complex quantities.
As for many non-Gaussian statistics, peak counts are highly sensitive We investigate in Sec. 4.4 the impact of the intrinsic alignment of
to the noise properties of the data. As such the simulations need to galaxies, where int is no longer chosen at random and instead corre-
reproduce exactly the position and shape noise of the real data, other- lates with the shape of dark matter haloes.
wise the calibration will be wrong. The solution, adopted in Liu et al.
(2015a), K16 and M18 is to overlay data and simulated light-cones,
and to construct mock surveys from the position and intrinsic shape of 3 THEORY AND METHODS
the former, and the convergence and shear of the latter.
Since the size of the full DES-Y1 footprint largely exceeds that of Since we validate our simulation suites with cosmic shear correlation
our individual light-cones, we connect the data and simulations with a functions and lensing power spectra, we begin this section with a re-
‘mosaic’ approach, where the DES-Y1 galaxy catalogues are divided view of the theoretical modelling and the measurement strategies re-
into smaller ‘tiles’11 that all fit inside 100 deg2 square areas. Each of lated to these quantities. We next move to the primary focus of this pa-
these tiles are then overlaid with simulated light-cones from which
lensing quantities are extracted. In total, 19 tiles are required to assem- 12 Note that the tiles are identical for all simulations (SLICS, cosmo-SLICS,
ble the full footprint with our mosaic, which is shown in Fig. 1. Every
SLICS-HR and Magneticum) since their light-cones all have the same opening
simulated light-cone from the Cosmology, Covariance or Systematics angle.
training sets is therefore replicated 19 times and associated with a full 13 In this process, we rotate both the celestial coordinate and the ellipticities of
realisation of the survey. every galaxy in the tile, to account for the modified distance to the South pole in
We emphasise that the simulated light-cones are discontinuous the new coordinate frame. The exact transformation uses the method presented
across tile boundaries, whereas data are not. To avoid significant cal- in the Appendix B of Xia et al. (2020), which rotates pairs of galaxies from
ibration biases caused by this difference, no measurement whatsoever any orientation on the sky onto the equator, placing one member at the origin.
In our case we instead map to the equator the straight line that bisects every
tile. Every tile has its unique rotation vector, which we also use to displace the
11 These tiles are sometimes called ‘patches’ in the literature, e.g. in K16. galaxies and to recompute their ellipticities (1/2 ) in this new coordinate frame.
© 2020 RAS, MNRAS 000, 1–26

per and describe our peak count statistics pipeline, detailing our treat- twice as small as the statistical error measured from the covariance
ment of the data, our approach to modelling the signal and estimating mocks (see below) and evenly scattered about the black points, vali-
the covariance matrix, and we finally describe our cosmological infer- dating our tiling method. We further verified that the difference on the
ence methods. inferred cosmology is negligible (see Sec. 5).
We also show, with the red squares in Fig. 2, the mean and ex-
pected 1σ error on the DES-Y1 data as estimated from the Covariance
3.1 ξ± statistics training set. The agreement between theory and simulations is excel-
Two-point correlation functions are well studied and have a key ad- lent at all scales, and the slight differences are well under the statistical
vantage over other measurement techniques: their modelling can be precision of the data. We can observe a slight loss of power at large
accurately related to the matter power spectrum, P(k, z), whose accu- angular scales in the ξ+ statistics, a finite box effect that we forward
racy is reaching the percent level far in the non-linear regime when model (see Appendix A1). For every simulated light-cone we gener-
calibrated with N-body simulations at small scales, and in absence ate a total of 10 realisations of the shape noise by rotating, as many
of baryonic physics (Heitmann et al. 2014; Euclid Collaboration: Kn- times, every galaxy in the catalogue, and recomputing new observed
abenhans et al. 2019). From this P(k, z), the lensing power spectrum ellipticities (with Eq. 1) and the correlation functions (Eq. 5). The red
between tomographic bins ‘i’ and ‘ j’ is computed in the Limber ap- squares in Fig. 2, as well as their associated error bars, correspond to
proximation as: one of these realisations; we observe no significant change in the other
Z χH i nine realisations, and recover the error bars reported by T17 to within
q (χ) q j (χ) ` + 1/2
5-15% over most angular scales, further demonstrating the robustness
C`i j = P , z(χ) dχ, (2)
0 χ2 χ of our training set. We do not expect a perfect match due to the slightly
where χH is the co-moving distance to the horizon, and the lensing different binning scheme.
kernels qi are computed from the redshift distributions n(z) as: Our cosmological analyses exclude the same angular scales as in
H 2 χ Z χH T18, removing the elements of the data vector where T18 conclude that
3 0 dz χ0 − χ 0
qi (χ) = Ωm ni (χ0 ) 0 dχ , (3) the uncertainty on the baryonic feedback and in the non-linear matter
2 c a(χ) χ dχ χ0 power spectrum is non-negligible. These scales are indicated by the
where c and H0 are the speed of light and the Hubble parameter, re- vertical lines in Fig. 2.
spectively. The cosmic shear correlation functions ξ±i j are computed The variation of ξ± with cosmology are well captured by cosmo-
from the C`i j as: SIS, and so are the responses to photometric redshift and shear cal-
Z ibration uncertainties (see Sec. 4). We therefore do not measure this
1 statistic in the Cosmology nor the Systematics training sets, and use
ξ±i j (ϑ) = C i j J0/4 (`ϑ) ` d`, (4)
2π ` ` instead the public modules provided in the latest cosmoSIS release to
with J0/4 (x) being Bessel functions of the first kind. Following T18, calculate these.
Eq. 4 is solved with the cosmological parameter estimation code cos-
moSIS14 (Zuntz et al. 2015), in which the matter power spectrum is
3.2 Shear peak count statistics
calculated by the Halofit model of Takahashi et al. (2012). The ξ±i j pre-
dictions for the DES-Y1 measurements are shown by the black lines in As mentioned in the introduction, the peak count statistic is a power-
Fig. 2 for all pairs of tomographic bins, at the SLICS input cosmology. ful alternative method to extract cosmological information from weak
The measurements of ξbi j from simulations and data are carried lensing data. It consists in measuring the ‘peak function’, i.e. the num-
±
out with Treecorr (Jarvis et al. 2004), a fast parallel tree-code that ber of lensing peaks as a function of their signal-to-noise, which is
computes shape correlations between pairs of galaxies ‘a, b’ separated very sensitive to cosmology and robust to systematics (see Zürcher
by an angle ϑ as: et al. 2020; Martinet et al. 2020, for recent comparisons with other
" # lensing probes).
Our measurement technique closely follows that described in K16
a,× a b,× b ∆ϑab

P i j i j
ab W W (θ )
a b a,t a b,t b (θ ) ± (θ ) (θ )
and M18, which we review here. Peaks are identified from local max-
ξ± (ϑ) =
bij
. (5)
ima in the signal-to-noise maps of the mass within apertures (Schnei-
P
ab Wa Wb S a S b
der 1996), Map (θ), searching for pixels with values higher than their 8
In the above expression, the sums are over all galaxies ‘a’ in tomo-
neighbours. This is one of many ways to estimate the projected mass
graphic bin i and galaxies ‘b’ in tomographic bin j; a,t
i
and a,×
i
are the
map from galaxy lensing catalogues, and was chosen primarily for its
tangential and cross components of the ellipticity of galaxy a in the di-
local response to data masking. This is to be contrasted with e.g. the
rection of galaxy b; Wa/b are weights attributed to individual galaxies,
Fourier methods of Kaiser & Squires (1993) in which masking intro-
which are set to unity in the Metacalibration shear inference method;
duces a complicated mode-mixing matrix that can affect all scales.
S a/b are the ‘shear response correction’ per object mentioned in Sec.
Other techniques such as Bayesian mass reconstruction (Price et al.
2.1 and provided in the DES-Y1 catalogue; ∆ϑab is the binning oper-
2020) or wavelets transforms (Leistedt et al. 2017) are also promising
ator, which is equal to unity if the angular separation between the two
and merit to be explored in the future (see as well Gatti et al. 2020, and
galaxies falls within the ϑ-bin, and zero otherwise. Our raw measure-
references therein).
ments are organised in 32 logarithmically-spaced ϑ-bins, in the range
From a lensing catalogue containing the position, ellipticity a
[0.5 − 475.5] arcmin, but not all angular scales are used in this work.
and shear response correction S a per galaxy, we construct an aperture
We present in Fig. 2 our measurements of ξbi j on the DES-Y1 data,
± mass map on a grid by summing15 over the tangential component of
showing with the black solid points the measurement on the full foot- the ellipticities from galaxies surrounding every pixel at coordinate θ,
print, and with the open magenta circles the measurements on the 19 weighted by an aperture filter Q. More precisely, we compute:
data tiles described in Sec. 2.6.2, which are combined with a weighted
1 X
mean using the Treecorr Npairs (ϑ) per tile as our weights. We see that Map (θ) = P a,t (θ, θa )Q(|θ − θa |, θap , xc ), (6)
the two results are similar, with differences that are everywhere at least ngal (θ) a S a a
14 bitbucket.org/joezuntz/cosmosis/wiki/Home 15 In practice, we use a link-list to loop only over nearby galaxies.
© 2020 RAS, MNRAS 000, 1–26

10 1 10 2 10 1 10 2 10 1 10 2 10 1 10 2
2
1
0
11 12 13 14
-1
4
4
2 2
0 0
44 22 23 24
4 4
2 2
0 33 33 34 0
34
4
4
2
2
0
24 23 22 44 0
2
1 DES-Y1 footprint
DES-Y1 tiles
0 SLICS
14 13 12 11 Theory
-1
1 2 1 2 1 2 1 2
10 10 10 10 10 10 10 10
Figure 2. Two-point correlation functions measured in the DES-Y1 data (filled and opened circles present measurements on the full survey footprint and from a
weighted mean over the tiles presented in Fig. 1, respectively) and in the SLICS simulations (red squares, with error bars showing the statistical error on the DES-Y1
data), compared to the analytical model (solid lines). The left- and right-hand side ladder plots present the ξ− and ξ+ statistics respectively, and the sub-panels in each
correspond to different combinations of tomographic bins. The vertical dotted lines indicate the angular scales excluded in the cosmological analysis, which match
those of T18.
where ngal (θ) is the galaxies density in the filter centred at θ, and θa is separated in tomographic bins (which we label 11, 22, 33 & 44 as
the position of galaxy a. The tangential ellipticity with respect to the for the 2PCF), and then from every combination of pairs of tomo-
aperture centre is computed as a,t (θ, θa ) = −[1 (θa ) cos(2φ(θ, θa )) + graphic catalogues (which we label 12, 13, 14 ... 34). As detailed in
2 (θa ) sin(2φ(θ, θa ))], where φ(θ, θa ) is the angle between both coordi- Martinet et al. (2020), analysing these ‘cross-tomographic’ catalogues
nates. Our filter Q(θ, θap , xc ), abridged to Q(θ) to shorten the notation, provides additional information that is not contained within the ‘auto-
is identical to that in Schirmer et al. (2007), which is optimal for de- tomographic’ case. They went further and also included triplets (123,
tecting haloes following an NFW profile (but faster than solving the 124, 134...) and quadruplets (1234), showing that these also contained
actual numerical NFW equation): additional information, but this gain is not as significant in our case,
tanh(x/xc ) −1 where the noise levels are much higher.
Q(x) = 1 + exp(6 − 150x) + exp(−47 + 50x) . (7)
x/xc
In the above expression, x = θ/θap , where θ is the distance to the filter
centre, and we adopt xc = 0.15 as in previous works, a choice that 3.2.1 Masking
maximises the sensitivity of the signal to the massive haloes, which
Weak lensing data are taken inside a survey footprint, and parts of
carry the majority of the cosmological information. The filter size of
the images are removed in order to mask out satellite tracks, bright
θap = 12.5 arcmin is adopted as in M18, however we also consider
stars, saturated foreground galaxies, etc. The effect of data masking
9.0 and 15.0 arcmin, and report results for these where appropriate.
on the aperture mass map can be significant: the signal and the noise
Hereafter, Eq. (6) defines the signal of our aperture mass map, which
are coherently diluted in apertures that strongly overlap with masked
we compute at every pixel location.
regions, generating regions where M is overly smooth. Therefore the
The variance about this map is calculated at every pixel location
survey mask must be included in the simulations and in the estimator
from:
such as to avoid biasing the statistics.
1 X
σ2ap (θ) = P 2 | a |2 Q2 (|θ − θa |), (8) If the masked pixels are known, this can be taken into account
2
2ngal (θ) a S a a by avoiding pixels for which e.g. more than half of the filter overlaps
with masked areas. Alternatively, one can examine the object density
where again the sum runs over all galaxies in the filter. Note that the in the aperture filter, ngal (θ), and require that it exceeds a fixed thresh-
magnitude of the measured galaxy ellipticities that enters this equa- old in order to down-weight or reject heavily masked apertures. In this
tion must also be calibrated by the shear response correction (see the method, pixels with little or no galaxies are treated as masked. We
Appendix A of Asgari et al. 2020b), hence the term [ a S a ]2 in the de- opted for the second method, setting the threshold to 1/π2 gal/arcmin2 ,
P
nominator. The signal-to-noise map from which peaks are identified, which directly identifies regions with very low galaxy counts. We fur-
M(θ), is computed by taking the ratio between Eq. (6) and the square ther augment the masking selection with an apodisation step that flags
root of Eq. (8) at every pixel location, e.g. as ‘also masked’ any pixel within a distance θap of a masked region
found in the first step.
M(θ) ≡ Map (θ)/σap (θ). (9)
Fig. 3 illustrates this procedure for one of the tiled catalogues,
Peaks catalogues are first constructed from the galaxy catalogues for an idealised noise-free case. For our fiducial choice of filter θap =
© 2020 RAS, MNRAS 000, 1–26

vestigations). 12 bins is also a good trade-off between our need to cap-

ture most cosmological information from Npeaks (S/N), while keeping
a small data vector for which the covariance matrix will be easier to
estimate. A number of recent studies (M18, Davies et al. 2019; Coul-
ton et al. 2020; Zürcher et al. 2020; Martinet et al. 2020; Davies et al.
2020) have shown that cosmological information is contained in peaks
of negative S/N or in lensing voids, however, as noted in Appendix B
of M18, the peaks with negative S/N strongly correlate with those of
positive S/N value and only marginally improve the constraints from
peak statistics in the case of Stage III surveys. We therefore focus only
on the positive peaks in this DES-Y1 analysis.
We show in Fig. 4 the peak function measured from the Cosmol-
ogy training set with θap = 12.5 arcmin, for all pair combinations of
the four redshift bins and colour-coded as a function of the input S 8 .
A pure noise case (Nnoise ), obtained from the average peak function af-
ter setting γ = 0 on 10 full survey realisations, has been subtracted to
highlight the cosmological variations. Off-diagonal panels present the
cross-tomographic measurements. The colour gradient is clearly visi-
ble in all tomographic bins; more precisely, all cosmologies present an
excess of large S/N peaks and a depletion of low S/N peaks com-
pared to pure noise. This is caused by the gravitational lensing sig-
nal, which create peaks and troughs in the Map map and smooths out
the smallest peaks. Importantly, these differences are accentuated for
high-S 8 cosmologies. Also shown with black squares are the measure-
ments from the DES-Y1 data, with error bars estimated from the Co-
variance training set. These demonstrate that most of the constraining
Figure 3. Example of the Map mass-reconstruction pipeline over one of our
power comes from the auto-tomographic bins 3 and 4 and from the
10 × 10 deg2 tiles. The larger panel on the bottom right presents the true κ
cross-tomographic bins. Some additional information is contained in
values at the position of the galaxies in this field, extracted from the cosmo-
SLICS model-00. The raw Map map is shown in the top left panel in the noise- the highest S/N peaks of the redshift bins 1-2, whereas the low S/N
free case. The number of galaxies in the filter (second panel) are then used to peaks of bin 2 mostly contribute noise.
construct a mask (third panel), which we apply on the raw Map maps (bottom
panel). The top right panel shows a zoom-in of the top left panel, highlighting
3.3 Analysis pipeline
the effect of masking on the raw reconstructed Map map.
In this analysis we extend multiple aspects of the K16 and M18
methodologies. Here is a summary of these improvements:
12.5 arcmin, we show in the upper two left panel the ‘raw’ Map (θ) map (i) We include a tomographic decomposition of the data, including
(e.g. before masking, computed directly from Eq. 6) as well as Ngal (θ). the cross-redshift pairs inspired by the method presented in Martinet
The masked regions are clearly visible in the latter but not so much in et al. (2020);
the former. A close inspection (top right panel) however reveals overly (ii) Our Cosmology training set (see Sec. 2.2) now includes four
smooth features in Map (θ), in regions where there are no galaxies (i.e. parameters (Ωm , σ8 , h and w0 ), and it would be straightforward to in-
in the blue regions of the Ngal (θ) map). The third left panel shows the crease that parameter list with additional training samples. Addition-
masked regions constructed from our pipeline, which is finally applied ally, the cosmo-SLICS simulations are more accurate than those of Di-
on the raw aperture map, resulting in the masked map shown on the etrich & Hartlap (2010), which were used in both K16 and M18: they
bottom panel. All choices of θap result in aperture maps that closely resolve smaller scales, and suffer less from finite box effects, having a
recover the true convergence (shown in the bottom right panel). volume almost 8 times larger;
It is clear from Fig. 3 that the unmasked area of our final maps is (iii) We deploy a fast emulator (see Sec. 3.5) that can model the
affected by the aperture filter size. Indeed, larger filters can be blind to signal at arbitrary cosmologies within the parameter volume included
small features in the mask, while the survey edges are more severly ex- in the training. In contrast with a likelihood interpolator, emulating the
cluded. This does not bias our cosmological inference since we apply data vector directly allows us to combine the summary statistics with
the same filter to the data and the simulations, but it does slightly af- other measurement methods such as the two-point correlation func-
fect the signal-to-noise of our measurement, which increases with the tions, to better include systematic uncertainties, and to easily interface
area of the survey. The net unmasked area in our final maps are (1426, with most MCMC samplers;
1408, 1366, 1327, 1284) deg 2 for θap = (6.0, 9.0, 12.5, 15.0, 18.0)’, (iv) We generate a Covariance training set from a larger ensem-
respectively. ble of independent survey realisations (see Sec. 2.3), and feed it into
a novel hybrid internal resampling technique that improves the accu-
racy and precision of lensing covariance matrices estimated from the
3.2.2 Peak function
suite (see Sec. 3.4). Moreover, the covariance training set is shown to
Peaks found in the (masked) M maps are counted and binned as a closely reproduce the published DES-Y1 cosmological constraints of
function of their pixel value, thereby measuring the peak function T18 when analysed with two-point correlation functions (see Table 5),
Npeaks (S/N). We use 12 bins covering the range 0 < S/N 6 4 in our thereby validating both the simulations themselves and the covariance
main data vector, which was found in K16 and M18 to avoid scales estimation pipeline. Our method is also compatible with joint-probe
where multiple systematics uncertainties such as the effects of baryon measurements;
feedback and intrinsic alignments of galaxies become large (we ex- (v) We construct a series of dedicated Systematics training sets
tend this range to higher S/N values in some of our systematics in- specifically tailored to our data, in which the most important cosmic
© 2020 RAS, MNRAS 000, 1–26

0 2 40 2 40 2 40 2 4
1000
-1000
11 12 13 14
1000
-1000
22 23 24
0
-1000
33 34 -2000
500
0
-500
44 -1000
S 8 = 0.61 0.68 0.76 0.83 0.90 -1500
Figure 4. Peak function Npeaks (S/N) in the DES-Y1 data (black squares) and simulations (coloured histograms), from which the expectation from pure shape noise
Nnoise (S/N) has been subtracted. The panels show different tomographic bin combinations, as labelled in their lower-right corners. The predictions are colour-coded
by their S 8 value. The DES-Y1 error bars are estimated from the Covariance training set.
shear-related systematics are infused. Specifically, we investigate the that these approaches do not make too much of an impact for cur-
impact of photometric redshift uncertainty (Sec. 4.1), of multiplica- rent cosmic shear data. Moreover, it is not our primary goal here to
tive shear calibration uncertainty (Sec. 4.2), of baryonic feedback (Sec. optimise these choices, as we are rather interested in establishing our
4.3), of possible intrinsic alignment of galaxies (Sec. 4.4), and of lim- simulation-based inference method as being robust, accurate and flexi-
its in the accuracy of the non-linear physics (Sec. 4.5). These tests ble. We therefore opted for an overall analysis pipeline that maximally
allow us to flag the elements of our data vector that are affected, and in resembles that of the fiducial DES-Y1 cosmic shear data, and leave
some case to model the impact. Following K16, we construct a linear some additional tuning for future work.
response model to the photometric redshift and multiplicative shear Aside from a different choice of n(z) calibration (see Sec. 2.1),
calibration uncertainties, but calibrate our models on a sample of 10 the key differences between our current pipeline and that of T18 are
deviations from the mean of the distribution (respectively ∆zi and ∆mi , the impossibility of ours to vary and marginalise over the other cosmo-
with i = 1..10) as opposed to one; logical parameters – the power spectrum tilt parameter ns , the baryon
density Ωb and the sum of neutrino mass mν . These would require
P
(vi) We implement the emulator, the covariance matrix and the lin-
ear systematic models within the cosmology inference code cosmoSIS more light-cone simulations such as the Mira-Titan (Heitmann et al.
(Zuntz et al. 2015), allowing us to carry out a joint MCMC analy- 2016) or the MassiveNuS (Liu et al. 2018), which are not folded in
sis based on peak statistics and shear two-point correlation functions, our training set at the moment, but form a natural extension to this
while coherently marginalising over the nuisance parameters that af- work. Also missing is a cosmology-dependent model for the effect of
fect both of these measurement methods. intrinsic alignment, which will be necessary in future analyses.
With this new pipeline, we are fully equipped to investigate the
impact of different measurement and modelling methods, of different
systematics mitigation strategies, but also of analyses choices related 3.4 Covariance matrix
to the likelihood sampling, such as the prior ranges, the specific com- The covariance matrix is a central ingredient to our cosmological infer-
bination of parameters to be sampled, or the manner in which maxi- ence as it describes the level of correlation between different elements
mum likelihoods and confidence intervals are reported. Indeed, it has of our data vector, and its inverse directly enters in the evaluation of the
been shown that these have a non-negligible impact on the final cos-
mological constraints, specifically in the context of weak lensing cos-
mic shear analyses (Joudaki et al. 2017; Chang et al. 2019; Joudaki
et al. 2020; Asgari et al. 2020a; Joachimi et al. 2020)16 . It turns out introduce biases when collapsing a high-dimensional hyper-volume into a one-
dimensional space. Instead, it is argued therein that a more accurate 1D infer-
ence is obtained by reporting the multivariate maximum a posteriori (MAP)
16 In particular, Joachimi et al. (2020) demonstrates that reporting the pro- distribution, along with a credible interval calculated using the projected joint
jected maximum likelihood value and the associated confidence interval can highest posterior density (PJ-HPD) of the full likelihood.
© 2020 RAS, MNRAS 000, 1–26

1
1 2 3 1 2 3 1 2 3 1 2 3
0.05
-0.05
0.05
0.5 0
-0.05
0.05
-0.05
0.05
0
0
-0.05
S8 = 0.61 0.68 0.76 0.83 0.90
Figure 6. Accuracy of the GPR emulator, computed with a leave-one-out cross-

validation test. The results are colour-coded with the input S 8 value of the re-
moved training point, and compared with the statistical uncertainty on the mea-
surement (shown with the black dashed lines). The black solid line shows the
-0.5 accuracy at the ΛCDM node, and the different panels show the ten combina-
Figure 5. This cross-correlation matrix highlights the correlations between the tions of tomographic bins.
different elements of the data vector. From left to right, the first 10 blocks show
the ξ+ tomographic measurements, followed by the 10 ξ− blocks, while the last
10 blocks show the tomographic peak count. Not all elements are used in the 3.5 Peak function emulator
analysis, see the main text for more details.
The peak count statistics measured from the Cosmology training set
(shown in Fig. 4 with the coloured histograms) is computed at 26
points in a wide 4-dimensional volume. From these we train a Gaus-
likelihood. In our analysis, it is estimated from the Covariance train- sian Process Regression (GPR) emulator that can model the peak func-
ing set, which is based on 124 independent light-cones, each replicated tion given an input set of cosmological parameters [Ωm , S 8 , h, w0 ] at
onto the 19 survey tiles such as to fully cover the DES-Y1 footprint. any point within the training volume. Directly adapted from the pub-
For each of these survey realisations, we further generate 10 shape lic cosmo-SLICS emulator17 described in Appendix A of HD19, we
noise realisations by randomly rotating the ellipticity measurements train our GPR emulator on the individual elements of the Npeaks data
from the data, which increases the number of pseudo-independent re- vector, first optimising the hyper-parameters from an MCMC analysis
alisations to Nsim = 1240 and is largely enough for the current analysis. that includes 200 training restarts, then ‘freezing’ the emulator once
The next step consists in combining the measurements obtained the best-fit solution has been found. As described in HD19, the train-
in the 19 different tiles into a final measurement of the covariance. ing can also involve a PCA decomposition and a measurement error;
This proves to be slightly tricky to achieve because of the way light- we include the former but find that the modelling is more accurate
cones are recycled over the survey. In the case of peak statistics, which without the latter.
are noise dominated, it is sufficient to concatenate the peak catalogues We evaluate the accuracy of the emulator from a leave-one-out
from the different tiles into a master survey peak function, then use cross-validation test: the emulator is trained on all but one of the train-
the different light-cones and shape realisations to evaluate the covari- ing nodes, then generates a prediction of the peak function at the re-
ance matrix (this method was used in both K16 and M18). Unfortu- moved cosmology, which is finally compared with the actual measure-
nately, this method cannot be directly transposed to the ξ± statistics ment. This test is performed for all nodes and provides an upper bound
because the measurement is signal-dominated; this means that in this on the interpolation error, since in this case the distance between the
case the sampling variance observed in every light-cone is coherently evaluation point and all other training nodes is significantly larger than
reproduced over all 19 tiles, resulting in an overall variance that is if all points had been present. Moreover, many of these points lie at
an order of magnitude too large. This can be avoided either by A) the edge of the training volume, hence removing them for this test
computing the covariance matrix in individual tiles, and combining requires the emulator to extrapolate from the other points, which is
them with an area-weighted average, or by B) mixing the light-cones significantly less accurate than the interpolation that is normally per-
at the survey construction stage, such that for each full survey assem- formed. As discussed in HD19, the node at the fiducial cosmology
bly, the 19 tiles are extracted from 19 different light-cones selected at was added by hand close to the centre of the wCDM Latin hypercube,
random. Both of these produce indistinguishable covariance matrices hence the cross-validation test performed at that single ΛCDM point
for ξ± , but the former is harder to combine with the peak statistics, for is more representative of the actual emulator’s accuracy.
which the covariance has to be computed on the full survey. We there- The results from this accuracy test are presented in Fig. 6, again
fore opted for method B) when estimating the joint covariance matrix colour-coded with S 8 . We achieve sub-percent interpolation accuracy
about [ξ±i j (ϑ); Npeaks
ij
(S/N)]. The net effect of including 10 shape noise over data points with S/N < 3, and for all points when testing the
realisations is to significantly lower the noise in the matrix, especially ΛCDM model (shown in thick black). We observe in some other mod-
over the terms where shape noise dominates. A similar technique is els a scatter of up to a few percent for peaks with S/N > 3, but this
applied in M18, who also find a negligible impact on the cosmological scatter over-estimates the true interpolation error for reasons explained
results from KiDS-450 data whether they average over 5 or 20 noise above. When compared to the statistical error on the DES-Y1 measure-
realisations. Fig. 5 shows the resultingpmatrix, normalised to unity of ment, the GPR emulator’s error is always subdominant (see the thick
the diagonal, e.g. r(x, y) = Cov(x, y)/ Cov(x, x) Cov(y, y). Whereas dashed lines in Fig. 6). We conclude from this that the accuracy of
the ξ+ block shows the highest level of correlation, the off-diagonal
blocks are mostly uncorrelated, which is promising for the prospect of
learning additional information from the joint analysis. 17 github.com/benjamingiblin/GPR Emulator/
© 2020 RAS, MNRAS 000, 1–26

our model is high enough, given our current data, and that it should over them. In particular, for the two-point correlation function sector,
introduce no noticeable bias in the cosmological inference. the photometric redshift errors ∆zi are included by shifting the tomo-
graphic redshift distributions, e.g. ni (z) → ni (z + ∆zi ), after which
new predictions for ξ±i j are computed from Eqs. (2 – 4). Errors on
3.6 Likelihood the shear calibration, ∆mi , are included directly on the statistics as
The GPR emulator is embedded within the cosmoSIS cosmological in- ξ±i j → ξ±i j (1 + ∆mi )(1 + ∆m j ). Finally, we include the two-parameters
ference package, which allows us to evaluate the likelihood at any cos- model of intrinsic alignments of galaxies that was used in T18. We
mology within the cosmo-SLICS training range, given measurements keep the T18 priors on the IA and shear calibration nuisance parame-
of ξ±i j (ϑ) and Npeaks
ij
(S/N) from the data, plus our joint covariance ma- ters, but use the J20 priors for the redshift uncertainty in our fiducial
trix. At the moment we can only provide predictions of the peak func- analyses. These are all summarised in Table 3. Inaccuracies at small
tion for the θap values on which the GPR was trained, but in the future scales due to uncertainty in the non-linear physics and in the baryonic
this could also be treated as a free parameter to be emulated, providing feedback are avoided with the angular scale cuts applied on ξ± .
even more flexibility to the prediction code and optimisation avenues. We expand the existing cosmoSIS infrastructure to include sys-
The predictions x at cosmology π are then compared with the data tematics models of the peak statistics based on our simulation training
d using a multivariate t-distribution likelihood following Sellentin & sets. Specifically, we added modules to include the effect of photo-
Heavens (2016): metric and shear calibration errors, which we parameterised by the
same ∆zi and ∆mi nuisance parameters as for the 2PCF, allowing us to
Nsim
marginalise coherently over these in a joint [ξ± ; Npeaks ] analysis. These
L(π|d) ∝ ln 1 + χ2 /(Nsim − 1) , with (10)
2 two models are detailed in Secs. 4.1 and 4.2. All other sources of un-
X certainty are identified and controlled by removing S/N > 4 bins,
χ2 = [x(π) − d]T Cov−1 [x(π) − d]. (11) which are identified from our Systematics training set as being con-
taminated beyond an acceptance threshold, or by verifying that they
This likelihood correctly takes into account the residual noise in the
do not impact the best-fit cosmological parameters. The following sub-
covariance matrix that stems from its sampling with a finite number of
sections detail our treatment of these sources of systematic uncertainty
simulations, and reduces to the standard multivariate Gaussian likeli-
for the peak count measurements; the reader hungry for results can
hood when Nsim → ∞. Since there are, at the very least, hundreds of
skip ahead to Section 5.
peaks in each of our bins, adopting this likelihood is justified.
For our first pipeline validation exercise, our choice of priors
matches that of T18 in our two-point statistics-only analysis (see Ta-
ble 3), allowing us to investigate the effect of replacing the fiducial 4.1 Photometric redshifts
(analytic) covariance matrix with our simulation-based matrix on the
As there is no analytic prescription to model the effect of photometric
parameter inference. At the same time this serves to validate our sim-
redshift uncertainty on the peak function, we investigate its impact di-
ulations.
rectly from the simulations: we infuse different shifts ∆zi in the galaxy
In our three fiducial analyses (2PCF, peaks and joint), the priors
distribution of the four tomographic bins (similar to the treatment in
reflect the parameter range probed by the cosmo-SLICS, and hence
the 2PCF), generate new mock light-cones and galaxy catalogues, and
we assign a flat prior on Ωm , σ8 , h and w0 (summarised in Table 3).
compute the peak statistics from these. Our approach is similar to the
All other parameters (i.e. the baryon density, the tilt in the primor-
linear model adopted by K16, which computes a two-point numerical
dial power spectrum and the sum of neutrino masses) are kept fixed to
derivative from simulations produced with a (single) shift in the mean
Ωb = 0.0473, ns = 0.969 and mν = 0.0eV, respectively. The lensing
P
n(z) by ∆z. Our model is slightly more sophisticated: we sample 10
constraints on these are very weak at the moment, hence we do not
∆zi values in every S/N bin, drawing numbers randomly from Gaus-
expect that holding them fixed should significantly affect our results.
sian distributions with means of zero and standard deviations σz given
We finally include the same 10 nuisance parameters as in T18: a
by J20 priors on the nuisance parameters (see Table 3). In the case of
shear calibration ∆mi and a photometric redshift calibration ∆zi per to-
cross-tomographic bins, we use the mean between the two σz values.
mographic bin, plus two parameters associated with the modelling of
For every shift, we generate five light-cones and use these to cover
the intrinsic alignments (IA) in the non-linear alignment model (Bridle
the full survey (with the tiling strategy described in Sec. 2.6.2); we
& King 2007, NLA). The latter two are not included in the peak count
next measure Npeaks (S/N) as a function of ∆zi , and fit a straight line
case for which we conduct instead a simulation-based assessment of
the impact of IA. We sample the likelihood with the MultiNest sam- through these 10 points in order to extract the numerical derivative
∂Npeaks /∂∆z for every S/N and tomographic bin. We carry out this
pler (Feroz et al. 2009), set with a tolerance parameter of 0.1 and an
calculation for the different aperture filters investigated in this paper,
efficiency of 0.3. We refer the interested reader to T18 and Krause
but also at two different cosmological models18 (model-FID and -00) in
et al. (2017) for more details about the DES-Y1 cosmology inference
order to assess the stability of the derivative with respect to cosmology.
pipeline.
The results are shown in the upper panels of Fig. 7.
We validate our implementation with a series of MCMC analy-
We note that the results for the two cosmological models are in
ses where the ‘data’ is taken from the mean measurement extracted
qualitative agreement, where the response of high (low) S/N peaks
from the Covariance training set, for which the cosmology and the
to an increased survey depth is positive (negative). This is caused by
systematic biases are known. We detail these results in Appendix A,
the fact that a greater depth increases the shear signal, which shifts the
and compare the 2PCF and the peaks wCDM performance on these
peak function towards higher S/N values. Some differences are ob-
mocks as well. We further validate our ΛCDM 2PCF segment both
served, the largest being at the low-S/N end of the derivative, with
against the T18 and J20 results in Sec. 5.1.
some also seen at the high end. With model-00 being quite distant
4 SYSTEMATICS
18 The cosmo-SLICS model-FID has Ωm = 0.2905, S 8 = 0.8231, h = 0.6898
The likelihood modules within cosmoSIS are equipped with an infras- and w0 = −1.0, while model-00 has Ωm = 0.3282, S 8 = 0.6984, h = 0.6766
tructure that allows us to define nuisance parameters and to marginalise and w0 = −1.2376.
© 2020 RAS, MNRAS 000, 1–26

1000 11 22 33 44 12 13 14 23 24 34
0
-1000
-2000
-3000
1000
0
-1000
-2000
-3000
0 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4
model-FID
model-00
Figure 7. Derivative of the peak function Npeaks with respect to shifts in the mean of the galaxy redshift distribution, ∆z (upper) and in the mean shear calibration ∆m
(lower). These derivatives are shown as a function of S/N bin, for every tomographic case and are computed from 10 deviation points (see main text for details). The
cosmo-SLICS model-FID is in red, model-00 in black, and the error bars are obtained from the scatter over 10 realisations of the full survey (see the footnote for their
values). Other filter sizes yield slightly different derivatives, but exhibit a similar level of agreement between the two cosmological models. The FID derivatives are
used in the cosmology pipeline to marginalise over these two nuisance parameters.
from model-FID – notably a 12% lower value of S 8 – we do not ex- both the shear estimate and the noise in a non-trivial manner (see Eqs.
pect the derivatives to be identical. Given the current statistical uncer- 6 and 8).
tainty of the measurements however, we ignore this uncertainty on the
derivative, but this will need to be revisited in the future. Similarly
to K16, we choose to model the redshift uncertainty by scaling the 4.3 Baryon feedback
measurement in each of the S/N bins with this linear model prior to
The uncertain impact of baryonic feedback on the peak count statistics
marginalising over ∆z in the likelihood
i j inference. Namely, we compute
has received an increasing degree of attention over the last few years
ij
Npeaks (∆zi j ) = Npeaks
ij
(∆z = 0) + ∂Npeaks /∂∆z ∆zi j , using the derivative
(Osato et al. 2015; Weiss et al. 2019; Coulton et al. 2020). The cur-
extracted from the model-FID cosmology, where ∆zi j = ∆zi + ∆z j /2. rent interpretation can be summarised as follows: radiative pressure
from sustained stellar winds, combined with supernovae explosions
and AGN activity combine to expel gas towards the outer regions of
4.2 Shear calibration the haloes. These mechanisms are maximally efficient on medium-size
(e.g. 1014 M ) clusters (e.g. McCarthy et al. 2017), since light haloes
The uncertainty in the shear calibration is forward-modelled with a generally do not host AGNs, while heavier haloes manage to keep the
similar method, except that no additional ray-tracing is required. In- material inside due to their deeper gravitational wells. This redistri-
stead we include a uniform (1 + ∆mi ) correction factor at the cata- bution of matter tends to reduce the number of high S/N peaks, and
logue level, which multiplies every observed ellipticities in Eqs. (6) possibly augment that of smaller S/N values, however the exact size
and (8). We next de-bias the peaks measurement with the original shear of this effect is highly uncertain and depends on the feedback model
response S a but deliberately ignore these additional ∆mi factors, re- adopted. Just as for cosmic shear two-point correlation functions, its
sulting in a net residual bias caused by an incorrect shear calibration, significance depends on the noise level of the data. We note that for
which is exactly what we wish to model. We repeat this process for 10 less massive haloes, stellar feedback is also important, however this
values of ∆mi sampled from the priors on the shear calibration errors occurs at significantly smaller scales not probed by our filter. Also,
(a Gaussian with width σm = 0.013, see Table 3), we measure the peak radiative cooling at high redshift should produce more concentrated
function for each of these samples on the full survey, average over 5 haloes and could enhance the lensing signal, acting in the opposite di-
light-cones, and fit a straight line to these points in every S/N bin rection of the AGN feedback.
and tomographic bin, allowing us toi j compute derivatives and model The impact on the lensing power-spectrum computed on a single
ij
Npeaks (∆mi j ) = Npeaks
ij
(∆m = 0) + ∂Npeaks /∂∆m ∆mi j . redshift slice (taken to be zs = 0.97 here) of the Magneticum reaches
The partial derivatives are calibrated this way for cosmological 15% at high-`, as seen in Fig. 8, an amplitude that is similar to those
models-FID and -00 and reported in the lower panels of Fig. 7. We of the BAHAMAS simulations mentioned in Sec. 2.5, and which are
observe a global agreement between the two cosmologies, similar in consistent with the PCA constraints.
shape to the effect of increasing the survey depth, but the amplitudes To investigate the relative impact of baryons on the peak function
of the derivative are in slight disagreement. Model-FID has a slightly specifically for our analysis, we tile the full DES-Y1 survey with the
larger derivative, which is linked to its higher S 8 value. We ignore Magneticum light-cones introduced in Sec. 2.5 either with the full hy-
these differences in this work due to the small size of the σm rela- drodynamical simulations or with the gravity-only solution, then eval-
tive to the statistical precision of the DES-Y1 data, and use solely the uate and compare the peak functions. We repeat this process and aver-
model-FID derivatives in the MCMC sampling. However, this could
be age the results over the 10 pseudo-independent light-cones, each fur-
easily addressed with our approach: the derivatives ∂Npeaks /∂∆z and ther sampled with 10 shape noise realisations. We also extend the data

∂Npeaks /∂∆m could be estimated at our 26 cosmological nodes, from vector to higher S/N values in order to verify where the baryons start
which we could train a Gaussian process emulator the same way we to play an important role. We additionally repeated the process in a
model our signal. This improved accuracy treatment of the derivatives non-tomographic set-up, where the catalogues from the four redshift
will likely need to be included in future analyses. bins are merged before producing the aperture mass map and counting
We finally note that in contrast with the 2PCF, where shifts in ∆m the peaks; this reduces the impact of the pure noise peaks to highlight
affect all scales equally, the derivative presented above exhibit a more subtle effects that occur in the underlying matter distribution.
complicated structure, caused by the fact that the m-calibration affects The results are shown in Fig. 9 for all tomographic bins, as well
© 2020 RAS, MNRAS 000, 1–26

1.4
measurements, although with some dispersion in the results. For ex-
ample, T18 finds a 2.5σ detection of the signal from the DES-Y1
1.2 data, with AIA = 1.3+0.5 −0.6 , the KV450 analysis by Hildebrandt et al.
(2018) found an AIA value consistent both with unity and with zero,
1 while Hikage et al. (2019) and Hamana et al. (2020) prefer a val-
ues of AIA = 0.38 ± 0.70 and 0.91+0.27 −0.32 from the power spectrum and
0.8 2PCF analyses of the HSC-Y1 data, respectively. Similar variations
of the best fit IA amplitudes are observed in the KiDS-1000 cosmic
2 3 4
10 10 10 shear analysis by Asgari et al. (2020a), who found AIA values between
+0.29
0.26+0.42
−0.33 and 0.97−0.38 depending on the estimator. These AIA measure-
ments are not expected to agree perfectly given some differences in the
Figure 8. Effect of baryon feedback on the convergence power spectrum mea-
sured from the Magneticum simulations, assuming a fixed redshift for the modelling of the IA signal, in the physical scales that are probed and
sources. in the selection of galaxies, but they clearly illustrate the complexity
of the problem at hand. It is also worth pointing out that the IA nui-
sance parameter marginalisation is degenerate with the photo-z errors
as for the case where no tomography is applied. We see that the effect as well (see, for example, the discussion in Appendix C of Heymans
is generally under a percent for S/N < 4, and in the no-tomographic et al. 2020; Wright et al. 2020). For our 2PCF analyses, we use the
case reaches five percent for 4 < S/N < 5. The statistical precision is same two-parameters NLA model as in T18, with priors on AIA and η
also reported on these plots, which shows that the impact of baryons listed in Table 3.
is everywhere sub-dominant compared to the uncertainty on the GPR The impact of IA on peaks statistics has not been well studied in
emulator. This reinforces our choice of selecting S/N 6 4 to ensure the literature so far, and a percent level calibration will require a sig-
a measurement mostly free of uncertainty related to baryons, although nificant level of development beyond what has been done so far. Nev-
this suggests that we could push the upper limit to higher S/N in the ertheless we present here a first step in this direction, with a measure-
tomographic case without much contamination, and that modelling the ment of the effect based on in-painting galaxies with intrinsic shapes
effect could be relatively straight-forward. This follows a logic similar determined by properties of the dark matter haloes contained inside
to Huang et al. (2020), where the impact of baryonic feedback on the the lensing light-cones. Again, the amplitude of the IA signal mea-
DES-Y1 2PCF is modelled and captured with a PCA decomposition, sured from peaks is not expected to be the same as for the two-point
allowing them to include smaller angular scales in the data vector and shear correlation function, largely because the physical scales and the
increase the constraints. number of galaxies involved in each estimator calculations are differ-
We use the ratios presented in Fig. 9 to construct a multiplicative ent.
correction factor that is optionally applied to the data vector during the Within the NLA model of Bridle & King (2007) with AIA = 1.0
cosmology inference pipeline, from which we can estimate the impact for example, intrinsic alignments can modify by up to 40% the ξ± sig-
of baryon feedback on the recovered best-fit parameters, as modelled nal in the DES-Y1 bin 1, 20% in bin 2, but only about 5% in bins 3
with the Magneticum baryon model. We note that a similar approach and 4. Even if the effective AIA increases with redshift (see the dis-
is adopted in the context of a Stage-IV survey in Weiss et al. (2019) cussion in Appendix A of Robertson et al. 2020), the lensing kernels
and in Martinet et al, 2020 (in prep.), where the increased galaxy den- of the GI and II terms are suppressed compared with the GG signal.
sity and overall statistical precision accentuates the bias caused by the Considering the IA model of Fortuna et al. (2020), we note that the IA
baryons. effect is more significant at small physical scales, which are only well
resolved at low redshift. For these reasons, we choose to only model
the peaks IA signal in the low redshift bins, more precisely on the bin
4.4 Galaxy intrinsic alignment combinations 1-1, 1-2 and 2-2. This likely leaves a residual, unmod-
The intrinsic alignment of galaxies is an astrophysical systematic sig- elled GI term present in the bins 1-3, 1-4, 2-3 and 2-4, and we estimate
nal that mimics weak lensing measurements, and arises from the fact its impact with a GI-suppressed analysis variation in which we remove
that the intrinsic orientation of galaxies is not exactly random (see Kirk all cross-tomographic bins. As it will become clear in Sec. 5.2, this is
et al. 2015; Kiessling et al. 2015, for a review). Indeed, it has been currently a limiting factor in our data analysis, which we will seek to
shown in multiple hydrodynamical simulations that the formation of improve with a better IA model in the future.
galaxies, and thus their final shape and alignment, is affected by their Our IA model is inspired from the methods of Heymans et al.
environment, notably by the neighbouring large scale structures (Chis- (2006) and Joachimi et al. (2013b), which assign an intrinsic ellipticity
ari et al. 2015; Codis et al. 2015) and tidal fields (Catelan et al. 2001), int to the galaxies based on the shape of their host dark matter halo (we
and by a complex relation with their host haloes (Chisari et al. 2017). summarise this method and detail our implementation in Appendix B).
The observed galaxy shape is therefore a combination of the intrinsic The model requires both light-cone haloes and in-pasted galaxies, two
(I) and the shear (G) term, which both contribute to the measured weak intensive post-processing steps that have not yet been completed on
lensing signal. the cosmo-SLICS. It is therefore not possible at the moment to explore
Intrinsic alignments have been directly measured and con- this in a cosmology-dependent manner. Instead, we use 26 light-cones
strained, notably in the COSMOS galaxies by Joachimi et al. (2013a), from the KiDS-HOD galaxy sample described in HD18, which have
who detect the signal for early-type (e.g. red galaxies) but hardly any these properties and have been downsampled to closely match the N(z)
signal for late-type galaxies. The WiggleZ blue galaxies were also in the four DES-Y1 tomographic bins. These also cover 100 deg2 each,
found to be consistent with no IA in Mandelbaum et al. (2011), while and since we are only interested in the relative effect, we do not tile
a significant IA signal was found in the BOSS LOWZ sample (Singh the full footprint and work instead directly on the light-cone galaxy
et al. 2015). Johnston et al. (2019) found similar results from the KiDS, samples. The effect of masking is hence not included, but it should
SDSS and GAMA surveys, with a null detection from the blue galax- equally impact the measurements with and without IA in these mocks.
ies and a 9σ detection for an IA signal for red galaxies, consistent with For every galaxy, the model outputs int ; this quantity is then in-
a value of AIA = 3.18+0.47
−0.46 when interpreted in the NLA model. serted in Eq. (1) alongside a randomly rotated version, from which we
The IA signal has also been indirectly inferred from cosmic shear compute observed ellipticities with or without IA. We finally run our
© 2020 RAS, MNRAS 000, 1–26

1.05
11 22 33 44 12 13 14 23 24 34 no-tomo
1
0.95
0.9
1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5
Figure 9. Ratio between the number of peaks measured in the Magneticum light-cones with and without including the baryonic physics. Results are shown for
different tomographic bins, and for an aperture filter size of θap = 12.5 arcmin; other filters show a similar relative effect. The dashed lines represent the statistical
precision, also plotted in Fig. 6.
0 1 2 3 0 1 2 3 0 2 4 0 2 4 0 2 4 0 2 4
1.2
0
1
11 12 13 14 0.8
-0.1
1.2
1
-0.2
22 23 24 0.8
0 1.2
1
-0.1 33 34 0.8
1.2
-0.2 1
44 0.8
Figure 10. Effect of intrinsic alignment on the peak count statistic, measured
in the lowest two redshift bins of our dedicated IA training set. The error bars Figure 11. Impact of the N-body force resolution on the peak statistics, mea-
show the error on the mean. Full details of this model are provided in Appendix sured from the ratio in S/N between the SLICS-HR and the covariance training
B. set. The dashed-lines show the statistical error of the measurement. Deviations
for S/N 6 4 peaks are generally under 0.2σ.
peak finder on these catalogues, count the peaks as for the other train-
ing sets, and compare the measurements in Fig. 10. We observe an im-
the amount of small-scale structures vary significantly for wave vec-
portant (10-15%) suppression of the number of peaks with S/N > 3,
tors larger than k = 1.0 h−1 Mpc. Fitting functions such as Halofit
which clearly exceeds the statistical uncertainty in our measurement
(Smith et al. 2003; Takahashi et al. 2012), simulation-based emulators
and therefore needs to be accounted for. Moreover, in the top two pan-
(Heitmann et al. 2014, 2016; Nishimichi et al. 2019) and halo model
els, peaks at all S/N values are suppressed by a few percent. We note
calculations (Mead et al. 2015, 2020) provide fast predictions for two-
that the results from Fig. 10 align well with those found in K16 (see
point function statistics, but these disagree sometimes at the 5-10%
their figure C3), even though a direct comparison is impossible due to
level, depending on the scales, redshifts and cosmological parameters.
differences in the source distribution of the samples. When examined
Recent efforts approach the 1% accuracy on the matter power spectrum
with 2PCF measurements, we find that our IA prescription is bracketed
(Euclid Collaboration: Knabenhans et al. 2019; Mead et al. 2020), at
by the NLA model with AIA ∈ [1.0 − 2.0], providing a consistent but
least for a subset of the cosmological parameter space.
slightly larger IA signal than that preferred by the data (see Appendix
It is possible to protect the shear 2PCF measurement against
B).
residual non-linear uncertainty by rejecting angular scales in ξ± (ϑ)
We recognise that our simple IA model is unlikely to represent
where differences between these models affect the cosmological infer-
accurately the real physical effect, and at the moment has no free pa-
ence beyond some threshold19 . This is one of the main drivers, along
rameter, which means that we cannot yet marginalise over different IA
with baryon feedback, for the choice of angular scales in the T18 2PCF
strengths like we do for 2PCF. Further development on the IA model
analysis (and hence ours).
will be required to achieve this in the future. However, we can get a
For weak lensing probes beyond two-point statistics that are cal-
sense of the impact of IA by infusing the relative effect observed in
ibrated directly on simulations however, one must additionally under-
Fig. 10 in the model returned by the emulator, and record the deviation
stand the impact of small-scale smoothing caused by limits in the the
from the no-IA case on the best-fit parameters. It turns out that this
mass/force resolution of the N-body code. To assess this, we rely on
results in a ∼0.1σ shift on the cosmological parameters, hence we do
the SLICS-HR simulations introduced in Sec. 2.4, in which the force
not include it in the fiducial peak count pipeline.
resolution has been significantly augmented, thereby resolving scales
almost an order of magnitude smaller. A comparison between the peak
4.5 N-body force resolution function of the SLICS-HR measured on 10 realisations of the full DES-
A large fraction of the weak lensing signal receives a contribution from

small scales, which are difficult to model accurately. Even in simpli- 19 There is a caveat to this argument, see the discussion in Asgari et al. (2020b)
fied gravity-only scenarios, different methods and codes to estimate about the small k-scale power leaking into ξ± (ϑ) to some level at all ϑ.
© 2020 RAS, MNRAS 000, 1–26

observed signal. Furthermore, these regions of high galaxy densities

1.2
will have a larger blending rate, meaning that source galaxies behind
1
clusters are more likely to be missed, while residual errors in the pho-
0.8 tometric redshift can wrongly assign cluster members to background
1.2
sources. When combined, these effects can result in a mis-calibration
1
between the data and the simulations, which can be corrected with a
0.8
1.2 ‘boost factor’ (Mandelbaum et al. 2005).
1 Boost factors for peak statistics can be evaluated in different
0.8 ways. K16 estimate the fractional over-crowding and over-blending
1.2 rates in peaks of different S/N from the Balrog catalogue (Suchyta
1 et al. 2016), a separate image simulation that matches the DES-SV
0.8 n(z) and blending properties. These rates are computed as a function
10 0 10 1 10 0 10 1 10 0 10 1 of distance to peaks centres, and a correction factor is used to cor-
rect the peak function found in their cosmological training set as a
function of S/N. They found that by restricting their measurement to
Figure 12. Profile of the excess galaxy clustering, as a function of angular S/N < 4, the impact is minimal (a shift in S 8 of about 0.01) and could
separation. The different columns show the profiles for three different S/N be neglected.
bins, while the rows present the results for the four tomographic bins. These Shan et al. (2018) instead use the Radovich et al. (2017) cluster
profiles (red lines) are averaged over 10 independent survey realisations and catalogue that overlaps with the KiDS-450 survey, and evaluated the
enter in the boost factor correction (Eq. 13, see main text). boost factor from the excess source density around these massive ob-
jects. They found that the contamination to the peak function reaches
27% for peaks with S/N = 5, however it caps at less than 6% for
1.05 S/N < 4. This effect is about twice the size of the baryonic feed-
11 22 33 44 back, but acting in the same direction, e.g. suppressing the number of
high S/N peaks. If overlooked, this could be mis-interpreted by the
1 simulations and lead to a best-fit inference with a S 8 that is too low.
We account for source-lens coupling by estimating the effect in
0.95 the DES-Y1 data and recalibrating our measurements, leaving the sim-
1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 ulations unchanged. We first extract the clustering of source galaxies
along the line-of-sight of peaks identified in the data, as a function of
their peak height. More precisely, we estimate widata j
(ϑ, S/N) for each
Figure 13. Impact of the boost factor on the peak statistics, measured from the combination of auto- and cross-tomographic bins ‘i j0 . These are next
ratio in Npeaks between the boost-corrected and the original peak count. The compared with a similar measurement carried on ten survey realisa-
dashed-lines show the statistical error of the measurement. Similar results are tions sampled from the Cosmology training set at the fiducial cosmol-
obtained for the cross-tomographic bins.
ogy, and the ratio of the two reveal clustering profiles:
 sim 
widata
j
(ϑ, S/N)  hNpeaks i 
ρ (ϑ, S/N) ≡ i j
ij
×  data  , (12)
Y1 footprint and that of our main Covariance training set is presented hwsim (ϑ, S/N)i Npeaks
in Fig. 11; it reveals slight differences for the S/N peaks we are prob-
which are shown in Fig. 12 for a sample of tomographic bins and S/N
ing, but strong deviations are observed for peaks with S/N > 4. These
bins. The brackets refer to the mean measurement in the above ex-
objects are rare, which explains the increased noise in the ratio to- data/sim
pression, and the right term involving Npeaks normalises the profiles.
wards S/N = 5, but the trend is clear: there is an overall shortcoming
It is clear from this figure that the largest peaks are generally more
of large S/N objects in the training set compared to the SLICS-HR,
severely affected by this boost factor correction, however the size of
which justifies our choice of restricting the data vector to S/N 6 4.
effect varies across redshift in a non-trivial manner. For example, the
Upon closer examination, the average size of the small deviations seen
fourth tomographic is less affected than the second or the third. In ab-
in that range correspond to no more than 20% of the statistical error,
sence of source-lens coupling these profiles would be flat. The excess
and are never higher than 0.5σstat . If we wanted to include these peaks
of galaxies in these profiles are for the most part cluster members;
in a future analysis, we would possibly need a new generation of train-
their shapes are therefore not sheared by the foreground matter over-
ing sets (for cosmology, and possibly covariance) with an increased
density and only dilute the lensing signal. We compensate for this by
mass resolution.
up-weighting the shear signal following the profile, which is most ef-
Just as for IA and baryons, we investigate the impact of mass res-
ficiently done by modifying the filter Q(r) in identified peaks as:
olution by extracting a correction factor from the black curves shown
in Fig. 11, which we consequently apply to the signal during the cos- Q(r) = Qorig (r) × ρi j (r, S/N), (13)
mology inference.
and re-evaluate the peak height with Eq. (7). Fig. 13 shows the ra-
tio between the corrected and original peak function in the four auto-
tomographic bins (similar results are obtained with the cross-bins).
4.6 Source-lens clustering
The effect if generally small, however it approaches the 1σ level in
One of the key differences between the mock DES-Y1 simulations some isolated data elements. The boost factor is included in our fidu-
from our Cosmology training set and the DES-Y1 data is the presence cial analysis, and has been applied to the data points shown in Fig.
of source-lens coupling. In real data, the source density is not homoge- 4. Following the method used in the previous sections, we isolate its
neous and in fact increases around foreground clusters. As explained impact on our cosmological inference by optionally removing this cor-
in Appendix C of K16, this introduces a coupling between the peak rection factor.
positions and the amplitude of the measured shear relative to the ex- Fig. 14 summarises all the correction factors we can include in
pected shear – the sources that are associated with the cluster dilute the our cosmology inference, from IA, mass resolution, baryons and no-
© 2020 RAS, MNRAS 000, 1–26

1.05 11 22 33 44 12 13 14 23 24 34
0.95
stat error baryons
0.9 IA finite mass
no-boost GPR
0.85
Figure 14. Impact on the peak statistics of various sources of systematic uncertainty (IA, mass resolution, baryons and boost), presented as a ratio between the
measurement on mocks with and without the effect. These are used as optional correction factor applied to the model in the cosmology inference, as described in Sec.
5.2. Also shown is the scatter in the GPR cross-validation test, as well as the statistical precision on the measurements.
Table 3. Priors used in the likelihood sampling. The ranges for the four cosmo- to the systematic effects described above. In this section we introduce
logical parameters are determined by the cosmo-SLICS simulations, the pho- these effects, and justify our choice to neglect them in this study.
tometric redshift ranges and priors are taken from the DIR errors found in J20, Our simulated light-cones are constructed with mass planes of
while those of the shear calibration and intrinsic alignments are taken from T18. constant thickness set to 256.5 h−1 Mpc (see Sec. 2.6.1), a choice which
Gaussian priors are characterised by a mean and a standard deviation (µ, σ). has an effect on the reconstructed lensing planes compared to a con-
struction made of hundreds of finer shells. This has been recently quan-
Parameter range prior tified in Zorrilla Matilla et al. (2020), where it is shown that the impact
Cosmology
on peaks with S/N 6 4 is below the one percent level, regardless of
Ωm [0.1, 0.55] Flat the thickness and of the source redshift.
σ8 [0.53, 1.3] Flat Correlations between the mass planes in our light-cones are ex-
h [0.6, 0.82] Flat plicitly suppressed by randomly shifting and rotating the mass sheets,
w0 [-2.0, -0.5] Flat which breaks the long line-of-sight correlations that exist in the data.
Nuisance
It was shown in Takahashi et al. (2017, see their Appendix B) that this
∆z1 × 102 [-10, 10] G(0, 0.8) affects the projected power-spectrum, reducing the power at intermedi-
∆z2 × 102 [-10, 10] G(0, 1.4) ate scales by a few percent on the sheets, however, the lensing kernels
∆z3 × 102 [-10, 10] G(0, 1.1) project an even larger volume and mixes these scales, which makes our
∆z4 × 102 [-10, 10] G(0, 0.9) measurements relatively immune to this.
∆mi × 102 [-10, 10] G(1.2, 1.3) The lensing plane construction has been carried out under the
IA Born approximation (see HD18), whereby light bundles record the
AIA [-5, 5] Flat convergence and shear along straight lines and ignore the deflection
η [-5, 5] Flat angles in this calculation. It has been shown (Hilbert et al. 2020) that
the difference between these methods induces variations smaller than
0.2% on the lensing power spectrum up to scales of ` = 2 × 104 . It is
boost. Our fiducial analysis includes none of them, since their overall therefore reasonable to expect that Born approximation plays a minor
impact is relatively minor and many of these effects act in opposite role in the peak statistics as well, however Castro et al. (2018) find that
directions. the impact on the PDF of the lensing maps is of a few percent. We
ignore this effect at the moment, but it will need to be revisited in the
future.
4.7 Cosmology inference Finite box effects are also known to plague the estimation of
2PCF and of their covariance matrix (Harnois-Déraps & van Waerbeke
We test our cosmology inference pipelines on mock data vectors taken 2015), being sensitive to Fourier modes larger than the simulation box.
to be the mean value from the Covariance training set, providing a This has an impact on the ξ± covariance matrix estimated from our
measurement that is almost noise-free. We present our results in Ap- Covariance training set, however it has a minor impact on the cos-
pendix A, notably in Fig. A1. This exercise reveals a high degree of mological contour, as shown by the good match between our analysis
similarity in the constraints between the two-point functions, the peak on the mocks and that of T18. Furthermore, since the peak statistics
statistics and the fiducial T18 analyses, despite major differences in measure quantities in local apertures, it is not sensitive to these large
our covariance matrix estimation techniques and in the observation scales, and hence are protected against this. Similar conclusions can
data vector (DES-Y1 data vs SLICS, 2PCF vs peaks). The best-fit be drawn regarding the incomplete account of the super-sample co-
cosmological parameters are also consistent with the input at the 1σ variance (Li et al. 2014, SSC hereafter) captured by our simulations:
level, with no noticeable shift between the probes, and the sizes of all HD19 found that the SLICS light-cones capture more than 75% of this
contours closely match that of the DES-Y1 analysis, two properties SSC term, yielding two-dimensional constraints on (Ωm , σ8 , h, w0 ) that
that respectively validate our cosmology calibration and our covari- match very well those of an analytical covariance matrix, with differ-
ance matrix. ences of less than 10% in the area of the 2σ regions. Given the sup-
pressed sensitivity of the peak statistics to these large-scale modes, we
expect the residual missing SSC contribution to play a minor role on
4.8 Others sources of uncertainty the full uncertainty, although this may need to be better quantified in
Our method relies on a certain number of well-justified approxima- the future.
tions, which could potentially contribute to the error budget in addition It has been recently shown that the depth variability across a lens-
© 2020 RAS, MNRAS 000, 1–26

Table 4. Properties of the different pipelines discussed in this paper. J20 CDM
T18 CDM
2PCF CDM
prior on Cov. n(z)
pipeline σ8 /As Matrix method
This work σ8 ∈ [0.53 − 1.3] SLICS DIR

0.9
DIR-wCDM As ∈ [0.5 − 5.0] × 10−9 SLICS DIR
0.8
T18 As ∈ [0.5 − 5.0] × 10−9 Analytic Stacked PDF
0.7
h
J20 ln 1010 As ∈ [1.5 − 5.0] Analytic DIR
0.6
1.2
Table 5. Cosmological pipeline comparison. The values used in the T18 com-
1.0
parison are taken from their Table 3, using a fixed neutrino mass density. Details
8
on the different pipelines are summarised in Table 4. Most posteriors on h and 0.8
w0 are prior-limited, so no constraints are reported here. The same applies to 0.6
Ωm in many cases. Tests on the mock data are presented in Appendix A.
0.8
pipeline S8 Ωm
0.7
S8
Peaks 0.780+0.019
−0.056
- 0.6
Fiducial 2PCF 0.753+0.043 +0.033
0.254−0.056
−0.043 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.6 0.8 1.0 1.2 0.6 0.7 0.8
Joint 0.766+0.033 -
−0.038 m h 8 S8
2PCF (T18, wCDM) 0.797+0.037 +0.079
0.290−0.051
−0.037
2PCF (DIR-wCDM) 0.752+0.042 +0.035
0.264−0.054 Figure 15. Constraints on the ΛCDM cosmological parameters inferred from
−0.037
2PCF (ΛCDM) 0.761+0.027
−0.027
+0.031
0.272−0.056 the 2PCF, obtained from the DES-Y1 cosmology inference pipeline, our
2PCF (T18, ΛCDM) 0.792+0.032
−0.021
+0.038
0.304−0.062 simulation-based covariance matrix and assuming the DIR n(z) (blue). These
2PCF (J20, ΛCDM) +0.036
0.765−0.031 +0.041
0.252−0.086 results are compared to the ΛCDM constraints from T18 (in red) based on the
Variations
Peaks (cross-tomo, with IA) 0.735+0.024
−0.032 - same modelling and likelihood sampling strategy, and to those of J20 (grey),
Peaks (cross-tomo, with baryons) 0.750+0.026
−0.031 - which also use the DIR redshift distribution but adopt different modelling, prior
Peaks (cross-tomo, with SLICS-HR) 0.734+0.025
−0.032 - ranges and likelihood sampling choices.
Peaks (cross-tomo, no-boost) 0.736+0.025
−0.032 -
Peaks (cross-tomo) 0.737+0.027
−0.031 -
Joint (cross-tomo) 0.743+0.024
−0.024 - 2PCF
T18 wCDM
Peaks (cross-tomo, no syst) 0.787+0.024
−0.024
+0.054
0.325−0.067 DIR-wCDM
Mocks Peaks (cross-tomo) 0.776+0.045
−0.045
+0.048
0.297−0.066
2PCF (FID) +0.042
0.772−0.042 +0.049
0.314−0.070
ing survey can impact the cosmic signal and variance (Heydenreich 0.8
S8
et al. 2020b; Joachimi et al. 2020). This will need to be the subject of
0.7
future investigations. Given the findings of Joachimi et al. (2020), we
expect the impact of unmodelled variable depth to be negligible, given 0.6
the statistical power of DES-Y1. 0.9
0.8
0.7
h
5 RESULTS 0.6
0.5
We present in this section the results from our cosmological inference 0.5
analyses, beginning with the 2PCF and the peak statistics pipelines. We
1.0
next discuss our tests on the importance of various systematic effects,
w0
1.5
before introducing the results from our joint [ξ± ; Npeaks ] analysis.
2.0
0.2 0.3 0.4 0.5 0.7 0.8 0.6 0.7 0.8 0.9 2.0 1.5 1.0 0.5
5.1 2PCF m S8 h w0
In this section we present the results from a series of cosmology infer- Figure 16. Constraints on the wCDM cosmological parameters inferred from
ence analyses carried out on the DES-Y1 data. We first report in Fig. the 2PCF with our fiducial pipeline (grey), from the T18 wCDM analysis (red)
15 the constraint on ΛCDM parameters obtained from our 2PCF mea- and from an intermediate pipeline, the DIR-wCDM, which uses the T18 pa-
surement, over-plotted with those from T1820 (red) and J20 (grey). Our rameter sampling and prior ranges on our tiled data with our simulation-based
constraints (in blue) assume the DES-DIR n(z), it uses our simulation- covariance matrix, assuming the DES-DIR n(z) (blue). Note that priors on h and
based covariance matrix, we marginalise over the 10 nuisance param- w0 are significantly tighter in our fiducial pipeline, and that our fiducial 2PCF
analysis hits the priors (shown with the dashed lines) on these two parameters.
eters listed in Table 3, but the cosmological sampling follows that of
T18. Another difference: as described in Sec. 3.1, our ξ± measurements
20 In this comparison, we used the values listed in their Table 3 using a fixed
neutrino mass density, which better match our fiducial pipeline.
© 2020 RAS, MNRAS 000, 1–26

are obtained from the weighted mean ξ± obtained over the 19 tiles. As
demonstrated in J20, the net effect of changing from the fiducial DES-
0.9 fiducial
Y1 n(z) to the DIR n(z) is to shift the amplitude of the modelled 2PCF
signal upwards, which translates into lower best-fit S 8 values. This can
IA
be seen by comparing the one-dimensional posteriors shown with the baryons
red and grey lines in the bottom right panel. When analysed this way, SLICS-HR
we obtain:
no-boost
S 82PCF,Λ = 0.761+0.027
−0.027 . 0.8 no-cross
S8
The priors and the parameter sampling in J20 are significantly differ-
ent from T18 and are responsible for some of the differences between
the blue and grey curves, notably the sharp edges in the h posterior,
and the more elongated contour in the [σ8 − Ωm ] plane. All parameter
constraints are summarised in Table 5. Replacing the ξ± data extracted
from the tiles with those measured on the full footprint results in a neg- 0.7
ligible change, with S 8 = 0.762 ± 0.026, thereby validating our mosaic
methodology. Additionally, the remarkable resemblance between our
confidence interval and those of T18 (the fractional errors on S 8 and
Ωm agree to within 0.002 and 0.005, respectively) speak of the accu- 0.2 0.3 0.4 0.5 0.6
racy of our simulation-based covariance matrix. The constraints on the m
nuisance parameters are mostly prior-dominated. We provide a more
complete comparison between T18, J20 and our 2PCF ΛCDM analy-
Figure 17. Impact of the different correction factors on the constraint from the
ses in Appendix C. peaks statistics. In most cases the likelihood hits the upper prior edge on Ωm
We next compare in Fig. 16 the constraints on the four wCDM within 2σ, as marked by the vertical line, which prevents us from reporting
parameters inferred from our 2PCF measurement (in grey), to the T18 constraints on that parameter.
wCDM results (in red). As explained previously, there are multiple
difference between these two pipelines, which we can dissect here.
We show (in blue) an intermediate pipeline, labelled the DIR-wCDM,
2PCF
which uses the T18 parameter sampling and prior ranges, but assume Joint
the DIR n(z), and uses our measurement on the tiled data and our Peaks
simulation-based covariance matrix. Table 4 summarises the differ-
ences between these pipelines. By construction, differences between
the blue and the grey contours are caused exclusively by the likelihood
sampling strategy: the former uses the T18 priors and samples As over 0.85
a flat prior, whereas the latter uses those listed in Table 3 and samples 0.80
0.75
S8
σ8 . In contrast, red and blue curves share the signal modelling as well 0.70
as the parameter sampling, but differ in the n(z) (which shifts down the 0.65
best S 8 and Ωm values, clearly visible in Fig. 16), and in the covariance
matrix (which weights slightly differently the various elements of the 0.8
compressed statistics and therefore affects the size and shape of the
0.7
h
contours).
In our fiducial analysis the MCMC sampler hits the priors on h 0.6
and w0 , which are limited by the range of values probed by our Cos-
mology training set. We note, however, that this should have no impact 0.5
on our analysis due to the low sensitivity of lensing to these particular 1.0
w0
parameters, and that we marginalise over these anyway. Consequently, 1.5

we are unable to report constraints on h and w0 in our 2PCF pipeline 2.0
with the current data. The best-fit parameters are reported in Table 5, 0.2 0.4 0.6 0.7 0.8 0.6 0.7 0.8 2 1
notably: m S8 h w0
S 82PCF,w = 0.753+0.043
−0.043 , Figure 18. Constraints on the four wCDM cosmological parameters obtained
from 2PCF (grey), peaks (blue) and from the joint analysis (red). The dashed
which is consistent with, but has a larger uncertainty than S 8T18,w =
lines indicate the prior ranges on Ωm , h and w0 , which are in most cases too
0.797+0.037
−0.037 reported in T18. The overall precision on the matter den- narrow to provide meaningful constraints on these parameters (see main text
sity is similar to that of the ΛCDM analysis (we measure Ωm = for exceptions).
0.254+0.033
−0.056 ), while the uncertainty on S 8 increases significantly, as ex-
pected from opening the parameter space. When adopting the DIR-
wCDM pipeline, we obtain:
size of the reported error bars (see the discussion on informative priors
S 8DIR−wCDM = 0.752+0.042
−0.037 , in Joachimi et al. 2020).
which best inferred value aligns with our 2PCF analysis, with error
bars slightly tighter. The fact that the T18 constraints (4.6% on S 8 )
5.2 Peaks
are tighter than these (5.3%), despite having larger uncertainty in their
redshift distributions, indicates that the T18 priors and sampling strat- We report in Fig. 17 the results of our peak count analyses for models
egy are informative about S 8 to some level, artificially decreasing the in which systematics are infused. All data presented from now on are
© 2020 RAS, MNRAS 000, 1–26

obtained from the 12.5 arcmin filter; we investigated other choices of case, the effect is marginal. We find
θap that yielded slightly less constraining results and hence we dropped +0.025
them from the analysis.
peaks
S 8,no−boost = 0.736−0.032 , (15)
Before examining differences between the various cases, we first a negligible change with respect to the fiducial pipeline.
note that, for all of them, the inferred matter density values are un- Overall, when including all redshift bins, we observe that the
expectedly high, with Ωm & 0.4 at the 1σ level. This is in tension baryons and IA cause the largest shifts. Excluding the cross-redshift
both with the 2PCF measurements and with external probes such as terms has the strongest impact, as it significantly offsets the best-fit
the Planck measurement of Ωm = 0.3153 ± 0.0073 (Planck Collab- value and degrades the constraining power. Nevertheless, the contours
oration et al. 2018, see also a quantitative assessment of these ten- are fully consistent with the fiducial case, and this choice is necessary
sions in Appendix A3). In the case where the cross-redshift bins are until an improved high-redshift IA model can be incorporated.
excluded however, the tension is significantly reduced. Interestingly, The fiducial constraints from the peak statistics are presented by
this is also the only case where the GI part of the intrinsic alignment the blue contours in Fig. 18 and compared to the 2PCF (in grey, same
signal is suppressed, which suggests that unmodelled IA systematics as in Fig. 16). The marginal constraints on S 8 are:
could be artificially pushing the likelihood analysis towards high val-
ues of the matter density. We therefore adopt a conservative approach
peaks
S 8,fid = 0.780+0.019
−0.056 .
and remove the cross-redshift bins from our fiducial analysis. We still These values are in excellent agreement with those reported by J20,
include them in one of our analysis variations, with the caveat that they T18 and with our 2PCF method, with marginal constraints on S 8 being
might be contaminated by an unaccounted secondary signal. consistent to within 1σ in all cases. More importantly, the peak count
In order to provide the most accurate account of the other system- analysis improves this fractional precision by 10% compared to our
atic effects, we include the cross-redshift bins in their measurements, fiducial 2PCF pipeline, providing a 4.8% measurement on S 8 . Includ-
shining light on their impact with the highest precision available. Start- ing the cross-tomographic bins results in S 8peaks = 0.737+0.027
−0.031 , which
ing with the IA, and recalling that our alignment model applies to to- improves the constraining power to a ∼3.9% measurement (see the
mographic bins 1-1, 1-2 and 2-2 only with no free parameter, we in- pipeline ‘cross-tomo’ in Table 5). Although these are excluded from
clude the systematic effect as an optional correction to the measured the fiducial analysis due to unaccounted for systematics, this demon-
data vector (see Sec. 4.4 and Fig. 14). The cosmology inference re- strates that even with the current noise level and in presence of sys-
sults, shown by the black dashed lines in Fig. 17, indicate that our tematics, some 20% additional information on S 8 can be recovered by
low-redshift IA model has a relatively mild impact, affecting the best- including these bins compared to the fiducial case, potentially bringing
fit S 8 value by less than 0.1σ. Notably, we have: the gain to 30% over the 2PCF analysis. This aligns with the findings
= 0.735+0.024 of Martinet et al. (2020), which observe a 50% gain in a systematics-
−0.032 ,
peaks
S 8,IA
free Stage-IV weak lensing survey setup, using five tomographic bins
which is almost unchanged compared to the baseline analysis: and including cross-terms for all possible combinations of slices (i.e.
larger than pairs). We finally note that the posterior on Ωm overlaps
peaks
S 8,cross−tomo = 0.737+0.027
−0.031 . with the upper prior edge within 2σ, while those on h and w0 are com-
We note that the best-fit Ωm increases by 0.02. This is consistent with pletely prior-dominated, hence we do not report constraints on these
the expectation that IA overall suppresses the lensing signal; therefore, parameters21 .
given a measurement (2PCF or peaks), the inferred S 8 and/or Ωm val-
ues increase with stronger IA model.
As mentioned earlier, the GI term, that we are currently unable 5.3 2PCF+Peaks
to model accurately, is most likely to impact the cross-tomographic The fact that both 2PCF and peaks count methods individually prefer
bin measurements. The fact that the best-fit S 8 shifts up by as much different values for the cosmological parameters is expected, given that
as 0.037 (and Ωm shifts down by 0.02) when removing the cross- they both probe slightly different physical scales, and react to the DES-
tomographic bins suggests that unmodelled systematics affect these Y1 shape noise in completely different ways. This is further supported
bins differently than the auto-tomographic bins. The GI contribution by our measurements on the mean of the SLICS mocks presented in
does exactly this, and a similar 1σ effect is detected in T18 (see their Appendix A, which correspond to noise-free data and show an excel-
Fig. 9). It is therefore plausible that unmodelled IA could be leading lent agreement in the best-fit values. We note that similar results have
to the observed preference for a high Ωm values when including the been observed in the literature already, notably when comparing the
cross-tomographic bins, however we postpone a more robust analysis inferred cosmology from the 2PCF and power spectra analyses of the
of IA to future work. HSC-Y1 data (Hamana et al. 2020; Hikage et al. 2019) and of the
Adopting a similar strategy, we apply the baryon correction fac- KiDS-1000 data (Asgari et al. 2020a). Some fluctuations are expected
tor (see Sec. 4.3, and shown in Fig. 14), and re-run our cosmological and statistically consistent, as demonstrated in Joachimi et al. (2020).
inference pipeline. The impact is stronger than the IA, yielding: We show with mock surveys in Appendix A3 that the observed dif-
ference of ∆Ωm = 0.2 is frequent, and that the tension between the
peaks
S 8,baryons = 0.750+0.026
−0.031 ,
two methods is low, enabling us to carry out a joint 2PCF+peaks like-
a 0.5σ shift in S 8 visualised by the red dashed lines in Fig. 17. When lihood analysis. For consistency, our fiducial case again excludes the
correcting the signal for a possible bias due to limits in mass resolution cross-tomographic redshift bins in the peak count data, and results in
of the Cosmology training sample, the values of S 8 is slightly reduced a best-fit value of:
as expected, with:
S 8joint = 0.766+0.033
−0.038 ,
peaks
S 8,SLICS−HR = 0.734+0.025
−0.032 . (14)
This demonstrates that our choice of S/N bins is unaffected by the 21 We expect this to change with the upcoming Stage-IV lensing surveys, as
known small-scales inaccuracies inherent to the N-body runs. Martinet et al. (2020) has shown that peak statistics could provide a 6% con-
Removing the boost factor correction can generally push the in- straint on Ωm and a 13% constraint on w0 from 100 deg2 of Euclid-like mocks
ferred cosmology to models with lower clustering, however in our built from the same SLICS and cosmo-SLICS suites.
© 2020 RAS, MNRAS 000, 1–26

with contours presented in red in Fig. 18. The best-fit value in this case measurements offer excellent fits. Similarly, the goodness-of-fit of the
is consistent with both the 2PCF and the peak statistics, with a ∼20% joint analysis with cross-terms varies between 1.2 and 1.7, depending
smaller fractional uncertainty on S 8 . The joint analysis achieves 4.6% on how it is evaluated. This only illustrates the complexity involved in
precision, fast approaching the 3.8% precision achieved by the DES- estimating a goodness-of-fit, and that the rule-of-thumb of χred ∼ 1 is
Y1 3×2-pts analysis (Abbott et al. 2018a). We note that the same 20% not always adequate nor sufficient to determine the true quality of a fit.
increase in precision is reported in M18 when comparing peaks and
2PCF. However, their treatment of the systematics is simpler in com-
parison to this study, which likely affects their real precision gain. The
analysis variant in which the cross-tomographic terms are included 6 CONCLUSIONS
yields a promising 3.2% measurement, with:
We analyse the DES-Y1 cosmic shear data with a cosmology infer-
joint
S 8,cross−tomo = 0.743+0.024
−0.024 , ence pipeline exclusively calibrated on suites of wCDM weak lensing
numerical N-body simulations. Our method is general and can be di-
however a better modelling of systematics is required before we can
rectly used with many non-Gaussian probes, however we opted for the
exploit this configuration.
tomographic lensing peak count statistics. Our pipeline interfaces with
We repeated the full cosmology joint-probe inference with two
the two-point functions via the public inference code cosmoSIS, which
other aperture filter sizes, and obtained similar constraining strength,
allows us to conduct a joint [ξ± ; Npeaks ] analysis. We model the peak
which leads us to the conclusion that given the current level of noise
statistics signal with the cosmo-SLICS, the covariance with the SLICS
in the data, the choice of aperture filter size does not make a large
and investigate the key cosmic shear systematics either by infusing
difference, at least for the scales we tested. It is also clear that under
them in our training set, or by extracting their impact from tailored
these conditions a multi-scale analysis would bring no additional in-
external mock data. Notably, the impact of baryons is assessed with
formation. M18 found a mild improvement when combining the 12.50
the Magneticum hydrodynamical simulations, the effect of finite parti-
and the 15.00 filters, but that was in the absence of tomography. In that
cle mass is quantified from high-resolution light-cones, while the im-
case, the noise levels are suppressed at the cost of losing the redshift
pact of intrinsic alignment is investigated by infusing intrinsic shapes
information.
to mock galaxies following a physical model based on the shape of
The dark energy equation of state can be constrained from the
dark matter haloes. Source-lens clustering is also studied by compar-
joint analysis, which has sufficient constraining power to not be prior-
ing galaxy excess in peaks of different height, and found to have a
dominated22 . We obtain:
negligible impact on our results given our peak selection criteria.
wjoint
0 > −1.55 We validate our method on mock data vectors and against the
DES-Y1 ξ± analyses of T18 and J20, which we reproduce well given
at 1σ, which is consistent with ΛCDM cosmology but not yet com- differences in our pipelines. We identify a residual unknown system-
petitive with T18, who find wT18 0 = −0.92+0.35
−0.40 . As shown in Martinet atic in the cross-tomographic redshift bins of the peaks data, which
et al. (2020), the cross-tomographic data points will greatly help with appears to shift the inferred cosmological parameters towards high Ωm
this measurement once properly calibrated. values. We identified the intrinsic alignment GI term as one possi-
The overall goodness-of-fit is evaluated with the reduced χ2 , de- ble and likely cause of this effect, and in absence of an accurate IA
fined as χred ≡ χ2 /νeff , νeff being the effective number of degrees of modelling, we decided to remove these redshift bins from our fiducial
freedom. This quantity is often estimated from the difference between analysis.
the number of data points and the number of free parameters in the We perform a joint MCMC analysis and set constraints on S 8 ,
model. However this does not exactly hold in the more accurate con- finding a ∼20% gain in precision compared to the correlation func-
text where these model parameters are not totally free but sampled tion analysis. The joint analysis yields a 4.6% measurement, with
from priors that often restrict the search, and sometimes degenerate. S 8joint = 0.766+0.033
−0.038 , approaching the 3.8% measurement reported by
For example, the KiDS-1000 cosmic shear analysis presented in Asgari the DES-Y1 3×2-pts analysis (Abbott et al. 2018a), and closely ap-
et al. (2020a) samples 12 parameters in their cosmological inference proaching the 2.9% S 8 measurement of Asgari et al. (2020a) with the
pipeline, but a careful study of the sampling reveals that the effective fourth KiDS data release. We show that after an improved calibration
number of free parameter is closer to 4.5 (Joachimi et al. 2020). of the cross-tomographic terms, the peak count method can achieve
In the current analysis, the size of the data vector varies from 48 3.9% on S 8 , and the joint analysis could reach 3.2%, one of the tight-
(peaks fiducial) to 377 (peaks + cross-tomography + 2PCF). The χ2 est measurements of this parameter so far. One possible caveat is that
for the three fiducial analysis pipelines, at their best-fit cosmologies, we include no IA modelling in our cross-tomographic analysis, and
are respectively χ2peaks = 71, χ22pcf = 209 and χ2joint = 463. If, follow- that its inclusion could possibly degrade our constraining power.
ing Joachimi et al. (2020), we also use 4.5 free parameters, this leads Our analysis pipeline builds from the infrastructure presented in
to goodness-of-fit values of χred = 1.60, 0.82 and 1.54, respectively. M18, and is inspired by many aspects of the K16 data analysis. We
The two values for the peak count and joint analysis data are compar- have significantly upgraded the underlying simulation support, we in-
atively large, indicating possible residual tensions that we could not clude a better treatment of the photometric and shear calibration un-
clearly identify, or a mis-estimation of the likelihood. For example, we certainties, we improve the inference method with a full integration
note that if we replace our t-distribution likelihood (Eq. 10) by a mul- within a MCMC sampler, and we demonstrate the robustness of our
tivariate Gaussian likelihood in which we debias the inverse covari- results to uncertainty on baryonic feedback, to intrinsic alignments,
ance matrix with the Hartlap et al. (2007) factor, we recover reduced to source-lens clustering and to limitations from the finite mass res-
χ2 values of 0.988, 0.498 and 0.928, where now the two problematic olution of the N-body runs, which are lacking in previous analyses.
Additionally, this paper is the first tomographic peak counts data anal-
22 ysis, and we confirm that including cross-redshift bins reinforces the
We adopt the criteria described in Appendix A of Asgari et al. (2020a) to
decide whether we can report a constrain: if the value of posterior at any edge
constraints, once secondary signals such as IA are properly calibrated.
of the prior is higher that 0.135 times its maximum value, the measurement is For the peaks statistics to remain competitive with the 2PCF, fur-
prior-dominated and the constraints should not be reported. Our measurement ther development will be required in the modelling of baryon feedback.
of w0 from the joint analysis marginally satisfies this criteria on the lower end For example, a 20% improvement on S 8 is obtained from the DES-
only, with posterior values of 0.134 at the lower prior edge. Y1 3x2pt data when modelling and marginalising over the impact of
© 2020 RAS, MNRAS 000, 1–26

baryons, as it enables to include additional small-scale elements in the a CNES Fellowship, TC is supported by the INFN INDARK PD51
ξ± data vector (Huang et al. 2020). A similar gain is observed in the grant and by the PRIN-MIUR 2015 W7KAWC grant, and KD by
cosmic shear-only DES-Y1 re-analysis by Asgari et al. (2020b), where the Deutsche Forschungsgemeinschaft (DFG, German Research
clean small scales are included via the COSEBIs estimators. Equiva- Foundation) under Germany’s Excellence Strategy – EXC-2094
lent procedures with peaks statistics must be investigated, which would – 3907833. BG acknowledges the support of the Royal Society
allow us to push back the S/N 6 4 limit. The ‘baryonification’ method through an Enhancement Award (RGF/EA/181006). We also ac-
described in Schneider et al. (2019) is one possible avenue to achieve knowledge support from the European Research Council under
this, as well as the direct training on a variety of hydrodynamical light- grant numbers 647112 (JHD, BG, CH & QX) and 770935 (HH),
cones, although the latter is a more expensive task and introduces sev- as well as the Deutsche Forschungsgemeinschaft (HH, Heisenberg
eral other uncertainties, e.g. how one changes the sub-grid physics for grant Hi 1495/5-1). CH further acknowledges support from the Max
different cosmologies. The goal here would be to construct a response Planck Society and the Alexander von Humboldt Foundation in the
model in order to include the effect in a MCMC inference analysis. framework of the Max Planck-Humboldt Research Award endowed
Intrinsic alignment of galaxies is the second topic where further by the Federal Ministry of Education and Research. Computations
development is critical, and improving the IA modelling is mandatory for the N-body simulations were enabled by Compute Ontario
for future analyses. The model we adopt here is too simple, and we (www.computeontario.ca), Westgrid (www.westgrid.ca) and Compute
demonstrate that IA likely impacts the inferred cosmology at the same Canada (www.computecanada.ca).
level as the baryons and potentially up to 1σ level, making it one of
the primary limiting factor in our analysis. A better approach would This project used public archival data from the Dark Energy
involve a suite of dedicated training samples where the modelling and Survey (DES). Funding for the DES Projects has been provided by the
the levels of IA can be adjusted, mimicking the varying (AIA , η) param- U.S. Department of Energy, the U.S. National Science Foundation, the
eters that control the amplitude of the secondary signal within the NLA Ministry of Science and Education of Spain, the Science and Technol-
model. One could also relate the simulated intrinsic galaxy shapes with ogy FacilitiesCouncil of the United Kingdom, the Higher Education
the halo-model based approach of Fortuna et al. (2020) directly from Funding Council for England, the National Center for Supercomput-
the HOD prescription. Even better, these approaches could be linked, ing Applications at the University of Illinois at Urbana-Champaign,
allowing us to marginalise coherently over a unique set of common the Kavli Institute of Cosmological Physics at the University of
astrophysical parameters. Chicago, the Center for Cosmology and Astro-Particle Physics at the
All of the improvements presented here further close the gap be- Ohio State University, the Mitchell Institute for Fundamental Physics
tween analytical and simulation-based cosmological inference tech- and Astronomy at Texas A&M University, Financiadora de Estudos
niques. While the former is preferred in cosmic shear analyses per- e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do
formed by the large weak lensing collaborations, the latter is required Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento
by most measurement methods based on non-Gaussian statistics, for Cientı́fico e Tecnológico and the Ministério da Ciência, Tecnologia
which a theoretical model of the signal and of the covariance matrix is e Inovação, the Deutsche Forschungsgemeinschaft, and the Collab-
generally not available. Given that the performance of non-Gaussian orating Institutions in the Dark Energy Survey. The Collaborating
estimators is now clearly shown to significantly exceed that of two- Institutions are Argonne National Laboratory, the University of
point correlation functions, and that this trend will intensify in future California at Santa Cruz, the University of Cambridge, Centro de In-
surveys (Li et al. 2019; Zürcher et al. 2020; Martinet et al. 2020), we vestigaciones Energéticas, Medioambientales y Tecnológicas-Madrid,
strongly advocate for a co-development of both approaches. the University of Chicago, University College London, the DES-
Most of the work that has been put into the development of our Brazil Consortium, the University of Edinburgh, the Eidgenössische
method can now be directly exploited with most of the alternative Technische Hochschule (ETH) Zürich, Fermi National Accelera-
statistics mentioned in the introduction: one only needs to measure tor Laboratory, the University of Illinois at Urbana-Champaign,
the desired statistic on our training sets (Cosmology, Covariance and the Institut de Ciències de l’Espai (IEEC/CSIC), the Institut de
Systematics), re-train the GPR emulator and adjust the linear models Fı́sica d’Altes Energies, Lawrence Berkeley National Laboratory,
that describe the photometric redshifts and the shear calibration. The the Ludwig-Maximilians Universität München and the associated
rest of the infrastructure is already in place, and we intend to explore Excellence Cluster Universe, the University of Michigan, the National
some of these alternative probes in the near future. For this reason we Optical Astronomy Observatory, the University of Nottingham, The
will make the mocks available upon request. We also intend to explore Ohio State University, the OzDES Membership Consortium, the Uni-
some extensions to the existing infrastructure, including for instance versity of Pennsylvania, the University of Portsmouth, SLAC National
a new Neutrino training set derived from the MassiveNuS simulations Accelerator Laboratory, Stanford University, the University of Sussex,
(Liu et al. 2018) and designed to measure the sum of the neutrino mass. and Texas A&M University. Based in part on observations at Cerro
Parallel pipelines tailored for the analysis of the KiDS-1000 and/or the Tololo Inter-American Observatory, National Optical Astronomy
HSC surveys could also be constructed directly, and we hope to see Observatory, which is operated by the Association of Universities for
these novel methods take an important role in the upcoming Stage-IV Research in Astronomy (AURA) under a cooperative agreement with
lensing analyses. the National Science Foundation.
All authors contributed to the development of this paper. JHD (lead)

conducted the analysis with significant scientific contribution from NM
ACKNOWLEDGEMENTS (co-lead); the other authors are listed in alphabetical order.
We would like to thank Marika Asgari for useful discussions and

comments, Carlo Giocoli for his contribution to an earlier ver-
sion of the Magneticum light-cone extractor, as well as Shahab
DATA AVAILABILITY
Joudaki making public the MCMC chains described in J20, which
can be obtained from https://github.com/sjoudaki/kidsdes/. JHD The SLICS numerical simulations can be found at
acknowledges support from an STFC Ernest Rutherford Fellowship http://slics.roe.ac.uk/, while the SLICS-HR, the cosmo-SLICS
(project reference ST/S004858/1). NM acknowledges support from and the Magneticum can be made available upon request. This
© 2020 RAS, MNRAS 000, 1–26

work also uses public DES-Y1 data, which can be found at Harnois-Déraps, J., & van Waerbeke, L. 2015, MNRAS, 450, 2857,
https://des.ncsa.illinois.edu/releases/dr1. 1406.0543
Hartlap, J., Simon, P., & Schneider, P. 2007, A&A, 464, 399,
arXiv:astro-ph/0608064
Heitmann, K. et al. 2016, ApJ, 820, 108
REFERENCES Heitmann, K., Lawrence, E., Kwan, J., Habib, S., & Higdon, D. 2014,
Abbott, T. M. C. et al. 2018a, Phys. Rev. D, 98, 043526, 1708.01530 ApJ, 780, 111, 1304.7849
——. 2018b, ApJS, 239, 18, 1801.03181 Heydenreich, S., Brück, B., & Harnois-Déraps, J. 2020a, arXiv e-
Alarcon, A. et al. 2020, arXiv e-prints, arXiv:2007.11132, prints, arXiv:2007.13724, 2007.13724
2007.11132 Heydenreich, S. et al. 2020b, A&A, 634, A104, 1910.11327
Alonso, D., Sanchez, J., Slosar, A., & LSST Dark Energy Science Heymans, C. et al. 2020, arXiv e-prints, arXiv:2007.15632,
Collaboration. 2019, MNRAS, 484, 4127, 1809.09603 2007.15632
Asgari, M. et al. 2020a, arXiv e-prints, arXiv:2007.15633, Heymans, C., White, M., Heavens, A., Vale, C., & van Waerbeke, L.
2007.15633 2006, MNRAS, 371, 750, arXiv:astro-ph/0604001
——. 2020b, A&A, 634, A127, 1910.05336 Heymans, C., et al. 2012, MNRAS, 427, 146
Becker, M. R. et al. 2016, Phys. Rev. D, 94, 022002 Hikage, C. et al. 2019, PASJ, 71, 43, 1809.09148
Benı́tez, N. 2000, ApJ, 536, 571 Hilbert, S. et al. 2020, MNRAS, 493, 305, 1910.10625
Bridle, S., & King, L. 2007, New Journal of Physics, 9, 444, Hildebrandt, H. et al. 2018, arXiv e-prints, arXiv:1812.06076,
0705.0166 1812.06076
Castro, T., Borgani, S., Dolag, K., Marra, V., Quartin, M., Saro, A., & ——. 2017, MNRAS, 465, 1454
Sefusatti, E. 2020, arXiv e-prints, arXiv:2009.01775, 2009.01775 Hirschmann, M., Dolag, K., Saro, A., Bachmann, L., Borgani, S., &
Castro, T., Quartin, M., Giocoli, C., Borgani, S., & Dolag, K. 2018, Burkert, A. 2014, MNRAS, 442, 2304, 1308.0333
MNRAS, 478, 1305, 1711.10017 Hoyle, B. et al. 2018, MNRAS, 478, 592, 1708.01532
Catelan, P., Kamionkowski, M., & Blandford, R. D. 2001, MNRAS, Huang, H.-J. et al. 2020, arXiv e-prints, arXiv:2007.15026,
320, L7, astro-ph/0005470 2007.15026
Chang, C. et al. 2013, MNRAS, 434, 2121 Izard, A., Fosalba, P., & Crocce, M. 2018, MNRAS, 473, 3051
——. 2019, MNRAS, 482, 3696, 1808.07335 Jarvis, M., Bernstein, G., & Jain, B. 2004, MNRAS, 352, 338
Cheng, S., Ting, Y.-S., Ménard, B., & Bruna, J. 2020, MNRAS, 499, Joachimi, B. et al. 2020, arXiv e-prints, arXiv:2007.01844,
5902, 2006.08561 2007.01844
Chisari, N. et al. 2015, MNRAS, 454, 2736, 1507.07843 Joachimi, B., Semboloni, E., Bett, P. E., Hartlap, J., Hilbert, S., Hoek-
Chisari, N. E. et al. 2017, MNRAS, 472, 1163, 1702.03913 stra, H., Schneider, P., & Schrabback, T. 2013a, MNRAS, 431, 477,
——. 2018, MNRAS, 480, 3962, 1801.08559 1203.6833
Codis, S. et al. 2015, MNRAS, 448, 3391, 1406.4668 Joachimi, B., Semboloni, E., Hilbert, S., Bett, P. E., Hartlap, J., Hoek-
Coulton, W. R., Liu, J., McCarthy, I. G., & Osato, K. 2020, MNRAS, stra, H., & Schneider, P. 2013b, MNRAS, 436, 819, 1305.5791
495, 2531, 1910.04171 Johnson, A. et al. 2017, MNRAS, 465, 4118
Davies, C. T., Cautun, M., Giblin, B., Li, B., Harnois-Déraps, J., & Johnston, H. et al. 2019, A&A, 624, A30, 1811.09598
Cai, Y.-C. 2020, arXiv e-prints, arXiv:2010.11954, 2010.11954 Joudaki, S. et al. 2020, A&A, 638, L1, 1906.09262
Davies, C. T., Cautun, M., & Li, B. 2019, MNRAS, 490, 4907, ——. 2017, MNRAS, 471, 1259, 1610.04606
1907.06657 Kacprzak, T. et al. 2016, MNRAS, 463, 3653
DES Collaboration et al. 2017, ArXiv e-prints, 1708.01530 Kaiser, N., & Squires, G. 1993, ApJ, 404, 441
Dietrich, J. P., & Hartlap, J. 2010, MNRAS, 402, 1049 Kiessling, A. et al. 2015, Space Sci. Rev., 193, 67, 1504.05546
Erben, T. et al. 2013, MNRAS, 433, 2545, 1210.8156 Kilbinger, M. 2015, Reports on Progress in Physics, 78, 086901
Euclid Collaboration: Knabenhans, M. et al. 2019, MNRAS, 484, Kilbinger, M. et al. 2017, MNRAS, 472, 2126, 1702.05301
5509, 1809.04695 Kilbinger, M., et al. 2013, MNRAS, 430, 2200
Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601, Kirk, D. et al. 2015, Space Sci. Rev., 193, 139, 1504.05465
0809.3437 Köhlinger, F. et al. 2017, MNRAS, 471, 4412, 1706.02892
Flaugher, B. et al. 2015, AJ, 150, 150, 1504.02900 Krause, E. et al. 2017, arXiv e-prints, arXiv:1706.09359, 1706.09359
Fluri, J., Kacprzak, T., Lucchi, A., Refregier, A., Amara, A., Hof- Laigle, C. et al. 2016, ApJS, 224, 24, 1604.02350
mann, T., & Schneider, A. 2019, Phys. Rev. D, 100, 063514, Leistedt, B., McEwen, J. D., Büttner, M., & Peiris, H. V. 2017, MN-
1906.03156 RAS, 466, 3728, 1605.01414
Fortuna, M. C., Hoekstra, H., Joachimi, B., Johnston, H., Chis- Li, Y., Hu, W., & Takada, M. 2014, Phys. Rev. D, 89, 083519
ari, N. E., Georgiou, C., & Mahony, C. 2020, arXiv e-prints, Li, Z., Liu, J., Zorrilla Matilla, J. M., & Coulton, W. R. 2019,
arXiv:2003.02700, 2003.02700 Phys. Rev. D, 99, 063527, 1810.01781
Fu, L. et al. 2014, MNRAS, 441, 2725, 1404.5469 Lima, M., Cunha, C. E., Oyaizu, H., Frieman, J., Lin, H., & Sheldon,
Gatti, M. et al. 2020, MNRAS, 498, 4060, 1911.05568 E. S. 2008, MNRAS, 390, 118, 0801.3822
Giblin, B. et al. 2018, MNRAS, 480, 5529, 1805.12084 Liu, J., Bird, S., Zorrilla Matilla, J. M., Hill, J. C., Haiman, Z., Mad-
Gruen, D., & Brimioulle, F. 2017, MNRAS, 468, 769, 1610.01160 havacheril, M. S., Petri, A., & Spergel, D. N. 2018, Journal of Cos-
Gruen, D. et al. 2018, Phys. Rev. D, 98, 023507, 1710.05045 mology and Astro-Particle Physics, 2018, 049, 1711.10524
Hamana, T. et al. 2020, PASJ, 72, 16, 1906.06041 Liu, J., Petri, A., Haiman, Z., Hui, L., Kratochvil, J. M., & May, M.
Harnois-Déraps, J. et al. 2018, MNRAS, 481, 1337 2015a, Phys. Rev. D, 91, 063507
Harnois-Déraps, J., Giblin, B., & Joachimi, B. 2019, A&A, 631, Liu, X. et al. 2015b, MNRAS, 450, 2888
A160, 1905.06454 Mandelbaum, R. 2018, ARA&A, 56, 393, 1710.03235
Harnois-Déraps, J., Pen, U.-L., Iliev, I. T., Merz, H., Emberson, J. D., Mandelbaum, R. et al. 2011, MNRAS, 410, 844
& Desjacques, V. 2013, MNRAS, 436, 540 Mandelbaum, R. et al. 2005, MNRAS, 361, 1287
© 2020 RAS, MNRAS 000, 1–26

Martinet, N., Harnois-Déraps, J., Jullo, E., & Schneider, P. 2020, 2PCF (Mocks)
T18 wCDM
arXiv e-prints, arXiv:2010.07376, 2010.07376 Peaks (Mocks)
Martinet, N. et al. 2018, MNRAS, 474, 712, 1709.07678
McCarthy, I. G., Schaye, J., Bird, S., & Le Brun, A. M. C. 2017,
MNRAS, 465, 2936, 1603.02702
Mead, A., Brieden, S., Tröster, T., & Heymans, C. 2020, arXiv e- 0.85
0.80
prints, arXiv:2009.01858, 2009.01858 0.75
S8
Mead, A. J., Peacock, J. A., Heymans, C., Joudaki, S., & Heavens, 0.70
A. F. 2015, MNRAS, 454, 1958 0.65
Morrison, C. B., Hildebrandt, H., Schmidt, S. J., Baldry, I. K., Bilicki, 0.80
M., Choi, A., Erben, T., & Schneider, P. 2017, MNRAS, 467, 3576, 0.75
0.70
h
1609.09085
Munshi, D., Namikawa, T., Kitching, T. D., McEwen, J. D., & 0.65
Bouchet, F. R. 2020, MNRAS, 2006.12832 0.6
Nishimichi, T. et al. 2019, ApJ, 884, 29, 1811.09504 1.0
Osato, K., Shirasaki, M., & Yoshida, N. 2015, ApJ, 806, 186,
w0
1.4
1501.02055
Petri, A., Liu, J., Haiman, Z., May, M., Hui, L., & Kratochvil, J. M. 1.8
0.2 0.3 0.4 0.5 0.7 0.8 0.65 0.70 0.75 1.5 1.0
2015, Phys. Rev. D, 91, 103511, 1503.06214 S8 h w0
m
Planck Collaboration et al. 2018, arXiv e-prints, arXiv:1807.06205,
1807.06205 Figure A1. Two-dimensional constraints on the wCDM cosmological param-
Price, M. A., McEwen, J. D., Pratley, L., & Kitching, T. D. 2020, eters (Ωm , w0 , S 8 , h) from 2PCF (blue) and peak count analyses (grey) of the
arXiv e-prints, arXiv:2004.07855, 2004.07855 mock DES-Y1 data corresponding to the mean SLICS values. The red contours
Radovich, M. et al. 2017, A&A, 598, A107, 1701.02954 show the fiducial T18 results wCDM for reference. The dashed lines indicate
Ragagnin, A., Dolag, K., Biffi, V., Cadolle Bel, M., Hammer, N. J., the input SLICS cosmology. The inferred S 8 value appears lower than the in-
Krukau, A., Petkova, M., & Steinborn, D. 2017, Astronomy and put value, but this is a plotting artefact; maximum a posterior (MAP) value is
S 8MAP = 0.793 and 0.785 for 2PCF and peaks respectively, which are very close
Computing, 20, 52, 1612.06380
to the input of 0.813 (see main text for more details).
Robertson, N., Alonso, D., Harnois-Déraps, J., et al. 2020, in prep.
Samuroff, S. et al. 2019, MNRAS, 489, 5453, 1811.06989
Schirmer, M., Erben, T., Hetterscheidt, M., & Schneider, P. 2007,
A&A, 462, 875, astro-ph/0607022 2013, MNRAS, 434, 1604, 1302.0183
Schneider, A., Teyssier, R., Stadel, J., Chisari, N. E., Le Brun, Zuntz, J. et al. 2015, Astronomy and Computing, 12, 45
A. M. C., Amara, A., & Refregier, A. 2019, JCAP, 2019, 020, ——. 2018, MNRAS, 481, 1149, 1708.01533
1810.08629 Zürcher, D., Fluri, J., Sgier, R., Kacprzak, T., & Refregier, A. 2020,
Schneider, P. 1996, MNRAS, 283, 837, astro-ph/9601039 arXiv e-prints, arXiv:2006.12506, 2006.12506
Schneider, P., van Waerbeke, L., Kilbinger, M., & Mellier, Y. 2002,
A&A, 396, 1, arXiv:astro-ph/0206182
Sellentin, E., & Heavens, A. F. 2016, MNRAS, 456, L132
APPENDIX A: VALIDATION OF THE COSMOLOGY
Shan, H. et al. 2018, MNRAS, 474, 1116, 1709.07651
INFERENCE PIPELINE
Sheldon, E. S., & Huff, E. M. 2017, ApJ, 841, 24
Singh, S., Mandelbaum, R., & More, S. 2015, MNRAS, 450, 2195, In this section we present a series of validation tests we performed on
1411.1755 our cosmology inference pipeline.
Smith, R. E. et al. 2003, MNRAS, 341, 1311
Springel, V. 2005, MNRAS, 364, 1105, arXiv:astro-ph/0505010
Suchyta, E. et al. 2016, MNRAS, 457, 786, 1507.08336 A1 Inference from 2PCF
Takahashi, R., Hamana, T., Shirasaki, M., Namikawa, T., Nishimichi, The ξ± measurements on the simulations, presented in Fig. 2, show an
T., Osato, K., & Shiroyama, K. 2017, ApJ, 850, 24 excellent agreement with the input cosmology, well within the statisti-
Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. cal precision of the DES-Y1 data. There is a small loss of power in the
2012, ApJ, 761, 152 simulations compared to the data that can be observed at large angular
Teklu, A. F., Remus, R.-S., Dolag, K., Beck, A. M., Burkert, A., scales, a known finite-box effect that is caused by the absence of den-
Schmidt, A. S., Schulze, F., & Steinborn, L. K. 2015, ApJ, 812, sity fluctuation modes larger than the simulated volume. The ξ± statis-
29, 1503.03501 tics measured from the SLICS are best modelled with a prediction that
Troxel, M. A. et al. 2018, Phys. Rev. D, 98, 043528, 1708.01538 include a minimum k-mode of kmin = 2π/505 hMpc−1 in the matter
van Uitert, E. et al. 2018, MNRAS, 476, 4662 power spectrum, which is readily implemented in cosmoSIS with the
van Waerbeke, L., et al. 2013, MNRAS, 433, 3373 kmin option. A correction scheme for the sample variance also exists
Weiss, A. J., Schneider, A., Sgier, R., Kacprzak, T., Amara, A., & (HD15), however the contribution to the total error budget from these
Refregier, A. 2019, JCAP, 2019, 011, 1905.11636 missing modes is negligible. We present in Fig. A1 the wCDM con-
Wright, A. H., Hildebrandt, H., van den Busch, J. L., Heymans, C., straints from our 2PCF inference pipeline, when analysing the mean
Joachimi, B., Kannawadi, A., & Kuijken, K. 2020, A&A, 640, L14, of the Covariance training set. For the exercise presented in this sec-
2005.04207 tion, the priors on nuisance parameters have all been centred on zero,
Xia, Q. et al. 2020, A&A, 633, A89, 1909.05852 since no systematic shifts are infused in the Covariance training set.
Zorrilla Matilla, J. M., Waterval, S., & Haiman, Z. 2020, AJ, 159, Our inference method produces contours that are highly similar
284, 1909.12345 to those obtained by the wCDM analysis of T18 (blue vs. red), with a
Zuntz, J., Kacprzak, T., Voigt, L., Hirsch, M., Rowe, B., & Bridle, S. best fit value that is consistent with the input cosmology, albeit with a
© 2020 RAS, MNRAS 000, 1–26

Peaks wCDM
T18 wCDM
2PCF wCDM
0.9
0.8
0.7
h
0.6
0.5
1.0
0.9
0.8
8
0.7
0.6
0.9
0.8
S8
0.7
0.6
0.5
1.0
w
1.5
2.0
0.05
m1
0.00
0.05
0.05
m2
0.00
0.05
0.05
m3
0.00
0.05
0.05
m4
0.00
0.05
2
0
AI
2
4
5
5
0.04
0.02
z1
0.00
0.02
0.04
0.04
0.02
0.00
z2
0.02
0.04
0.02
z3
0.00
0.02
0.03
0.00
z4
0.03
0.06
0.2 0.4 0.6 0.8 0.6 0.8 1.0 0.6 0.7 0.8 0.9 2.0 1.5 1.0 0.5 0.05 0.00 0.05 0.05 0.00 0.05 0.05 0.00 0.05 0.05 0.00 0.05 4 2 0 2 5 0 5 0.03 0.00 0.03 0.04 0.00 0.02 0.02 0.05 0.00
m h 8 S8 w m1 m2 m3 m4 AI z1 z2 z3 z4
Figure A2. Pipeline validation: Two-dimensional constraints on the wCDM cosmological parameters (Ωm , h, σ8 , S 8 , w0 ) and on the 10 nuisance parameters (4 × ∆zi ,
4 × ∆mi , AIA and η). Mean and covariance are obtained from the SLICS 2PCF (blue) and peaks (grey), and compared to the T18 constraints (red). The thin dashed
lines indicate the input values in the simulations, which are well recovered by the blue and grey contours. Note that the T18 priors differ, most notably in h, ∆mi and
∆zi .
slightly lower S 8 value. This shift is caused by our parameter sampling, nuisance parameters used in this analysis, which are associated with
more precisely by the upper limit on w0 . A closer look on the [w0 − S 8 ] photometric uncertainty (4 × ∆zi ), shear calibration (4 × ∆mi ) and in-
panel reveals that these two quantities are degenerate, and that the w0 trinsic alignments of galaxies (AIA and η). The blue contours are ob-
posterior is limited by the prior, especially on the upper bound. Had tained from our 2PCF inference pipeline, executed on the mean 2PCF
we access to higher values, the contours would extend further in the extracted from the SLICS and using the SLICS covariance matrix. We
upper-right direction, bringing centroid to higher S 8 values. This is observe that the constraints on the four cosmological parameters are
further supported by the fact that the maximum a posterior (MAP) well centred on the input cosmology, and the sizes of the contours are
value is located at S 8MAP = 0.793, which is very close to the input of fully consistent with those of T18.
0.813.
We also note that our pipeline accurately recovers the input (null)
We next present in Fig. A2 the two-dimensional marginal con- amplitude of the intrinsic alignment in the SLICS, which also shows a
straints on the wCDM parameters (Ωm , h, σ8 , S 8 , w0 ) and on the 10 degeneracy with S 8 , whereas the data prefers values closer to AI = 1,
© 2020 RAS, MNRAS 000, 1–26

0 2 40 2 40 2 40 2 4 1 10
1000
0 5
-1000 0
11 12 13 14 0.9
1000 0 0.2 0.4
0
-1000
22 23 24
0
0.8
-1000
33 34 -2000
500
0 0.7
-500
44 -1000
-1500
S8 = 0.61 0.68 0.76 0.83 0.90
0.6 2PCF
Peaks
Figure A3. Same as Fig. 4, but the black points represent the mean over the full
SLICS sample instead of the data. The error bars are showing the error on the 0.1 0.2 0.3 0.4 0.5 0.6
mean.
Figure A4. Scatter in the best-fit parameters inferred from 50 simulations. The
as reported in T18 and seen in the red contours. The fact that the poste- thin black lines link the solutions from peaks statistics (blue) and 2PCF (orange)
rior of the η parameter is almost identical in the data and in the IA-free analyses associated with the same survey realisation. The distribution of the
simulation indicates that the evidence for a redshift evolution in the length of these lines is shown in the inset.
data is not strong. In fact, it instead suggests that the measurement of
this parameter is sensitive to something else, present in both simula-
tions and data (see the discussion on degeneracies between IA param-
eters and photometric redshift nuisance parameters in Appendix C of We report on Fig. A4 the best-fit parameters for the 50 peaks
Heymans et al. 2020). All other nuisance parameters are prior domi- analyses (in blue, including the cross-tomographic terms) and 2PCF
nated, and we observe a clear difference between our nuisance priors analyses (in orange), and further link the pairs associated with the
and those of T18, especially for the photometric bias; the DIR method same simulations. The inset shows the distribution of the difference
has a smaller error, which reflects here in smaller contours for all ∆zi . in Ωm , which extends beyond 0.3. The value measured from the data is
∆Ωm = 0.21, which is indicated by the vertical dotted line. In this met-
ric, the difference in matter density inferred by both probes is common
A2 Inference from peaks and likely, given the noise levels present in the DES-Y1 data. Repeat-
ing the same exercise for S 8 reveals a smaller scatter, with ∆S 8 6 0.2;
We next run our peak statistics wCDM inference pipeline on the same
the difference of 0.2 observed in the data is therefore a rare event. Ex-
data vector as in Appendix A1, e.g. mean of the Covariance train-
cluding the cross-redshift bins significantly increases the width and the
ing set measurements. The mean peak function is presented by the
contours, which making the observed ∆S 8 = 0.3 more likely.
black symbols in Fig. A3, showing the noise-free signal which lies
well within the range covered by the Cosmology training set. Again,
in this validation exercise, the priors on the nuisance parameters are
centred on zero, and we marginalise over the 8 nuisance parameters
(4 × ∆zi , 4 × ∆mi , the two IA parameters are not included in the peaks APPENDIX B: INFUSION OF INTRINSIC ALIGNMENT
analysis). We next feed this data vector into our cosmoSIS pipeline, and In this paper we estimate the impact of IA on peak count statistics by
confirm in Fig. A1 that we recover the same cosmology as the 2PCF infusing an alignment between mock galaxies and properties of their
analysis. Once again, the S 8 is slightly lower than expected, which can host dark matter haloes, based on the method described in Heymans
be explained by the asymmetric sampling of w0 that causes a bias in the et al. (2006) and Joachimi et al. (2013b). Since most of the simu-
inferred contour plots. The maximum a posterior value is well within lated light-cones used in this paper do not use dark matter haloes to
1σ, with S 8MAP = 0.785. We also recover similar contours on Ωm and assign galaxy positions, we instead have to rely on a separate simula-
S 8 from peaks and 2PCF. Results are summarised in Table 5. tion suite. We use for this task the KiDS-HOD mock data described
in Harnois-Déraps et al. (2018), in which dark matter haloes from the
light-cones are populated with an HOD based on a conditional lumi-
A3 Towards a joint analysis
nosity function, then sub-sampled to match the galaxy redshift distri-
Peak statistics and 2PCF are affected differently by the noise in the bution of the KiDS-450 data up to z = 1.5. Since these mocks are
data, and even though they converge in the noise-free limit (as shown slightly denser than the DES-Y1 data, we further downsample them
in Fig. A1), it is expected that a given noise realisation will scatter the such as to closely match the DES tomographic bins. These DES-HOD
best-fit cosmological parameters inferred by the two pipelines. We esti- mocks lack some of the very low redshift galaxies at z < 0.2 and all
mate the level of scatter by comparing the best-fit parameters obtained galaxies with z > 1.5, but these are all down-weighted by the lens-
at the maximum likelihood returned by 50 survey realisations taken ing kernel, leaving the expected lensing signal almost unchanged. We
from the Covariance training set. Following the recommendations of match the mean redshift in each bins to within 0.04, and the galaxy
Joachimi et al. (2020), we improve the accuracy of the solution by re- density to within 0.07 gal/armin2 , which is sufficient for this exercise.
peating the process multiple times by varying the starting points in the A number of physical models are presented in Joachimi et al.
[σ8 − Ωm ] plane and recording only the solution with the lowest χ2 . In (2013b), and in this first study we opted for the simplest possible case:
this exercise all nuisance parameters are held to zero, and we include we assume that all central galaxies are early-type red galaxies, and that
the cross-redshift bins in the peak count data to exacerbate the effect. all satellites are late-type blue galaxy. This is of course inaccurate, but
© 2020 RAS, MNRAS 000, 1–26

provides a good staring point upon which we can build and improve
-4
the model in future work. 10 11 12
The ellipticity of the central galaxies is known to correlate with

the shape of the host dark matter halos, in a complicated way that
depends on galaxy type, redshift, and possibly merger history (for a 10 -6
review see Kiessling et al. 2015). The galaxy catalogues that we con-
struct from the DES-HOD mocks contain the inertia matrix of the host
dark matter haloes, from which one can compute the eigenvalues and 10 -8
10
0
10
1 10 -4 22
eigenvectors (ωµ , sµ ). The projected ellipses reconstructed from these
are described by the symmetric tensor W (Joachimi et al. 2013a):
Theory Mocks
3 10 -6
X s⊥,µ s⊥,µT AIA=0.0 noIA
W −1 = , (B1)
ωµ2
AIA=1.0-2.0 IA
µ=1 II II
-GI -GI -8
10
where s⊥,µ is the eigenvector projected along the line of sight, and 10
0
10
1
the semi-major axes are given by 1/ωµ . Halo ellipticities h can be

p
Figure B1. Effect of the intrinsic alignment of galaxy on ξ+ . The upper and
obtained from:
lower thin black lines show the theoretical predictions based on the linear non-
W1,1 − W2,2 linear model of Bridle & King (2007) with AIA = 1.0 and 2.0, respectively,
h,1 = √ (B2)
W1,1 + W2,2 + 2 det(W) while the thick black line presents the predictions without IA. The solid red
2W1,2 and the dashed blue lines show the ‘II’ and ‘−GI’ contributions. The symbols
h,2 = √ . (B3) represent the measurements on the dedicated IA mocks. The thin red symbols
W1,1 + W2,2 + 2 det(W) present the absolute value of the ‘II’ measurements, which becomes negative
Once these are determined, we opted for a 100% alignment between at small scales. Three panels show the combinations between the two lowest
the halo and the central galaxy ellipticities. This is likely to slightly redshift bins, as indicated in the top right corner. Combinations involving higher
redshifts show a similar agreement between mocks and models, albeit the effect
over-estimate the effect, a deliberate choice that we make when devel-
is milder.
oping a relatively conservative approach. The absence of scatter be-
tween the two, combined with the approximation that all centrals are
early-type galaxies both act as to maximise the IA signal in our model.
Hildebrandt et al. (2018). These differences are responsible for a loss
We also find that this model does not work well at high redshift, since
of precision on the S 8 constraints: while T18 finds S 8T18,Λ = 0.792+0.032
−0.021 ,
the haloes are not fully relaxed, and their shapes are less well mod-
J20 reports S 8J20 = 0.763+0.037
−0.031 , e.g. a ∼ 1σ shift compared to the orig-
elled. We therefore concentrate on the lower redshift bins only.
inal DES-Y1 results, and an error almost 30% larger. As explained in
We make a second important approximation by treating all satel-
J20, the shift in S 8 is largely driven by differences in the n(z) estima-
lite galaxies as blue, late-type, and assign them no intrinsic alignment,
tions.
consistent with recent findings (Mandelbaum et al. 2011; Samuroff
et al. 2019; Johnston et al. 2019). We could assign the halo elliptic-
ities to the centrals directly, but doing so strongly biases the ellipticity This paper has been typeset from a TEX/ LATEX file prepared by the
distribution to lower values compared with the data. Instead we keep author.
the ellipticities drawn from the Gaussian distribution, and rotate them
until they align with h . Once this is done, we use Eq. (1) to shear
these intrinsic ellipticities, and compute ξ+11 , ξ+12 and ξ+22 with given
either from the IA model described above or from the no-IA case. We
additionally compute the same 2PCF statistics from the pure intrinsic
shapes int , thereby estimating the II term, as well as the combina-
tion h int γi to compute the GI term. Results are presented in Fig. B1,
and compared with the theoretical predictions with AIA = 1.0 and 2.0.
We see in all three tomographic bin combinations that the mocks re-
produce reasonably well the IA, noIA and GI models, however the II
terms remains very noisy.
The IA model infused in these mocks is not accurate enough to
be used for signal calibration, but is adequate for diagnostic tests such
as those for which they are designed here. Future developments with
more flexible options regarding early-types/late-types separation, in-
clusion of satellite alignment, and possibly different N-body runs, will
be the subject of future work.
APPENDIX C: COMPARISON BETWEEN ΛCDM ANALYSES

Differences in signal modelling and in the likelihood sampling strategy
are responsible for the broader range of accepted Ωm and σ8 values in
J20, compared to T18. Notably, they replaced the halofit predictions
of the non-linear matter spectrum by the halo-based model HMcode
(Mead et al. 2015). Furthermore, the sampling overthe amplitude pa-
rameter As is replaced by a sampling over ln 1010 As , the sum over the
neutrino mass is fixed, and the ranges of priors are changed to that of
© 2020 RAS, MNRAS 000, 1–26

Cosmic Shear Cosmology Beyond 2-Point Statistics: A Combined Peak Count and Correlation Function Analysis of DES-Y1

Uploaded by

Document Informationclick to expand document information

Document Informationclick to expand document information

Copyright:

Available Formats

Cosmic Shear Cosmology Beyond 2-Point Statistics: A Combined Peak Count and Correlation Function Analysis of DES-Y1

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Cosmic Shear Cosmology Beyond 2-Point Statistics: A Combined Peak Count and Correlation Function Analysis of DES-Y1

Uploaded by

Copyright:

Available Formats

Mon. Not. R. Astron. Soc. 000, 1–26 (2020) Printed 7 December 2020 (MN LATEX style file v2.

Cosmic Shear Cosmology Beyond 2-Point Statistics: A Combined Peak

Joachim Harnois-Déraps1,2? , Nicolas Martinet3 , Tiago Castro4,5,6,7 , Klaus Dolag8,9 ,

Accepted XXX. Received YYY; in original form ZZZ

© 2020 RAS, MNRAS 000, 1–26

IM3SHAPE (Zuntz et al. 2013) that were both fully implemented

tomo ZB range No. of objects neff σ hzDIR i

© 2020 RAS, MNRAS 000, 1–26

2.4 Systematics training set: mass resolution

© 2020 RAS, MNRAS 000, 1–26

2.6 Simulation post-processing

© 2020 RAS, MNRAS 000, 1–26

14 bitbucket.org/joezuntz/cosmosis/wiki/Home 15 In practice, we use a link-list to loop only over nearby galaxies.

© 2020 RAS, MNRAS 000, 1–26

© 2020 RAS, MNRAS 000, 1–26

vestigations). 12 bins is also a good trade-off between our need to cap-

© 2020 RAS, MNRAS 000, 1–26

© 2020 RAS, MNRAS 000, 1–26

Figure 6. Accuracy of the GPR emulator, computed with a leave-one-out cross-

© 2020 RAS, MNRAS 000, 1–26

© 2020 RAS, MNRAS 000, 1–26

© 2020 RAS, MNRAS 000, 1–26

© 2020 RAS, MNRAS 000, 1–26

A large fraction of the weak lensing signal receives a contribution from

© 2020 RAS, MNRAS 000, 1–26

observed signal. Furthermore, these regions of high galaxy densities

© 2020 RAS, MNRAS 000, 1–26

© 2020 RAS, MNRAS 000, 1–26

This work σ8 ∈ [0.53 − 1.3] SLICS DIR

© 2020 RAS, MNRAS 000, 1–26

parameters, and that we marginalise over these anyway. Consequently, 1.5

© 2020 RAS, MNRAS 000, 1–26

© 2020 RAS, MNRAS 000, 1–26

© 2020 RAS, MNRAS 000, 1–26

All authors contributed to the development of this paper. JHD (lead)

We would like to thank Marika Asgari for useful discussions and

© 2020 RAS, MNRAS 000, 1–26

© 2020 RAS, MNRAS 000, 1–26

© 2020 RAS, MNRAS 000, 1–26

© 2020 RAS, MNRAS 000, 1–26

© 2020 RAS, MNRAS 000, 1–26

The ellipticity of the central galaxies is known to correlate with

the semi-major axes are given by 1/ωµ . Halo ellipticities  h can be

APPENDIX C: COMPARISON BETWEEN ΛCDM ANALYSES

© 2020 RAS, MNRAS 000, 1–26

You might also like

tomo ZB range No. of objects neff σ hzDIR i

the semi-major axes are given by 1/ωµ . Halo ellipticities h can be