UNIONS: The impact of systematic errors on weak-lensing peak counts

UNIONS is an ongoing deep photometric multi-band survey of the Northern sky. As part of UNIONS, CFIS provides r-band data which we use to study weak-lensing peak counts for cosmological inference. We assess systematic effects for weak-lensing peak counts and their impact on cosmological parameters for the UNIONS survey. In particular, we present results on local calibration, metacalibration shear bias, baryonic feedback, the source galaxy redshift estimate, intrinsic alignment, and the cluster member dilution. For each uncertainty and systematic effect, we describe our mitigation scheme and the impact on cosmological parameter constraints. We obtain constraints on cosmological parameters from MCMC using CFIS data and MassiveNuS N-body simulations as a model for peak counts statistics. Depending on the calibration (local versus global, and the inclusion of the residual multiplicative shear bias), the mean matter density parameter $\Omega_m$ can shift up to $-0.024$ ($-0.5\sigma$). We also see that including baryonic corrections can shift $\Omega_m$ by $+0.027$ ($+0.5 \sigma$) with respect to the DM-only simulations. Reducing the impact of the intrinsic alignment and cluster member dilution through signal-to-noise cuts can lead to a shift in $\Omega_m$ of $+0.027$ ($+0.5 \sigma$). Finally, with a mean redshift uncertainty of $\Delta \bar{z} = 0.03$, we see that the shift of $\Omega_m$ ($+0.001$ which corresponds to $+0.02 \sigma$) is not significant. This paper investigates for the first time with UNIONS weak-lensing data and peak counts the impact of systematic effects. The value of $\Omega_m$ is the most impacted and can shift up to $\sim 0.03$ which corresponds to $0.5\sigma$ depending on the choices for each systematics. We expect constraints to become more reliable with future (larger) data catalogues, for which the current pipeline will provide a starting point.


Introduction
Weak gravitational lensing has been used as a cosmological probe in recent years with great success, for example with Dark Energy Survey 1 (DES), Kilo-Degree Survey 2 (KiDS), Hyper Suprime-Cam 3 (HSC), and Canada France Hawaii Lensing Survey 4 (CFHTLens). It corresponds to the small distortions we observe in the images of background sources (such as highredshift galaxies) due to the deflection of photons as they pass et al. 2019; Ajani et al. 2020, and references therein), the starlet 1 norm (Ajani et al. 2021), the scattering transform (Cheng et al. 2020), wavelet phase harmonic statistics (Allys et al. 2020), and machine learning-based methods (Fluri et al. 2018;Shirasaki et al. 2021;Fluri et al. 2022, among others), have been introduced to account for non-Gaussian information.
In this work we choose to focus on peak counts extracted from real data, using the pipeline developed to perform the study presented in Ajani et al. (2020), which had previously only been tested on simulations. Peaks in weak-lensing convergence maps are tracers of overdense regions. They are the local maxima defined as a pixel that is larger than all eight of its neighbors. The peak function -that is, the number of peaks as a function of peak height (in a convergence map) or signal-to-noise ratio (S/N)-depends on the nonlinear and non-Gaussian part of the LSS. This higher-order weak-lensing statistic can be used to constrain cosmological parameters. Peak counts are complementary to second-order shear statistics (Jain & Waerbeke 2000), and by combining both parameters degeneracies can be removed (Dietrich & Hartlap 2010).
Weak-lensing peaks are an indirect tracer of dark-matter halos: large peaks are strongly correlated with massive halos, whereas lower-amplitude peaks are generally created by multiple smaller halos along the line of sight (Yang et al. 2011). Lowamplitude peaks can also be caused by mass outside dark-matter halos, or by galaxy shape noise (Liu & Haiman 2016;Yang et al. 2011).
Explicit expressions or complete theoretical predictions of peak counts are still an active area of research. It is, however, possible to generate weak-lensing simulations densely sampled in cosmological parameter space in order to interpolate them and use the interpolation as a prediction. The advantage of simulations is the possibility to incorporate the exact survey mask and shape noise. For example, Dietrich & Hartlap (2010) created a set of N-body simulations in the (Ω m , σ 8 ) plane for 158 cosmologies, and Liu et al. (2015b) and Kacprzak et al. (2016) used ray-tracing N-body simulations in the (Ω m , σ 8 , w) plane for 91 cosmologies. Here we use the MassiveNuS simulations (Liu et al. 2018) to predict peak counts as in Li et al. (2019) and Ajani et al. (2020). These simulations are described in more detail in Sect. 3.
The first cosmological constraints from peak counts were obtained on real data by Liu et al. (2015c) using the Canada-France-Hawaii Telescope (CFHT) Stripe 82 Survey, and by Liu et al. (2015b), using CFHTLenS data. DES Science Verification data have been analyzed by Kacprzak et al. (2016). KiDS 450 deg 2 data have been studied by Shan et al. (2017) and Martinet et al. (2017). The first tomographic analysis was performed by Harnois-Déraps et al. (2021a) for the DES-Y1 data release. Recently, Zürcher et al. (2022) analyzed the DES-Y3 data release using an emulator approach. These analyses all use peak counts and complement second-order statistics analyses, constraining cosmological parameters.
Weak-lensing observables have to be corrected for systematic effects, which can have an observational and astrophysical origin and can induce biases into the cosmological constraints if not properly taken into account. These artifacts can easily be created by the atmosphere, the telescope, and the detector, and during the data analysis. Astrophysical correlations such as intrinsic galaxy alignments add to the lensing correlations in a nontrivial way. Furthermore, to be able to interpret weak-lensing observables in a cosmological context, the physics of small scales needs to be reliably estimated.
Extensive studies of many of these systematics for secondorder statistics exist, for example Massey et al. (2012), Troxel, M.A. & Ishak, Mustapha (2015), Hildebrandt et al. (2016a), and Mandelbaum (2018). For current and future surveys, precise requirements on the instrument and data processing are routinely derived for those statistics. However, how to properly include systematic effects in the context of higher-order statistics is an ongoing research topic (Kacprzak et al. 2016;Coulton et al. 2020;Zürcher et al. 2021;Harnois-Déraps et al. 2021b;Pyne & Joachimi 2021).
This paper is the first to study several weak-lensing systematics and uncertainty, their effect on peak counts, and the resulting constraints on cosmological parameters from the Ultraviolet Near-Infrared Optical Northern Survey (UNIONS)/Canada-France Imaging Survey (CFIS) galaxy survey. In addition, we develop a novel method for locally calibrating the measured shear, including shear and selection biases. We investigate mitigation schemes for these effects and quantify biases in cosmological parameters.
This paper is organized as follows: In Sects. 2 and 3 we describe the data and simulations used in this work, respectively. In Sect. 4 we describe the measurement of weak-lensing peak counts. In Sect. 5 we introduce local shear calibration and compare it to the standard global calibration. In Sect. 6 we present results for each of the studied biases and uncertainties and discuss mitigation methods. Finally, we draw our conclusions in Sect. 7.

Image processing
We used weak-lensing data from UNIONS 5 . This ongoing survey will provide 4, 800 deg 2 of multiband photometric images in the Northern Hemisphere. Founded in 2018, UNIONS is a collaboration of several groups and surveys, each providing data in different bands. These participating surveys are CFIS, the Panoramic Survey Telescope And Rapid Response System (Pan-STARRS), Wide Imaging with Subaru HSC of the Euclid Sky (WISHES), and the Waterloo Hawai'i IfA G-band Survey (WHIGS).
For the CFIS part of the survey r-and u-band images are taken at the CFHT with MegaCAM, a wide-field optical imaging facility. The CFIS r-band images have a median seeing of 0.65 arcsec reflecting the extremely stable atmosphere at Manua Kea together with the excellent CFHT optical system. Each observed sky location is covered by at least three single exposures, which are dithered by one-third of the focal plane, or 0.33 • . The exposure time varies between 100 and 300 seconds, where smaller exposure times are chosen in better observing conditions. This survey strategy provides images of a very homogeneous depth. In this work, we only use r-band data from CFIS for our weaklensing peak count study. In particular, the weak-lensing data from P3, a patch of size 34.7 × 17.7 deg 2 is analyzed. This patch overlaps with CFHTLenS, which is used to infer the redshift distribution (see Sect. 2.3). These data have been processed and validated (see Guinot et al. 2022 for details). The preprocessing of the CFHT data consists of a calibration step with the MegaPipe pipeline (Gwyn 2008). The single exposures are first calibrated before they are combined with SWARP 6 to build stacked images. MegaPipe provides an astrometric calibration of the survey using Gaia Data Release 2 (Brown et al. 2018), and photometric number density of galaxies n gal 7 arcmin −2 pixel size A pix 0.4 2 arcmin 2 /px 2 global ellipticity dispersion σ e 0.44 size of the field 34.7 × 17.7 deg 2 calibration of the r-band data relative to the Pan-STARRS1 survey (Chambers et al. 2016). Both astrometric and photometric calibration are excellent so we do not expect calibration errors to significantly impact weak-lensing measurements (Guinot et al. 2022).

The weak-lensing catalog
Since the stacks have a larger S/N than the single-exposure images, we use the former to detect galaxy candidates. Due to the large dithers, the point-spread function (PSF) on the stacks is very inhomogeneous and discontinuous. For this reason, we detect stars and construct the PSF model on single exposures. Galaxy shapes are then obtained using the multi-epoch modelfitting method ngmix (Sheldon & Huff 2017) 7 . Those measurements are calibrated with metacalibration (Huff & Mandelbaum 2017) to provide shear estimates. The creation and validation of the shear catalog are fully described in Guinot et al. (2022). We note that this is a preliminary version, called "version 0," of the catalog where the source density is conservative and does not reflect future versions of the CFIS shear data. Specifics of the shear catalog are presented in Table 1.

Redshift distribution of source galaxies
Multiband data of UNIONS are still sparse, and photometric redshifts have not been obtained yet. Using the r-band data only, redshifts have been obtained for a population of galaxies with two different methods, as follows. For both methods, CFIS galaxies were matched to the deeper CFHTLenS on the 50 deg 2 W3 field. First, Guinot et al. (2022) approximated the CFIS redshift distribution as the histogram of the best-fit photometric redshifts from CFHTLenS of that matched subsample. Photo-z's for CFHTLenS had been obtained in Hildebrandt et al. (2012) from u, g, i, r, z multiband data, calibrated with various spectroscopic deep data sets. Second, employing the direct calibration technique, Spitzer et al. (2022) re-weighted the Deep Extragalactic Evolutionary Probe 2 (DEEP2) spectroscopic sample (Newman et al. 2013) in a 5D space spanned by the u, g, i, r, z photometric bands, to match the density of the matched subsample in that space. The re-weighted DEEP2 spectroscopic redshift distribution is an estimate of the CFIS n(z). This distribution was fit in Spitzer et al. (2022) by an analytical function with two components. The first, exponential, component accounts for the bulk of the distribution, whereas the second, Gaussian, term, models the tail at z > 2. However, the addition of this second term somewhat overestimates the reweighted DEEP2 redshift distribution between z = 1.5 and 2.
The mean redshift is obtained asz = z max 0 dz z n(z), where the integral over the normalized redshift distribution n(z) is carried out nominally up to the limiting redshift of the survey, z max . Both redshift distributions have a significantly non-vanishing probability at high redshifts z 2. This is most likely not a physical feature, since we do not expect a large number of galaxies in the CFIS sample at those redshifts. Our best estimate ofz is obtained by integrating over the CFHTLenS-matched n(z) with a limit of z max = 2, resulting in z = 0.65. Table 2 presents alternative estimates with varying n(z) and z max . Three further, reasonable combinations yield a slightly higherz of 0.68. We use this alternative value in Sect. 6.3 to test the impact of the estimated redshift uncertainty.

Simulations
To get the predictions for the summary statistics that we use to perform cosmological inference, we employed the MassiveNuS simulations, a suite of cosmological dark-matter-only N-body simulations that explore different cosmologies including massive neutrinos in the range m ν = 0 -0.62 eV. The simulations have a 512 Mpc h −1 box size with 1024 3 cold dark matter particles. The pixel size is 0.4 arcmin. The implementation is performed using a modified version of the public tree-Particle Mesh (tree-PM) code Gadget2 8 with a neutrino patch, describing the effect of massive neutrinos on the growth of structures up to k = 10h Mpc −1 . A complete description of the implementation and the products is provided in Liu et al. (2018). The cosmological parameters vary across the simulations within the range M ν ∈ [0, 0.62], Ω m ∈ [0.18, 0.42], and A s ∈ [1.29, 2.91] × 10 −9 . We thus worked on the constraints on these three cosmological parameters, which are well sampled by the simulations. They include the effects of radiation on the background expansion and the impact of massive neutrinos is included with a linearresponds method: neutrinos are evolved perturbatively, while their clustering is caused by the nonlinear dark-matter evolution.
The simulations assume a flat universe with Hubble constant H 0 = 70 km s −1 Mpc −1 . The primordial power-spectrum scalar index is n s = 0.97, the baryon density Ω b = 0.046, and the darkenergy equation of state w = −1. A fiducial cosmology is set to [M ν , Ω m , 10 9 × A s ] = [0.1, 0.3, 2.1]. Simulations are available for 101 cosmologies, with 10000 realizations for each cosmology, obtained by randomly rotating and shifting the lensing potential planes.
Following Ajani et al. (2020), the peak counts were computed for each of the MassiveNuS cosmologies from the simulated convergence maps, averaged over the 10000 realizations for each model. A model with massless neutrinos corresponding to [M ν , Ω m , 10 9 × A s ] = [0.0, 0.3, 2.1] is also provided, and we used it to compute the covariance matrix.

Effective redshift distribution
For each of the 101 cosmological models, five source redshifts, z s = {0.5, 1.0, 1.5, 2.0, 2.5}, are present. To match the simulations to the observed redshift distribution of CFIS, we made the following approximations. We matched the mean redshift between MassiveNuS and CFIS (Sect. 2.3) and neglected the shape of the redshift distribution. For that, we interpolated the two convergence maps closest in redshift toz and obtained the new effective map as κz = κ z=0.5 λ + κ z=1 (1 − λ). (1) This effectively defines the new redshift distribution by a weighted sum of two Dirac delta distributions, To match our best estimate,z = 0.65, we set λ = 0.7. As anticipated in Sect. 2.3, we also employedz = 0.68 in order to estimate the uncertainty related to the redshift estimation. To do so, we used the relation in Eq. 1 to the new effective maps at κz =0.65 and κ z=1 and imposedz = 0.68. This results in an interpolation parameter λ = 0.91.

Noise
To include the CFIS shape noise in the simulations that we employ to perform inference, we add Gaussian noise. First, we compute the global ellipticity dispersion σ 2 e . The ellipticity noise per smoothed pixel is then where n gal is the galaxy number density and A pix is the pixel size.
For CFIS data, we used the value listed in Table 1. We applied the Gaussian noise to the convergence maps because we do not have direct access to the simulated shear.

Weak-lensing peak counts
We employed weak-lensing peak counts as summary statistics for our analysis, measured from the weak-lensing convergence maps.

Convergence maps
The reduced shear g is the main observable in weak gravitational lensing by galaxies, it is estimated from galaxy ellipticities. We can define it as where γ is the shear and κ the convergence. For the data, the Eand B-mode convergence maps are built through Kaiser-Squires inversion (Kaiser & Squires 1993) from the reduced shear provided by the weak-lensing catalog. The CFIS-P3 map is quite large for a projection but peaks are local, and any distortion will rather affect large scales. In our analysis we worked with the E-mode of the convergence map as the B-mode contain mostly noise (Guinot et al. 2022). Hereafter, when we speak about convergence, we mean the E-mode of the convergence map. The CFIS-P3 convergence map is shown in Fig. 1. Every pixel has a size of 0.4 arcmin, as in the simulations. The 13 black squares are boxes of 512 × 512 pixels, which correspond to the size of the simulation convergence maps. These boxes are placed such that they do not overlap with larger masked or missing areas. For the simulations, we employed the already existing maps from the MassiveNus 9 suite, obtained with the LensTools 10 (Petri 2016) ray-tracing package. We mimicked the CFIS shape noise by adding the noise introduced in Sect. 3.2 to the maps coming from the simulations, and we smoothed them with a Gaussian kernel of width σ = 2 arcmin. Li et al. (2019) find that for the same simulations, using peak counts, the optimal smoothing to get tighter constraints is 2 arcmin. As a first work, we chose this smoothing scale to be consistent with the simulations. Moreover, Hildebrandt et al. (2016b) find that a smoothing scale of 2 arcmins is a good choice for KiDS data, which have the same number of galaxies per pixel as us. For future work, it can be interesting to adapt simulations to the data to determine the optimal smoothing scale for future releases. We thus computed signal-to-noise ratio (S/N) maps, where the S/N is defined as the noisy convergence map, smoothed with a Gaussian filter over the standard deviation of the noise defined in Eq. (3). Then we computed the peak counts with the lenspack 11 python package on the S/N maps collecting the local maxima, namely computing the pixels with higher values with respect to their neighboring pixels. In our analysis, we consider linearly spaced bins in the range S/N= [−2, 6]. The peak counts distribution used for parameter inference (Sect. 6) corresponds to the mean of the peak counts from the 13 patches. This corresponds to an area of 13 × 12.25 deg 2 ≈ 160 deg 2 . Figure 2 presents the S/N map and the peak counts histogram from one patch of the CFIS-P3 data.

Modeling and parameter inference
We modeled the peak function with the MassiveNuS N-body and ray-tracing simulations (Liu et al. 2018). In these simulations the matter density parameter, Ω m , the primordial powerspectrum normalization amplitude, A s , and the total mass of neutrinos, M ν = m ν , were varied (see Sect. 3 for more details). We constructed a likelihood function as follows. First, the theoretical model was obtained by numerically computing the peak function for each simulated parameter combination. These functions were then interpolated to arbitrary parameters within the parameter boundaries using a Gaussian process. The error of prediction was always below the CFIS statistical error and on the order of a few per cent. The covariance matrix was computed from the variations of realizations at the fiducial model and is thus assumed to be parameter-independent. As explained in the previous section, for the data, the peaks were computed as the mean of the peaks on 13 mask-free patches; thus, the re-scaling factor of the covariance when we infer parameters from the data is 1/13. The likelihood was taken as a multivariate Gaussian as a function of the data vector. As prior on the parameters, we used a flat prior for the three parameters: M ν ∈ [0.06, 0.62], Ω m ∈ [0.18, 0.42] and A s ∈ [1.29, 2.91] × 10 9 . The prior of Ω m and A s are the bounds of the simulation. For M ν , 0.06 is the minimum from oscillation experiments (Patrignani et al. 2016) and 0.62 is the upper bound of the simulations. We explored the posterior distribution with Monte Carlo Markov chain (MCMC) sampling using the python package emcee. Specifically, we employed 250 walkers initialized in a tiny Gaussian ball of radius 10 −3 around the fiducial The squares indicate the regions free of large masks, which were used to compute the peak count. The total peak count is the mean of the peaks over the 13 patches.

Local shear calibration
In this section, we describe and measure a spatially varying calibration of the estimated shear. For the multiplicative shear bias, including both shear and selection biases, we use the technique of "metacalibration" (Huff & Mandelbaum 2017), which we briefly describe in the following subsection.

Metacalibration
This method consists in measuring the response matrix R of a shape measurement algorithm to a shear artificially applied to an image. The i th component of the observed ellipticity of a galaxy, ε i , is an estimator of the galaxy shear, ε i = g obs i , however, not an unbiased estimator in general. The linearized relation between true and estimated shear for a given galaxy can be written as where the response matrix can be described as R i j = ∂g obs i /∂g true i and c i is the additive bias. The trace of the shear response matrix is also parameterized as tr (R) = 2(1 + m) where m is the multiplicative shear bias. We numerically compute the shear response matrix by replacing the derivative with finite differences. For this, we create four new images for each detected galaxy. These images are deconvolved with the model PSF at that position, then sheared in two directions for both shear components, and re-convolved with a circularized PSF that is slightly larger than the original one. We note that a fifth image is created, which is unsheared but has the same re-convolved PSF. This image is used to measure shapes that are consistent with the metacalibration correction.
The response matrix has two additive components, the shear response matrix and the selection response matrix, both of which can be computed via metacalibration, as follows. First, the shear response matrix is computed for each galaxy individually. Denoting with g obs,±, j i the observed shear of the galaxy image sheared by ±∆g in the j direction, we get Article number, page 5 of 19 A&A proofs: manuscript no. aanda We used ∆g = 0.01. Since this is a very noisy measurement, often resulting in singular matrices, we compute the mean value over all galaxies, Second, the selection response matrix quantifies selection biases originating from correlations between shear and selectiondependent image properties. For example, if a cut is applied on the S/N of the galaxy sample and the S/N varies with shear, the effective sample shear after the cut is modified, and no longer representative of the underlying population shear, which introduces a bias.
The selection matrix cannot be obtained for a single galaxy, but only for a sample of galaxies. For that, we first apply our selection S to the four samples of sheared galaxies. This provides us with a selection mask S ±, j for each of the four cases. These four masks are now in turn applied to the fifth sample without added shear. We compute the mean observed shear for each of the four masked samples, denoted as g obs i S ±, j . Any difference in these mean values is due to shear-dependent selection criteria, and can be used to define the selection response matrix as The total mean response matrix is The estimated shears are calibrated by matrix multiplication with R −1 , to obtain g true i for all i following Eq. (5). In the following sections of the paper, we will need the use the diagonal and off-diagonal terms of the R shear and R selection matrix. They are defined as

Global values of the multiplicative shear bias and additive bias
Usually, the metacalibration method is applied to the entire catalog, providing us with global values for the shear calibration quantities. We reproduce those results for CFIS-P3 here. The response matrices for CFIS-P3 are The additive shear bias components are The errors are computed via jackknife resampling of all galaxies.

Local calibration
The multiplicative shear bias depends on quantities that vary spatially, such as the PSF properties (Paulin-Henriksson, S. et al. 2008), galaxy size (for example Spindler 2018; Kuchner et al. 2017) and magnitude (Miller et al. 2013), or the local galaxy density (Hoekstra et al. 2017). These spatial variations may be correlated with shear: in some cases, they both vary with the LSS environment, such as galaxy density and shear. In other cases, residual errors may create cross-correlations, for example, an imperfect PSF calibration influences both calibration and shear.
The understanding and mitigation of spatial patterns in shear calibration is an active field of research for future weak-lensing surveys (Kitching et al. 2021;Cragg et al. 2022). Here, we investigate local shear calibration, the dependence of calibration on observed quantities, and the impact on cosmological parameters.
In the context of weak-lensing mass maps as tracers of the LSS, Van Waerbeke et al. (2013) carry out a local calibration in pixels of 1 arcmin size. Their multiplicative shear bias is computed as a smooth, two-parameter fitting function from image simulations. This is in contrast to our case of metacalibration, which provides very noisy calibration estimates for individual galaxies. These estimates need to be averaged over substantially larger areas to reduce the uncertainty such that the corresponding response matrices are numerically stable enough for inversion.
In this work, we employed a series of square patches with a size of d s = 4, 2, 1, 0.5, and 0.25 degrees over which the shear calibration estimates are averaged in turn. To easily divide the observed sky area into an integer number of square-sized subpatches, we first projected the data into a Cartesian plane. Next, we extended this area via zero padding to size N x × N y deg 2 , such that N x and N y are multiples of the largest sub-patch size d s, max = 4 deg.
We computed local versions of R shear (7) and R selection (8) by carrying out the averages per sub-patch. If the number of galaxies in a sub-patch was smaller than a threshold n gal, 0 , we replaced the local response matrices in that sub-patch by their global mean, to avoid numerical instabilities. We chose n gal, 0 =n gal /2. This occurs mainly at the edges of the field where sub-patches overlap with the area of zero padding. We calibrated the estimated shear of each galaxy in a given sub-patch by the total local response matrix, in analogy to (9). Similarly, we computed the local shear bias c as the average per sub-patch.
We computed the error of the local response matrices by creating jackknife resamples from the smallest sub-patch size d s, min = 0.25 deg. Each sub-patch of size d s > d s, min was thus split into (d s /d s, min ) 2 jackknife subsamples; for the smallest subpatch size, we could not compute the error. The additive bias is easily obtained locally by computing jackknife errors over galaxies per sub-patch.

Variation of the shear and selection response
We show the spatial variation of R shear and R selection in subpatches of different size in Fig. 3. From top to bottom the shear and selection responses are computed on sub-patches with increasing size of d s = 0.25, 0.5, 1, 2 and 4 degrees. The smaller the sub-patch size, the more the response shows variations corresponding to the local spatial environment. The larger the subpatch size, the more the calibration approximates the global cal-ibration. The errors of the shear response are small, even in the case of the 0.5 degree sub-patch size. Concerning the selection matrix, we can see that on 0.5 deg 2 , the errors have the same order of magnitude as R selection . On 1 deg 2 , the errors are smaller than R selection indicating that the computation of the local selection calibration is precise when done on 1 deg 2 or more. To emphasize these conclusions, the 1D distributions of R shear and R selection are plotted in Fig. A.1.

Variation of the additive bias
The additive bias and its standard deviation are shown in Fig. 4. The standard deviation is as large as the value itself for local calibration on small scales. Except for the largest sub-patch sizes, the values of c 1 and c 2 per pixel are consistent with zero. This reflects the large fluctuations and randomness of the additive bias on small scales. To emphasize these conclusions, the 1D distributions of c 1 and c 2 are plotted in Fig. A.2.

Relative errors
We computed the errors of the response matrices at different subpatch sizes, averaged over all sub-patches (blue lines), compared to the relative errors of the global calibration (black lines). The results are shown in Fig. 5.
The relative error of the shear matrix is below 1% in all cases but always above the relative error of the global calibration, meaning that there are fluctuations at all the studied scales, but these are small. The fractional error of the selection matrix is higher. It reaches 17% at small sub-patch sizes but decreases up to 3% at 4 deg 2 . For the shear and selection matrix, we see that when we go on a smaller calibration size, the relative errors asymptotically approach the one of the global calibration.

Parameter correlation matrix
The calculation of the local shear bias allows us to further explore possible origins of shear bias, and their influence on other quantities obtained from the galaxy sample. To that end, we computed the correlation matrix between different quantities, for which we used patches of size d s = 1 deg. We also considered correlations between quantities other than shear bias. The correlation matrix is calculated using the pandas DataFrame.corr function 12 ; it uses the Pearson method, which computes the standard correlation coefficient. The correlation matrix between the different quantities listed in Table 3 is shown in Fig. 6. To better see the correlation terms by terms, another correlation matrix is shown in Fig. B.1.
First, we see that the E-and B-mode smoothed convergence, κ E,sm and κ B,sm , are not strongly correlated. This indicates that there is no significant systematic effect that mixes both modes. Neither mode is strongly correlated to other quantities. Also, the number of peaks n peak is not correlated to observational effects.
Further, there is no visible correlation between the PSF size and the additive bias. This is evidence for the correct estimation of PSF size since a bias in the PSF size typically leads to an additive shear bias. We also note the absence of a correlation between the PSF ellipticity and the shear selection matrix. This gives us confidence that the PSF deconvolution in the metacalibration process works well. There is, however, a 20 -30% negative correlation between the selection response elements and the Table 3. Symbols and description of quantities used in the crosscorrelation matrix (Fig. 6).

Symbol
Description κ E,sm , κ B,sm Smoothed E-and B-mode convergence, respectively n gal number of galaxies per pixel n peak number of peaks |e psf | modulus of the point spread function T gal galaxy size n epoch number of single-exposure epochs r r-band galaxy magnitude SNR gal galaxy signal-to-noise ratio T psf psf size w galaxy weight R shear diag , R shear off-diag Average of the diagonal and off-diagonal shear response matrix element, respectively, Eq.
Average of the diagonal and off-diagonal selection response matrix element, respectively, Eq. (8) c Average of the additive bias element size of the PSF. A plausible explanation is that in areas with larger seeing, fewer small and faint galaxies make it into the weak-lensing sample. This leads to a stronger (more negative) selection bias, which is reflected in the anticorrelation. A negative correlation can be seen between the number of epochs, n epoch , and the PSF ellipticity modulus |e PSF |. Since the former is the mean over all contributing epochs, a reason for this anticorrelation might be the PSF ellipticity gets more circular when averaged over more independent observations. There is a negative correlation of the r-band magnitude with galaxy weight w, S/N, and size T gal , as well as a positive correlation with n epoch . This reflects the expectation that faint galaxies have a lower weight and S/N, and are easier observed with more exposures.

Impact of uncertainty and systematic effects on cosmological parameters
This section explores different systematic effects and uncertainties, and their impact on cosmological parameter constraints. We quantify the potential biases on cosmology from the spatially varying shear calibration (Sect. 6.1), the redshift uncertainty (Sect. 6.3), a residual shear bias (Sect. 6.2), baryonic feedback (Sect. 6.4), intrinsic galaxy alignment and cluster member dilution (Sect. 6.5. All these effects are studied jointly in Sect. 6.6. For each effect, we will show constraints on cosmological parameters as compared to those obtained with the "ideal" case. This case, which is represented in blue in all the following figures represents the constraints obtained when we use the data calibrated globally without residual shear bias, a mean redshift ofz = 0.65, without baryonic correction, cluster member dilution or intrinsic alignment. This assumes an unrealistic bestcase scenario of vanishing spatial variation of shear calibration and residual shear bias, and no biases from baryons, intrinsic alignment or cluster members. When the parameters are wellconstrained, we specify the shift compared to the ideal case, which is the case for Ω m . The parameters M ν and the lower bound for A s are in most cases not well constrained within the prior.

Local calibration
To study the impact on cosmology of spatially varying shear biases, we performed a local shear calibration using different scales, following the method explained in Sect. 5.3 for the multiplicative shear bias. For the additive bias, we used the value obtained with the global calibration because the local calibration value is very noisy. The result is shown in Fig. 7. We see that a local calibration will always shift Ω m to a lower value except for the calibration on 0.5 deg 2 . These variations are all within the statistical error bars. There is no systematic variation that becomes evident when going to smaller calibration sizes. For A s and M ν , the trend is not clear. Ω m shifts by −0.015, which corresponds to −0.3σ (i.e., 0.3 times the statistical uncertainty) when the calibration goes from global to local on 4 deg 2 for example. From this work, it is not clear which calibration size is the best one because, depending on the size of the calibration, we will capture more or less the local effect. Future work has to be done to determine the criteria to choose the size of the calibration. In this work, when we need a local calibration, we use the one done on 1 deg 2 because the relative error of the selection matrix (which is the limiting one) is below 10%.

Residual multiplicative shear bias
The shear bias m computed earlier with metacalibration is not perfect. A residual bias ∆m = m metacal − m true remains and gen-erally has to be quantified with image simulations. We used the value ∆m = 0.007 based on the results from Guinot et al. (2022), who found ∆m = 0.007 ± 0.03, estimated using CFISlike image simulations of isolated galaxies and ignoring the effect of blending. In comparison, other surveys also find residual biases at around the per cent level or below; for example, MacCrann et al. (2021) and Gatti et al. (2021) state a value of ∆m = −0.0208 ± 0.0012 for DES-Y3.
We modeled the effect of the residual multiplicative shear bias by adding the residual bias ∆m to the response matrix in the local case. This corresponds to a conservative, worst-case scenario where the residual bias is constant in space. Thus, Eq. (9) is modified to where I is the identity matrix. The impact of this bias on the cosmological constraints is shown in Fig. 8. In blue, the calibration is global whereas it is done on 1 deg 2 for the green and red cases. Moreover, the red case includes the multiplicative shear bias. Including the residual shear bias changes the values of the parameters in the same direction as the local calibration compared to the global calibration. Ω m shifts by −0.024, which corresponds to −0.5σ with respect to the global calibration. A positive ∆m effectively reduces the estimated shear, resulting on average in a smaller number of peaks, and thus smaller clustering parameters. The 1D marginalized contours for A s shift to the right, contrary  to what we would expect. However, the joint Ω m -A s 2D contours clearly reflect the smaller calibrated shear and shift toward smaller clustering amplitudes.

Redshift uncertainty
As discussed in Sect. 2.3,z = 0.65 is our best estimate of the redshift of the data observed from CFIS. This estimate has however a level of uncertainty. To quantify the impact of this uncertainty on our results, we compared cosmological constraints using the simulated convergence maps interpolated toz = 0.65 and z = 0.68, which corresponds to the two mean redshift estimates discussed in Sect. 2.3. The result is shown in Fig. 9. Fitting the data with the model at a higher mean redshift only slightly shifts the posterior: Ω m shifts by +0.001, which corresponds to +0.02σ when the redshift goes from 0.65 to 0.68. These shifts are well within the statistical uncertainties of the two parameters. Such a shift is, however, expected to have a significant influence on analyses using a larger survey area and/or redshift tomography. We also performed a test using only simulations, where the data vector was the mean of the fiducial simulation at z = 0.65 or at z = 0.68. The result is presented in Appendix C. We note, however, that when using simulations only, the shift appears to be A&A proofs: manuscript no. aanda  Table 3. The colors indicate the amplitude of the correlation, ranging between −1 and 1. All quantities are mean values.
larger, this may indicate some residual systematics that impacts the results on the data shown in this section.

Baryonic feedback
The impact of baryons on the LSS is important for cosmological analyses with weak-lensing peak counts (Osato et al. 2015;Harnois-Déraps et al. 2021a;Coulton et al. 2020). The redistribution of matter due to baryonic processes tends to reduce the number of high S/N peaks and augment that of smaller S/N values. We used the results from Coulton et al. (2020) who model the fractional difference of the number of peaks, ∆N peaks /N peaks , for three different baryon physics scenarios based on the BAHAMAS hydrodynamical simulations (McCarthy et al. 2016) compared to MassiveNuS dark-matter-only simulations (Liu et al. 2018). This assumes that baryonic processes are independent of the underlying cosmology. The three scenarios are denoted as LowAGN, Fiducial, and HighAGN, where the amount of bary- Fig. 7. 1D and 2D marginal posteriors for CFIS-P3 using different calibration sizes in the metacalibration. The 2D inner and outer contours show the 68% and 95.5% credible region, respectively. The case when the calibration is done globally (blue) is compared to the case where the calibration is done on 0.5 (black), 1 (green), 2 (red), and 4 deg 2 (pink). Ω m shifts by −0.015, which corresponds to −0.3σ when the calibration goes from global to local on 4 deg 2 . Fig. 8. 1D and 2D marginal posteriors for CFIS-P3, with the residual multiplicative shear bias added. The 2D inner and outer contours show the 68% and 95.5% credible region, respectively. Ω m shifts by −0.01, which corresponds to −0.2σ when the calibration goes from global (blue contour) to local on 1 deg 2 (green contour). When ∆m is added (red contour), Ω m shifts by −0.024, which corresponds to −0.5σ with respect to the global calibration.  onic feedback increases in that order. The fiducial baryonic correction comes from the BAHAMAS simulations with a feedback model designed to best match the observations. The LowAGN and HighAGN corrections are simulations where the active galactic nucleus (AGN) heating is lowered or raised by 0.2 dex, respectively, and thus the simulations skirt the lower and upper bounds of the observed gas fraction. We modified the model predictions of peak counts by multiplying the number of peaks from the MassiveNuS (dark-matteronly) simulations by N peaks /N peaks,DM , shown in Fig. 10. Fig. 11. 1D and 2D marginal posteriors for CFIS-P3, with the peak count predictions corrected with different baryonic feedback. The 2D inner and outer contours show the 68% and 95.5% credible region, respectively. The case without baryonic correction (blue) is compared to the three baryonic scenarios: LowAGN (green), Fiducial (red), and HighAGN (pink). Ω m shifts by +0.027, which corresponds to +0.5σ when the fiducial goes from dark-matter-only to LowAGN correction.
We can see how the constraints evolve when we use real CFIS-P3 data and corrected simulations. To see the effect of the baryonic correction only, we use the global calibration, simulations at z = 0.65, no intrinsic alignment or cluster member dilution correction. We do not know how the baryonic feedback influences the observed data thus we correct the MassiveNuS simulations with the three different baryonic feedback corrections. In Fig. 11 we show the results of these corrections (LowAGN in green, Fiducial in red, and HighAGN in pink) compare to the case of the dark-matter-only simulation (blue).
Baryonic feedback has a significant impact on the best-fit parameters. Ω m and A s are shifted to higher values with increasing baryonic modifications. For example, Ω m shifts by +0.027, which corresponds to +0.5σ when the fiducial goes from darkmatter-only to LowAGN correction. Li et al. (2019) found that it is the high S/N peaks that dominate the constraints and in Fig. 10 we see that the number of high peaks is suppressed by the baryonic feedback. This results in a shift toward higher values in Ω m , compensating for the suppression of peaks for S/N>3. We also see that the strength of the feedback is not as important as just taking it into account.

Further systematic effects
At least two further effects can affect constraints: intrinsic alignment and cluster member dilution. As we briefly explain below, both effects can be suppressed, in first approximation, following Harnois-Déraps et al. (2021b), with a cut in S/N.

Intrinsic alignment
Due to the radial alignment of satellite galaxies within darkmatter halos, galaxies are not randomly oriented. Thus, their shape and alignments are affected by their environment and tidal fields. The intrinsic alignment has two components: the intrinsicintrinsic correlations caused by the alignment of galaxies that are physically linked together and the gravitational-intrinsic correlations, which is the alignment of halo galaxies.

Cluster member dilution
The source density is not homogeneous and increases around foreground clusters. Around a cluster at a given redshift, there are more galaxies but, as we do not know their redshift, they are included in the signal but may not lensed. This effect leads to a coupling between the peak positions and the amplitude of the measured shear relative to the expected shear (Kacprzak et al. 2016). Moreover, these regions of clusters have a larger blending rate, and thus galaxies behind clusters are more likely to be missed. These effects can result in a miscalibration between data and simulations and lead to a reduction in the mean shear signal (as we count more galaxies for the signal). An analysis of cluster member dilution in Subaru HSC weak-lensing mass maps is available, for example, in Oguri et al. (2021). -Déraps et al. (2021b) find that the effect of intrinsic alignment is small for peaks with S/N < 3 and that the effect of the cluster member dilution is small for S/N < 4. The local calibration might capture part of the cluster dilution effect, but as we have no way of knowing it without simulations, we chose to cut the S/N range to be conservative. Thus, as a first approximation of how these effects impact the cosmological parameters, we computed the constraints using only the range for peak counts of S/N < 3. This selection should minimize the impact of both systematic contaminations. We tested both cases (−2 < S/N < 3 and −2 < S/N < 6) on mock peaks from MassiveNuS and find consistent results. The comparison of the resulting constraints obtained with S/N < 3 (red) or S/N < 6 (blue) is shown in Fig. 12. Constraints are slightly larger when we cut the S/N range, which is expected. For both M ν and A s parameters, the change is small.

Harnois
The high-S/N peaks are affected by intrinsic alignment and cluster member dilution, which leads to a reduction in peak S/N (Harnois-Déraps et al. 2021b). When using the full S/N range and not accounting for those two effects, the reduction in peak S/N results in a lower inferred clustering amplitude and a lower Ω m . This can be seen in the Ω m -A s 2D posterior distribution. We also perform a test using simulations only, with the range for peak counts of S/N < 3. As in previous tests, the data vector is the mean of the fiducial simulation at z = 0.65 and the result is presented in Appendix C.

Parameter constraints that combine all systematic effects
Our conservative model, which is the most suitable to represent the data, combines different mitigation schemes, and uses the following parameters and settings: The mean redshift is set tō z = 0.65. We calibrate the shear locally at a scale of 1 deg 2 . This accounts for spatial variations of shear bias and is the smallest . Ω m shifts by +0.027, which corresponds to +0.5σ when we cut the high S/N peaks.
scale for which the estimate of the selection is not dominated by noise. The residual shear bias is set to the value estimated from image simulations, ∆m = 0.007. Baryonic feedback is accounted for by using the "fiducial baryon" case of Coulton et al. (2020). The data vector is composed of peak counts with −2 < S/N < 3 to minimize the systematic errors from intrinsic alignment and cluster member dilution.
In Fig. 13 we show constraints using this conservative model, and compare those to the ones obtained under the ideal (optimistic) case described at the beginning of this section. Ω m shifts by +0.008, which corresponds to +0.2σ when we go from the ideal to the conservative model. As expected, the conservative model results in wider constraints. In Fig. 14, we can see the marginalized distributions for the 68% confidence interval of the different parameters depending on the different cases we have tested.

Conclusions
This study highlights the importance of properly accounting for systematics, such as the local calibration, residual multiplicative shear bias, intrinsic alignment and cluster member dilution, redshift uncertainties, and baryonic corrections, in the context of a peak-count cosmological analysis with CFIS data. We performed a likelihood analysis on the sum of neutrino masses, M ν , the matter density parameter, Ω m , and the amplitude of the primordial power spectrum, A s . For this purpose, we used the CFIS-P3 catalog and the MassiveNuS N-body simulations because they have the advantage of providing us with a large set of cosmological models and different tomographic bins. First, we obtained constraints with simulations only to validate the methodology used to model the redshift uncertainty, baryonic feedback, in- Fig. 13. 1D and 2D marginal posteriors for CFIS-P3 using the ideal or conservative model. The 2D inner and outer contours show the 68% and 95.5% credible region, respectively. The ideal one (blue) is obtained with global calibration, no baryonic correction, no intrinsic alignment, no boost factor, and on −2 < S/N < 6. The conservative model (red) is obtained with local calibration on 1 deg 2 , with residual bias, modeled under the fiducial baryonic correction, and on −2 < S/N < 3. Both models are computed atz = 0.65. The 2D contours show the 95.5%, and the 1D filled area corresponds to the constraints within the 1 sigma confidence level. Ω m shifts by +0.008, which corresponds to +0.2σ between the reference (blue) and the fiducial model (red). trinsic alignment, and cluster member dilution. We then used the CFIS-P3 data to quantify the effect of systematics on cosmological constraints. We remark that Ω m is the parameter that is constrained the most and is the most impacted by the systematics. As a summary plot, Fig. 14, which shows the marginalized distributions for the 68% confidence interval of the parameters per studied case. On the one hand, the baryonic corrections, the intrinsic alignment, and the cluster member dilution shift this parameter to higher values. On the other hand, a local calibration and a multiplicative shear bias ∆m = 0.007 shift Ω m to smaller values. Concerning A s , the baryonic correction and the addition of ∆m shifts this parameter to higher values, whereas cutting the high S/N peaks shifts A s to a lower value. For M ν , the parameter is not constrained enough to allow for any conclusions to be drawn. More specifically, concerning the shear calibration, it is probable that a local one is preferred over a global one, as long as the error and standard deviation are small, because it accounts for the local effects of the catalog. Nevertheless, some work has to be done to determine which calibration size is the better one. We notice that using a calibration on 1 deg 2 with a value of the residual multiplicative shear bias ∆m = 0.007 shifts Ω m by −0.024 (−0.5σ) compared to the case in which a global calibration is applied. Choosing a reasonable calibration size and having a robust estimate of ∆m is important for having a better estimate of the constraints. We note that the residual bias used here computed from image simulations neglected effects such as galaxy image blending. A more realistic estimate of ∆m is necessary for future analyses of CFIS weak-lensing data.
Concerning the uncertainty on the mean redshift estimate, for a bias of ∆z = 0.03, as considered in this study, the impact on the constraints is small enough to not be considered: the shift in Ω m is only +0.001 (+0.02σ). Nevertheless, we are aware that further work is needed to have a more complete description of the redshift distribution and obtain more accurate constraints. Using simulations only, the shift in Ω m is larger, 0.7σ. This indicates that other effects in the data might lead to an underestimation of the actual shift.
To account for baryonic corrections, we considered three different flavors of baryonic feedback, labeled in this study as HighAGN, Fiducial, and LowAGN correction. We conclude that the specific flavor of baryonic feedback is less important than the difference between baryonic feedback and pure dark-matter models. Applying the LowAGN baryonic correction shifts Ω m by +0.027 (+0.5σ), showing the importance of modeling the baryonic feedback. Using a prediction of peak counts based on hydrodynamical simulations to model the baryonic feedback can help in getting more accurate constraints.
Then, we see that when we minimize the effect of the cluster member dilution and intrinsic alignment by only considering peaks with S/N < 3, the posterior on Ω m is subject to an offset of 0.027 (0.5σ) toward higher values. Using simulations that model these effects will significantly improve the constraints since a cut in the S/N range will no longer be necessary.
Finally, we computed a conservative model where we used the parameters and settings that are most representative of the the data. This gives larger constraints on M ν , Ω m , and A s .
We noticed that the value of Ω m in particular can shift a lot due to different systematics. Our aim in this paper is not to estimate the final reference cosmological parameters, but rather to investigate, for the first time with UNIONS and peak counts, the impact of some of the different systematics at play for this survey. Among all the systematics considered in this paper, we have shown that the one with the highest impact on Ω m is related to how baryonic corrections are implemented (Fig. 11); the choice of cut on the S/N (Fig. 12) also causes a substantial shift. Other systematics not considered here may further shift the final parameters. It is necessary to further pursue this effort to include other systematic effects, as well as to extend the way we incorporate the baryonic feedback in the analysis and the choice of the calibration size. The constraints will be more accurate when simulations are able to also model intrinsic alignment and cluster member dilution; hydrodynamical simulations that include baryonic feedback could further improve the robustness of the results. In the future, having access to simulations built for UNIONS will allow us to make more precise investigations of various sources of bias. We also expect constraints to become more reliable with future (larger) data catalogs, such as the full UNIONS data set, for which the current pipeline will provide a starting point.
To have a better visualization of the variation of the response matrix and the additive bias, we compute the histogram of their distribution in Figs. A.1 and A.2, respectively. In every histogram, we add as a black line the value found by the global calibration. From those histograms, we can clearly see that when the pixels are smaller the dispersion is larger, which is because the number of galaxies in one patch is lower and therefore there is more noise. In all cases, we see that the values are spread around the mean one.
We confirm that the values are spread around the mean one, with larger dispersion at smaller scale of calibration.

Appendix B: Correlation matrix
The correlation matrix presented in Fig. B.1 is presented to see the correlation between the individual elements of the matrices.
We can observe a strong correlation ( 60 -80%) between R shear 11 and R shear 22 , which is expected because if the shear comports a bias, both component will be impacted. A correlation of 20 -40% between R selection 11 and R selection 22 is seen for the same reason as the correlation between R shear 11 and R shear 22 . The R selection 11 and R selection 22 are anticorrelated ( 20 -40%) with the size of the psf T psf because if the psf is larger, we will miss more objects. Thus, the correlation effect is stronger. 1D and 2D marginal posteriors using as data the mean fiducial from the simulations at different redshifts. The 2D inner and outer contours show the 68% and 95.5% credible region, respectively. Ω m shifts by +0.026, which corresponds to +0.7σ when the redshift goes from 0.65 (blue) to 0.68 (red).

Appendix C.2: Intrinsic alignment and cluster member dilution
The method used here to decrease the impact of intrinsic alignment and cluster member dilution is to cut the high S/N range. In Fig. C.2, we show the constraints when we use a S/N range between −2 and 6 (blue) or a reduced S/N range between −2 and 3 (red). When the reduced range is used, the constraints are larger, which is expected because there is less information. is compared to the case where we use −2 < S/N < 3 to minimize the effect of residual bias (red). The constraints are larger when we cut the high S/N peaks.