UNIONS-3500 Weak Lensing: III. 2D Cosmological Constraints in Configuration Space

Introduction

Since its first detection two and a half decades ago (Bacon et al. 2000; Kaiser et al. 2000; Van Waerbeke et al. 2000; Wittman et al. 2000), weak gravitational lensing has emerged as a powerful tool to probe the Universe, particularly the properties of its most elusive components: dark energy and dark matter. Weak lensing exploits the fact that along the line of sight, the observed shapes of background galaxies are distorted—or ‘sheared’—as their light is deflected by the intervening gravitational potential of foreground matter distributions. By measuring the shapes of these sheared galaxies, or more specifically, the correlations between them, we are able to trace the underlying dark matter distribution, more commonly known as the large-scale structure (LSS) of the Universe.

An additional advantage of weak lensing as a cosmological probe is that it offers insight into the evolution of the LSS at low redshifts and on relatively small scales, thus making it highly complementary to other cosmological probes such as the cosmic microwave background (CMB) and baryonic acoustic oscillations (BAO). Weak lensing best constrains the combination of $\sigma_8$, the amplitude of clustering at a scale of $8\,h^{-1}\,\mathrm{Mpc}$, and the present-day total matter density $\Omega_\mathrm{m}$, often expressed as a single combined parameter $S_8 \equiv\sigma_8\sqrt{\Omega_\mathrm{m}/0.3}$.

There have been several galaxy surveys dedicated to weak lensing science, with Stage III surveys (see Albrecht et al. 2006, for a definition of the different weak-lensing survey ‘stages’) such as the Kilo-Degree Survey (KiDS; Wright et al. 2024) and the Dark Energy Survey (DES; Bechtol et al. 2025) publishing their last data releases, as well as the Hyper-Suprime Cam survey releasing their latest data set comprising three years’ worth of observations (HSC Y3; Li et al. [78]). In recent years, low-redshift weak lensing and galaxy clustering analyses have consistently favoured a lower value of $S_{8}$ than that inferred from high-redshift CMB measurements (e.g. Planck Collaboration [91]), with discrepancies at the level of $1$–$3 \sigma$ (see for example Di Valentino et al. 2025, for an in-depth review). Whether this tension signals new physics or simply reflects residual systematics, the current cosmological landscape appears poised for an answer. Most recently, the KiDS-Legacy cosmic shear analysis (Wright et al. 2025a), which incorporates improved modelling of systematic effects, together with a re-analysis of HSC Y3 data using clustering-redshift techniques (Choppin de Janvry et al. 2025), seem to suggest a possible resolution.

Looking ahead, Stage IV surveys, including Euclid (Euclid Collaboration: Mellier et al. [31]) and the Vera C. Rubin Observatory’s Legacy Survey of Space and Time (LSST; Ivezić et al. 2019) have begun operations and are expected to publish their first data releases in the coming years. With this new generation of experiments, sub-percent precision on key cosmological measurements becomes achievable, making the coming decade an exciting one for both the quantity and quality of observational data.

Into this landscape comes the Ultraviolet Near Infrared Optical Northern Survey (UNIONS; Gwyn et al. 2025), an ongoing photometric Stage III survey targeting up to 6250 deg$^{2}$ of the northern sky. In particular, UNIONS $r$-band imaging with the Canada-France Hawai’i Telescope (CFHT) enables high-quality weak lensing measurements through its excellent seeing.

This paper, along with four other accompanying papers, constitutes the first major UNIONS cosmic shear analysis, based on 3500 deg$^{2}$ of data accumulated over more than five years of coordinated observations from three wide-field telescopes in Hawai’i. We present configuration space cosmological constraints from the largest northern hemisphere weak lensing data set to date, which we name UNIONS-3500, hence offering a fitting complement to preceding Stage III analyses ([1, 2, 79] (Dalal et al. 2023; Anbajagane et al. 2025; Wright et al. 2025a)), which have been largely derived from observations in the southern sky.

This paper is structured as follows: Sect. 2 provides a brief overview of the UNIONS-3500 galaxy catalogue. In Sect. 3, we outline the theoretical modelling of the cosmic shear two-point correlation function (2PCF), the central probe used in this work, along with the additional components required for a complete cosmic shear analysis, including systematic error modelling, the covariance matrix and redshift distribution estimation, and a brief outline of external data sets included with the analysis. In Sect. 4, we describe our inference pipeline and detail how we take into account systematic effects and define our scale cuts. In Sect. 5 we present our results in terms of the marginalised cosmological parameter constraints, across various analysis setups and data set combinations. Finally, we offer our conclusions in Sect. 6.

This paper forms part of a larger release: Table 1 provides a summary of the accompanying papers, which detail catalogue construction (Paper I; Hervas Peters et al. 2026); B-mode validation (Paper II; Daley et al. 2026); cosmological constraints in configuration space (Paper III; this work); cosmological constraints in harmonic space (Paper IV; Guerrini et al. 2026), and lastly image simulations and shear calibration (Paper V; Hervas Peters et al. 2026b).

Paper Index	Author	Title
I	Hervas Peters et al. (2026a)	Weak lensing catalogues
II	Daley et al. (2026)	$B$-mode validation and comparison
III	This work	Cosmological constraints in configuration space
IV	Guerrini et al. (2026)	Cosmological constraints in harmonic space
V	Hervas Peters et al. (2026b)	Simulations and validation

The UNIONS data set

UNIONS combines multi-band photometric images from multiple telescopes. The Canada-France Imaging Survey (CFIS) provides $u$- and $r$-band images from CFHT. The shape measurement of galaxies relies on high-quality images taken in the $r$ band, which benefit from exquisite seeing of $\sim$0.7 arcsec, making it ideal for weak lensing science. The Panoramic Survey Telescope and Rapid Response System (Pan-STARRS) provides imaging in the $i$- and $z$-bands. The Subaru telescope also takes images in the $z$-band within the framework of the WISHES (Wide Imaging with Subaru HSC of the Euclid Sky) programme, and in the $g$-band through the WHIGS (Waterloo-Hawai’i IfA $g$-band Survey) programme. UNIONS provides wide-field multiband photometry over the optical bands, contributing to Euclid’s photometric redshifts in the northern sky. An extension to UNIONS is currently collecting data in all five bands down to 15° in declination. For a full review of the survey and its strategies, see Gwyn et al. (2025).

UNIONS-3500 weak lensing catalogue

This paper presents the first cosmological parameter estimation using the UNIONS-3500 2D weak lensing catalogue. Since this galaxy sample comprises images collected specifically in the $r$ band, we do not yet possess colour information to construct multiple tomographic bins; this analysis therefore relies on a single, two-dimensional redshift distribution (see Sect. 3.3 for more details).

The shape measurement was performed with ShapePipe (Farrens et al. 2022; Guinot et al. 2022), which also incorporates the galaxy point-spread function (PSF) modelling based on PSFex (Bertin 2011). Our lensing sample gathers data collected until the end of 2022, and is composed of over 611 million galaxies totalling an area of 3500 deg$^{2}$, corresponding to 2894 deg$^{2}$ of effective area after masking. The shapes of the galaxies were measured using ngmix (Sheldon 2015) and the calibration of the galaxy ellipticities was performed using Metacalibration (Huff & Mandelbaum 2017; Sheldon & Huff 2017). This gives an effective galaxy number density of 4.96 arcmin$^{-2}$ and per-component shape noise of $\sigma_{\mathrm{e}}=0.27$. Calibrated shapes using Metacalibration are saved but undergo a second correction step to remove residual per-component PSF leakage, an empirical process adapted from Li et al. [80]. The catalogue contains the spin-2 ellipticities for each galaxy, $e_{1}$ and $e_{2}$, corrected for the per-component additive biases $c_{1}$ and $c_{2}$. We also include their non-leakage corrected counterparts. Specific details on the PSF fitting and validation, shape measurement and leakage correction methods can be found in Paper I.

Additionally, different size cuts and masking schemes were explored before converging on the one employed to create the fiducial UNIONS-3500 shear catalogue. In its current version, a size cut of $r_{\mathrm{ngal}}/r_{\mathrm{PSF}}>0.7$ was applied, where $r_{\mathrm{h}}$ is defined as the half-light radius. The effective mask is an amalgamation of flags generated during MegaPipe image processing (bad pixels, chip defects), at ShapePipe processing time (foreground objects, blended galaxy images), and in post-processing (cosmic rays, exposure coverage); see Paper I for more details. In total, this nominally excludes objects detected as blends, regions around Messier and NGC objects, bright stars, stellar halos, diffraction spikes, defects, and cosmic rays. We also explored a more comprehensive masking of halo-like emissions around faint and bright stars, but $B$-mode null tests showed this additional masking introduced spurious systematic contamination; ultimately, it was not employed. For an in-depth study of the impact of size cuts and masking schemes, we refer the reader to Paper I and Paper II.

Modelling

In this section, we detail the modelling of the cosmic shear 2-point correlation function (2PCF), including validation tests for systematic effects, the covariance matrix, calibration of the shear magnification bias, and estimation of the source redshift distribution.

Cosmic shear 2-point correlation function

In the regime of weak lensing, the magnitude of distortion of the galaxy images is too small to be detected individually. Hence, an ensemble of galaxies is typically analysed together by computing the correlation of the measured shapes between pairs of galaxies. There exist several summary statistics that capture this information, most notably the two-point correlation functions in configuration space ($\theta$), and the angular power spectra $C_\ell$ in Fourier or harmonic space. In principle, these two statistics capture similar information on the sky, albeit with different sensitivities to scales and masking effects. For comprehensive reviews on weak lensing theory and methodologies, see for example Bartelmann & Schneider (2001) and Kilbinger (2015). For results on the cosmic shear analysis of UNIONS-3500 in Fourier space, we refer the reader to Paper IV.

The cosmic shear 2PCF is the expectation value of the product of the tangential and cross components of the galaxy shear signal, $\gamma_t$ and $\gamma_\times$ respectively Kaiser (1992),

which is a function of $\theta$, the angle of separation between a pair of galaxies on the sky. The $\xi_{\pm}(\theta)$ functions can also be derived from their Fourier space counterpart by a Hankel transform, such that, assuming the flat-sky approximations, it is given by Kaiser (1998)

where $J_n$ is the $n$th-order Bessel function of the first kind, with $n=0$ for $\xi_+$ and $n=4$ for $\xi_-$. Note that we have dropped the notation for tomographic indices $i$ and $j$ since we are only considering a 2D analysis. The cosmic shear power spectrum $C_\ell$ is given by the integral of the product of the lensing efficiency as a function of the comoving distance $q(\chi)$, and the nonlinear matter power spectrum $P_{\mathrm{NL}}(k,z(\chi))$,

where under the Limber approximation $k=\frac{\ell+1/2}{\chi}$, $\chi_{\mathrm{H}}$ is the comoving distance to the horizon, and

Here $H_0=100h\ \mathrm{km\,s^{-1}Mpc^{-1}}$ is the Hubble constant and $n(\chi)$ is the galaxy redshift distribution.

The cosmic shear observable is the observed ellipticity of a galaxy $e^{\mathrm{obs}}$, which is the sum of its intrinsic ellipticity $e^{\mathrm{s}}$ and shear $\gamma$ (boldface denotes spin-2 quantities):

Following the notation of Schneider et al. (2002), we define the angular bin width as $\Delta\theta$ and $\Delta\theta(\phi)=1$ for $\theta-\Delta\theta/2\leq\phi\leq\theta+\Delta\theta/2$ and zero elsewhere. The 2PCF estimator for a galaxy pair $i,j$ at separation $|\theta_i-\theta_j|$ is

where $w$ is the galaxy weight, computed during the shape measurement step, and $N_{\mathrm{p}}^{ij}(\theta)$ is the effective number of galaxy pairs in the galaxy bin.

In Figure 1 we present the $\hat{\xi}_{\pm}(\theta)$ data vectors computed by TreeCorr Jarvis et al. (2004), where we have binned the separation angle $\theta$ into 20 logarithmically-spaced bins over 1–250 arcmin. We also include their associated uncertainties given by the diagonals of the covariance matrix, whose calculation is detailed in Sect. 3.5.

Vector (PDF)

Systematic effects and validation tests

We assess the robustness of our data vectors to systematic measurement effects, focusing on PSF modelling, configuration-space $E/B$-mode decomposition, and the Complete Orthogonal Sets of $E/B$-mode Integrals (COSEBIs; Schneider et al. 2010; Asgari et al. 2012). These tests also drove the evolution of the catalogue version that led to the fiducial sample used here.

PSF systematic effects

The shear signal of a galaxy is estimated from its observed ellipticity. However, this measurement is often noisy, since it is contaminated by the PSF of the instrument. PSF systematics and mismodelling can contribute to the observed galaxy ellipticity in the form of an additive bias factor, $e^{\mathrm{sys}}$, such that

Quantifying this bias is essential for recovering an unbiased shear signal. A commonly used model for PSF leakage is given by Paulin-Henriksson et al. 2008; Jarvis et al. 2016

where $\alpha_{\mathrm{PSF}}$, $\beta_{\mathrm{PSF}}$ and $\eta_{\mathrm{PSF}}$ are constant free parameters, $e^{\mathrm{p}}$ is the ellipticity of the PSF model, $\delta e^{\mathrm{p}}=e^*-e^{\mathrm{p}}$ is the PSF ellipticity residual and $\delta T^{\mathrm{p}}=e^*(T^*-T^{\mathrm{p}})/T^*$ is the PSF size residual. We use the superscript ‘p’ to refer to the properties of the PSF (size or ellipticity), and ‘*’ to denote those of point-like sources (also referred to as stars) which, by convention, have zero size before convolution. with the PSF. Residual terms can only be evaluated at the positions of the stars.

The PSF-PSF and galaxy-PSF 2PCFs can be used to derive the amplitude of the leakage bias in the measured galaxy-galaxy correlation function. This is done by estimating the $\rho$ statistics (introduced by Rowe (2010); Jarvis et al. (2016)), representing PSF-PSF correlations, and $\tau$ statistics (introduced by Hamana et al. (2020); Gatti et al. (2021); Giblin et al. (2021)), quantifying galaxy-PSF correlations, directly from the data. They are respectively given by (where we have dropped the ‘PSF’ subscript in $\alpha$, $\beta$ and $\eta$ for clarity of notation)

where we only consider the ‘+’ component of each field, since the ‘−’ component tends to be noise-dominated and does not significantly increase cosmological information.

The second equality in Eq. (10) shows how $\alpha$, $\beta$, and $\eta$ are estimated from the $\rho(\theta)$ and $\tau(\theta)$ statistics. Further details can be found in Guerrini et al. (2025). In Figure 2, we plot the $\tau_{0,+}(\theta)$ and $\tau_{2,+}(\theta)$ statistics estimated from the catalogue, along with their best-fit signal recovered from the inference.

Vector (PDF)

Subsequently, we can compute the 2PCF estimator of the PSF systematics as

such that the total estimated 2PCF signal is a sum of the true shear signal and a contribution from the PSF systematics,

Section 4 of Paper I details the measurement of the $\xi^{\mathrm{sys}}_{\pm}(\theta)$ signal, where a non-negligible leakage bias at large scales was detected for both the $\xi_{\pm}$ data vectors, despite a noticeable improvement after applying an empirical leakage correction. We therefore modelled PSF systematics by additionally sampling the $\alpha$ and $\beta$ parameters at the inference step (see Sect. 4.3 for details), fixing $\eta=0$ since Guerrini et al. (2025) demonstrated that it is strongly correlated with $\alpha$, and that fixing $\eta$ does not significantly modify the estimate of the leakage bias. We also used this systematic test to inform our scale cuts: requiring that PSF systematics contribute less than 10% of the total signal gave an upper scale cut of 83 arcmin for both $\xi_{\pm}$.

E/B mode validation tests

Gravitational lensing produces a shear field that is curl-free to leading order; higher-order effects and intrinsic alignments can produce $B$-modes but remain below Stage III sensitivity (Schneider et al. 2022). Significant $B$ modes are therefore evidence of residual systematics, and we use them as one of the determining factors informing our scale cuts. For a detailed comparison of configuration- and harmonic-space $B$ modes across catalogue variants (changes in PSF size cuts, masking, etc.), see Paper II; here we summarise the results for the fiducial catalogue used for the cosmic shear analysis.

We computed two $B$-mode statistics: pure-mode correlation functions $\xi^{E/B}_{\pm}(\theta)$ in configuration space (Schneider et al. 2022); and the first six COSEBIs modes $B_n$. For each statistic, we calculate a probability to exceed (PTE) and require PTE $>0.05$ to pass the null test. Preliminary blinded analyses with the initial version of the catalogue passed configuration-space $B$-mode null tests but failed in harmonic space. We therefore introduced a more conservative size cut, removing approximately 25 per cent of galaxies and reducing $n_{\mathrm{eff}}$ from 6.48 to 4.96 arcmin$^{-2}$ and the per-component shape noise from 0.28 to 0.27. The stricter cut improved PSF systematics statistics, and the harmonic-space $B$-mode null tests passed.

Over the full angular range, COSEBIs show significant $B$-modes with a PTE of $1.37 \times10^{-5}$, while $\xi^{B}_{\pm}$ passes with a PTE of 0.40. After restricting to scales of 12–83 arcmin for both $\xi_{\pm}$ signals, both null tests pass with PTEs of 0.78 and 0.18 respectively. These results then informed our fiducial scale cuts (see Sect. 4.6).

Redshift distribution

We estimate the redshift distribution of our weak lensing source sample based on the colour-redshift relation. Our shear catalogue used for this analysis, however, is not yet fully covered by UNIONS multi-band photometry. In order to still employ the colour-based approach for the redshift calibration, we exploit the spatial overlap between the UNIONS $r$-band data and the Canada–France–Hawai‘i Telescope Lensing Survey (CFHTLenS; Heymans et al. 2012; Erben et al. 2013) W3 field, which covers 44.2 deg$^2$. CFHTLenS provides significantly deeper $ugriz$ photometry (Hildebrandt et al. 2012) than UNIONS, such that essentially all detected sources have CFHTLenS counterparts. We therefore cross-match these sources to CFHTLenS and adopt the associated $ugriz$ magnitudes. The matched sample is assumed to trace the same underlying colour–redshift distribution as the full UNIONS population.

To calibrate the redshift distribution, we assemble a spectroscopic sample ($\sim65\,500$ galaxies) that occupies the same photometric space as the matched UNIONS-CFHTLenS sources. This sample combines data from three deep spectroscopic surveys: the DEEP2 Galaxy Redshift Survey (Newman et al. 2013), the VIMOS VLT Deep Survey (VVDS; Le Fèvre et al. 2005), and the Vimos Public Redshift Survey (VIPERS; Scodeggio et al. 2018). All of these surveys were observed with the CFHTLenS $ugriz$ filters, enabling direct comparisons in colour–magnitude space. After applying standard survey-specific quality cuts, we obtain a clean and representative set of galaxies for the calibration of the redshift distribution.

Using the multi-band photometry of the spectroscopic sample, we train a self-organising map (SOM; Kohonen 1982), which arranges galaxies based on their positions in the multi-dimensional magnitude space (Masters et al. 2015; Wright et al. 2020). The initial SOM is defined on a 101 $\times$ 101 cell grid, providing a fine-grained tiling of colour–magnitude space. For robust statistical sampling, we subsequently hierarchically cluster the SOM into approximately 5000 effective resolution elements. This preserves the structure of the photometric manifold while ensuring that each region contains a sufficient number of spectroscopic objects for more robust mean statistics.

We then populate the SOM with the UNIONS sources by assigning each galaxy to its best-matching SOM cell based on its $ugriz$ photometry. Each UNIONS galaxy carries two lensing-related weights: (i) the galaxy weight $w$ (as explained in Paper I); and (ii) an additional weighting based on the average shear response $\langle R_{\gamma}\rangle$.

For the average shear response, we use the values that result from binning the whole UNIONS-3500 sample in a two-dimensional space defined by the signal-to-noise ratio (SNR) and the size ratio parameter $r_{\mathrm{h}}/r_{\mathrm{PSF}}$ (shown in Figure 3 of Paper I). Each galaxy $j$ is assigned the average shear response value of its bin $\langle R_{\gamma}\rangle_{\mathrm{bin}(j)}$, which we use as an additional multiplicative weight.

To account for potential selection biases in the spectroscopic sample, we adopt the “prior volume weighting” scheme described in Wright et al. (2025b). This correction addresses the fact that the probability distribution of redshift at a fixed colour may differ between the spectroscopic and wide-field samples due to differing selection functions or spectroscopic success rates (e.g., colour-redshift degeneracies, Hartley et al. 2020). Following the methodology of Wright et al. (2025b), we define a corrective weight for each spectroscopic galaxy $k$ as the ratio of the redshift probability density functions (PDFs) of a simulated wide-field sample ($P_{\mathrm{w}}$) and the calibration sample ($P_{\mathrm{c}}$):

The wide-field PDF, $P_{\mathrm{w}}(z)$, is derived from an analytic model for a magnitude-limited sample, $N(z,m)$, which was calibrated using SURFS-Shark light cones [73] (Elahi et al. 2018). For our UNIONS sample, we approximate this selection using an $r$-band magnitude-limited model ($20 \le r \le23.75$) to mimic the effective depth and the impact of the lensing shape weights. This prior volume weight ensures that even before the SOM cell assignment, the spectroscopic galaxies are adjusted to better represent the expected underlying redshift baseline of the UNIONS survey.

Combining these weighting schemes, the SOM weight for cell $i$ is defined as

where the numerator sums the weighted counts of UNIONS galaxies $j$ assigned to cell $i$, and the denominator is the sum of the prior volume weights $w_k^{\mathrm{prior}}$ of spectroscopic calibration galaxies $k$ in that cell (Wright et al. 2025b). This definition ensures that the spectroscopic sample is reweighted to match the effective distribution of sources contributing to the shear signal.

Finally, the redshift distribution of the UNIONS weak-lensing sample is obtained by applying the SOM weights to the spectroscopic redshift distributions in each cell:

where $n_i^{\mathrm{spec}}(z)$ is the redshift distribution of spectroscopic galaxies located in SOM cell $i$. To obtain a statistically robust estimate of the final redshift distribution, we conduct a bootstrap resampling over the spatially binned spectroscopic calibration sample. We generate 1000 realisations of $n(z)$ and adopt the mean of these realisations as our final estimate. The final calibrated $n(z)$ distribution is shown in Figure 3.

Shear multiplicative bias

Errors in the shear measurement step potentially introduce systematic biases, such that the true ellipticity of a galaxy is given by

where $\mathbf{c}$ is the additive component, corrected for during catalogue processing, and $m$ is the multiplicative bias, for which we detail its method of estimation in this section.

In Stage III ground-based surveys, blended objects, which are galaxies that receive flux from a close angular neighbour and cannot be spatially resolved, have been determined to be the leading contributor of multiplicative bias ([80] (MacCrann et al. 2021)). In Paper V, the image simulations and shear calibration efforts are detailed in full. In short, two suites of images were run: one with

galaxies placed on a grid and one with realistic galaxy placements obtained from an $N$-body simulation. Each is run with four shears ($\gamma_1 = \pm0.025$, $\gamma_2 = \pm0.025$) to test how well the input shear is recovered.

The image simulations revealed three important sources of multiplicative bias that proved challenging to disentangle. Firstly, a bug in the shear measurement pipeline led to the WCS offsets between different epochs not being correctly propagated, giving rise to a $\Delta m$ contribution of approximately $-0.02$. Secondly, a selection effect was introduced when SExtractor Flags = 0 were selected, since objects with Flags 1 and 2 are, on average, more elliptical. While most selection effects are calibrated through Metacalibration, this one was not captured, as the detection was not repeated on the differently sheared branches. This also led to a multiplicative bias $\Delta m \approx-0.02$. Finally, physically blended objects estimated from simulated images with galaxies placed at realistic positions based on an $N$-body simulation gave a multiplicative bias $\Delta m \approx-0.017$.

Jointly, all effects were cumulatively responsible for a multiplicative bias $m = -0.057 \pm0.0047$. We conservatively multiply the uncertainty by a factor of 3 to account for the fact that these strong effects might not be exactly captured by the simulated images. This leads to a total multiplicative bias $m = -0.057 \pm0.014$.

Covariance modelling

We calculate the covariance matrix between angular bins by taking into account the Gaussian (G), connected non-Gaussian (nG), and super sample covariance (SSC) contributions. In a tomographic analysis, the total covariance matrix of two angular power spectra $C_{\ell_1}^{ij}$ and $C_{\ell_2}^{kl}$ in Fourier space is given by

Following the notation of Hu & Jain (2004), the Gaussian covariance term of any galaxy $g$ or shear $\gamma$ field can be expressed as

where $A,B,C,D$ represents either $g$ or $\gamma$ for the observed fraction of the sky $f_{\mathrm{sky}}$, and noise power spectra $N_{ij}$. In the case of cosmic shear, $A=B=C=D=\gamma$, hence it is simply $N_{\ell}^{\gamma\gamma,ij}=\frac{\sigma_{\mathrm{e}}^2}{2\bar{n}_i}\delta_{ij}^{K}$. Here $\bar{n}_i$ is the mean galaxy number density of bin $i$ and $\sigma_{\mathrm{e}}$ the shape noise, while $\delta_{ij}^{K}$ denotes the Kronecker delta. However, transitioning to small scales, the cosmic shear field becomes significantly non-Gaussian, leading to extra contributions arising from the connected four-point function of these fields. For the full expression of the connected non-Gaussian covariance terms, see for example Cooray & Sheth (2002) and Takada & Jain (2009).

Lastly, the SSC term, which captures the uncertainty due to a change in the background density at modes larger than the survey area, can be modelled following Takada & Hu (2013):

where $\sigma_b^2$ is the covariance of the background density field $\delta_b$ within the survey window. In the case of a cosmic shear-only analysis, $q_i^{A,B,C,D}(\chi)$ is the cosmic shear lensing efficiency given in Eq. (4), and $\sigma_b$ is given by Lacasa & Rosenfeld (2016)

where $j_0$ is the spherical Bessel function of zeroth order, and $P_{\mathrm{L}}(k,z)$ refers to the linear matter power spectrum.

The covariance matrix for the 2PCF in configuration space follows from a Hankel transformation,

where once again $n=0$ for $\xi_+$ and $n=4$ for $\xi_-$. We calculate our covariance matrix based on the above equations, using the CosmoCov software (Krause & Eifler 2017; Fang et al. 2020), where we have adopted the fiducial input cosmology following table 1 of Reischke et al. (2025). We plot the full correlation matrix (where $\operatorname{Corr}_{ij}=\operatorname{Cov}_{ij}/\sqrt{\operatorname{Cov}_{ii}\times\operatorname{Cov}_{jj}}$) in Figure 4.

Effect of masking on covariance estimation

As was investigated in Troxel et al. (2018) and Friedrich et al. (2021), masking can result in a noncontiguous and patchy survey footprint, thus inducing additional uncertainties in the covariance matrix of the correlation functions. Following Schneider et al. (2002), the pure shape noise term of the covariance matrix in configuration space (referred to as $N_\gamma^{ij}$ in the previous section) is given by

where $N_{\mathrm{p}}^{ij}(\theta)$ is as defined in Eq. (6). When boundary and masking effects are neglected, it can be approximated by

where $A$ is the survey area. On the other hand, a more accurate estimation of $N_{\mathrm{p}}^{ij}$ can be obtained by taking into account the survey mask, such that

where the mask power spectrum $w_{\mathrm{mask}}(\theta)$, normalised to the survey area, modifies $N_{\mathrm{p}}^{ij}$.

We examine the effect of the UNIONS mask on the covariance matrix, and present the full ratio of the masked to unmasked covariance matrix, including the cross-covariances and non-Gaussian and SSC contributions, in Figure 5. We find that overall, including masking effects increases the diagonal terms of the auto-covariance (i.e. the noise terms as defined in Eq. (22)), with a peak in the ratio occurring at approximately 1 arcmin for $\xi_+$ and 20 arcmin for $\xi_-$. This ratio decreases with increasing $\theta$, as the angular separation of the measurement becomes larger than the masking scale. The overall impact of masking is larger for the uncertainties in $\xi_-$ than $\xi_+$. Our results are consistent with the trend and values found in Troxel et al. (2018), as shown in the bottom panel of their Fig. 1. On the other hand, the ratios of the off-diagonal and cross-covariance terms are reduced (except at large scales). This is because the off-diagonal terms include an SSC contribution, which is suppressed under masking. This is apparent in the integrand for $\sigma_b^2$, where the survey window is rescaled by the mask power spectrum, thereby downweighting the contribution of modes suppressed by the survey mask to the overall background variance.

Vector (PDF)

We adopt the masked covariance matrix as our fiducial choice, and also present results with the unmasked covariance to quantify its impact on our analysis.

External data sets

Here we summarise the external data sets incorporated into our inference runs. These include highly complementary measurements from the CMB, which provide information at high redshifts, as well as BAO data, which offer information on the expansion history at low redshifts. Combined with cosmic shear, they give a more complete description of the Universe at both the background and perturbation levels, tightening overall cosmological constraints.

Cosmic Microwave Background

The CMB, relic radiation from free-streaming photons that decoupled from baryons at redshifts $z\simeq1100$, provides highly precise constraints on early-Universe physics. We include data from the 2018 data release of the Planck mission (Planck Collaboration [90]), specifically those of the $TT$, $TE$ and $EE$ auto- and cross-spectra, spanning a multipole range of $2\leq\ell\leq2508$ for $TT$, and $2\leq\ell\leq1996$ for $TE$ and $EE$. To avoid having to sample over the large number of CMB nuisance parameters, we employ instead the plik_lite likelihood (Prince & Dunkley 2019), which assumes a Gaussian likelihood and only marginalises over the nuisance parameter $A_{\mathrm{lens}}$, the correction amplitude due to lensing effects. Furthermore, we do not include the lensing power spectrum likelihood, as we currently lack a prescription for the cross-correlation in the covariance matrix between cosmic shear and lensing data.

Baryonic Acoustic Oscillations

We include observations of the BAO scale from the second data release of the Dark Energy Spectroscopic Instrument (DESI DR2; Abdul Karim et al. 2025a,b), in the form of measurements of the transverse and comoving BAO distance ($D_{\mathrm{M}}$ and $D_{\mathrm{H}}$ respectively). They span a redshift range of $0.1 < z < 3$, and represent the most comprehensive and precise measurement of the BAO scale to date. These results attracted significant attention, as they favour a cosmological model with a varying dark energy equation of state over a cosmological constant $\Lambda$. However, we shall only consider a $\Lambda$CDM model in all our analyses, for which BAO can lend competitive constraining power on $\Omega_{\mathrm{m}}$. We employ the full range of BAO measurements taken from the Bright Galaxy Survey (BGS), Luminous Red Galaxies (LRG), Emission Line Galaxies (ELG), quasars (QSO) and Lyman-$\alpha$ forest (Ly$\alpha$) samples.

Inference pipeline setup

We use Bayesian inference to derive constraints on the cosmological parameters of the concordance $\Lambda$CDM model from the data vectors, covariance matrix, and redshift distribution described above.

Bayesian inference relies on Bayes’ Theorem, which gives the probability distribution of the parameters $\theta$ given a model $M$ and the observed data $d$. This probability distribution $P(\theta|d,M)$, known as the posterior distribution, is defined as

where $\Pi(\theta|M)$ is the prior distribution, quantifying our initial knowledge of the distribution of $\theta$, $Z(d|M)$ is the evidence, which gives the probability of observing the data given $M$, and $\mathcal{L}(d|\theta,M)$ is known as the likelihood function, which is the probability of observing the data $d$ given the model $M$ with parameter values $\theta$.

We adopt a Gaussian likelihood:

where $d$ is the data vector, $T(\theta)$ is the theory vector derived from the model, and $C$ is the covariance matrix of the data.

Likelihood inference

We use the Einstein–Boltzmann solver CAMB (Lewis et al. 2000) to calculate cosmological quantities and the linear matter power spectrum, and the CosmoSIS (Zuntz et al. 2015) pipeline for parameter inference. Our baseline sampler is Polychord (Handley et al. 2015a,b), following many other Stage III survey analyses ([1, 79] (Dalal et al. 2023; Anbajagane et al. 2025)). Polychord produces reliable posterior estimates at reasonable computational cost (Lemos et al. 2023), though it becomes expensive as the parameter space grows. We therefore also run Nautilus (Lange 2023), which leverages deep learning to improve sampling efficiency (Euclid Collaboration: Cañas-Herrera et al. [30]; Wright et al. (2025a)), as a consistency check.

We adopt prior ranges for the cosmological parameters as presented in Table 2 of Abbott et al. [2], except that we additionally impose a prior on the baryonic energy density parameter $\omega_{\mathrm{b}} \equiv\Omega_{\mathrm{b}}h^2$, which was derived from Big Bang nucleosynthesis (BBN) constraints. Computed from the PRyMordial code (Burns et al. 2024; Schönberger et al. 2024), it gives $\omega_{\mathrm{b}} = 0.02218 \pm0.00055$ assuming a $\Lambda$CDM model. In total, the cosmological parameters being sampled are $\{\omega_{\mathrm{c}},\omega_{\mathrm{b}},H_0,n_{\mathrm{s}},S_8\}$: the energy density of cold dark matter; the energy density of baryons; the present-day value of the Hubble parameter; the spectral index; and the amplitude of clustering $S_8 \equiv\sigma_8\sqrt{\Omega_{\mathrm{m}}/0.3}$.² Since the optical depth of reionisation $\tau_{\mathrm{reio}}$ is not well constrained with cosmic shear data, we fix it at the Planck 2018 best-fit value of $\tau_{\mathrm{reio}} = 0.0544$ (Planck Collaboration [91]). However, when including CMB data, we do not impose the BBN prior, and additionally sample $\tau_{\mathrm{reio}}$. We also assume a flat Universe and one massive and two massless neutrino species, with a total mass of $\sum m_\nu= 0.06$ eV. The following subsections describe the modelling of intrinsic alignment, the nonlinear matter power spectrum, and PSF systematics. Table 2 lists the full set of priors for the cosmological and nuisance parameters that we sample in our analysis.

Parameter	Prior
Cosmology
$\omega_{\mathrm{b}}$ [no CMB]	$\mathcal{N}(0.02218, 0.00055)$
$\omega_{\mathrm{b}}$ [with CMB]	$\mathcal{U}(0.0190, 0.0260)$
$\omega_{\mathrm{c}}$	$\mathcal{U}(0.051, 0.255)$
$h$	$\mathcal{U}(0.64, 0.82)$
$n_{\mathrm{s}}$	$\mathcal{U}(0.84, 1.10)$
$S_8$	$\mathcal{U}(0.1, 1.3)$
$\log_{10}(T_{\mathrm{AGN}}/\mathrm{K})$	$\mathcal{U}(7.3, 8.0)$
$\tau_{\mathrm{reio}}$ [with CMB]	$\mathcal{U}(0.01, 0.80)$
Nuisance
$A_{\mathrm{IA}}$	$\mathcal{N}(0.83, 0.78)$
$\alpha_{\mathrm{PSF}}$ (leakage corrected)	$\mathcal{N}(0.005, 0.0022)$
$\beta_{\mathrm{PSF}}$ (leakage corrected)	$\mathcal{N}(0.810, 0.1148)$
$\alpha_{\mathrm{PSF}}$ (non-leakage corrected)	$\mathcal{N}(0.022, 0.0026)$
$\beta_{\mathrm{PSF}}$ (non-leakage corrected)	$\mathcal{N}(0.798, 0.1150)$
$\Delta z$	$\mathcal{N}(-0.003, 0.018)$
$m$	$\mathcal{N}(-0.057, 0.014)$

When reporting the best-fit parameter values for all parameters except $S_8$ and $\Omega_{\mathrm{m}}$, we quote their 1D marginalised mode, or marginalised maximum a posteriori (MAP), and their 68% confidence interval (CI), as in recent analyses such as Dalal et al. (2023); Li et al. [79]. This is obtained by evaluating the peak of the 1D posterior distribution after marginalising over all remaining parameters, which we estimate using the in-built Gaussian KDE smoothing kernel of getdist (Lewis 2025). Other commonly adopted statistics include the weighted mean of the marginalised posterior (e.g. DES Collaboration et al. [22]), while others report the MAP value with a 68% confidence region defined by the projected joint highest posterior density (PJ-HPD; Joachimi et al. (2021)), as in Wright et al. (2025a). Each approach has distinct advantages and limitations (see Abbott et al. [2], for a detailed discussion). In particular, the MAP and its associated PJ-HPD interval can be challenging to estimate robustly from nested sampling chains, due to the relatively sparse sampling of high-posterior regions (Muir et al. (2020)), and can be unstable in the presence of strong parameter degeneracies. Conversely, statistics based on marginalisation, such as the weighted mean and marginalised mode, are susceptible to volume projection effects, especially for non-Gaussian posteriors. We thus attempt to balance a trade-off by choosing instead to report the 2D marginalised mode for $S_8$ and $\Omega_{\mathrm{m}}$, the two parameters which are the most important and also most sensitive to projection effects in cosmic shear analyses (see Sect. A). By marginalising over both parameters simultaneously, we mitigate the worst of the impact due to projection effects when marginalising in one dimension. For completeness, we also quote the weighted mean of $S_8$ in Table B1.

Analysis choice	$S_8$ 2D MAP	$S_8$ Weighted mean	$\xi_+$	$\xi_-$	$\tau_{0,2}$	Total
Fiducial	$0.831^{+0.067}_{-0.078}$	$0.829$	$\chi^2/n=4.29$, PTE $0.745$	$\chi^2/n=5.27$, PTE $0.627$	$\chi^2/n=9.50/14$, PTE $0.798$	$\chi^2/n=65.77/40$, PTE $0.00629$; $\chi^2/n=75.33/54$, PTE $0.0291$
Including $\xi_+$ small scales	$0.886^{+0.049}_{-0.104}$	$0.860$	$\chi^2/n=16.46/10$, PTE $0.0872$	$\chi^2/n=4.60/7$, PTE $0.708$	$\chi^2/n=22.85/17$, PTE $0.154$	$\chi^2/n=65.77/40$, PTE $0.00628$; $\chi^2/n=86.83/57$, PTE $0.00663$
Flat $\alpha_{\mathrm{PSF}}$ and $\beta_{\mathrm{PSF}}$ priors	$0.825^{+0.071}_{-0.075}$	$0.830$	$\chi^2/n=4.66/7$, PTE $0.70$	$\chi^2/n=5.47/7$, PTE $0.603$	$\chi^2/n=9.99/14$, PTE $0.763$	$\chi^2/n=65.78/40$, PTE $0.00627$; $\chi^2/n=75.90/54$, PTE $0.0263$
No $\xi_{\mathrm{sys}}^{\pm}$	$0.824^{+0.063}_{-0.068}$	$0.826$	$\chi^2/n=4.62/7$, PTE $0.706$	$\chi^2/n=5.27/7$, PTE $0.626$	$\chi^2/n=9.87/14$, PTE $0.771$	—; $\chi^2/n=9.87/14$, PTE $0.771$
No object-wise leakage correction	$0.818\pm0.073$	$0.820$	$\chi^2/n=4.44/7$, PTE $0.728$	$\chi^2/n=5.39/7$, PTE $0.613$	$\chi^2/n=9.61/14$, PTE $0.790$	$\chi^2/n=67.41/40$, PTE $0.0043$; $\chi^2/n=77.24/54$, PTE $0.0207$
Flat $\Delta z$ prior	$0.808^{+0.121}_{-0.093}$	$0.844$	$\chi^2/n=9.81/7$, PTE $0.200$	$\chi^2/n=17.67/7$, PTE $0.0136$	$\chi^2/n=24.28/14$, PTE $0.0424$	$\chi^2/n=65.77/40$, PTE $0.00628$; $\chi^2/n=93.25/54$, PTE $0.000727$
$\Delta z$	$0.830^{+0.058}_{-0.071}$	$0.825$	$\chi^2/n=5.17/7$, PTE $0.639$	$\chi^2/n=6.35/7$, PTE $0.499$	$\chi^2/n=11.17/14$, PTE $0.672$	$\chi^2/n=65.76/40$, PTE $0.00629$; $\chi^2/n=77.28/54$, PTE $0.0205$
Flat $A_{\mathrm{IA}}$ prior	$0.895^{+0.172}_{-0.134}$	$0.918$	$\chi^2/n=9.13/7$, PTE $0.244$	$\chi^2/n=14.67/7$, PTE $0.0404$	$\chi^2/n=21.31/14$, PTE $0.094$	$\chi^2/n=65.77/40$, PTE $0.00628$; $\chi^2/n=89.57/54$, PTE $0.00168$
$A_{\mathrm{IA}}$	$0.770^{+0.029}_{-0.059}$	$0.752$	$\chi^2/n=4.26/7$, PTE $0.750$	$\chi^2/n=5.27/7$, PTE $0.627$	$\chi^2/n=9.45/14$, PTE $0.801$	$\chi^2/n=65.77/40$, PTE $0.00629$; $\chi^2/n=75.30/54$, PTE $0.0293$
$m$ bias	$0.809^{+0.050}_{-0.079}$	$0.796$	$\chi^2/n=4.69/7$, PTE $0.697$	$\chi^2/n=5.04/7$, PTE $0.656$	$\chi^2/n=9.62/14$, PTE $0.789$	$\chi^2/n=65.76/40$, PTE $0.00629$; $\chi^2/n=75.49/54$, PTE $0.0283$
Unmasked covmat	$0.822^{+0.072}_{-0.070}$	$0.823$	$\chi^2/n=4.78/7$, PTE $0.686$	$\chi^2/n=6.93/7$, PTE $0.436$	$\chi^2/n=11.39/14$, PTE $0.655$	$\chi^2/n=65.76/40$, PTE $0.0063$; $\chi^2/n=77.48/54$, PTE $0.0198$
Halofit	$0.800^{+0.068}_{-0.074}$	$0.806$	$\chi^2/n=4.58/7$, PTE $0.711$	$\chi^2/n=4.83/7$, PTE $0.680$	$\chi^2/n=9.40/14$, PTE $0.804$	$\chi^2/n=65.77/40$, PTE $0.00628$; $\chi^2/n=75.19/54$, PTE $0.0299$
HMCode no baryons	$0.813^{+0.072}_{-0.068}$	$0.822$	$\chi^2/n=4.84/7$, PTE $0.680$	$\chi^2/n=5.78/7$, PTE $0.566$	$\chi^2/n=10.46/14$, PTE $0.728$	$\chi^2/n=65.77/40$, PTE $0.00628$; $\chi^2/n=76.38/54$, PTE $0.0242$
Nautilus sampler	$0.826^{+0.068}_{-0.069}$	$0.826$	$\chi^2/n=4.57/7$, PTE $0.713$	$\chi^2/n=5.83/7$, PTE $0.559$	$\chi^2/n=10.28/14$, PTE $0.742$	$\chi^2/n=65.76/40$, PTE $0.0063$; $\chi^2/n=76.16/54$, PTE $0.0251$
UNIONS-3500 + CMB	$0.855^{+0.020}_{-0.017}$	$0.857$	$\chi^2/n=4.33/7$, PTE $0.741$	$\chi^2/n=4.92/7$, PTE $0.67$	$\chi^2/n=9.15/14$, PTE $0.822$	$\chi^2/n=65.76/40$, PTE $0.0063$; $\chi^2/n=75.01/54$, PTE $0.031$
UNIONS-3500 + CMB + BAO	$0.849^{+0.015}_{-0.017}$	$0.849$	$\chi^2/n=4.40/7$, PTE $0.733$	$\chi^2/n=4.99/7$, PTE $0.661$	$\chi^2/n=9.22/14$, PTE $0.817$	$\chi^2/n=65.77/40$, PTE $0.00628$; $\chi^2/n=75.16/54$, PTE $0.030$

Intrinsic alignment

Multiple strategies to disentangle intrinsic alignment (IA) from cosmic shear require accurate estimates of galaxy redshifts (Joachimi & Schneider 2010). Recent cosmic shear measurements have been somewhat successful in directly estimating the intrinsic alignment parameters when jointly fitting the shear-shear ($\gamma\gamma$), shear-intrinsic ($\gamma I$) and intrinsic-intrinsic (II) correlations to the signal (Asgari et al. 2021; Secco et al. 2022). This success can be traced back to the very different ways these effects act with respect to galaxy separation: small angular bins are more affected by II contributions, while bins with larger separation have a strong $\gamma I$ contribution.

In Fourier space, the noiseless cosmic shear angular power spectrum decomposes as

Models exist which give an approximate theoretical prescription for the latter two terms, with one of the widely employed models being the redshift-dependent nonlinear linear alignment model (NLA; Hirata & Seljak 2004; Bridle & King 2007). Here, the 3D power spectra $P^{\gamma I}(k,z)$ and $P^{II}(k,z)$ are expressed as a rescaling of the nonlinear matter power spectrum

where $A_{\mathrm{IA}}$ is a dimensionless IA amplitude parameter to be marginalised over, $\rho_{\mathrm{crit}}$ is the critical energy density, $D(z)$ is the growth factor and $C_1=5\times10^{-14}h^{-2}M_{\odot}^{-1}\mathrm{Mpc}^{3}$.

In this analysis, the use of a single redshift bin prohibits us from efficiently constraining $A_{\mathrm{IA}}$, necessitating the imposition of a well-motivated prior. Since multiple direct intrinsic alignment measurement models have been proposed (Joachimi et al. 2011; Mandelbaum et al. 2011; Singh et al. 2015; Johnston et al. 2019; Fortuna et al. 2021b; Samuroff et al. 2023; Navarro-Girones et al. 2025), we choose to rely on these, further ignoring the impact of any higher-order terms, such as those proposed in the tidal alignment and tidal torquing model (TATT; Blazek et al. 2019). That strategy has also been adopted in previous works ([81] (Fortuna et al. 2021a; Wright et al. 2025a)).

We use a data-driven approach to estimate $A_{\mathrm{IA}}$, making use of the DES Y3 C3HLenS BPZ (Benítez 2000) catalogue from which the $n(z)$ was derived. We then adopt the commonly employed method of splitting intrinsic alignment contributions coming from red and blue galaxies (Krause et al. 2016), motivated by the distinction between pressure-supported and rotationally supported galaxies, respectively. The parametrisation is given by

where $f_{\mathrm{r}}$ indicates the fraction of red galaxies, and $f_{\mathrm{b}}$ indicates the fraction of blue galaxies within our sample. For simplicity, we do not introduce any luminosity or redshift dependence, hence $A_{\mathrm{IA}}$ is a constant parameter. To separate blue from red galaxies, we choose to assign all galaxies with $T_{\mathrm{B}}<1.9$ as red galaxies and treat the rest as blue. The determination of the $T_{\mathrm{B}}$ index is done through BPZ, where it utilises six model templates, arranged by increasing star formation activity, and performs linear interpolation between neighbouring templates to find the best-matching spectral energy distribution. The output parameter, $T_{\mathrm{B}}$, represents the selected best-fit template, or more precisely the blend of adjacent templates, in steps of 0.1. From there, we obtain estimates for $f_{\mathrm{r}}=0.245$ and $f_{\mathrm{b}}=0.745$.

For the values $A_{\mathrm{IA,r}}$ and $A_{\mathrm{IA,b}}$, we use prior values derived from direct measurements. These rely on accurate redshift estimation, either photometric or spectroscopic, to compute the projected shape-density correlation function $w_{\mathrm{g+}}$ with limited separation. For the blue galaxy sample, we adopt the result from Johnston et al. (2019), which fitted the available measurements to obtain $A_{\mathrm{IA,b}}=0.21\pm0.37$. While other works using this model have set the intrinsic alignment of blue galaxies to 0 (consistent with the data), we prefer a more conservative approach and adopt this observationally grounded prior with a non-negligible Gaussian width. For an estimate of $A_{\mathrm{IA,r}}$, we sample from the posterior distribution of the double power-law fit done on $A_{\mathrm{IA}}(L/L_0)$ at the mean luminosity value of red galaxies measured on W3, following Hervas Peters et al. (2025). The luminosities are obtained from LePhare (Arnouts & Ilbert 2011), where the mean value for UNIONS is found to be $(L/L_0)=-0.77$. Sampling the posterior gives $A_{\mathrm{IA,r}}=2.75\pm0.49$, a value in accordance with most KiDS estimates, but not with fiducial DES values (Abbott et al. [2]). We then add the Gaussian errors in quadrature to arrive at a prior of $A_{\mathrm{IA}}=0.83\pm0.39$. We have tested the robustness of this prior against changes in $T_{\mathrm{B}}$, finding only small differences.

This method is mainly limited by two factors. First, the galaxies for which reliable spectroscopic or photometric redshifts are available are not representative of either the blue or red population, since they mostly need to be bright to have reliable photometric redshifts or are targeted by spectroscopic surveys. Second, these various methods are all obtained with different lensing surveys. It has been shown that the shape algorithm used (Singh & Mandelbaum 2016) and the band the galaxy has been measured in Georgiou et al. (2019) impact the measured amplitude of the intrinsic alignment signal. These two reasons motivate the doubling of the width of the Gaussian prior to ±0.78 instead. For comparison, we also conduct a run assuming an uninformative flat prior on $A_{\mathrm{IA}}$, and additionally assume the absence of any IA effects (setting $A_{\mathrm{IA}} = 0$).

PSF systematics

To mitigate systematic effects from PSF mismodelling, we employ two safeguards: firstly, we apply scale cuts to ensure that the additive contribution from the leakage signal, $\xi^{\mathrm{sys}}_{\pm}$, remains below 10% of the total measured signal. Additionally, we have chosen to jointly fit the observed correlation function signal $\xi^{\mathrm{obs}}_{\pm}$ with the $\tau$ statistics as introduced in Sect. 3.2.1. Specifically, we sample $\alpha$ and $\beta$ as nuisance parameters and compute the theoretical $\tau$ statistics according to Eq. (10), given the $\rho(\theta)$’s estimated from the data. Since we set $\eta= 0$, this would require only the computation of $\tau_0(\theta)$ and $\tau_2(\theta)$, where we only consider the $\tau_+$ signal. The $\tau$ statistic likelihood is thus calculated as

where the subscripts ‘d’ and ‘t’ denote the data and theory vectors respectively. $\mathbf{C}_{\tau}$ is a semi-analytical covariance, presented in Figure 4 concatenated with the shear-shear covariance; see Guerrini et al. (2025) for its full derivation. Cross-correlations between the shear-shear 2PCF and the $\tau$ statistics are neglected for simplicity.

Additionally, we include the PSF systematic contribution as an additive bias to the cosmological signal, by calculating the theoretical $\xi^{\mathrm{sys}}_{\pm}$ following Eq. (11), given the sampled values of $\alpha$ and $\beta$. The total log-likelihood being minimised is thus

This joint fit accounts for the additive leakage bias and marginalises over PSF systematic uncertainty. Other Stage III surveys employ a similar approach (see for example Li et al. [79]; Zhang et al. (2023)).

For the fiducial analysis, we use Gaussian priors on $\alpha$ and $\beta$ rather than flat priors. For data vectors with significant PSF systematics, unconstrained $\alpha$ and $\beta$ might risk biased posteriors. Specifically, if these parameters are left unconstrained (i.e., with wide priors) during the main cosmological inference step, they may fit the cosmological signal $\xi^{\mathrm{obs}}_{\pm}$ rather than the systematic $\xi^{\mathrm{sys}}_{\pm}$ (Eq. (12)), since non-trivial degeneracies exist between leakage and cosmological parameters. Consequently, the resulting posteriors on $\alpha$ and $\beta$ would be inconsistent with those obtained instead from a direct $\rho/\tau$ fit.

To obtain the priors on $\alpha$ and $\beta$, we first run a separate ‘PSF systematics inference step’ whereby we solely fit $\alpha$ and $\beta$ based on the $\rho$ and $\tau$ statistics estimated from the data, per Eq. (10). This is the methodology described above. We use a vanilla emcee sampler to minimise $\chi^2_{\tau_{0,2}}$, and the resulting posteriors of $\alpha$ and $\beta$ serve as priors in the cosmological inference step (see Table 2 for their values). We refer the reader to Guerrini et al. (2025) for a detailed description of the methodology of this PSF systematics inference step.

We also conduct a PSF inference run where we use the catalogue containing the set of galaxy ellipticities that have not been corrected for PSF leakage, i.e., not corrected object-wise (see Paper I for details). These ellipticities give different values of $\tau(\theta)$, $\rho(\theta)$, and subsequently different posteriors for $\alpha$ and $\beta$. Figure 6 shows the marginalised posteriors of $\alpha$ and $\beta$ for both the object-wise leakage-corrected and uncorrected ellipticities. The non-leakage-corrected case (in pink) shows a larger mean $\alpha$, as expected from the greater residual PSF leakage when this correction is not applied. On the other hand, $\beta$ remains largely constant, showing that it mainly accounts for the PSF model ellipticity, $e^{\mathrm{p}}$. Nonetheless, leakage not corrected at the object level can be absorbed into a larger $\alpha$, and thus a larger $\xi^{\mathrm{sys}}_{\pm}$ signal, giving consistent cosmological posteriors with the leakage-corrected case.

Vector (PDF)

For comparison, we shall also conduct a cosmological inference run where we adopt flat priors on $\alpha$ and $\beta$ (where $\alpha\in[-0.1, 0.1]$ and $\beta\in[-2.0, 2.0]$), as well as a setup where we use the non-leakage corrected set of ellipticities, to assess the impact of the PSF inference step on the resultant estimate of $\xi^{\mathrm{sys}}_{\pm}$.

Redshift calibration bias estimation

Biases in the calibration of the redshift distribution $n(z)$ propagate into cosmological parameter inference. We therefore marginalise over $\Delta z$, the bias in the mean redshift. Here we describe how we obtain an informed prior on this parameter.

Because the true $n(z)$ of the data is unknown, this bias can only be assessed using simulations, where the true redshifts are available. By constructing mock catalogues that mimic the photometric properties, selection functions, and redshift-calibration steps of the real UNIONS analysis, we directly compare the recovered SOM-based $n(z)$ to the true underlying distribution, obtaining estimates of the systematic shift $\Delta z$.

We construct our mock catalogues from the MICE2 galaxy catalogue, generated from the MICE $N$-body simulation (Fosalba et al. 2015b), which adopts a flat $\Lambda$CDM cosmology and evolves a large-volume dark-matter distribution with sufficient resolution for weak-lensing applications. Dark-matter haloes are identified using a friends-of-friends algorithm (Crocce et al. 2015), and galaxies are populated up to $z=1.4$ through a hybrid scheme combining halo abundance matching with a halo occupation distribution model, calibrated to reproduce the observed luminosity function and galaxy clustering (Carretero et al. 2015). To account for empirical luminosity evolution, we apply the redshift-dependent magnitude-evolution correction of Fosalba et al. (2015a).

To perform the same redshift calibration analysis as in the real UNIONS data, we require that the MICE2 galaxies have realistic multi-band photometry. In the observational data, the UNIONS $r$-band sources are matched to CFHTLenS to obtain five-band $ugriz$ magnitudes (Hildebrandt et al. 2012), which are used for the SOM-based redshift calibration. For the mock catalogues, we therefore replicate this photometric setup by generating CFHTLenS-like noisy fluxes for every simulated galaxy. We adopt a noise-modelling framework inspired by van den Busch et al. [119], adapted to the depth and filter properties of the CFHTLenS imaging (Erben et al. 2013). Through this procedure, we obtain noisy fluxes depending on galaxy size, seeing variations and the intrinsic brightness of each source.

In a further step, we employ a kNN-based matching procedure following the approach presented in Wright et al. (2025b) to ensure that the simulated catalogues replicate the photometric and spectroscopic selection functions of the UNIONS analysis. The UNIONS weak-lensing sample is matched to deeper CFHTLenS $ugriz$ photometry in the W3 field. To reproduce this in the mocks, we match the noisy MICE2 galaxies to the real UNIONS–CFHTLenS sources using kNN-matching in colour–magnitude space. For each real galaxy, we select the nearest simulated object and transfer its multi-band properties to construct a photometric sample that follows the joint distributions of colour and magnitude observed in the data. We also transfer the weights $w$ and the average shear response $\langle R_\gamma\rangle$ from the real UNIONS shear catalogue to the simulations, since we require these in the SOM-based redshift calibration (as described in Sect. 3.3).

To replicate the spectroscopic calibration set consisting of DEEP2, VVDS, and VIPERS, we create a second mock sample using kNN-matching between these surveys and the noisy MICE2 galaxies. Matching in slices of redshift in colour–magnitude space ensures that the mock calibration sample follows the same selection function as the combined real spectroscopic data set.

The mock photometric and spectroscopic samples are processed with the same SOM-based redshift calibration pipeline that is used for the real UNIONS data (Sect. 3.3). The SOM is trained on the spectroscopic-like mock sample, and the photometric mocks are projected onto the trained map to yield a recovered $n(z)$, weighted by both the shape and shear response weights. From the bootstrap resampling, we recover a redshift distribution, which we compare to the true redshift distribution of the simulated sources. This provides a direct measure of the redshift-calibration bias, which we find to be $\Delta z=-0.003$.

The width of the prior on $\Delta z$ is derived directly from the data: we use the standard deviation of the bootstrap realisations of the SOM-based $n(z)$ obtained by resampling the spectroscopic calibration sample. This bootstrap-based estimate of the uncertainty is fully data-driven, and we verified with the mock catalogues that it is consistent with the scatter expected from the true redshift-calibration error budget. Moreover, we include a systematic error from using several mock patches and by changing the parameters of the SOM (e.g., hierarchical clustering). We therefore estimate the total uncertainty to be $\pm0.018$. For completeness, we also consider two additional setups where we impose a flat, uninformative prior on $\Delta z\in[-0.1,0.1]$, and assume no uncertainty (fixing $\Delta z=0$).

Nonlinear matter power spectrum

Cosmic shear probes small scales of the matter power spectrum, requiring accurate modelling of the nonlinear regime. Following the methodology of Doux et al. (2021), we derive the contribution of the nonlinear matter power spectrum at each wavenumber $k$ to the 2PCF signal at different $\theta$ angles, by calculating the integral given in Eq. (2). We adopt the Planck 2018 best-fit cosmology (Planck Collaboration [91]) as our fiducial model, and employ the Limber approximation $k=(\ell+1/2)/\chi(z)$, to compute, for a given $k_{\max}$, the corresponding radial comoving distance up to which we should integrate Eq. (3). Figure 7 presents a heatmap of the ratio of the total 2PCF signal as a function of $k$ and $\theta$ scale, for both $\xi_+$ and $\xi_-$. With the conservative criterion that scales beyond $k_{\max}=3\,h\,\mathrm{Mpc}^{-1}$ contribute less than 10% of the signal, our $n(z)$ implies minimum angular scales of 1.5 arcmin for $\xi_+$ and 11 arcmin for $\xi_-$. However, as described in Sect. 4.6, by considering the effects of $B$-mode contamination, we instead adopt a fiducial angular lower bound of $\theta=12$ arcmin, which would then imply that scales beyond $k_{\max}=0.43\,h\,\mathrm{Mpc}^{-1}$ do not contribute to more than 10% of the $\xi_+$ signal, and scales beyond $k_{\max}=2.85\,h\,\mathrm{Mpc}^{-1}$ for the $\xi_-$ signal. These scale cuts are safely within the range of accuracy of the nonlinear prescriptions adopted below. We mark out these boundaries in Figure 7 as well.

Vector (PDF)

For our fiducial analysis, we use the most updated version of HMCode2020 (Mead et al. 2021), which employs the halo model formalism as well as baryonic feedback modelling to calculate the nonlinear matter power spectrum, with up to 2.5% accuracy at scales $k\leq10\,h\,\mathrm{Mpc}^{-1}$ when compared to $N$-body simulations. It has also been widely adopted in other Stage III cosmic shear analyses (Hikage et al. 2019; Gatti et al. 2021; Wright et al. 2025a).

When modelling the small-scale baryonic feedback arising from active galactic nuclei (AGN), an additional parameter $\log(T_{\mathrm{AGN}})$ is introduced, which quantifies the ‘temperature’ or strength of this feedback. Calibrations with $N$-body simulations estimate it to be around 7.6–8.0. In our fiducial analysis, we marginalise over this parameter, while also including a test where we do not account for baryonic feedback to quantify its impact on the clustering cosmological parameters $\sigma_8$ and $S_8$. We also run a test without baryonic feedback (disabling the ‘feedback’ option in HMCode2020) and one using the nonlinear prescription of Halofit (Takahashi et al. 2012) instead, which is based on a halo-fitting model, to quantify the sensitivity of our data to a difference in nonlinear modelling.

Scale cuts

Our baseline scale cuts were determined by three criteria: PSF systematics, $B$-mode analyses, and potential mismodelling of the nonlinear power spectrum at small scales, as described in Sects. 3.2.1, 3.2.2, and 4.5, respectively. By setting the requirements that the real space and COSEBIs $B$-mode null tests are passed with a PTE of at least 0.05, and that the PSF systematics signal $\xi^{\pm}$ does not contribute more than 10% of the total signal, we obtain $\theta\in[5,83]$ arcmin for $\xi_{+}$ and $\theta\in[12,83]$ for $\xi_{-}$.

To assess the sensitivity of our scale cut choices to cosmology, we conducted blinded inference analyses by successively increasing the values of the lower scale cut of $\xi_{+}$, with $\theta_{\min}$ varying from 3.98 to 12.0 arcmin. The results showed a nontrivial $S_{8}$ dependence on scale, with a systematic downward drift in $S_{8}$ and an improved $\chi^{2}$ value with increasing $\theta_{\min}$. We therefore remove angular bins where the induced shift in the mean $S_{8}$ exceeds $0.2\,\sigma$. This behaviour eventually stabilises, going from the ninth to the tenth angular bin (representing a $\theta_{\min}$ increase from 9.10 arcmin to 12.0 arcmin), where we also find the largest $\chi^{2}$ improvement of 0.7. This could be correlated with the exclusion of the localised feature around 10 arcmin of the $\xi_{+}$ function (see Figure 1), coincident with the 9.4 arcmin size of the MegaCam CCD. This feature persisted with similar amplitude across different variations of PSF size cuts and masking schemes, appearing in both the cosmological and systematic $\xi_{+}$ signals. Additionally, we found that scale cut combinations of $\theta\in[5,83]$ arcmin and

$\theta\in[12,83]$ arcmin pass the COSEBIs $B$-mode null test, but not for the particular case of $\theta\in[9,83]$ arcmin (see Paper II for details). The fiducial scale cut we therefore adopt is $\theta\in[12,83]$ arcmin for both the $\xi_{\pm}$ correlation functions. On the other hand, we utilise the full angular range to fit the $\tau(\theta)$ statistics. This amounts to a total of 54 data points that we fit to: 7 each from $\xi_{\pm}$, and 20 each from $\tau_{0,2}$.

For completeness, we also present inference results where we allow for small-scale contributions in $\xi_{+}$, i.e. when $\theta\in[5,83]$ arcmin.

Pipeline validation with mock catalogues

With our inference pipeline set up, we validate it by running it on 341 random galaxy mock catalogues generated with GLASS (Tessore et al. 2023). The catalogues were created using the same mask, effective number density, and $\sigma_{\epsilon}$ as the data, with a Planck 2018 best-fit fiducial cosmology. In addition, one of the blinded $n(z)$ distributions was chosen at random. We refer the reader to Appendix A for a validation of the mock galaxy survey properties. The mock catalogues are not contaminated by PSF systematics, therefore the $\tau$ statistics are consistent with zero. We nonetheless account for PSF systematics uncertainty by sampling $\alpha$ and $\beta$ with priors obtained from the $\rho$ and $\tau$ statistics of the mock data vectors. Both parameters recover a zero mean to within $1\,\sigma$.

We measure the cosmic shear and PSF correlation functions, compute their respective covariance matrices, and run both the harmonic and configuration space inference pipelines on each mock catalogue, employing the parameters and priors as in Table 2, and the same fiducial scale cuts as detailed in Sect. 4.6. This test acts as a criterion for unblinding (see Sect. 4.8). First, we assess whether we recover the fiducial cosmology within a reasonable level of uncertainty. Additionally, we run both our pipelines on the same set of mocks to ascertain that the average $S_{8}$ mean obtained from both analyses differs by less than $0.5\,\sigma$.

We plot in Figure 8 the histogram of the difference in the 341 best-fit values of $S_{8}$ obtained by both pipelines. We see that when the weighted mean is used to determine the best-fit, we obtain an offset from 0 (albeit still within $1\sigma$), while this discrepancy vanishes when using the 2D marginalised mode to estimate the best-fit instead (albeit with slightly more spread). This demonstrates the robustness of both our pipelines which, despite having been developed independently, are able to deliver internally consistent results. We refer the reader to Appendix A for further details.

Vector (PDF)

catalogue in configuration space, which are given by

Blinding Strategy

We conduct a blinded analysis to remain agnostic to the results and to mitigate potential confirmation biases in our analysis choices. The blinding procedure was implemented by an external collaborator, who generated three versions of the $n(z)$ distributions by applying random shifts to two of them. These shifts induce corresponding offsets in the marginalised $S_8$ posterior distributions. We then computed covariance matrices and performed cosmological inference for all three $n(z)$ variants.

Because the data vector itself was not blinded, we were able to directly apply the tests described in Sect. 3.2 as diagnostics of potential systematic contamination. While the covariance matrix used to compute the $B$-mode PTEs varies with each $n(z)$ shift, we verified that the relative differences between them remain below 10%. The resulting variations in PTE do not affect our conclusions regarding unblinding criteria or scale cuts. We therefore computed PTEs using all three covariance matrices and adopted the most conservative values when defining our scale cuts.

The majority of this manuscript was prepared prior to unblinding and underwent internal review within the broader UNIONS collaboration, including the external blinding coordinator. Cosmological inference was carried out for all three blinded data vectors, alongside variations in modelling (including PSF systematics, nuisance parameter priors, and nonlinear prescriptions) and analysis setup, as described in the preceding sections. We verified that differences in the inferred cosmological posteriors arising from these choices were consistent and exhibited similar trends across all three blinded analyses.

In addition, as described in the previous section, we imposed a further validation criterion requiring consistency between the Fourier- and configuration-space pipelines when applied independently to mock catalogues (Sect. 4.7), specifically in terms of the recovered $S_8$ constraints. Three preliminary versions of the results section and figures were sketched out in advance, such that unblinding required only selecting the appropriate version and finalising the manuscript.

For transparency, we note that the primary authors of the cosmology analysis papers were inadvertently unblinded during the final stages of the project. However, this knowledge did not influence any analysis choices or changes in the inference pipeline. Furthermore, a bug in the code to estimate the mean $\Delta z$ was discovered only after unblinding, which necessitated a rerun of the inference with the correct $\Delta z$ prior. That being said, the post-unblinding rerun did not affect the pipeline validation, since the mock catalogues were run with $\Delta z = 0$. Furthermore, since no analysis-setup changes were made after unblinding, we ensured that the rerun did not introduce any confirmation bias.

Results

Fiducial analysis

We report our best-fit 2D marginalised mode for the $S_8$ parameter, and the matter density $\Omega_{\mathrm{m}}$, as well as the 1D marginalised mode of the amplitude of clustering $\sigma_8$, derived from the UNIONS-3500

catalogue in configuration space, which are given by

Figure 9 presents the whisker plot of the three parameters for all of our analysis runs, as well as the fiducial results obtained in Paper IV, and comparison to external survey results. In Figure 10 we plot the 2D marginalised contours of our fiducial results in configuration and harmonic space, and compare it with those reported in external surveys.

Vector (PDF)

Vector (PDF)

Additionally, our results are consistent at the $1\sigma$ level with cosmic shear analyses from other Stage III galaxy surveys, such as DES Y6 ([22]), KiDS-Legacy (Wright et al. 2025a), and HSC-Y3 ([79] (Dalal et al. 2023)), albeit with a tendency toward systematically higher values. Our measurement provides an independent constraint from the northern hemisphere, among the first of its kind, thus offering a valuable cross-check of other results derived primarily from the southern hemisphere. Our $S_8$ constraint is also fully consistent with that reported by [91], and therefore does not indicate any significant tension between early- and late-time measurements of the clustering amplitude. We anticipate that future tomographic analyses will substantially improve the constraining power of our catalogue, enabling us to become competitive with other Stage III surveys.

We estimate the goodness-of-fit of our result by comparing it to those previously derived from our mocks. In Figure 11 we plot the histogram of mock $\chi^2$ values and fit a $\chi^2$ distribution, where we find that the effective degree of freedom is 51.9, while the $\chi^2$ value of our fiducial analysis is 75.33. This is a more robust method of obtaining the degrees of freedom than simply taking the difference between the number of data points and the number of sampled parameters, since it is often the case in cosmic shear experiments that parameters are prior-dominated and highly degenerate.

Vector (PDF)

We obtain a PTE of 0.0147, which is mainly driven by the fit to the $\tau_{0,2}$ statistics. When focusing on the pure $\xi_+$ data vector, we find a good fit of the $\xi_+$ signal with a PTE of 0.39. In Table B1, we present the best-fit values of $S_8$ and their 68% CI, quoting both the 2D marginalised mode and the weighted mean, together with the $\chi^2$ and PTE for each inference setup considered.

Consistency with harmonic space analysis

We find that our results are internally consistent between both summary statistics, where we obtain an $S_8$ value that is $\sim0.6\sigma$ higher in harmonic space. This difference arises from the distinct range of scales probed by the two approaches.

When expressed in $k$-space, we see that a substantial range of the $\ell$ scales employed in the harmonic space analysis correspond to angular scales of $\lesssim5$ arcmin, which were excluded from the configuration space analysis based on systematics-motivated scale cuts (see Sect. 4.6). Notably, our data exhibit a trend of stronger clustering on these small scales, evident from Figure 1, which then leads to a higher $S_8$ best-fit found in the pseudo-$C_\ell$ analysis. We investigate this in more detail in Sect. 5.2.4, where we include more data points at small scales by varying scale cuts.

Robustness to modelling choices

We now consider the impact of various systematic effects and analysis choices, as described in Sect. 4. In particular, we highlight the most important ones, namely the impact of prior choices on nuisance parameters describing intrinsic alignments, PSF systematics, and redshift calibration uncertainties, as well as choices related to scale cuts and nonlinear modelling.

Intrinsic alignment modelling

In Figure 12 we present the 2D marginalised contours of $S_8$, $\Omega_{\mathrm{m}}$ and $A_{\mathrm{IA}}$ for different IA setups: our fiducial choice of imposing a Gaussian prior on its amplitude $A_{\mathrm{IA}}$, using a flat prior, and not including IA effects entirely (i.e. setting $A_{\mathrm{IA}} = 0$). We see that employing an uninformative flat prior has the largest impact on $S_8$ and $A_{\mathrm{IA}}$ whereby their constraints are greatly broadened, although their best-fit values remain consistent: we obtain $S_8 = 0.895^{+0.172}_{-0.133}$ and $A_{\mathrm{IA}} = 2.61^{+1.55}_{-1.96}$, in the case of a flat $A_{\mathrm{IA}}$ prior, compared to $S_8 = 0.831^{+0.067}_{-0.078}$ and $A_{\mathrm{IA}} = 1.13^{+0.688}_{-0.798}$ in the fiducial setup with a Gaussian prior.

Vector (PDF)

This is due to the strong degeneracy between $A_{\mathrm{IA}}$ and $S_8$, a well-known phenomenon in cosmic shear analyses whereby both parameters work in opposing manners to impact the shear signal. Since an increase in $A_{\mathrm{IA}}$ and thus the IA power spectra $C_{\ell}^{\mathrm{I}\gamma}$ and $C_{\ell}^{\mathrm{II}}$ contribute negatively to the overall shear signal, an enhanced amplitude of clustering would then be required to compensate for it, leading to a strong linear relationship between these two parameters. This is evident in the contour plot of Figure 12. Nonetheless, we are able to recover consistent results from both flat and Gaussian priors, demonstrating the robustness of our technique in estimating the uncertainty of $A_{\mathrm{IA}}$ based on galaxy populations (Sect. 4.2).

On the other hand, not accounting for IA naturally gives the tightest constraints, at the expense of a slight bias: we obtain a lower value of $S_8 = 0.770^{+0.029}_{-0.059}$, which is the largest shift in $S_8$ from the fiducial amongst all the inference tests conducted. Moreover, from Table B1 we see that this modelling choice returns the largest $\chi^2$ and a very low PTE value, thus signifying that the model is a sub-optimal fit to the data. We anticipate that with future tomographic analyses, additional redshift information can break the degeneracy between $S_8$ and $A_{\mathrm{IA}}$, which would also imply larger biases should the effects of IA not be modelled accurately.

Redshift calibration

We subsequently consider the impact of our choice of priors for $\Delta z$, the uncertainty of the shift in $n(z)$ obtained during redshift estimation and calibration. We present in Figure 13 the 2D marginalised contours of $S_8$, $\Omega_{\mathrm{m}}$ and $\Delta z$ where, analogous to the previous section, we compare the results from three different analysis choices: the fiducial case of imposing a Gaussian prior, using a flat prior, and fixing $\Delta z = 0$. We see that our results on $S_8$ remain largely consistent in all cases, save for a significant widening of the $\Delta z$ posterior when an uninformative prior is employed: $\Delta z = -0.0450^{+0.077}_{-0.036}$ compared to the fiducial result of $\Delta z = -0.0043^{+0.016}_{-0.018}$, which equates to a three times increase in the size of error bars. Since the estimated mean value of $z$ is already very small (on the order of $10^{-3}$), fixing it at zero does not produce an obvious shift in $S_8$ results. However, this would not be the case should the magnitude of $\Delta z$ be larger.

A degeneracy between $S_8$ and $\Delta z$ is also evident from the contour plot, although in a different direction to that of $A_{\mathrm{IA}}$; once again, this has to do with how they impact the amplitude of the $\xi_{\pm}$ signal. A larger $\Delta z$ would imply a galaxy source distribution with higher $n(z)$ and thus stronger lensing efficiency, as evident from Eq. (4), which would then require a lower $S_8$ to counteract. Inclusion of tomographic information would then be able to tighten constraints on these parameters. Furthermore, we see that even when an uncertainty in the mean redshift shift is not taken into account, we recover consistent results, reaffirming the accuracy and robustness of our SOM calibration methodology.

PSF systematics

We quantify the impact of PSF systematics within our data set, and assess our ability to mitigate it. In Figure 14 we plot the 2D marginalised distributions of $S_8$, $\Omega_{\mathrm{m}}$, and the leakage parameters $\alpha_{\mathrm{PSF}}$ and $\beta_{\mathrm{PSF}}$ for four different cases: when imposing Gaussian priors on $\alpha_{\mathrm{PSF}}$ and $\beta_{\mathrm{PSF}}$, when using flat priors instead, when not accounting for $\xi_{\pm}^{\mathrm{sys}}$ (i.e. not doing a joint fit with the $\tau_{0,2}$ statistics), and when using instead the catalogue derived from the set of ellipticities that were not corrected for per-object PSF leakage. We see that our results for the cosmological parameters are almost perfectly consistent. Employing flat priors on $\alpha$ and $\beta$ works only to slightly enlarge their

Vector (PDF)

resultant posterior distributions, demonstrating that even without the additional step of first fitting for $\alpha$ and $\beta$ (as outlined in Sect. 4.3), we are still able to accurately recover their posteriors.

Moreover, comparing the pink and orange contours in the $\Omega_{\mathrm{m}} - S_8$ plane of Figure 14, we find that the impact of $\xi^{\mathrm{sys}}$ is negligible, as evidenced by the recovery of consistent $S_8$ posteriors even when it is not jointly modelled in the inference step. This indicates that, at the current level of data precision, adopting a threshold of 10% for the maximum allowed contribution of $\xi^{\mathrm{sys}}$ to the total signal, as was used to define our scale cuts, is sufficient to mitigate any bias arising from PSF contamination. Furthermore, even when an object-wise leakage correction is not applied (grey contours), it can still be accounted for by using the corresponding set of priors for $\alpha$ and $\beta$ instead. This exercise highlights how PSF leakage correction is not just assumed to work at catalogue level, but is also propagated into the cosmological inference with a non-negligible impact.

Scale cuts

We consider a setup where we include more contribution from small scales, imposing a scale cut of $\theta\in[5,83]$ arcmin for the $\xi_{+}$ data vector. This gives an increased $S_8 = 0.886^{+0.049}_{-0.104}$, which nonetheless remains consistent with the fiducial result, and does not seem to significantly tighten the contours in the $\Omega_{\mathrm{m}} - S_8$ plane (see Figure B1).

Vector (PDF)

Not surprisingly, we also find that the $S_8$ results obtained when including smaller angular scales become more consistent with the fiducial harmonic space results (where $S_8 = 0.891^{+0.057}_{-0.084}$), and the small-scale best-fit $\xi_{+}$ vector in Figure 15 moves towards that derived from the harmonic space analysis. This further reinforces the observation made in Sect. 5.1.1: our data exhibit enhanced clustering on small scales, which drives the higher best-fit value of $S_8$ as the lower scale cut is reduced from 12 arcmin to 5 arcmin. When expressed in terms of $k$-space, this smaller scale cut combination has a larger overlap with the $\ell$ range probed by the pseudo-$C_\ell$ estimator, leading to more consistent $S_8$ values between the two approaches.

Vector (PDF)

However, when naively computing the reduced $\chi^2$ of the pure $\xi_{\pm}$ data vector, where $\chi^2_{\mathrm{red}}=\chi^2/n$ and $n$ is the number of data points, we obtain $\chi^2_{\mathrm{red}}=1.34$, which gives a worse fit than in the fiducial case (where $\chi^2_{\mathrm{red}}=0.68$), most likely owing to the anomalously high data point at 10 arcmin.

Nonlinear scales

Finally, we consider the impact of modelling choice for the nonlinear power spectrum. We present the best-fit data vectors for all three cases considered in Figure B2: the fiducial case of employing HMCode2020 and including the effects of baryonic feedback, HMCode2020 without baryonic feedback, and using Halofit instead. We see that not including the effects of baryons leads to the lowest amplitude of clustering, while Halofit gives largely consistent results with our fiducial analysis.

Vector (PDF)Vector (PDF)

From Figure 9, we do not see appreciable changes in the best-fit values of $S_8$, $\sigma_8$ or $\Omega_{\mathrm{m}}$. This is generally to be expected, since our adopted scale cuts exclude most of the small scales where nonlinear modelling makes a significant impact (in our case, 90% of the signal receives contributions from $k$-scales larger than $0.43\,h\,\mathrm{Mpc}^{-1}$). Additionally, using a single tomographic bin results in a broad projection along the line of sight, mixing contributions from a wide range of redshifts and effectively diluting sensitivity to small-scale information. We thus expect to unlock more constraining power from the small scales, and better constrain the amplitude of baryonic feedback log($T_{\mathrm{AGN}}$) with future tomographic analyses and more precise data in future releases.

Combination with external data sets

Finally, we consider combined constraints by incorporating external data, including CMB from Planck 2018 and BAO from DESI DR2, which we present in Figure 16. We find that, for the UNIONS-3500 + Planck combination, the constraints on $S_8$ are largely driven by CMB data. In particular, the orange contour reduces to the blue one, which are themselves nearly identical to the green contour obtained using Planck-only data (with the plik_lite likelihood). We note that our Planck-only posteriors do not exactly match the fiducial results published in Planck Collaboration [91], due to differences in sampler, likelihood, parameters and priors employed in our inference pipeline. Including BAO then reduces constraints on $\Omega_{\mathrm{m}}$ to that favoured by the data, while $\sigma_8$ stays approximately constant since BAO in itself is insensitive towards the amplitude of clustering.

Vector (PDF)

Conclusion

We have conducted the first 2D cosmic shear analysis and presented cosmological constraints in configuration space with the UNIONS-3500 photometric weak lensing catalogue, the largest of its kind in the northern hemisphere to date. We report a fiducial result of $S_8 = 0.831^{+0.067}_{-0.078}$ and $\Omega_{\mathrm{m}} = 0.265^{+0.130}_{-0.075}$ with a flat $\Lambda$CDM model. This is $\sim0.6\,\sigma$ lower compared to the analysis in harmonic space. This difference is primarily driven by the difference in $k$-space sampled by the two summary statistics, which then responded differently to the scale-dependent features in our data. Furthermore, our results are consistent with those from previous Stage III weak lensing surveys, as well as with early-time constraints from Planck CMB data.

We have also outlined our methodology to derive these constraints, ensuring their robustness and ability to mitigate contamination due to systematic effects. We adopted $B$-mode-informed scale cuts, and ensured that PSF systematics is adequately accounted for by jointly fitting the systematic and cosmological signals, on top of employing an additional safeguard of first sampling the leakage parameters $\alpha_{\mathrm{PSF}}$ and $\beta_{\mathrm{PSF}}$ to obtain more data-informed priors in the subsequent inference step. We estimated the redshift distribution using SOMs and validate our catalogue with dedicated simulations, which are also used to assess shear calibration biases. We then derived data-driven priors on the nuisance parameters quantifying the magnification bias, the redshift uncertainty and the amplitude of intrinsic alignment, thereby maximising the constraining power of our analysis. Finally, both harmonic and configuration space pipelines were simultaneously validated using mock galaxy catalogues, enabling a direct assessment of internal consistency before unblinding; to our knowledge, this represents a novel approach in cosmic shear analyses.

Our final results demonstrate the reliability of our analytical methodology across a wide range of modelling assumptions and analysis configurations. In particular, we find that the inferred cosmological constraints are most sensitive to intrinsic alignment modelling and uncertainties in the source redshift distribution. Adopting a flat prior on $A_{\mathrm{IA}}$ increases the uncertainty on $S_8$ by nearly a factor of two, while fixing $A_{\mathrm{IA}}=0$ produces a shift of $-0.71\sigma$ in the recovered $S_8$ value, reflecting the well-known degeneracy between intrinsic alignments and the clustering amplitude. Similarly, adopting a flat prior on the redshift shift parameter $\Delta z$ produces the second-largest increase in the $S_8$ uncertainty among the systematic tests considered. Taken together, these findings might suggest that the current constraining power of our data is driven not only by the statistical precision of the shear measurements, but also by nuisance parameter degeneracies and conservative systematic treatments.

Yet it is important to note that the priors adopted for these nuisance parameters were intentionally chosen to be conservative relative to their statistical uncertainties, in order to avoid underestimating the impact of systematic effects. This then opens up possibilities for future improvement going forward, both in terms of strengthening the systematic characterisations of redshift distributions, shear biases, PSF systematics and intrinsic alignment models, as well as including more data with future releases and from external surveys.

In this regard, we have already shown that combining UNIONS with external probes such as Planck CMB measurements and DESI BAO data substantially tightens constraints on both $S_8$ and $\Omega_{\mathrm{m}}$. A natural next step will therefore be to conduct a full tomographic 3 × 2 point analysis incorporating spectroscopic galaxy clustering and galaxy–galaxy lensing measurements using DESI lens galaxies, enabling us to fully exploit the substantial overlap between the two surveys across the northern sky.

As one of the final Stage III weak lensing surveys, UNIONS provides a valuable complement to preceding cosmic shear measurements, offering an independent validation of results with comparable robustness in methodology. Looking ahead, the techniques and validation strategies developed in this work can help pave the way for upcoming Stage IV surveys, whose unprecedented statistical power and improved control of systematics promise to usher in a new and exciting era of weak lensing science.

Acknowledgements

We would like to thank our external blinding coordinator, Koen Kuijken. We would also like to thank Axel Guinot and Douglas Scott for their comments and helpful feedback. LWK Goh thanks the University of Edinburgh School of Physics and Astronomy for a postdoctoral Fellowship. HH is supported by a DFG Heisenberg grant (Hi 1495/5-1), the DFG Collaborative Research Center SFB1491, an ERC Consolidator Grant (No. 770935), and the DLR project 50QE2305. MJH and LVW acknowledge support from NSERC through a Discovery Grant. AHW is supported by the Deutsches Zentrum für Luft- und Raumfahrt (DLR), under project 50QE2305, made possible by the Bundesministerium für Wirtschaft und Klimaschutz, and acknowledges funding from the German Science Foundation DFG, via the Collaborative Research Center SFB1491 "Cosmic Interacting Matters - From Source to Signal". This work was made possible by utilising the CANDIDE cluster at the Institut d’Astrophysique de Paris. The cluster was funded through grants from the PNCG, CNES, DIM-ACAV, the Euclid Consortium, and the Danish National Research Foundation Cosmic Dawn Center (DNRF140). The authors acknowledge the use of the Canadian Advanced Network for Astronomy Research (CANFAR) Science Platform operated by the Canadian Astronomical Data Centre (CADC) and the Digital Research Alliance of Canada (DRAC), with support from the National Research Council of Canada (NRC), the Canadian Space Agency (CSA), CANARIE, and the Canada Foundation for Innovation (CFI). It is maintained by Stephane Rouberol. We are honoured and grateful for the opportunity of observing the Universe from Maunakea and Haleakala, which both have cultural, historical and natural significance in Hawai’i. This work is based on data obtained as part of the Canada-France Imaging Survey, a CFHT large program of the National Research Council of Canada and the French Centre National de la Recherche Scientifique. Based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA Saclay, at the Canada-France-Hawai’i Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the Institut National des Science de l’Univers (INSU) of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. This research is based in part on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. Pan-STARRS is a project of the Institute for Astronomy of the University of Hawai’i, and is supported by the NASA SSO Near Earth Observation Program under grants 80NSSC18K0971, NNX14AM74G, NNX12AR65G, NNX13AQ47G, NNX08AR22G, 80NSSC21K1572 and by the State of Hawai’i.

Data Availability

A subset of the raw data underlying this article is publicly available via the Canadian Astronomical Data Centre at http://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/en/megapipe/. The remaining raw data and all processed data are available to members of the Canadian and French communities via reasonable requests to the principal investigators of the Canada-France Imaging Survey, Alan McConnachie and Jean-Charles Cuillandre. All inference chains and software used to produce the results will be publicly available to the international community upon paper acceptance.

as expected. We do not expect systematic effects to contribute to the mock signal, such that the nuisance parameters listed in Table 2 should be consistent with zero. Nonetheless, we sample over them, keeping their Gaussian priors with the same standard deviations but centred over a zero mean.

We also plot the histogram of the 2D marginalised modes of $\Omega_{\mathrm{m}}$ and $S_8$ in Figure A2, specifically those obtained in the configuration space mock pipeline analysis.

Vector (PDF)

Validation with GLASS mocks

We present validation tests using the GLASS galaxy mock catalogues. The redshift distribution of the mocks agrees with the data (Figure A1),

Vector (PDF)

Additional inference results

We present additional results from our inference.

This paper has been typeset from a TEX/LATEX file prepared by the author.

Vector (PDF)