SBi3PCF: cosmology with Dark Energy Survey Year 3 cosmic shear

Prepared for submission to JCAP

Introduction

Cosmic shear, the coherent distortion imprinted in the observed shapes of background galaxies by the weak gravitational lensing (WL) effect of the foreground matter distribution [1, 2, 3], has become one of the most powerful probes of the large-scale structure (LSS) of the Universe. By analysing statistical correlations in the cosmic shear field, we can probe the distribution and the properties of the matter density field and constrain key cosmological parameters, such as the matter density $\Omega_{\mathrm{m}}$, the amplitude of matter density fluctuations $\sigma_{8}$ (often expressed through the structure growth parameter $S_{8} = \sigma_{8}\sqrt{\Omega_{\mathrm{m}}/0.3}$, and the dark energy equation of state parameter $w_{0}$.

Current Stage-III WL surveys, such as the Dark Energy Survey (DES) [4, 5, 6], the Kilo-Degree Survey (KiDS) [7, 8], and the Hyper Suprime-Cam (HSC) survey [9], have placed powerful constraints on these parameters, primarily by measuring two-point statistics (2-point correlation functions $\xi_{\pm}$ (2PCFs) and power spectra) of the shear field [10, 11, 12, 13, 14, 15, 16]. However, $\xi_{\pm}$ only probes the information contained in the second moments of the density field and is not sensitive to the higher-order moments of the non-Gaussian late-time matter distribution [17]. To extract this additional information, one has to turn to higher-order statistics (HOS), which have already been applied to weak lensing survey data, for example through third-order shear statistics measured on CFHTLenS [18], DES [19, 20, 21], KiDS [22, 23, 24] and HSC [25] data, moments and cumulative distribution functions of the DES Y3 mass maps [26, 27], peak counts on DES [28, 29], KiDS [30, 31] and HSC [32] data, the convergence PDF on HSC [33, 34] data, persistent homology on DES [35] and the scattering transform on HSC data [36]. Beyond these, a range of further HOS is currently being developed, with the Euclid papers [37, 38] giving an overview, including the weak lensing convergence 1-point PDF [39, 40, 41, 42, 43], critical points of the convergence field [44], the parity-odd lensing trispectrum [45] and machine-learning-derived summary statistics [46, 47, 48].

A practical way of accessing higher-order information is through the integrated 3-point correlation function (i3PCF) [49, 50, 51], which correlates $\xi_{\pm}$ measured locally in small patches of the sky with the aperture mass at the patch centres. It is mostly sensitive to squeezed triangle configurations of the lensing bispectrum and is measured directly from the cosmic shear field, without the need to reconstruct a convergence map, or to apply a full 3PCF estimator. In contrast to many other HOS, it can also be modelled analytically and interpreted physically using the response approach to perturbation theory [52, 53].

Extracting the cosmological information contained in such HOS such as the i3PCF, however, requires a likelihood function, which in general is not known analytically. In simulation-based inference (SBI), the likelihood of the summary statistic is instead learned from a suite of forward-modelled simulations that include the survey and noise properties, and systematic effects using generative neural networks such as normalising flows [54, 55]. This avoids approximations such as the widely used Gaussian likelihood, which might not hold for HOS. SBI has been applied to cosmic shear data from DES SV [56] and KiDS-1000 [57, 58], and in a number of DES Y3 analyses using wavelet harmonics, scattering transforms and moments of mass maps [59, 60], persistent homology [61], neural compression at the map level [62, 63], and combinations of the power spectrum with neural network summaries [64].

In the first paper of this series [65] (hereafter Paper 1), we developed and validated SBi3PCF, a complete SBI framework for a joint analysis of the cosmic shear 2PCFs $\xi_{\pm}$ and the i3PCFs $\zeta_{\pm}$. Its forward model, built on the CosmoGridV1 N-body simulations [66], generates realistic DES Y3-like cosmic shear mocks including the survey mask, shape noise, and a range of crucial systematic effects such as intrinsic alignments and baryonic feedback. There, we have shown on simulated measurements that the pipeline produces well-calibrated posteriors, recovers the true cosmological parameters without bias, and that adding $\zeta_{\pm}$ to a standard $\xi_{\pm}$ analysis significantly improves the cosmological constraining power.

In this work, we apply this framework to the DES Y3 data and present the first cosmological analysis of $\zeta_{\pm}$ in weak lensing survey data. We measure $\xi_{\pm}$ and $\zeta_{\pm}$ from the DES Y3 weak lensing shape catalogue [6] and infer posteriors in the $w$CDM parameter space spanned by CosmoGridV1, for each statistic individually as well as for their joint analysis. To avoid confirmation bias, we carry out the analysis blind. This means that all analysis choices are frozen, and that a series of tests (coverage tests, tests of the accuracy of the learned likelihoods and internal consistency tests on the blinded measurement) has to be passed before we look at the cosmological constraints [67]. We then present our constraints on $\Omega_m$, $\sigma_8$, $S_8$ and $w_0$ and compare them with the other DES Y3 analyses of the same source catalogue, the latest DES Y6 $3 \times2$pt results [68], and also with external data: KiDS-1000 [69], KiDS-Legacy [13, 70], Planck [71] and DESI [72]. Aside from the constraints themselves, this analysis prepares $\zeta_{\pm}$ for the beyond-two-point SBI analyses of the upcoming Stage-IV surveys Euclid [73] and Rubin LSST [74], where the statistical power of the data will require the treatment of systematic effects such as intrinsic alignments [75], baryonic feedback [76, 77] and photometric redshift outliers [78] in the forward model (for a broader overview of these systematic effects, see Section 2.2.1 of Ref. [79]).

We structure this paper as follows. In § 2, we describe the DES Y3 data, the measurement of the correlation functions and the simulation-based forward model. § 3 briefly summarises our SBI pipeline, for which we refer the reader to Paper 1 for details. We present the pre-unblinding tests in § 4, our cosmological constraints in § 5 and the comparison with external data in § 6, and conclude in § 7.

Data and simulations

The forward model and the measurements are identical to Paper 1, so that we only repeat the essentials here. For further details, we refer the reader to Sections 2 and 3 of Paper 1 [65].

DES Y3

The Dark Energy Survey (DES) spans $\sim5000\,\mathrm{deg}^2$ in the grizY filters, observed with the Dark Energy Camera at the 4m Blanco telescope [4]. This work uses data from the first three years of observations (Y3) with the original data products, i.e. the Gold catalogue [5], the weak lensing shape catalogue [6], the PSF modelling [80], the image simulations used for the shear calibration [81] and the SOMPZ redshift calibration [82]. The only DES Y3 data product that we do not use is the shear-ratio (SR) likelihood [83], as we have neither a lens sample nor an analytic likelihood into which it could enter.

The DES Y3 weak lensing shape catalogue [6] uses METACALIBRATION-based [84, 85] shear estimates for 100,204,026 galaxies over an effective area of 4143 $\mathrm{deg}^2$, with a weighted effective number density of 5.59 $\mathrm{arcmin}^{-2}$ in the definition of Ref. [86]. The galaxies are divided into the four DES Y3 tomographic source bins [82], with the per-bin number densities, mean redshifts and ellipticity dispersions of Table I of Ref. [12], as in Paper 1. The redshift distribution of each bin was calibrated via SOMPZ [82] and is shown in Figure 2 of Paper 1. The Gaussian priors on the per-bin shear calibration biases $m_i$ listed in Table 1 are taken from the image simulations of [81]. They include a 2–3% residual bias due to blending, which enters through the non-zero means of the priors and is marginalised over in the analysis.

We have pixelised the catalogue onto HEALPix [87] maps with NSIDE = 512, corresponding to a pixel scale of 6.9′ (pixel area of 47 $\mathrm{arcmin}^2$). For each tomographic bin, we subtract the mean shear and apply the METACALIBRATION response, after which the value of a pixel is the weighted combination of the ellipticities of all galaxies falling within it, with the pixel weight given by the sum of the individual galaxy weights. Pixels outside the DES Y3 footprint are set to have zero weights, such that the mask is applied to the data and the mocks in the same way. Since every mock of the forward model carries the same per-pixel weights (§ 2.3), the observed data enter the measurement in exactly the same pixelised form as the mocks.

We have blinded our analysis at the catalogue level, following the DES Y3 scheme [6]. All galaxy ellipticities were rescaled by an unknown multiplicative factor, and every subsequent step (the pixelisation, the measurement, the compression and the inference) was run on this blinded catalogue. We did not apply any additional parameter shifts to the posteriors [67]. Instead, all axis values and numerical summaries of the blinded posteriors were hidden, and these posteriors were never shown together with literature, fiducial or other chains. This differs from the other DES Y3 SBI analyses, which blinded at the posterior level, either by not looking at the posteriors of the data [60, 61] or by shifting them by an unknown offset [63, 64], and is closest to the DES Y3 $\xi_\pm+ \langle M_{\mathrm{ap}}^3\rangle$ analysis [21], which also measured its data vector on a catalogue with rescaled ellipticities. Because the blinding acts on the catalogue, the measured data vector enters every test of § 4 unchanged. All analysis choices were frozen before these tests, and the unblinded catalogue was processed only after all of them had been passed. All figures and numbers of § 2–§ 4 were first produced with the blinded measurement and replaced by the unblinded ones afterwards, and after unblinding no changes were made to the pipeline, the analysis choices, the scale cuts or the validation criteria.

Correlation Functions

We use the cosmic shear 2PCFs $\xi_\pm$ and the i3PCFs $\zeta_\pm$, which correlate the local $\xi_\pm$ measured in small, roughly degree-scale sky patches with the aperture mass $M_{\mathrm{ap}}$ at the patch centres [49, 50, 51]. $\xi_\pm$ are estimated from the tangential and cross components of the pixel shears relative to the line connecting a pixel pair, as the weighted average over all pixel pairs within a separation bin, with the pixel weights of § 2.1. $M_{\mathrm{ap}}$ is estimated as the weighted sum of the tangential shear of all pixels around the patch centre with the compensated filter of Refs. [88, 89, 90] and aperture radius $\theta_{\mathrm{ap}}$, and $\zeta_\pm$ as the mean product of the local $\xi_\pm$ and aperture masses over all patches, minus the product of their averages over the footprint. The estimators are given in Eqs. (3.4)–(3.8) of Paper 1 and are applied in the same form to the pixelised DES Y3 maps and to every mock.

The measurement setup is identical to Paper 1. We measure the global $\xi_\pm$ in 8 logarithmically spaced bins in the angular separation $\alpha$ between 15′ and 250′. The lower limit is about twice the pixel scale at NSIDE = 512 and therefore excludes the small scales that would be most affected by baryonic feedback, while the upper limit matches the DES Y3 $\xi_\pm$ setup [10, 11]. For $\zeta_\pm$ we use five aperture radii $\theta_{\mathrm{ap}} \in\{50', 70', 90', 110', 130'\}$. The local $\xi_\pm$ are measured in the same binning as the global ones, truncated at $\alpha_{\max} \leq2\theta_{\mathrm{ap}} - 5'$, which leaves 5, 6, 6, 7 and 8 angular bins for the five apertures. The patch centres are placed on a HEALPix grid at NSIDE = 64 for $\theta_{\mathrm{ap}} = 50'$ and 70′ and at NSIDE = 32 for the larger apertures, of which we keep those centres for which both the patch (radius $\theta_{\mathrm{ap}}$) and the region over which the filter is evaluated (radius $5\theta_{\mathrm{ap}}$) are at most 50% masked, which leaves between 5763 (50′) and 846 (130′) overlapping patches. The patch centres are identical for the DES Y3 measurement and for every mock footprint.

For both statistics we measure the auto- and cross-correlations of all tomographic bin combinations. For $\xi_\pm$ these are the 10 pairs $(i,j)$ with $i \leq j$, and for $\zeta_\pm$ the 20 combinations $(i,j,k)$ with $i \le j \le k$, where $M_{\rm ap}$ is measured in bin $i$ and the local $\xi_{\pm}$ between bins $j$ and $k$. The $\xi_{\pm}$ data vector thus has 160 elements (10 tomographic bin combinations $\times$ 2 components $\times$ 8 angular bins), the $\zeta_{\pm}$ data vector 1280 (20 bin combinations $\times$ 2 components $\times$ 32 angular bins summed over the apertures), and the joint $\xi_{\pm}+\zeta_{\pm}$ data vector a total length of 1440.

All measurements, on the DES Y3 data and on all mocks, are performed with the same pipeline as in Paper 1, i.e. the global $\xi_{\pm}$ with TreeCorr [89, 91] on the pixelised maps and $\zeta_{\pm}$ with our GPU code PixPoint (formerly CosmoFuse; Appendix § 9 of Paper 1). As the same estimators are applied to identically pixelised maps, the measured data vector is directly comparable to mock data vectors produced by our forward model.

The measured $\xi_{\pm}$ and $\zeta_{\pm}$ are shown in Appendix § 11 together with the posterior predictive means of the analyses of § 5. With the covariance estimated from the 16,000 fiducial footprints of the forward model (§ 2.3), the signal-to-noise ratio, S/N $=\sqrt{x^{\mathsf{T}}C^{-1}x-N_{\rm d}}$, which removes the noise contribution of the $N_{\rm d}$ data points [19], is 31.0 for $\xi_{\pm}$, 12.8 for $\zeta_{\pm}$ and 33.0 for the joint data vector.

Finally, we do not perform any additional PSF or B-mode null tests on $\zeta_{\pm}$, since these tests have already been carried out at the level of the DES Y3 shape catalogue [6, 11, 12, 80], which we use as released. $M_{\rm ap}$ and the local $\xi_{\pm}$ are linear and quadratic, respectively, in the validated pixelised shear field as the global $\xi_{\pm}$, so that PSF residuals and B-modes at the level allowed by the catalogue tests, well below the statistical uncertainty of the DES Y3 $\xi_{\pm}$, are also expected to be negligible for the shape-noise-dominated $M_{\rm ap}$ measurements. Hence, we do not measure the B-mode combinations, i.e. $M_{\rm ap}$ from $\gamma_{\times}$ and $\xi_{\times}$.

Forward Model

The forward model was developed and validated in Paper 1, which is described in detail in its § 3. It is based on the CosmoGridV1 particle shells, from which we create weak lensing mocks with the DES Y3 mask, per-pixel shape noise from the observed catalogue and systematic effects added at the map level. In the following, we summarise the essential parts of the forward model for the data analysis and refer the reader to Paper 1 for the details.

The forward model varies 17 parameters: six $w$CDM parameters, two baryonification [92, 93] parameters capturing baryonic feedback, one intrinsic alignment amplitude (in the NLA model [94, 95]) and four parameters each for the shear calibration bias and the mean photometric redshift shift. Table 1 shows the prior ranges used in this work, which are the same as in Paper 1, together with the fiducial parameter values. For the cosmological, baryonification and intrinsic alignment parameters we use uniform priors. Apart from these, there are further restrictions in the $\Omega_{\rm m}$–$\sigma_8$ and $\Omega_{\rm m}$–$w_0$ planes, which are inherited from the CosmoGridV1 grid (see § 2 of Ref. [66]). The restriction in the $\Omega_{\rm m}$–$\sigma_8$ plane was chosen to concentrate on the simulations along the known degeneracy direction, and the latter in the $\Omega_{\rm m}$–$w_0$ plane removes those phantom DE cosmologies ($w_0 < -1$) for which the relativistic treatment of the simulations breaks down before $z=0$ (see Appendix § 9 of Ref. [57]). On the shear calibration biases $m_i$ and the photometric redshift shifts $\Delta z_i$ we use the DES Y3 Gaussian priors following Refs. [10, 11]. As in Paper 1, $A_{\rm IA}$, $m_i$ and $\Delta z_i$ are drawn from their priors once per full-sky map and shared by its four footprints, with $A_{\rm IA}$ drawn from a Sobol sequence [96] and $m_i$ and $\Delta z_i$ from their Gaussian priors. The fiducial values of Table 1 are those of the CosmoGridV1 fiducial cosmology, at which the covariance and the parameter derivatives entering the compression (§ 3) are evaluated. These priors are used for every inference in § 5. For a full corner plot including a chain sampled from the prior alone, see Figure 11 in Appendix § 10.

Parameter	Fiducial	Prior
$\Omega_m$	0.26	$\in[\ 0.10,\ 0.50\ ]$
$\Omega_b$	0.0493	$\in[\ 0.03,\ 0.06\ ]$
$H_0$	67.3	$\in[\ 64.0,\ 82.0\ ]$
$n_s$	0.9649	$\in[\ 0.87,\ 1.07\ ]$
$\sigma_8$	0.84	$\in[\ 0.40,\ 1.40\ ]$
$w_0$	$-1.0$	$\in[\ -2.00,\ -0.33\ ]$
$\log M_c^0$	13.82	$\in[\ 12.0,\ 15.0\ ]$
$\nu$	0.0	$\in[\ -2.00,\ 2.00\ ]$
$A_{\mathrm{IA}}$	0.0	$\in[\ -5.00,\ 5.00\ ]$
$m_1$	$-0.0063$	$\mathcal{N}(-0.0063,\ 0.0091)$
$m_2$	$-0.0198$	$\mathcal{N}(-0.0198,\ 0.0078)$
$m_3$	$-0.0241$	$\mathcal{N}(-0.0241,\ 0.0076)$
$m_4$	$-0.0369$	$\mathcal{N}(-0.0369,\ 0.0076)$
$\Delta z_1$	0.0	$\mathcal{N}(0.0,\ 0.018)$
$\Delta z_2$	0.0	$\mathcal{N}(0.0,\ 0.015)$
$\Delta z_3$	0.0	$\mathcal{N}(0.0,\ 0.011)$
$\Delta z_4$	0.0	$\mathcal{N}(0.0,\ 0.017)$

The CosmoGridV1 [66] suite of N-body simulations was run with PKDGRAV3 [97] and is provided as full-sky HEALPIX particle shells at NSIDE $=512$, with a mass resolution that is sufficient for the statistics and scales ($\ge15'$) considered here (see the discussion in § 3.1 of Paper 1) and three degenerate neutrino species with $\sum m_\nu=0.06\,\mathrm{eV}$. Throughout this work, we use the set of shells in which baryonic feedback has been included at the map level via the shell baryonification method [66, 92, 93] with the $(M_c,\nu)$ of each grid cosmology. Here, $M_c$ is the characteristic halo mass below which feedback has pushed a large fraction of the gas out of the haloes, flattening their gas profiles [93], on which we place a log-uniform prior, and $\nu$ describes its redshift evolution as $M_c(z)=M_c^0(1+z)^\nu$ [57, 66]. The simulation set consists of a grid of 2500 flat $w$CDM cosmologies with 7 realisations each, made up of two Sobol-sampled [96] sets of 1250 cosmologies that span the full prior box of Table 1 and a narrower central region, and of a fiducial set of 200 realisations at the fiducial cosmology and at a $\pm$ step in each of the 6 cosmological and 2 baryonification parameters. The test cosmologies of § 3, which are excluded from training, are taken from the full-box set.

From the particle shells, we create the weak lensing maps with UFalcon2 [98, 99, 100] as described in § 3.2 of Paper 1: the shells are weighted by the lensing kernel of each tomographic bin in the Born approximation, the resulting convergence maps are converted to shear maps through a full-sky spin-2 harmonic transform, akin to Kaiser-Squires reconstruction [101], and the DES Y3 mask is applied. By rotating each full-sky map, four non-overlapping footprints are extracted, leading to 2500 cosmologies $\times$ 7 realisations $\times$ 4 footprints $=70{,}000$ independent mocks for training. From the fiducial maps, 200 realisations $\times$ 20 random orientations $\times$ 4 footprints = 16,000 fiducial footprints are obtained, from which the covariance entering the compression (§ 3) is estimated. To each footprint, we then add per-pixel shape noise, obtained by randomly rotating the ellipticities of the galaxies in the observed DES Y3 catalogue and summing them up within each pixel, so that every mock carries the same per-pixel weights as the data, as in the other DES Y3 map-based SBI analyses [59, 62, 63]. The systematic effects are added at the map level as described in Section 2.4 of Paper 1: intrinsic alignments with the nonlinear alignment (NLA) model [94, 95] with a free amplitude $A_{\mathrm{IA}}$, without redshift evolution ($\eta= 0$) and with $C_1 = 5 \times10^{-14}h^{-2}M_{\odot}^{-1}\mathrm{Mpc}^{3}$ [94, 102], i.e. in the same normalisation as the DES Y3 NLA amplitude $a_1$ at the pivot redshift, the photometric redshift uncertainty as a coherent shift $n_i(z) \to n_i(z + \Delta z_i)$ of each bin’s SOMPZ distribution when projecting the shells, and the shear calibration bias as a multiplicative factor $(1 + m_i)$ per bin on the shear maps. Reduced shear and source galaxy clustering are not part of the forward model, as they were shown to change the data vector by far less than the statistical uncertainty in § 5.2 of Paper 1.

Simulation-based inference

We infer the cosmological parameters by learning the likelihood of the compressed data vector from the forward-modelled mocks, i.e. by neural likelihood estimation (NLE) [55, 103], as developed and validated in Paper 1. The pipeline consists of three steps, which are described in detail in § 4 of Paper 1: (i) a linear MOPED [104, 105] compression of the 1440-element data vector to 17 coefficients (§ 3.1), (ii) a normalising flow trained on the parameters and compressed data vectors of the training set to estimate the likelihood $p(x \mid\theta)$ (§ 3.2), and (iii) MCMC sampling of the posterior $p(\theta\mid x_{\mathrm{obs}}) \propto p(x_{\mathrm{obs}} \mid\theta)\,\pi(\theta)$ with the flat prior $\pi(\theta)$ of Table 1, for which we use our own GPU implementation of the affine-invariant sampler [106] (§ 3.3). The SBI pipeline does not assume a Gaussian likelihood or an analytical model of the data vector. The only place where the fiducial covariance is used is in the construction of the compression vectors (§ 3.1).

Compared to Paper 1, the pipeline is now trained and tuned separately for every analysis variation, i.e. $\xi_{\pm}$, $\zeta_{\pm}$, their combination and the data-vector variations of § 4, and then applied to the measured (blinded) DES Y3 data vector. We further use neural spline flows (NSFs) instead of masked autoregressive flows (MAFs), a different hyperparameter search space (see Appendix § 12), the binomial instead of the bootstrap uncertainty band for the coverage test (see Appendix § 9.§ 1), and a faster MCMC sampler that runs many chains at once on a GPU. On a single NVIDIA A100 GPU, one trial of the hyperparameter search (§ 3.2) takes about 10 minutes, so that its 300 trials amount to 2–3 GPU days per analysis variation. With the batched sampler, the 1,250 test and 800 fiducial posteriors of an analysis variation are sampled in about 3 and 2 hours, respectively, and the 10 million-sample DES Y3 chain in about 25 minutes.

MOPED Compression

To reduce the dimensionality of our data vectors, we employ the MOPED compression [104, 105], a linear transformation $y_{\alpha} = b_{\alpha}^{T}x$ that maximises the Fisher information of the data vector $x$ at the fiducial cosmology. It reduces the 1440-element joint data vector (160 for $\xi_{\pm}$ alone, 1280 for $\zeta_{\pm}$ alone) to 17 coefficients, one per forward-model parameter: the 6 cosmological parameters, $M_{\mathrm{c}}$, $\nu$, $A_{\mathrm{IA}}$, the four $m_i$ and the four $\Delta z_i$. The compression vectors $b_\alpha\propto C^{-1}\partial\mu/\partial\theta_\alpha$, where $\mu(\theta)$ is the mean data vector at the parameters $\theta$, are Gram-Schmidt orthogonalised in the order of Table 1, so that the coefficients are uncorrelated at the fiducial cosmology (Eqs. 4.4–4.5 of Paper 1). The covariance $C$ is estimated from the 16,000 fiducial footprints (§ 2.3) and its inverse carries the Hartlap correction factor [107] (for 16,000 footprints and 1440 elements the correction factor is 0.91). While this constant factor is not relevant for the MOPED coefficients, we still include it for the later goodness-of-fit tests in § 4. The parameter derivatives $\partial\mu/\partial\theta_\alpha$ with respect to the cosmological and baryonification parameters come from the dedicated CosmoGridV1 derivative simulations around the fiducial cosmology, as central finite differences of the noiseless mean over the 200 realisations of the $\pm$ simulations, and those with respect to $A_{\mathrm{IA}}$, $m_i$ and $\Delta z_i$ from the fiducial maps reprocessed with the respective parameter shifted.¹ We compute a separate set of compression vectors for each analysis variation defined in § 4. For the variations that remove a tomographic bin, the derivative with respect to its shear calibration bias vanishes, so that the compressed data vector has 16 instead of 17 coefficients.

We note that the compression is performed at the fiducial parameter set. It is lossless only for a Gaussian likelihood with a parameter-independent covariance and a data vector that depends linearly on the parameters. Away from the fiducial parameters, some information is lost, but no bias is introduced, since the likelihood of the compressed data vector is learned directly (see Paper 1).

Neural Density Estimation

To learn the likelihood $p(x \mid\theta)$ of the MOPED coefficients, we use neural spline flows (NSFs) [108], autoregressive flows whose transformations are monotonic rational-quadratic splines. They replace the masked autoregressive flows of Paper 1 and are implemented in zuko [109], which we use through the sbi package [110]. Both the parameters and the compressed data vectors are standardised with the mean and standard deviation of the simulation set before training. The base distribution of the flow is a standard normal, and the loss function is the negative mean log-likelihood of the training batch. As in Paper 1, we only provide the 9 parameters of interest, $\Omega_{\mathrm{m}}$, $\Omega_{\mathrm{b}}$, $H_0$, $n_{\mathrm{s}}$, $\sigma_8$, $w_0$, $A_{\mathrm{IA}}$, $M_{\mathrm{c}}$ and $\nu$, as labels. The remaining 8 systematic parameters ($m_i$ and $\Delta z_i$) are thereby implicitly marginalised over while learning the likelihood [111].

Of the 70,000 forward-modelled mocks, 60,000 are used for training, 7,500 for validation, and 2,500 are excluded from training (one footprint of one realisation per cosmology) as test set. The 1,250 of these belonging to the full-box Sobol set (§ 2.3) form the test set of § 4. We note that training on a finite number of simulations adds noise to the learned likelihood, similar to the Dodelson-Schneider effect [112] for covariance matrices estimated from simulations. Unlike for a Gaussian likelihood, however, this additional uncertainty propagates into the width of the posteriors [113].

The hyperparameters of the flows are tuned with optuna [114] using a TPESampler [115] for 300 trials, separately for every analysis variation, and we select the model with the highest likelihood on the validation set, i.e. the lowest “validation loss”. The search ranges, the fixed training settings and the selected hyperparameters of every analysis variation are listed in Appendix § 12. We use a single flow per analysis variation rather than an ensemble, as this already yields well-calibrated posteriors in the coverage tests (Appendix § 9.§ 1). As these test the calibration on average over the test set, we have also applied the five next-best trials of the joint analysis to the unblinded DES Y3 data vector. In units of the 68% half-width of the respective joint posterior, the posterior medians shift by at most 0.26 in $\Omega_{\mathrm{m}}$, 0.16 in $S_{8}$ and 0.14 in $w_{0}$, which is small compared to the statistical uncertainty.

Sampling posteriors

We use a flat prior over the parameter space spanned by CosmoGridV1 (Table 1), including the restrictions in the $\Omega_{\mathrm{m}}$–$\sigma_{8}$ and $\Omega_{\mathrm{m}}$–$w_{0}$ planes, which are enforced as a hard boundary with a log-probability of $-\infty$ outside the prior. The posteriors are then sampled by MCMC with the learned likelihood. For this, we have written a dedicated GPU implementation of the affine-invariant ensemble sampler with the stretch move ($a = 2$) [106] in PyTorch [116], which evaluates the likelihood of all walkers at once and runs up to 200 chains simultaneously on a single GPU. This replaces the emcee [117] sampling of Paper 1 and makes the thousands of posteriors needed for the tests of § 4 computationally feasible.

Each chain has 1,000 walkers, which are initialised in a Gaussian ball of 0.3 prior standard deviations around the CosmoGridV1 fiducial parameters for the DES Y3 measurement, or around the true parameters for the simulated measurements. For the simulated measurements (test and fiducial), each walker runs for 5,000 steps, of which the first 2,000 are discarded as burn-in and the remainder is thinned by 3, giving $10^{6}$ samples per posterior. For the DES Y3 measurement, each walker runs for 70,000 steps with 20,000 steps of burn-in and a thinning of 5, giving $10^{7}$ samples per posterior. With an integrated autocorrelation time of $\tau\approx390$ steps per walker for the DES Y3 joint chain, this corresponds to about $1.3 \times10^{5}$ independent samples. The simulated-measurement chains retain $\approx8\tau$ per walker, which is sufficient for the medians and 68% intervals used in the validation tests.

From the chains, we derive the parameter $S_{8} = \sigma_{8}\sqrt{\Omega_{\mathrm{m}}/0.3}$ and quote posterior medians with 68% credible intervals. In addition, we determine the maximum of the posterior (MAP), i.e. the best-fitting parameters, of every chain by maximising the log-posterior, the flow log-likelihood plus the prior, with the Adam optimiser [118] in PyTorch: starting from the 16 chain samples with the highest posterior probability, we take 300 steps with a learning rate of $10^{-2}$ and keep the best end point. The MAP is used for the goodness-of-fit tests of § 4 and quoted in § 5. The same pipeline, applied to the test simulations, is used for the pre-unblinding tests described in the next section.

Pre-unblinding tests

Before unblinding the cosmological constraints, the analysis had to pass a two-stage series of tests. First, the statistical calibration of the posteriors and the accuracy of the learned likelihoods were tested on simulated measurements. Second, the internal consistency between analysis variations, their goodness of fit, and posterior predictive distribution (PPD) tests [119] were evaluated on the blinded measurement. The analysis variations are the three main analyses, i.e. $\xi_{\pm}$, $\zeta_{\pm}$ and their combination, and five variations of the joint $\xi_{\pm} + \zeta_{\pm}$ data vector, in which either one tomographic bin (‘no BIN1’ to ‘no BIN4’) or the two smallest angular bins (‘no small scales’) are removed. All criteria were fixed before assessing the blinded results, and every test was passed before unblinding. In this section, we briefly summarise the tests and their outcome. Their full description, together with all figures and tables, is given in Appendix § 9, where we show the unblinded measurement and quote the corresponding blinded results in the text. After unblinding none of the test outcomes changed. The unblinding criteria, in the order in which they were checked, and their results are:

1. The coverage of all eight inference pipelines tested with TARP [120] lies within the $2\sigma$ binomial band (Appendix § 9.§ 1) with all pipelines passing the coverage tests. The largest deviation of the empirical coverage from the credibility level $\alpha$ is 0.007, compared to a $2\sigma$ band of 0.028 for which $\alpha= 0.5$.

2. The learned likelihoods are consistent with the $\chi^2$ expectation on the test set, verified by computing the squared Mahalanobis distance of each test data vector with respect to its own learned likelihood (this is a qualitative criterion, Appendix § 9.§ 2). While the test shows that the learned likelihoods are neither systematically too wide nor too narrow, the exaggerated tail reveals a (small) number of poorly fit cosmologies at the corners of the prior. The DES Y3 measurement itself lies at the 56th (joint), 60th ($\xi_{\pm}$) and 73rd ($\zeta_{\pm}$) percentile of $\chi^2$, i.e. the learned likelihoods describe the data measurement to be indistinguishable from a typical mock.

3. The blinded compressed data vector lies within the range spanned by the simulations and its log-posterior at the MAP falls inside the test-set distribution (qualitative, Appendix § 9.§ 4). The measurement lies at the 28th (joint), 20th ($\xi_{\pm}$) and 24th ($\zeta_{\pm}$) percentile of the test set, and the MOPED coefficients for $\Omega_m$, $\sigma_8$ and $w_0$ deviate from the mean of the learned likelihood at the MAP by at most $0.13\sigma$. The only coefficients deviating by more than $2\sigma$ are the redshift evolution $\nu$ of the baryonification (+2.9$\sigma$) and the shear calibration bias $m_3$ (−2.1$\sigma$).

4. The posteriors of the blinded measurement are consistent between the data-vector variations (qualitative, without a numerical threshold, Appendix § 9.§ 3). The largest shift of a posterior median relative to the joint analysis is the upward $S_8$ shift of ‘no BIN3’ (+2.2$\sigma$), where $\sigma$ denotes the 68% half-width of the joint posterior, followed by ‘no small scales’ (+1.3$\sigma$). All other shifts in the cosmological parameters are at most 1.1$\sigma$.

5. The calibrated PPD $p$-value of the test statistic $T^{\mathrm{cos}}$ on the cosmological MOPED coefficients, i.e. those of $\Omega_m$, $\sigma_8$ and $w_0$, satisfies $\tilde{p}(T^{\mathrm{cos}}) > 0.01$ for every PPD test, following the criterion used by DES [119] (Appendix § 9.§ 5). We find $\tilde{p}(T^{\mathrm{cos}}) = 0.08$–0.89 for the five data-vector variations, the $\xi_{\pm}$-only and $\zeta_{\pm}$-only analyses and the two cross-prediction tests between $\xi_{\pm}$ and $\zeta_{\pm}$, with the lowest value again for ‘no BIN3’. The parameters inferred from any subpart of the data vector therefore correctly predict the removed part of the measurement in $\Omega_m$, $\sigma_8$ and $w_0$.

Two of these results are relevant for the following sections. First, the same PPD statistic evaluated on the full compressed data vector (i.e. of all 17 MOPED coefficients), $\tilde{p}(T)$, is more marginal, with 0.021 for the goodness of fit of the joint analysis and 0.001 for ‘no BIN3’. This tension is concentrated in the $\nu$ and $m_3$ coefficients named above, while the cosmological coefficients show no tension. Our consistency criterion excludes these two directions, since the learned likelihood is not reliable in them: the parameter value of $m_3$ is never passed as input to the flow during training, and the data vector is barely sensitive to $\nu$, so that the likelihood in this direction is prior-dominated (Appendix § 9.§ 5). Second, the upward $S_8$ shift of ‘no BIN3’, the −2.1$\sigma$ deviation of $m_3$ and the lowest $\tilde{p}(T^{\mathrm{cos}})$ all indicate that the amplitude of the third tomographic bin is slightly lower than what the joint model predicts, which pulls $S_8$ down when this bin is included. Similarly, removing the two smallest angular scales, $15'$ and $30'$, raises $S_8$ by $1.3\sigma$ and lowers $M_{\rm c}$ by $0.7\sigma$, meaning these scales prefer a lower amplitude and stronger baryonic feedback. The same effect can be seen on the fiducial simulations, where the baryonic model is correct by construction and the shifts in $S_8$ and $M_{\rm c}$ are anti-correlated. Both of these shifts in $S_8$ and $w_0$ move along the $S_8$–$w_0$ degeneracy seen in Figure 1.

Cosmological constraints

In this section, we present the unblinded DES Y3 constraints of our three analyses on the cosmological parameters (§ 5.1) and on the systematic parameters (§ 5.2), and quantify the improvement of the joint analysis over $\xi_{\pm}$ alone (Section 5.3).

Constraints on $\Omega_m$, $S_8$ and $w_0$

Figure 1 shows the marginalised DES Y3 posteriors of $\Omega_m$, $S_8$ and $w_0$ for our three main analyses, and Table 2 lists the corresponding medians and 68% credible intervals of $\Omega_m$, $S_8$, $\sigma_8$ and $w_0$. The full posterior of all sampled parameters is shown in Appendix § 10. The joint $\xi_{\pm}+\zeta_{\pm}$ analysis clearly tightens the constraints compared to $\xi_{\pm}$ alone, most of all for $\sigma_8$ and in the $\Omega_m$–$S_8$ plane. The $\zeta_{\pm}$-only posterior is broader and fully consistent with the other two analyses though its degeneracy direction in the $\Omega_m$–$S_8$ plane differs from that of the $\xi_{\pm}$ posterior. $\xi_{\pm}$ constrains $S_8$ itself, while $\zeta_{\pm}$ constrains another combination of $\sigma_8$ and $\Omega_m$ ($\sigma_8\Omega_m^\alpha$ with $\alpha> 0.5$). This is in line with theoretical expectations, as at leading order $\zeta_{\pm}$ is proportional to $\sigma_8^2$ on large scales, while $\xi_{\pm}$ is sensitive to $\sigma_8^2$ [17, 49], and our findings in previous simulated analyses (Section 5.4 of Paper 1). The joint posterior lies at the intersection of the two, which is where the improvement in $\Omega_m$ and in the figure of merit come from. Adding $\zeta_{\pm}$ narrows and shifts the contours, as listed in Table 2: compared to $\xi_{\pm}$ alone, $S_8$ decreases by $0.6\sigma$, $w_0$ by $0.5\sigma$, and $\Omega_m$ increases by $0.4\sigma$.

Analysis	$\Omega_m$	$S_8$	$\sigma_8$	$w_0$	$\chi^2_{\rm MAP}/N_{\rm d}$
$\xi_{\pm}$ ($N_{\rm d}=160$)	$0.259^{+0.045}_{-0.042}$	$0.798^{+0.029}_{-0.027}$	$0.861^{+0.089}_{-0.074}$	$-0.93^{+0.25}_{-0.23}$	$1.11$ (54th)
$\zeta_{\pm}$ ($N_{\rm d}=1280$)	$0.341^{+0.105}_{-0.121}$	$0.780^{+0.035}_{-0.043}$	$0.740^{+0.138}_{-0.087}$	$-1.05^{+0.38}_{-0.35}$	$1.10$ (70th)
$\xi_{\pm}+\zeta_{\pm}$ ($N_{\rm d}=1440$)	$0.272^{+0.038}_{-0.033}$	$0.784^{+0.024}_{-0.021}$	$0.823^{+0.058}_{-0.052}$	$-1.04^{+0.25}_{-0.19}$	$1.11$ (69th)

The best-fitting parameters, i.e. the MAP in the 9-dimensional parameter space obtained by maximising the log-posterior (§ 3), at which the goodness-of-fit tests in § 4 evaluate the model, are listed in Table 10 in Appendix § 10. The MAP is a point estimate in the full parameter space and not a marginalised quantity, and the prior-dominated parameters lie at the prior boundaries, so that the MAP in these directions is not determined by the data. For the joint analysis, the MAP of $w_0$ lies $0.6\sigma$ below the median, which illustrates the projection effects we expect for posteriors that lie along the curved $\Omega_m$–$w_0$ degeneracy and are cut by the prior boundary: the medians in Table 2 are affected by the volume of the posterior, while the MAP is not.

Parameter	$\xi_{\pm}$	$\zeta_{\pm}$	$\xi_{\pm}+\zeta_{\pm}$
$\Omega_m$	0.269	0.322	0.243
$\Omega_b$	0.041	0.060	0.054
$H_0$	79.6	65.8	81.2
$n_s$	0.895	1.063	1.068
$\sigma_8$	0.869	0.748	0.871
$w_0$	$-0.72$	$-1.42$	$-1.18$
$A_{\mathrm{IA}}$	4.02	$-1.11$	2.45
$\log M_c$	15.00	14.76	14.55
$\nu$	0.57	0.04	1.99
$S_8$	0.823	0.775	0.784

Astrophysical systematic parameters

For the intrinsic alignment amplitude, we find $A_{\mathrm{IA}} = 2.7^{+1.4}_{-1.9}$ for the joint analysis, consistent with the $\xi_\pm$ analysis, while it is unconstrained for $\zeta_\pm$ alone. This is larger than the DES Y3 cosmic shear result [11] and than the values of the map-based SBI analyses [60, 61, 63], which lie between 0.2 and 0.6 (see § 6). In our case, the constraint on $A_{\mathrm{IA}}$ is driven by the amplitude of the low-redshift bins of $\xi_\pm$ (the ‘no BIN1’ variation loses the constraint entirely, see Table 7). The high $A_{\mathrm{IA}}$ has only a small effect on the cosmological parameters: $A_{\mathrm{IA}}$ and $S_8$ are uncorrelated in the joint posterior ($r = 0.15$), and removing the first tomographic bin loses almost all constraining power on $A_{\mathrm{IA}}$ while only slightly raising $S_8$, with a scatter that is comparable to leaving out any of the other tomographic bins (Table 8). A large $A_{\mathrm{IA}}$ could also point to a strongly aligned galaxy population within the DES Y3 source sample, which a colour/galaxy-property-split analysis [121, 122] could test. For baryonic feedback, we find $\log M_c = 13.9^{+0.8}_{-1.1}$, with the median in the upper half of the prior, i.e. only a mild preference for the strong-feedback half of the prior (Figure 11), while $\nu$ remains unconstrained. This is expected, since the scales $\ge15'$ contain little information on baryonic feedback.

Analysis	$\Omega_m$	$S_8$	$w_0$	$A_{\mathrm{IA}}$
$\xi_{\pm}+\zeta_{\pm}$	$0.272^{+0.038}_{-0.033}$	$0.784^{+0.024}_{-0.021}$	$-1.04^{+0.25}_{-0.19}$	$2.7^{+1.4}_{-1.9}$
no BIN1	$0.264^{+0.041}_{-0.033}$	$0.800^{+0.025}_{-0.023}$	$-0.94^{+0.24}_{-0.22}$	$-0.1^{+3.4}_{-3.3}$
no BIN2	$0.263^{+0.040}_{-0.032}$	$0.773^{+0.024}_{-0.020}$	$-1.09^{+0.22}_{-0.16}$	$2.9^{+1.3}_{-2.0}$
no BIN3	$0.254^{+0.041}_{-0.036}$	$0.834^{+0.029}_{-0.030}$	$-0.80^{+0.23}_{-0.25}$	$2.9^{+1.4}_{-2.2}$
no BIN4	$0.294^{+0.052}_{-0.044}$	$0.790^{+0.028}_{-0.026}$	$-0.94^{+0.29}_{-0.26}$	$2.9^{+1.4}_{-2.1}$
no small scales	$0.260^{+0.044}_{-0.034}$	$0.813^{+0.029}_{-0.027}$	$-0.90^{+0.26}_{-0.24}$	$2.0^{+1.8}_{-2.5}$

Variation	$\Omega_m$	$S_8$	$w_0$	$A_{\mathrm{IA}}$	$M_c$
no BIN1	$-0.23$	$+0.73$	$+0.47$	$-1.67$	$+0.02$
no BIN2	$-0.25$	$-0.47$	$-0.19$	$+0.10$	$+0.25$
no BIN3	$-0.50$	$+2.25$	$+1.13$	$+0.14$	$-0.89$
no BIN4	$+0.60$	$+0.28$	$+0.49$	$+0.10$	$-0.15$
no small scales	$-0.34$	$+1.30$	$+0.69$	$-0.40$	$-0.71$

Improvement from adding $\zeta_{\pm}$

Table 3 quantifies the improvement of the joint analysis over $\xi_{\pm}$ alone for the DES Y3 data vector, for the individual parameters and for the figure of merit (FoM). As in Paper 1, the improvement for an individual parameter is largest in $\sigma_8$ and smallest in $w_0$. This improvement for the DES Y3 data is typical: the FoM improvements lie between the 65th and the 77th percentile of the same quantity for the 800 fiducial noise realisations and for the 1250 test cases. We note that it is obtained on comparatively large angular scales, as the global and local $\xi_{\pm}$ start at 15′ and the smallest aperture radius is 50′, on which the forward model only weakly depends on baryonic feedback (see § 6).

	$\Omega_m$	$S_8$	$\sigma_8$	$w_0$	$\mathrm{FoM}(\Omega_m,S_8)$	$\mathrm{FoM}(\Omega_m,S_8,w_0)$	$\mathrm{FoM}(\Omega_m,\sigma_8,w_0)$
Improvement	22%	26%	49%	10%	54%	59%	99%

We note that all factors in Table 3 are ratios for the single DES Y3 data realisation. The median improvements over the simulated sets are a different quantity and consistent with Paper 1 [65], considering that retraining and retuning for every analysis variation (§ 3) leads to minor changes. We also obtain the precision on the DES Y3 data that we expected from Paper 1: the widths of the constraints differ by less than 20% from those of the simulated inference on a fiducial mock (Section 5.4 of Paper 1).² This agreement between the data and the simulated expectation of Paper 1, together with the analytical understanding of the i3PCF and its parameter dependence [49, 51, 52], gives us physical confidence in the SBI results beyond the statistical tests alone.

Comparison with external data

Figure 2 compares the joint, the $\xi_{\pm}$-only and the $\zeta_{\pm}$-only constraints with published $w$CDM results from the same source catalogue and from independent data. The constraints of every row are listed in Table 4, recomputed from the public chains where available and taken from the publications otherwise. Table 5 summarises the data, the likelihood, the treatment of intrinsic alignments and baryons, the amplitude and neutrino priors and the $w_0$ prior of every row, so that the differences in the priors and models between the analyses can be read from there.

Analysis	$\Omega_m$	$S_8$	$w_0$
$\xi_{\pm} + \zeta_{\pm}$ (this work)	$0.272^{+0.038}_{-0.033}$	$0.784^{+0.024}_{-0.021}$	$-1.04^{+0.25}_{-0.19}$
DES Y3 cosmic shear [11, 12]	$0.314^{+0.066}_{-0.056}$	$0.731^{+0.039}_{-0.029}$	$-1.37^{+0.46}_{-0.39}$
DES Y3 $3\times2$pt [10]	$0.350^{+0.039}_{-0.038}$	$0.776^{+0.024}_{-0.026}$	$-0.95^{+0.23}_{-0.31}$
DES Y3 $2+3$pt shear* [21]	$0.242^{+0.038}_{-0.037}$	$0.749^{+0.027}_{-0.026}$	$-1.39 \pm0.31$
DES Y3 moments+wavelets* [60]	$0.273 \pm0.029$	$0.817 \pm0.021$	$< -0.72$
DES Y3 map-level CNN* [63]	$0.283^{+0.020}_{-0.027}$	$0.804^{+0.025}_{-0.017}$	$< -0.80$
DES Y3 persistent homology* [61]	$0.304 \pm0.037$	$0.821 \pm0.018$	$< -0.75$
DES Y3 hybrid statistics* [64]	$0.325 \pm0.024$	$0.808 \pm0.017$	$< -0.766$
DES Y6 $3\times2$pt* [68]	$0.325^{+0.032}_{-0.035}$	$0.781^{+0.021}_{-0.020}$	$-1.12^{+0.26}_{-0.20}$
KiDS-1000 cosmic shear [69]	$0.34^{+0.10}_{-0.10}$	$0.743^{+0.049}_{-0.052}$	$-1.02^{+0.45}_{-0.81}$
KiDS-1000 $3\times2$pt [69]	$0.302^{+0.017}_{-0.018}$	$0.760^{+0.020}_{-0.021}$	$-1.03^{+0.10}_{-0.13}$
KiDS-Legacy cosmic shear* [70]	$0.332 \pm0.059$	$0.809 \pm0.041$	$-1.10^{+0.45}_{-0.47}$
Planck 2018 + BAO [71]	$0.304^{+0.011}_{-0.012}$	$0.826^{+0.010}_{-0.011}$	$-1.04 \pm0.05$
DESI DR2 BAO+CMB+SNe [72]	$0.305 \pm0.005$	$0.814^{+0.008}_{-0.007}$	$-0.995 \pm0.023$

Analysis	Data	Likelihood	IA	Baryons	$P(k)$ Ampl.	$\sum m_\nu$	$w_0$ prior
This work	$\xi_{\pm}, \zeta_{\pm}$	SBI, CosmoGrid	NLA (1)	$(M_c,\nu)$	flat	$\sigma_8$ fixed	$[-2,-1/3]^\dagger$
cosmic shear	$\xi_{\pm} + \mathrm{SR}$	analytic	TATT (5)	scale cuts	flat	$A_s$	marg.	$[-2,-1/3]$
$3\times2$pt	$3\times2\mathrm{pt} + \mathrm{SR}$	analytic	TATT (5)	scale cuts	flat	$A_s$	marg.	$[-2,-1/3]$
$2+3$pt shear*	$\xi_{\pm}, \langle M^3_{\mathrm{ap}}\rangle$	analytic	as DES Y3	scale cuts	flat	$A_s$	marg.	$[-2,-1/3]$
moments+wavelets*	map statistics	SBI, GS	NLA (2)	scale cuts	flat	$\sigma_8$	marg.	$(-1,-1/3)$
map-level CNN*	$C_\ell+ \mathrm{CNN}$	SBI, GS	NLA (2)	scale cuts	flat	$\sigma_8$	marg.	$(-1,-1/3)$
persistent homology*	Betti numbers	SBI, GS	NLA (2)	scale cuts	flat	$\sigma_8$	marg.	$(-1,-1/3)$
hybrid statistics*	$C_\ell+ \mathrm{neural}$	SBI, GS	NLA (2)	scale cuts	flat	$\sigma_8$	marg.	$(-1,-1/3)$
Y6 $3\times2$pt*	$3\times2\mathrm{pt}$	analytic	TATT (5)	scale cuts	flat	$A_s$	marg.	$[-2,-1/3]$
KiDS-1000	band powers	analytic	NLA (1)	HMCode16	flat	$S_8$	fixed	$[-3,-1/3]$
KiDS-1000 $3\times2$pt	+ clustering	analytic	NLA (1)	HMCode16	flat	$S_8$	fixed	$[-3,-1/3]$
KiDS-Legacy*	COSEBIs	analytic	mass-dep.	HMCode20	flat	$S_8$	marg.	$[-3,0]$
Planck 2018 + BAO	CMB, BAO	analytic	–	–	flat	$\ln A_s$	fixed	$[-3,1]$
DESI DR2	BAO, CMB, SNe	analytic	–	–	flat	$\ln A_s$	fixed	$[-3,1]$

Published results

The DES block contains analyses of the same Y3 source catalogue with different summary statistics and models. Any differences between these rows are therefore due to the methods and not due to the data. The DES Y6 $3\times2$pt result [68] is an update of the Y3 results and is not independent of them. It uses a new source catalogue and a new lens sample, and we quote its fiducial result. The DES Y3 $\xi_{\pm} + \langle M^3_{\mathrm{ap}}\rangle$ analysis [20, 21] is the closest analogue of our data vector, as it combines second- and third-order shear statistics, but using an analytical model and an explicit likelihood of Ref. [123]. Its $w_0$ prior follows the DES Y3 $\xi_{\pm}$ analyses and allows $w_0 < -1$, which gives a two-sided constraint ($w_0 = -1.39 \pm0.31$). The four map-based SBI analyses [60, 61, 63, 64], on the other hand, all use the Gower Street simulation suite, a $\nu w$CDM parameter space and a non-phantom $w_0$ prior that covers half our range. Because of this, they only quote upper limits (from $< -0.72$ to $< -0.8$, compatible with our result).

Our improvement over $\xi_{\pm}$ (Table 3) is obtained on comparatively large angular scales, as the global and local $\xi_{\pm}$ start at $15'$ and the smallest aperture radius is $50'$. The map-based DES Y3 higher-order analyses instead include much smaller scales in their forward model, working on HEALPIX maps with a pixel scale of $6.9'$ up to $\ell\approx1000$ [60, 61, 63, 64] or, with a smoothing that depends on the tomographic bin, up to $\ell\approx600$–$1400$ [62]. For this, they rely on dedicated robustness tests against baryonic feedback or, in Ref. [62], on the marginalisation over the baryonification parameters. In Paper 1 we instead chose a larger scale cut that excludes the scales most affected by baryonic feedback and resolution effects from the simulations, at the cost of some non-Gaussian information. Our improvement over $\xi_{\pm}$ is therefore not due to additional small-scale information, but comes from the third-order information on scales where the forward model mildly depends on baryonic feedback, which can also be seen in the weak constraint on $M_{\mathrm{c}}$ (§ 5.2).

The independent block contains results from the KiDS-1000 beyond-$\Lambda$CDM analysis [69], which uses band powers instead of COSEBIs, both for cosmic shear alone and for 3$\times$2pt with BOSS and 2dFLenS clustering, and KiDS-Legacy from the extended-cosmology analysis [70] of the fiducial COSEBIs setup [13]. We show Planck 2018 [71] together with BAO, because the CMB alone cannot constrain $w_0$ due to the geometric degeneracy. DESI DR2 BAO [72] is shown in combination with the CMB (Planck PR4 CamSpec and CMB lensing) and the Pantheon+ supernovae sample that does not include DES data.

Tension with independent data

We quantify the tension with the independent (external) data with the parameter-shift probability introduced in [124, 125]. It is based on the distribution of the parameter difference between two independent posteriors,

which is integrated above its zero-shift iso-contour,

and converted to an effective significance $n_\sigma= \sqrt{2}\,\operatorname{Erf}^{-1}(\Delta)$. We compute it with the kernel density estimate of tensiometer [125, 126] on our chains and each public chain, both in the $(\Omega_{\mathrm{m}}, S_8)$ plane and in the full $(\Omega_{\mathrm{m}}, S_8, w_0)$ space, and list the resulting shifts for the three main analyses in Table 6. For the joint analysis, the shifts are below $1\sigma$ for both KiDS-1000 results and reach $2.4\sigma$ for Planck 2018 + BAO and $1.6\sigma$ for DESI DR2 in this three-dimensional parameter space. These larger shifts come from the $(\Omega_{\mathrm{m}}, w_0)$ plane and not from the $w_0$ marginal, for which the simple marginal posterior estimate described below gives $\leq0.2\sigma$. The $\xi_{\pm}$-only posterior, which is broad and centred at higher $\Omega_{\mathrm{m}}$, is compatible with all of them, while the $\xi_{\pm}$-only posterior shows a similar $\approx2\sigma$ three-dimensional shift with respect to Planck 2018 + BAO as the joint analysis. For KiDS-Legacy, which has no public chain, we use the simple marginal posterior estimate $|\Delta\theta|/\sqrt{\sigma_1^2+\sigma_2^2}$ [127], with the marginal means and symmetrised 68% intervals, as published for KiDS-Legacy, and find shifts below $1\sigma$ (Table 6). For results with a public chain, we use the medians instead.

External data	$\xi_\pm+ \zeta_\pm$ 2D	$\xi_\pm+ \zeta_\pm$ 3D	$\xi_\pm$ 2D	$\xi_\pm$ 3D	$\zeta_\pm$ 2D	$\zeta_\pm$ 3D
KiDS-1000 cosmic shear	$0.4\sigma$	$0.3\sigma$	$0.6\sigma$	$0.4\sigma$	$0.2\sigma$	$0.1\sigma$
KiDS-1000 $3\times2\mathrm{pt}$	$0.8\sigma$	$0.7\sigma$	$1.0\sigma$	$0.8\sigma$	$0.4\sigma$	$0.1\sigma$
Planck 2018 + BAO	$1.2\sigma$	$2.4\sigma$	$0.9\sigma$	$2.2\sigma$	$1.3\sigma$	$1.1\sigma$
DESI DR2	$1.0\sigma$	$1.6\sigma$	$0.7\sigma$	$1.6\sigma$	$0.9\sigma$	$0.6\sigma$

Overall, we do not find a significant tension with any of the independent results. The $S_8$ of the joint analysis lies below the Planck 2018 + BAO and DESI DR2 values, which tends in the direction of the low value of $S_8$ found in weak lensing surveys compared to CMB analyses. However, we note that this should only be interpreted bearing in mind the different priors and parameter spaces used in the analyses, in particular the $w_0$ ranges.

Comparison with DES Y3 analyses

For the other DES analyses, we do not compute a tension metric, as these analyses use the same source catalogue and are therefore correlated with our results. Using the simple estimate above, the $S_8$ of the joint analysis lies within $0.3\sigma$ of DES Y3 $3\times2\mathrm{pt}$, $1.3\sigma$ above DES Y3 cosmic shear and $0.6$–$1.2\sigma$ below the four map-based SBI analyses, which find $S_8=0.80$–$0.82$.

Part of this difference can be explained by their non-phantom $w_0$ prior: in our posterior, $S_8$ and $\Omega_m$ are positively correlated with $w_0$ (Figure 1), and restricting it to $w_0>-1$ gives posterior means of $S_8=0.803$ and $\Omega_m=0.290$. The other analysis combining second- and third-order shear statistics, $\xi_\pm+\langle M_{\mathrm{ap}}^3\rangle$ [21], which also allows $w_0<-1$, lies within $0.6\sigma$ in $\Omega_m$ and $1.0\sigma$ in $S_8$.

Figures 3 and 4 show the $\xi_\pm$-only posterior against the public DES Y3 cosmic shear $w\mathrm{CDM}$ chain, and the joint analysis against the DES Y3 $3\times2\mathrm{pt}$ chain and Planck 2018 + BAO, in the $\Omega_m$–$S_8$–$w_0$ space. Figure 3 is a direct comparison of our $\xi_\pm$-only analysis with the official DES Y3 cosmic shear result, which uses the same statistic and the same catalogue, so that the differences between the two are only in the modelling (Table 5) and in the MOPED compression. Our $\xi_\pm$-only posterior is tighter, which is mainly due to the smaller scale cuts: the DES Y3 analysis removes the scales of $\xi_\pm$ on which baryonic feedback affects the analytical model, and for $\xi_-$ these reach well above our lower limit of $15'$, while in our case the feedback is part of the forward model and marginalised over. The fixed neutrino mass and the one-parameter IA model also contribute to the tighter posterior. The IA amplitude of our joint analysis, $A_{\mathrm{IA}}=2.7^{+1.4}_{-1.9}$ (§ 5.2), is $1.2\sigma$ larger than the DES Y3 cosmic shear result $a_1=0.34^{+0.25}_{-0.23}$ [11], obtained with the same NLA model without redshift evolution or shear ratios (their $a_1$ corresponds to our $A_{\mathrm{IA}}$), and also larger than the values of the map-based SBI analyses $(0.37^{+0.17}_{-0.29}$ [60], $0.21^{+0.17}_{-0.25}$ from the additional $A_{\mathrm{IA}}$ compression of Ref. [63], and $0.64 \pm0.35$ [61]), which use the NLA model with a free redshift evolution. Compared to all of these, our $A_{\rm IA}$ is, however, considerably less constrained, and the difference may be caused by the different priors, scale cuts and treatments of baryonic feedback. The $\xi_{\pm}$-only posterior lies at the high-$S_8$ edge of the 68% contour of DES Y3 cosmic shear, whose broad tail towards low $S_8$ causes the $1.3\sigma$ difference quoted above. In Figure 4, the joint analysis overlaps with the $3\times2\mathrm{pt}$ contour in $S_8$ and lies at lower $\Omega_{\rm m}$, and Planck 2018 + BAO lies in the high-$\Omega_{\rm m}$, high-$S_8$ corner of the joint 95% contour, which is in the $(\Omega_{\rm m}, w_0)$ direction that is responsible for the $2.4\sigma$ shift. As for $A_{\rm IA}$, the differences in $\Omega_{\rm m}$, $S_8$ and $w_0$ from the DES Y3 cosmic shear and $3\times2\mathrm{pt}$ analyses may be due to modelling choices such as the different priors, scale cuts and treatments of baryonic feedback. In particular, $w_0$ is only weakly constrained and correlated with $\Omega_{\rm m}$ and $S_8$, so that differences in the $w_0$ prior volume, including the restriction in the $\Omega_{\rm m}$–$w_0$ plane (§ 2.3), can shift the marginalised constraints through projection effects.

Conclusion

Maximising the scientific return from Stage-IV cosmic shear surveys like Euclid and LSST will require the analysis of higher-order statistics alongside $\xi_{\pm}$. However, before such statistics can be applied at Stage-IV precision, it is important to demonstrate that they work on Stage-III survey data. This requires a realistic forward model that includes the relevant systematic effects, together with validation tests on blinded data.

In this work, we have presented the first cosmological analysis of the integrated 3-point correlation function on survey data. We have used SBi3PCF, the simulation-based inference framework developed and validated in Paper 1, for a blinded analysis of the DES Y3 weak lensing shape catalogue and obtained $w$CDM posteriors for the cosmic shear statistics $\xi_{\pm}$, $\zeta_{\pm}$ and their combination. Our main findings are as follows:

1. We report the first detection of $\zeta_{\pm}$ on survey data, with a signal-to-noise ratio of S/N = 12.8 for $\zeta_{\pm}$ alone and of 33.0 for the joint $\xi_{\pm} + \zeta_{\pm}$ data vector (§ 2).

2. All pre-unblinding tests, i.e. coverage, accuracy of the learned likelihoods, internal consistency and goodness of fit of the compressed and the uncompressed data vector, were passed on the blinded measurement. The internal consistency tests between parts of the data vector and between the two statistics are calibrated on simulations instead of relying on raw $p$-values (§ 4).

3. For the joint $\xi_{\pm} + \zeta_{\pm}$ analysis, we find $\Omega_{\mathrm{m}} = 0.272^{+0.038}_{-0.033}$, $S_{8} = 0.784^{+0.024}_{-0.021}$ and $w_{0} = -1.04^{+0.25}_{-0.19}$. We also find the intrinsic alignment amplitude $A_{\mathrm{IA}} = 2.7^{+1.4}_{-1.9}$ and a weak constraint on baryonic feedback, $\log M_{\mathrm{c}} = 13.9^{+0.8}_{-1.1}$. Adding $\zeta_{\pm}$ improves the figure of merit in the $\Omega_{\mathrm{m}}$–$S_{8}$–$w_{0}$ parameter space by a factor of 1.59 over $\xi_{\pm}$ alone (Table 3), in line with the expectation from the simulated measurements of Paper 1. We note that it is obtained using only angular scales $\ge15'$ and aperture radii $\ge50'$, meaning that we do not rely on the small scales that the map-based higher-order analyses use (§ 5).

4. Our constraints agree with the other DES Y3 analyses of the same source catalogue and with the DES Y6 update, the $\xi_{\pm}$-only analysis is consistent with the official DES Y3 cosmic shear result in a direct comparison, and we find no significant tension with independent data, with parameter shifts of at most $2.4\sigma$ for Planck 2018 + BAO (§ 6).

The present analysis has some limitations. The forward model is based on a single simulation suite, and we could not test the recovery of the input cosmology on simulations from an independent N-body code. The prior is inherited from CosmoGridV1 and spans only $w$CDM at a fixed total neutrino mass $\sum m_{\nu}$, so that there is no corresponding $\Lambda$CDM analysis. Intrinsic alignments are only infused with the NLA model, with no redshift dependence ($\eta= 0$), but we find no impact of $A_{\mathrm{IA}}$ on cosmology. The MOPED compression vectors are fixed at one fiducial cosmology, which can lead to a loss of information, but not to a bias (§ 3). Since the fiducial values lie within the 68% credible intervals of the joint analysis, we expect this loss to be small. The $\nu$ and $m_{3}$ directions of the learned likelihood are not reliable and were not used for the consistency criterion in § 4, while the cosmological coefficients show no tension. Finally, source clustering and reduced shear were left out of the forward model based on the tests of Paper 1, which found their effect negligible.

We have thus shown that $\zeta_{\pm}$ can be analysed on Stage-III survey data with a fully simulation-based approach, an important step towards Stage-IV surveys like Euclid and Rubin LSST. At their precision, the map-level forward model will need tidal-field-based alignment models [75], the inclusion of photometric redshift outliers [78] and a treatment of baryonic feedback beyond two parameters [77], as well as larger simulation suites or emulators of the data vector to reduce the noise in the learned likelihood due to the finite number of simulations [113]. These extensions can build on the pipeline presented here, as can joint analyses with other higher-order statistics of the shear maps or with statistics of photometric galaxy clustering maps. Together with the comparison with the analytic-likelihood $\zeta_{\pm}$ analysis of the same data (Halder et al. in preparation), they can contribute to the growing efforts in the weak lensing community towards robust HOS analyses that more fully exploit the non-Gaussian information encoded in the large-scale structure of the Universe.

Acknowledgments

We would like to thank Lina Castiblanco, Oliver Friedrich, Marco Gatti, Daniel Gruen, Joachim Harnois-Dér̨aps, Jed Homer, Mike Jarvis, Tomasz Kacprzak, Kai Lehman, Nils Schöneberg, Arne Thomsen, Simone Vinciguerra, and especially the late Stella Seitz for fruitful discussions at various stages of this project. We are grateful to Tomasz Kacprzak and the CosmoGridV1 team for making the simulations publicly available. We acknowledge support from the Excellence Cluster ORIGINS, which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094-390783311. Some of the numerical calculations have been carried out on the ORIGINS computing facilities of the Computational Center for Particle and Astrophysics (C2PAP). DG and CU were supported by the European Union (ERC StG, LSS_BeyondAverage, 101075919). DG acknowledges support from the German Academic Exchange Service (DAAD) through a Research Scholarship for Doctoral Candidates, and thanks the Institute of Astronomy, University of Cambridge, for hosting him as a visiting student. AH has been supported by funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programmes (grant agreement no. 101018897 CosmicExplorer). ZG acknowledges support via the funding from the U.S. Department of Energy Cosmic Frontier program, grant DE-SC002599.

This project used public archival data from the Dark Energy Survey (DES). Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana–Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey. The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas–Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l’Espai (IEEC/CSIC), the Institut de Física d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, The Ohio State University, the OzDES Membership Consortium, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, and Texas A&M University. Based in part on observations at Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.

During the preparation of this work, the authors used Claude Code and several Claude models (Anthropic, 2026 versions) to assist with writing analysis and plotting code and to improve the language and readability of parts of the text. All AI-assisted code was reviewed and validated by the authors, and all figures are produced directly from the measured data and simulations. The authors reviewed and edited all AI-assisted text and take full responsibility for the content of the manuscript.

This work made use of numpy [128], scipy [129], matplotlib [130], ChainConsumer [131], PyTorch [116], sbi [110], zuko [109], optuna [114], TARP [120] and TreeCorr [89, 91].

Data availability: The full analysis code is publicly available at SBi3PCFII, and the data products required to reproduce every figure and quoted number of this paper are archived on Zenodo (https://doi.org/10.5281/zenodo.22018495). The training and covariance sets can be reconstructed from the public CosmoGridV1 simulations [66], with the weak lensing maps built using UFalcon2 [98, 99, 100], the $\xi_{\pm}$ measurements performed with TreeCorr [89, 91] and the $\zeta_{\pm}$ measurements with PixPoint, following the procedure of Paper 1. These sets can also be shared by the authors upon request.

Pre-unblinding tests

In this appendix, we describe the pre-unblinding tests that are summarised in § 4. The first two tests are carried out on the simulated measurements of the test set, the remaining three on the blinded data vector, and all of them had to be passed before unblinding. Where the DES Y3 measurement enters, the figures show the unblinded measurement and the corresponding blinded results are quoted in the text.

Coverage tests

We verify that the posteriors of all eight inference pipelines used in this work are statistically well calibrated with the TARP [120] coverage test, as in Paper 1. The test is evaluated on $N = 1250$ test simulations that are distributed following the prior and were excluded from training (Paper 1 used 500), with the reference points of TARP drawn uniformly over the prior box. The empirical coverage probability $\operatorname{ECP}(\alpha)$ is the fraction of the test inferences for which the true parameter values are enclosed within the $\alpha$-level credible interval, and the binomial band $\pm\sqrt{\alpha(1-\alpha)/N}$ around the diagonal gives its expected scatter for $N$ simulations.

Figure 5 shows that all curves lie within the 1 (2) $\sigma$ binomial band around the one-to-one line, so that the posteriors of the main models and of the data-vector variations are statistically well calibrated. The largest deviation is $\lvert\operatorname{ECP}(\alpha)-\alpha\rvert= 0.007$ for ‘no BIN3’, compared to a $2\sigma$ binomial band of 0.028 at $\alpha= 0.5$, and the mean deviations lie between $-0.003$ and $+0.002$, i.e. the curves scatter around the diagonal without a systematic offset.

We note that passing such a coverage test is a necessary condition for accurate posteriors, but not a sufficient one. The equivalence between exact expected coverage and posterior accuracy shown in Ref. [120] only holds in the limit of infinitely many test simulations and every possible reference point, so that a passed test on a finite test set limits, but does not exclude, a residual miscalibration of the posteriors.

Accuracy of the learned likelihoods

While the coverage test checks the posteriors, we also test the learned likelihoods directly. For this, we compute the squared Mahalanobis distance of each test data vector with respect to its own learned likelihood $p(x \mid\theta_{\mathrm{true}})$, which for a Gaussian likelihood is exactly the $\chi^{2}$, and compare its distribution over the test set with the $\chi^{2}$ distribution expected for the $d=17$ elements of the compressed data vector. This checks whether the widths of the learned likelihoods are correct. The distance is evaluated in the space of the standardised MOPED coefficients, with the mean and covariance of each likelihood estimated from 1000 draws of the flow and the inverse covariance corrected with the Hartlap factor [107]. Figure 6 shows the resulting distributions for the three main analyses. The medians per degree of freedom are 0.93–0.94, which shows that the learned likelihoods are not systematically too wide or too narrow, while the means (1.18–1.28) are raised by a tail of poorly fit cosmologies: 7–9% of the test cases lie above the 95% quantile of $\chi^{2}(17)$, and the 9–12 cases with $\chi^{2}>100 (outside of the plotted range) all lie at the corners of the prior, where the training set is sparsest.

We compute the same statistic for the DES Y3 measurement with respect to the learned likelihood at its MAP and find 17.1 (joint), 17.8 ($\xi_{\pm}$) and 20.0 ($\zeta_{\pm}$) for 17 degrees of freedom, corresponding to the 56th, 60th and 73rd percentile of $\chi^2(17)$. For the blinded data vector, the values were 20.5, 21.1 and 24.7. The learned likelihoods thus describe the DES Y3 measurement as well as a typical mock.

Internal consistency tests

Figure 7 shows the DES Y3 posteriors for the three main analyses and the five variations of the joint $\xi_{\pm} + \zeta_{\pm}$ data vector, in which either one tomographic bin or the two smallest angular bins are removed. All variations are consistent with the joint posterior. To quantify this, Table 8 lists the shifts of the posterior medians of the five variations relative to the joint analysis in units of the 68% half-width of the joint posterior,

which is the metric of Section 5.3 of Paper 1 with the joint median in place of the true value. Positive values lie above the joint median. The largest shift is the upward $S_8$ shift of ‘no BIN3’ (+2.2), followed by ‘no small scales’ (+1.3); all other shifts in the cosmological parameters are at most 1.1, and the largest shift in a systematic parameter is −1.7 in $A_{\mathrm{IA}}$ for ‘no BIN1’. We note that these should not be read as significances, since removing part of the data changes both the position and the width of the posterior. The quantitative consistency test is the posterior predictive test of Appendix § 9.§ 5.

Table 7 lists the marginalised constraints shown in Figure 7. Relative to the joint analysis, the variations widen the constraints by at most a factor of 1.35 in $\Omega_{\mathrm{m}}$ and 1.26 in $w_0$ (‘no BIN4’) and 1.33 in $S_8$ (‘no BIN3’), and ‘no BIN1’ loses the constraint on $A_{\mathrm{IA}}$ (×2.0, with a median compatible with zero), as the lowest-redshift bin carries most of the intrinsic alignment signal. On the blinded measurement, we found the same pattern of shifts (‘no BIN3’ +2.29 and ‘no small scales’ +1.28 in $S_8$, $A_{\mathrm{IA}}$ −1.66 for ‘no BIN1’), so that the outcome of this test is unchanged by unblinding.

The shifts of the variations are also consistent with the other tests in this appendix. The upward shift in $S_8$ for ‘no BIN3’ matches the −2.1$\sigma$ deviation of the $m_3$ coefficient in Figure 9 and the lowest $\tilde{p}(T)$ of the posterior predictive tests below (Table 9), all indicating that the amplitude of the third tomographic bin is slightly lower than what the joint model predicts. Similarly, removing the two smallest angular bins raises $S_8$ (+1.3) and lowers $M_{\mathrm{c}}$ (−0.7), i.e. the scales between 15′ and 30′ prefer a slightly lower amplitude, as one would expect from baryonic suppression. As discussed in § 4, both shifts move along the $S_8$–$w_0$ degeneracy seen in Figure 1, and we consider them to be a single shift rather than two separate ones.

Conditioning test	$\tilde{p}(T^{\mathrm{cos}})$	$\tilde{p}(T)$
no BIN1	0.391	0.018
no BIN2	0.741	0.024
no BIN3	0.084	0.001
no BIN4	0.889	0.043
no small scales	0.253	0.015
$\xi_{\pm}$	0.651	0.064
$\zeta_{\pm}$	0.695	0.394
$\zeta_{\pm}\mid\xi_{\pm}$	0.511	0.068
$\xi_{\pm}\mid\zeta_{\pm}$	0.235	0.438
$\xi_{\pm}+\zeta_{\pm}$	0.964	0.021

Likelihood of the observed data vector

Figure 8 compares the log-posterior at the MAP of the DES Y3 measurement with the distribution of the same quantity over the test set, for which the MAP of each of the 1250 test simulations is evaluated from its own posterior in the same way as for the measurement (§ 3). Since the prior is flat within its support, the differences within a panel are differences in the log-likelihood. The measurement lies at the 28th (joint), 20th ($\xi_{\pm}$) and 24th ($\zeta_{\pm}$) percentile of the test-set distributions, and the blinded data vector at the slightly lower 23rd, 16th and 21st percentile. The observed data vector is therefore a typical realisation of the forward model, and we find no indication that the forward model does not describe the data.

Figure 9 repeats this check for each of the nine MOPED coefficients corresponding to the model parameters, with the learned likelihood at the MAP represented by 1000 draws of the joint flow. We can see that the coefficients for $\Omega_m$, $\sigma_8$ and $w_0$ deviate from the mean of the learned likelihood by at most $0.13\sigma$ and are therefore typical draws from it. All other coefficients, including the eight shear calibration and photo-$z$ coefficients not shown in the figure, lie within $1.3\sigma$, with the exception of the redshift evolution $\nu$ of the baryonification ($+2.9\sigma$) and the shear calibration bias $m_3$ ($-2.1\sigma$), to which we come back in the posterior predictive tests below.

Posterior predictive consistency tests

Complementary to the comparison of the posteriors in Figure 7, we test the internal consistency of the data-vector variations, and of the two statistics with each other, with posterior predictive distribution (PPD) tests following the DES methodology of Ref. [119]. The idea behind these tests is to run the inference with part of the data removed, to predict the removed part from the resulting posterior, and to then compare this prediction to the actual measurement using a test statistic. The raw PPD $p$-values obtained in this way are not uniformly distributed, but skewed towards low values if the two parts of the data vector constrain different directions in parameter space. Following Ref. [119], we therefore calibrate them on simulated realisations, which are consistent by construction: the calibrated $\tilde{p}$ is the fraction of the 800 fiducial noise realisations that have a larger test statistic than the measurement, and as in DES, we require $\tilde{p} > 0.01$ for consistency.

Ref. [119] uses an analytical Gaussian likelihood, draws realisations of the removed data from the PPD and compares them to the measurement with a $\chi^2$-like distance. Here, the likelihood is only available as the trained flow, so we use the likelihood of the actual measurement as the test statistic instead,

which is the joint flow evaluated at the full MOPED data vector and averaged over 2000 thinned chain samples $\tilde{\theta}$ from the posterior $\theta_{\mathrm{post}}$ of the analysis with the reduced data vector. With this choice, we need neither realisations drawn from the PPD nor a distance measure in data space. As the normalisation of the flow is arbitrary, the absolute values of $T$ have no meaning on their own, but it cancels out in the calibration, and we only ever compare the calibrated $\tilde{p}$. We also note that the full data vector still contains the part of the data on which the posterior was conditioned, which biases the test toward agreement. This bias is present in the simulated measurements in the same way and is therefore also removed by the calibration. We compute $T$ for the five data-vector variations, for the two single-statistic analyses (the joint flow with the $\xi_{\pm}$-only or the $\zeta_{\pm}$-only posterior, with the other statistic as the removed data), and for the joint analysis itself, which amounts to a goodness-of-fit test.

For $\xi_{\pm}$ and $\zeta_{\pm}$, we additionally perform two cross-prediction tests in which the data used for the posterior and the data at which the likelihood is evaluated do not overlap, which is the direct analogue of the consistency tests in Ref. [119]. Here, the flow of statistic $A$ is evaluated at its own data vector using parameters that were inferred from the other statistic $B$ alone,

with $\tilde{\theta}$ drawn from the posterior of statistic $B$. In this case, there is no bias towards agreement, since the data at which the likelihood is evaluated never enters the posterior, and the correlated noise of the two statistics, which are measured on the same maps, is again absorbed by the calibration. These tests differ from the single-statistic tests above only in how the removed statistic enters: there as one part of the jointly compressed MOPED coefficients, here through a separate compression.

Since the marginal distributions of the flow over a subset of the coefficients are not available, we also evaluate a Gaussian version of the statistic for the three coefficients corresponding to $\Omega_m$, $\sigma_8$ and $w_0$, with the mean and covariance estimated from one draw of the flow per chain sample. We denote this statistic by $T^{\mathrm{cos}}$ and show its calibration in Figure 10. It directly tests whether all parts of the data vector prefer the same values for the parameters that we report.

All tests on the cosmological coefficients are passed (Table 9) with $\tilde{p}(T^{\mathrm{cos}}) = 0.08$–$0.89$. The value is lowest for ‘no BIN3’, where it is still eight times the threshold. The tests conditioning the full data vector on $\xi_{\pm}$ and $\zeta_{\pm}$ give 0.65 and 0.70 respectively, and the separate cross-predictions give 0.51 and 0.24. We can therefore conclude that the parameters inferred from any of the subparts of the data vector correctly predict the removed part of the measurement in $\Omega_m$, $\sigma_8$ and $w_0$.

For the statistic $T$ on the full data vector, the results are more marginal, with the lowest $\tilde{p}(T) = 0.001$ for ‘no BIN3’ again, and even the goodness-of-fit test of the joint analysis gives only 0.021, which is lower than the 28th percentile of the MAP log-posterior in Figure 8 and the 56th percentile of the $\chi^2$ of Appendix § 9.§ 2. There are two reasons for that. Firstly, the fiducial realisations used for the calibration do not contain the poorly fit cosmologies at the corners of the prior, which widen the test-set distributions. And secondly, $T$ is evaluated with the full density fitted by the flow and not with its Gaussian approximation, in which a deviation in a single direction is diluted over the 17 degrees of freedom. Looking at individual coefficients, we find that the tension is concentrated in the $+2.9\sigma$ deviation in $\nu$ and the $-2.1\sigma$ deviation in $m_3$ of Figure 9, while the cosmological coefficients show no tension. However, we do not base the consistency test on these two directions, as the learned likelihood is not reliable in them: during training we never pass the parameter value of $m_3$ to the flow, meaning it cannot learn the dependence of its corresponding MOPED coefficient, and since the data vector is barely sensitive to $\nu$, the likelihood in this direction is dominated by the prior. In these directions, the calibrated test does therefore not probe the internal consistency, but only how close the (noise-dominated) scatter of the observed data vector lies to the centre of the prior, and we base the consistency criterion on $T^{\mathrm{cos}}$.

Comparing the two statistics further supports this interpretation. All tests that are conditioned on data containing $\xi_{\pm}$ have a marginal $\tilde{p}(T)$ (0.001–0.07), whereas the two tests conditioned on $\zeta_{\pm}$ have the highest values of all tests (0.39 and 0.44), as one would expect if the $\nu/m_3$ tension enters through the $\xi_{\pm}$ part of the data vector. We note though, that the $\zeta_{\pm}$-only posterior is much wider, and that conditioning on it is therefore a weaker test. The more meaningful direction is predicting $\zeta_{\pm}$ from $\xi_{\pm}$, for which we find $\tilde{p}(T^{\mathrm{cos}}) = 0.51$ and $\tilde{p}(T) = 0.068$. On the blinded data vector, the results for $T^{\mathrm{cos}}$ were the same, with all tests passing in both cases, while the values of $\tilde{p}(T)$ were lower throughout (e.g. ‘no BIN1’ 0.005 and ‘no small scales’ 0.006), meaning the $\nu/m_3$ tension is slightly weaker in the unblinded data vector, which is consistent with the MAP percentiles above.

Full Posterior

Figure 11 shows the complete 1D and 2D marginalised posterior distributions of all nine sampled parameters and the derived $S_8$ for the $\xi_{\pm}$-only (blue), the $\zeta_{\pm}$-only (orange) and the joint (green) analysis, complementing the results for the main cosmological parameters of § 5. The regions excluded by the prior are shown in grey; they include the restrictions in the $\Omega_m$–$\sigma_8$ and $\Omega_m$–$w_0$ planes of § 2.3. The 1D prior marginals on the diagonal are obtained from a chain sampling the prior and rescaled to the panel height. Table 10 lists the best-fitting parameters of the three analyses, i.e. the MAP in the 9-dimensional parameter space (§ 3), at which the goodness-of-fit tests of § 4 evaluate the forward model.

To quantify which parameters are constrained by the data, we compute the ratio of the 68% widths of the posterior and of the prior chain. $\Omega_{\mathrm{b}}$, $H_0$, $n_{\mathrm{s}}$, $M_{\mathrm{c}}$ and $\nu$ are prior-dominated in all three analyses (ratios $\ge0.90$, with the median of $M_{\mathrm{c}}$ in the upper half of the prior), and as in Paper 1, the data vector adds little information on these parameters. In the joint analysis, $A_{\mathrm{IA}}$ (0.49) and $w_0$ (0.57) are partially constrained, while the well-constrained parameters $\Omega_{\mathrm{m}}$, $\sigma_8$ and $S_8$ (0.30, 0.19 and 0.11) are reported in § 5. For $\xi_{\pm}$ alone, the ratios are slightly larger (0.37, 0.29 and 0.14 for $\Omega_{\mathrm{m}}$, $\sigma_8$ and $S_8$, 0.63 for $w_0$ and 0.52 for $A_{\mathrm{IA}}$), while $\zeta_{\pm}$ on its own only constrains $\sigma_8$ (0.40) and $S_8$ (0.19) and leaves $\Omega_{\mathrm{m}}$, $w_0$ and $A_{\mathrm{IA}}$ at the prior width.

Uncompressed DES Y3 measurements

Figures 12, 13 and 14 show the measured DES Y3 $\xi_{\pm}$ and $\zeta_{\pm}$ entering the compression for all tomographic bin combinations, where for the $\zeta_{\pm}$ combination $(i,j,k)$ $M_{\mathrm{ap}}$ is measured in bin $i$ and the local $\xi_{\pm}$ between bins $j$ and $k$. We show $\zeta_{\pm}$ for the aperture radius $\theta_{\mathrm{ap}}=90'$ , as the other radii look very similar; their dependence on the aperture radius is shown in Figure 15 for the combination $(3,3,4)$, which, together with $(2,2,2)$, is one of the two combinations with the highest signal-to-noise ratio of all twenty, with the number of angular bins growing with the radius as $\alpha_{\max}\leq2\theta_{\mathrm{ap}}-5'$. The measurements are scaled by the angular separation, the error bars show the scatter of the 16,000 fiducial footprints (the diagonal of the covariance), and the two smallest angular bins of every correlation, which are removed in the ‘no small scales’ variation, are shaded. In Figure 12, the posterior predictive means of the joint and the $\xi_{\pm}$-only analysis are nearly indistinguishable, which is the small difference between the two analyses of Table 2 seen in data space.

As we do not have an emulator, we use the posterior predictive mean of the forward-modelled mocks as the model prediction in data space, i.e. the mean of the 67,500 training and validation data vectors weighted by a Gaussian kernel in parameter space that is centred on the MAP and has the marginal posterior widths of the respective analysis as its widths. The effective number of mocks contributing to this mean, $N_{\mathrm{eff}=(\sum_k w_k)^2/\sum_k w_k^2}$ with the kernel weights $w_k$ [132], is 70 for the joint analysis and 790 and 1627 for the $\xi_{\pm}$-only and $\zeta_{\pm}$-only posteriors, as their wider kernels include more of the training set. The goodness of fit, $\chi^2=(x-\mu)^T C^{-1}(x-\mu)$ of the measurement with respect to this prediction, uses the inverse of the fiducial covariance of that analysis corrected by the Hartlap factor [107], and the percentiles in Table 2 are computed from the same statistic for the MAP cosmology of each test simulation. For the blinded data vector, we found $\chi^2=1748$ ($N_d=1440$, 78th percentile) for the joint analysis, 190/160 (67th) for $\xi_{\pm}$ and 1548/1280 (80th) for $\zeta_{\pm}$, which is slightly worse than for the unblinded measurement (Table 2), similar to the MAP percentiles of Appendix § 9.§ 4.

The signal-to-noise ratios of the individual parts of the data vector are 19.2 for $\xi_{+}$, 29.7 for $\xi_{-}$, 7.7 for $\zeta_{+}$ and 11.8 for $\zeta_{-}$. $\xi_{-}$ has the higher signal-to-noise ratio of the two, because at a given separation it probes smaller scales than $\xi_{+}$, whose covariance on 15′–250′ is dominated by cosmic variance.

Machine learning parameters

All density estimators are trained with the Adam optimiser [118] and a batch size of 250. The learning rate is multiplied by a tuned decay factor whenever the validation loss has not improved for 5 epochs, and training is stopped once the validation loss has not improved for 20 epochs, after which the model of the best epoch is kept. For every analysis variation, the hyperparameters of the NSF are tuned with optuna [114] over the search space listed in Table 11 (the number of transformations and of spline bins sampled uniformly over the integer ranges, the learning rate log-uniformly, the remaining parameters from the listed discrete choices), using a TPESampler [115] with 100 random start-up trials and 300 trials in total, and selecting the model with the highest likelihood on the validation set. The resulting hyperparameters of the three main analyses and the five data-vector variations of § 4 are listed in Table 12.

Hyperparameter	Search range
Transforms	$4-12$
Spline bins	$5-50$
Blocks	$\{2,3,4\}$
Features	$\{32,64,128,256,512\}$
Learning rate	$[10^{-5},10^{-3}]$, log-uniform
Decay factor	$\{0.1,\sqrt{0.1}\}$

Analysis	Transforms	Spline bins	Blocks	Features	Learning rate	Decay factor
$\xi_{\pm}$	4	5	2	32	$9.2 \times10^{-4}$	0.1
$\zeta_{\pm}$	4	7	3	32	$9.0 \times10^{-4}$	0.1
$\xi_{\pm} + \zeta_{\pm}$	4	5	2	32	$9.7 \times10^{-4}$	$\sqrt{0.1}$
no small scales	4	6	2	32	$9.9 \times10^{-4}$	0.1
no BIN1	4	5	2	32	$9.4 \times10^{-4}$	$\sqrt{0.1}$
no BIN2	6	5	2	32	$8.0 \times10^{-4}$	0.1
no BIN3	4	8	4	64	$4.1 \times10^{-4}$	0.1
no BIN4	4	5	2	32	$8.8 \times10^{-4}$	$\sqrt{0.1}$