<figure><img src="image_25.png" /><figcaption>**Fig. 23.** Stray light map and survey fields. The coloured squares show the stray light level in VIS dark exposures as a function of spacecraft orientation angles. The log-scaled greyscale map shows the density of fields in the latest survey configuration including calibrations. The survey minimises stray light over the EWS and EDS, with the majority of the observations to be taken at $AA = -4.5$. The CPC fields are NISP-specific and include higher SAA positions (above the jagged black line); while NISP is not affected by stray light, parallel VIS observations must still be taken.</figcaption></figure>

The halo catalogue and the set of 2D dark-matter maps are the main inputs for the Flagship mock galaxy catalogue. A detailed description of the catalogue production is given in Euclid Collaboration: Castander et al. (2025), which we summarise here. Galaxies were generated following a combination of halo-occupation distribution (HOD) and abundance matching (AM) techniques. Following the HOD prescription, halos were populated with central and satellite galaxies. Each halo contains a central galaxy and a number of satellites that depends on the halo mass. The halo occupation was chosen to reproduce observational constraints of galaxy clustering in the local Universe (Zehavi et al. 2011).

The luminosities of the central galaxies were assigned by performing abundance matching between the halo mass function of the simulation halo catalogue and the galaxy luminosity function (LF). We took as the reference luminosity function a parameterised prescription that tries to fit the observed LF throughout the redshift range of the simulation. We then applied a 15% scatter to the resulting luminosities. The application of scatter is necessary to reproduce the galaxy clustering dependence on luminosity. The satellite luminosities were assigned assuming a universal Schechter LF for satellites in which the characteristic luminosity depends on the central luminosity in a way that ensures that the global luminosity function agrees with observations. Galaxies were split into three colour types, namely red, green and blue, and the central and satellite galaxies in each group were distributed to match the observed clustering as a function of colour by Zehavi et al. (2011). The radial positions of the satellites within their halos follow Navarro–Frenk–White

<figure><img src="image_26.png" /><figcaption>**Fig. 24.** Reach in ecliptic longitude around transit for the leading side of the survey, for SAA $\in$ [87°, 104°]. Shown in red is the reach for the originally planned symmetric AA range. To minimise the stray light in VIS, the range was shifted to $AA \in$ [-8.4°, -3.0°] with much reduced visibility (grey) that would not permit the completion of the survey. By allowing the fields to rotate by up to 3° with respect to the tessellation, a much larger area of the sky becomes accessible (blue). For observations in the trailing side, the areas must be rotated by 180° around the origin.</figcaption></figure>

(NFW) profiles (Navarro et al. 1997) with each colour type having its own concentration: green and blue galaxies are distributed using a concentration that is respectively one-half and one-quarter times that of the red galaxies, consistent what Collister & Lahav (2005) found for blue versus red galaxies in the nearby Universe. This colour segregation is assumed to hold at all redshifts (see Sect. 10.6 for a discussion of this assumption).

To determine the simulated redshifts, we assumed that the central galaxy is at rest in its halo. The satellite velocities were drawn using formula derived by solving the Jeans equation of local dynamical equilibrium for each type, assuming velocity anisotropy profiles consistent with what was measured from the kinematics of low-redshift regular clusters (Mamon et al. 2019), with green and especially blue galaxies having more radial orbits around their host haloes. Redshifts of the galaxies were then obtained by projection of these 3D velocities along the LoS.

SEDs were assigned to each galaxy with a procedure that aims to mimic the observed colour distributions as a function of redshift. The resulting SEDs are a linear combination of the ones presented in Ilbert et al. (2009). The stellar masses were computed from the galaxy luminosities and the mass-to-light ratios of the SEDs. The star-formation rates were computed from the ultraviolet luminosity of the SED. The luminosity of the $H\alpha$ line was computed from the star-formation rate following the Kennicutt recipe (Kennicutt 1998a). The galaxy clustering measurements of Euclid rely on the detection of this line in the galaxy spectra (see Sect. 2.1). We therefore want to simulate their distribution to our best current knowledge. Consequently, we then refined the $H\alpha$ luminosities to match the models of Pozzetti et al. (2016), using abundance-matching techniques. The luminosities of the other main emission lines were assigned
<figure><img src="image_27.png" /><figcaption>**Fig. 25.** EWS coverage and colour-coded yearly progress in an all-sky Mollweide projection. The blue borders enclose the 16 000 deg² RoI that contains the 13 416 deg² observed sky of the EWS. Small dark regions within the EWS are masks for stars brighter than 4 AB mag.</figcaption></figure>

using observed relations, taking the Hα line as reference. The shapes and sizes of the galaxies were assigned following relations based on HST observations (Miller et al. 2013; Dimauro et al. 2018). The observed fluxes of each galaxy were computed by integrating the SED with the filter transmission for the surveys that are expected to be used to obtain photometric redshifts.

As for galaxy lensing properties, we followed the ‘onion universe’ approach presented in Fosalba et al. (2008) and Fosalba et al. (2015), to compute all-sky weak lensing observables (convergence, shear, and deflection) within the Born approximation. The latter agrees within 0.1% out to a multipole $\ell = 10^4$ with the much more complex and central processing unit (CPU) time-consuming ray-tracing technique (Hilbert et al. 2020). Figure 26 shows the convergence (colour-coded) and the overlaid shear field (sticks indicating the amplitude and direction of the shear) for a source redshift $z = 1$ in a patch of approximately 50 deg$^2$ of the EFS.

## 6.2. Simulated Euclid data

Pixel image simulations serve as test data for the development of the data-processing pipeline and computing infrastructure, and to validate the stringent requirements regarding performance and data quality (e.g. Euclid Collaboration: Gabarra et al. 2023). In this respect, several end-to-end science performance verification (SPV) tests have been performed throughout the mission preparation, which allowed for the reproduction of certain instrumental issues to assess their impact and guide decision making (see Sect. 8.1 for more details). Image simulations have also been essential to develop alternative solutions to critical problems that were discovered during on-ground tests, such as the non-conformity of one of the three red grisms for the NISP instrument (Euclid Collaboration: Scaramella et al. 2022; Euclid Collaboration: Jahnke et al. 2025); or to evaluate the impact of unexpected features discovered in-flight, such as the

contamination of VIS images by X-ray photons during solar flares (see Sect. 5.3 and Fig. 21). Finally, image simulations are necessary to calibrate intrinsic biases related to the methods used to measure the shapes of the galaxies in the cosmic shear analysis. A large volume of very accurate and representative simulations is required to determine these biases to the required level of precision (e.g. Hoekstra et al. 2017).

To ensure a common input for the instruments, the starting point is a ‘true universe’ catalogue that contains all input sources and their corresponding parameters, spectra, and shapes. These parameters are based on the output from the EFS (Sect. 6.1), which provides both spectra (continuum and emission lines) and morphological parameters (bulge and disc models). The stars are simulated using a hybrid catalogue, using actual stars from Gaia DR3 (Gaia Collaboration 2023b) at the bright end ($G < 18.5$), merged with the fainter end of deep stellar population simulations using the Besançon model (Robin et al. 2003; Czekaj et al. 2014; Lagarde et al. 2021), including binaries, down to $I_E = 26$. All objects are simulated using common coherent libraries to ensure consistent results between the simulated VIS imaging, NISP photometry and spectroscopy, and ground-based imaging. In addition to individual sources, the simulations include zodiacal light, diffuse stray light from stars beyond the simulated pointing, and thermal irradiance caused by the heat of the various elements of the telescope and instruments. The instrument models and the reference survey characteristics are drawn from the central mission database. This database (versioned and controlled by dedicated change control boards) provides the instrument simulators with all the parameters necessary to simulate the numerous instrumental features of each simulation channel.

To create simulated VIS observations and the associated calibration frames, we have developed the Euclid VIS simulator (ELViS). It includes an accurate emulation of the optical response, which is based on the complex PSF modelling tools that are described in Sect. 7.6.4. Although challenging, ELViS
<figure><img src="image_28.png" /><figcaption>**Fig. 26.** The image on the left shows the lensing convergence for sources with $z_s = 1$ for a simulated patch of sky covering 50 deg$^2$. A zoom-in of the central square degree is shown on the right, with the sticks indicating the direction and amplitude of the corresponding shear. The colour bar of the convergence field displays values within the range ±3σ, where σ is the rms value of the full-sky map. The stick at the bottom of the zoom-in image shows a reference amplitude for the shear sticks overlaid on that area of the mass map.</figcaption></figure>

can capture the complexity of the chromatic PSF and its spatial variations. The sources are projected on the simulated mosaic of 36 CCDs, and an extensive list of instrumental signatures can be included, such as bias, pixel response non-uniformity, cosmetic defects, saturation, bleeding, shutter movement, and ghosts from the dichroic plate. Particularly important for Euclid are simulations of the imperfect charge transfer during the readout. The simulated readout electronics include a nonlinear response, saturation of the analogue-to-digital converter, and electric cross-talk, with parameters determined during the on-ground test campaigns.

To capture the dual use of the NISP instrument, two separate codes are used, Imagen and TIPS, to simulate the photometric and spectroscopic channels, respectively. As is done for VIS, both science and calibration images are produced. As the photometric and spectroscopic channel share the same optical path and detector units, the NISP-P simulator, Imagen and NISP-S simulator, TIPS employ common background, PSF, optics, detector, and electronics models. Consistency across the two channels is particularly important because the spectroscopic analysis relies on measurements from the imaging data. The simulations start from the same list of astronomical sources that is used by ELViS. The simulated NISP PSFs contain the same types of effects that have been included for VIS instrument, but with a simplified module to capture the variation of the PSF (given the less strict requirements, we opted to use a tabulated PSF as opposed to recomputing it for each object). Each different multiple accumulated sampling readout (MACC mode) is simulated for photometric and spectrometric images. At the detector level, variable QE is simulated using QE estimates obtained from

on-ground tests, while the pixel-response non-uniformity (PRNU) is taken from the on-ground flatfield data. To capture biases introduced during readout, nonlinearity and gain are simulated together with the MACC readout modes. On the spectroscopic side (TIPS simulator), the slitless light dispersion is handled by aXeSIM<fnref n="16" /> (Kümmel et al. 2009), developed by the Space Telescope Science Institute. The simulation of the slitless spectra requires a complete characterisation of the trace dispersion, with sensitivity, diffraction coefficients, grism tilts, vignetting, and PSFs at each dispersion order (expressed as a Taylor series expansion). The dispersion order of interest is the first order, where the main spectrum can be recovered. However, the zeroth and second orders of dispersion need to be simulated as well. The readout electronics are shared with the photometric channel and the simulations are therefore also similar (except for the MACC modes, which are photometry and spectroscopy specific).

We also produce pixel simulations of the external ground-based surveys that are used to characterise the SEDs of the stars and galaxies and to determine photometric redshifts. Here, it is important to capture the key characteristics of the ground-based data, summarised in Sect. 4.4, while providing realistic images that can be used to examine the performance of the various processing steps. With the exception of LSST, the surveys have already collected data, so that realistic prescriptions for their main features can be readily implemented. For instance, the simulations include realistic values for the background levels, zero points, filter transmissions, bias levels, flat-field characteristics, detector defects, observed cosmic rays, and PSFs for each of the

<footnote n="16">¹⁶ http://axe-info.stsci.edu/axesim/</footnote>
simulated surveys. Importantly, the simulated ground-based data use the same input catalogues as the simulations of the *Euclid* instruments, so that all the detection and cataloguing steps can be tested, as well as the determination of photometric redshifts. For a more detailed description of the pixel simulations, the reader is referred to Euclid Collaboration: Serrano et al. (2024).

## 7. Survey data products

*Euclid* provides high-quality optical and NIR imaging, as well as slitless NIR spectroscopy over a large fraction of the extragalactic sky. The processing of these data and supporting ground-based observations is performed by the SGS. In this section, we present the most salient parts of the pipeline. The key steps in the processing of the *Euclid* optical and NIR imaging data are highlighted in Sects. 7.1 and 7.2, respectively, while the processing of the supporting ground-based imaging is summarised in Sect. 7.3. These data form the basis for the creation of the object catalogues, described in Sect. 7.4. The processing of the slitless spectroscopy and the subsequent redshift determination are discussed in Sect. 7.5.

In Sect. 7.6 we describe how these results are used to derive the key ingredients for the 3×2pt analysis. The science-ready data products, summarised in Sect. 7.7, include the various two-point statistics, as well as weak lensing convergence maps and catalogues of clusters of galaxies. Although many of these data products are excellent starting points for myriad scientific investigations, several high-level data products will be released as well, with a particular focus on the interpretation of the primary probes. These data, as well as a range of calibrated data products, will be released to the scientific community in a number of data releases (Sect. 7.8).

### 7.1. Processing of VIS imaging data

The high resolution imaging data provided by the VIS instrument is the starting point for the weak lensing measurements. To reach the main objectives of *Euclid*, the shapes of about 1.5 billion galaxies need to be measured with unprecedented accuracy. To this end, a range of instrumental effects need to be carefully accounted for (Massey et al. 2013; Cropper et al. 2013; Euclid Collaboration: Paykari et al. 2020). This involves fully characterising the performance of the instrument through the processing and analysis of an extensive amount of calibration data. The subsequent processing of the raw VIS data relies on these calibration products to correct all relevant instrumental effects to enable robust shape measurements.

Before launch, each pipeline processing element was tested and validated using the simulations described in Sect. 6.2. These simulations enabled an assessment of how well instrument models (based on pre-launch observations) could be derived using realistic observing sequences. They also allowed an evaluation of how well measurements (photometry, astrometry, and object shape) could be made after the application of these models. These findings are being updated using the PV observations and in-flight calibrations.

At its most basic level, the processing of the VIS data can be divided into three categories: the calibration pipeline that generates or updates the VIS calibration models, which by definition do not alter the input pixels; the science pipeline that alters the input pixels by applying the calibration models; and the validation pipeline designed to assess the performance of the data processing. Because the requirements for weak lensing shape

measurement are so strict, the adopted approach is to apply the minimal number of pixel-level corrections, and, if feasible, to provide the information to reverse the correction if necessary (or to provide images for which a given correction has not been applied).

In detail, the pipeline that processes the VIS data calibrates and corrects a comprehensive set of pixel-level effects, which would otherwise bias the shear estimation (Cropper et al. 2013). These include electronic offset (bias), dark current, PRNU, detector-chain nonlinearity, brighter-fatter effect (Antilogus et al. 2014), charge-transfer inefficiency (Israel et al. 2015), illumination correction, and the flagging of cold, hot, and saturated pixels, as well as optical and electronic ghosts and cosmic rays. Moreover, the astrometric solution, required to be better than 30 mas, and photometric solution are computed, using Gaia DR3 (Gaia Collaboration 2023b) as the reference catalogue. Apart from robust shape measurements, these processing steps enable relative photometry measurement with an accuracy better than 1%. This implies that the collective contribution of residuals from the detection chain, small- and large-scale flat fielding, source extraction, and background and scattered light correction, is smaller than this target.

The main outputs of this step in the processing are the calibrated individual exposures, as well as stacks and their associated source catalogues. However, important supporting data are also provided, such as the calibration models, the background and flag maps, the distortion model, and a first estimate of the PSF model that is used for multi-band photometry during the generation of the main survey catalogue (Sect. 7.4). For reference, Fig. 27 displays a 2.5 × 2.0 view of a single raw frame (left) and the corresponding processed stack (right) of the *Euclid* self-calibration field. The left panel demonstrates the need for the robust detection of cosmic rays, whereas the deep image on the right shows the potential of *Euclid* to study low-surface-brightness features around galaxies.

### 7.2. Processing of NIR imaging data

The interpretation of the observed weak lensing signal requires accurate estimates of photometric redshifts (Sect. 2). As shown by Abdalla et al. (2008), the NIR images from the NISP instrument in the $Y_E$, $J_E$, and $H_E$ bands provide key information to improve the photometric redshift precision (see also Sect. 7.6.1). To reach this goal, NISP images require an exquisite calibration that takes into account all instrumental effects and a possible time variation of the telescope's throughput (Sect. 5.1). Most notably, the relative photometric calibration needs to be better than 1.5% for the entire magnitude range down to 24.0 AB. This requirement has been the main driver for the design of the processing pipeline.

The generation of the calibration products needed for the reduction of NISP images is carried out through a number of dedicated pipelines running on well-defined blocks of calibration observations. This is the case, for instance, for bad pixel identification, dark current, detector nonlinearity, PRNU, large-scale illumination correction, geometric distortions, persistence image creation, and relative and absolute photometric calibration factors.

The reduction of the NIR exposures can be divided into three main steps: the preprocessing, where all detector-related effects common to both NISP photometric and spectroscopic exposures are accounted for; the calibration part, which includes PSF estimation, computation of the astrometric solution based on VIS
<figure><img src="image_29.png" /><figcaption>**Fig. 27.** VIS view of a 2.5 × 2.0 wide area of Euclid's self-calibration field (Sect. 4.2.3) taken during the PV phase. Left: an unprocessed single exposure, where cosmic rays are clearly visible. Right: a VIS-processed stack using 42 exposures, or about 10 times the exposure time of the EWS. The ability of Euclid to reveal low-surface-brightness features is evident.</figcaption></figure>

catalogues, application of flat-fielding and photometric calibration based on the precomputed products, and source catalogue extraction from individual dithered observations; and lastly the stacking of frames for each observation sequence along with the generation of their PSF and source catalogues. Similarly to the VIS processing, we keep the amount of pixel-level operations to a minimum, while providing at the same time all relevant information about modelling and corrections, so that they may be reversed if needed or tailored for specific purposes.

For each of the three NIR filters, the output of the NIR processing consist of calibrated individual dithered observations and stacks with their associated source catalogues. These data products are provided to the subsequent processing step for the creation of the main survey catalogue (Sect. 7.4), along with all ancillary information, including root mean square (RMS), data quality flags, background, and PSF images.

Prior to launch, the design and implementation of the various processing steps have been extensively tested through a number of simulation campaigns with increasing realism and coverage of instrumental features, mostly based on modelling from ground calibration campaigns and in some cases by design. Further improvements will be made using the calibrations obtained during the PV phase (Sect. 4.3.1). Figure 28 gives an impression of the quality of data. It shows a small area in the Euclid self-calibration field (Sect. 4.2.3).

### 7.3. Processing of ground-based imaging data

The Euclid imaging data are complemented by supporting ground-based observations (Sect. 4.4) to improve the precision of the photometric redshifts and to determine the SED-weighted PSF model for the galaxies used in the weak lensing analysis (Sect. 7.6.3). The challenge is to uniformly calibrate this vast and heterogeneous external data set with extraordinary accuracy: the EWS is composed of three surveys (UNIONS, DES and the LSST) employing five telescopes with differing *u*, *g*, *r*, *i*, and *z* filters. The EDS and EAFs have primarily been observed using HSC, but there are additional datasets available featuring, for example intermediate band filters, that are being included to aid in the testing and validation of the photometric redshift

estimation and the VIS PSF modelling. Hence, the data are complex, with a wide range of single exposure depths and seeing conditions (from roughly 0.5 to 1.5), sampled with over 450 different detectors. The resulting data set will be large, ultimately comprising roughly 1.5 PB of raw (uncompressed) science observations. Crucially, to ensure robust photometric redshifts, the processing of these data should result in 1% single-filter photometric homogeneity over the entire EWS, and subpercent colour homogeneity over the VIS FoV (Eriksen & Hoekstra 2018).

For the cataloguing process described in Sect. 7.4, the required external data product inputs are coadded ground-based images. These coadded images are supplemented with associated PSF models for all objects that are identified in the VIS- and NIR-based detection process. These coadded images are built using a single pipeline whose software components were initially developed as part of the DES data management system (Mohr et al. 2008) and have been tuned and validated to serve our needs and then integrated into the SGS code base. The pixel coaddition is based on the widely used SWarp code (Bertin et al. 2002; Bertin 2010). We developed a dedicated code to model the PSF to take into account the position-dependent PSF models of all input images.

To facilitate the processing of the individual ground-based exposures, a common data model is enforced across all external ground-based data sets. This common input data product–termed a single-epoch frame (SEF)– consists of a detrended and astrometrically and photometrically calibrated single CCD image, the associated position-dependent PSF model and an associated catalogue that includes, at a minimum, the sky positions and PSF-fitted magnitudes of the brighter, unresolved sources. Typically, these SEFs are created using output data products from the external surveys. For example, the ensemble of *i*-band SEFs from Pan-STARRS is prepared using the software within the Pan-STARRS collaboration (Magnier et al. 2020; Waters et al. 2020). Similarly the *g*- and *z*-band SEFs from WHIGS and WISHES are produced using output data products from HSCpipe (Bosch et al. 2018, 2019), which is software developed by the HSC and LSST teams. The *u*- and *r*-band SEFs from CFHT are created using the MegaPipe software developed by the CFHT team (Gwyn 2008).
<figure><img src="image_30.png" /><figcaption>**Fig. 28.** False-colour NISP image of a 4.5 × 3.0 area of *Euclid*'s self-calibration field (Sect. 4.2.3). Filters $Y_E$, $J_E$, and $H_E$ are shown in blue, green, and red, respectively. The depth is that of the EDS (Sect. 4.2), about 26.4 AB mag per band. The bright star has 11.5 AB mag, showcasing *Euclid*'s excellent performance for in-field stray light suppression. Field rotation between observations is evident from the diffraction spikes.</figcaption></figure>

For the EDS data, we receive the ground-based data from the Cosmic Dawn team (Euclid Collaboration: McPartland et al., in prep.; Sect. 4.4) in the case of Euclid Deep Field North (EDF-N) and Euclid Deep Field Fornax (EDF-F), while we directly downloaded public data from the HSC archive for other EAFs such as COSMOS and produce the data ourselves in the case of Euclid Deep Field South (EDF-S). For the publicly available DES data (Abbott et al. 2018b) and other Dark Energy Camera (DECam) (Flaugher et al. 2015) data, we produce and calibrate the SEFs ourselves using extended versions of pipelines originally developed for DES (Mohr et al. 2012; Desai et al. 2012) that include also image masking (Desai et al. 2016) and have been reorganised and rewritten to simplify the large-scale processing required. Finally, thanks to a collaboration agreement, the data products from the LSST are provided through a joint Euclid-Rubin working group (see Guy et al. 2022, for details).

To enable an accurate and uniform photometric calibration of the ground-based external data across the full extragalactic sky that will be observed by *Euclid*, we leverage *Gaia* (*Gaia Collaboration 2016*) photometry and spectroscopy that are stable across the sky with a systematic uncertainty of ~2 mmag (*Gaia Collaboration 2023b*). Initial testing of statistical transformations from *Gaia G*, *BP* and *RP* to each of the external *griz* bands in DES demonstrated a high level of consistency between the DES DR2 (Abbott et al. 2021) calibration and the *Gaia*-based predictions (George et al. 2020). Further tests with the UNIONS dataset provided indications that the external photometry predictions from statistical transformations

based on *Gaia G*, *BP* and *RP* photometry exhibited improved stability across the sky and improved internal self-consistency in comparison to the original UNIONS calibration. This demonstrated the promise of adopting *Gaia* data as a basis for calibrating the heterogeneous ground-based data sets to ensure consistent photometric redshifts and stellar SEDs across the sky.

A challenge in employing the statistical transformation function from *Gaia G*, *BP* and *RP* (George et al. 2020) to the external data bands is that systematic changes in the stellar populations over the *Euclid* sky could bias the *Gaia*-based predictions, introducing errors larger than the advertised ~2 mmag systematic uncertainties in *Gaia* photometry. Therefore, with the availability of the *Gaia* DR3 (*Gaia Collaboration 2023b*) *BP* and *RP* calibrated spectra, we have transitioned to using *Gaia* spectra and the associated synthetic magnitudes to calibrate the external data. We calculated *Gaia* synthetic magnitudes using the *Gaia* calibrated spectra together with the appropriate bandpass for each external survey band and camera combination. Under the assumption that the external data bandpasses are robust, these *Gaia* synthetic magnitudes are highly accurate, enabling us to meet the *Euclid* requirements. In the case that there are significant residual errors in the external data bandpasses, these errors would enter both the ground-based photometric calibration and the SED fitting being performed for photometric redshifts and stellar SED constraints. In addition, we expect to be able to use the *Gaia* spectra to aid in characterising bandpass variations across the focal plane in the ensemble of
cameras being used to obtain the ground-based data needed for *Euclid*.

For the UNIONS and DES external data sets there are enough *Gaia* synthetic magnitude constraints per SEF to enable a photometric calibration of individual SEFs that approaches the 2 mmag systematic floor of the *Gaia* mission. For the EDS and EAF data sets, the situation is more complicated, because the typical HSC integrations are longer, pushing the saturation limits in these SEFs to fainter magnitudes and reducing the overlap between the *Gaia*-based synthetic magnitudes and the external data sets. To address this we incorporate also relative photometric constraints between overlapping pairs of SEFs. This allows us to combine the direct photometric constraints from *Gaia* synthetic magnitudes across an ensemble of co-located SEFs, delivering a comparable level of photometric calibration accuracy in the EDS and EAF datasets as we achieve in the EWS.

## 7.4. Catalogue creation

The imaging data form the basis of the object catalogue that is used for most of the subsequent analyses. To avoid multiple entries for the same object, the survey area is divided into predefined tiles. Each tile consists of a rectangular extended area and a core area that is defined by a set of HEALPix indices (Górski et al. 2005). All data needed for the detection and photometry are generated to cover the extended tile area, which overlaps with neighbouring tiles. Only the objects in the core area, which is unique for each tile, are actually selected for insertion into the object catalogue (Kümmel et al. 2022). This procedure avoids multiple detections of identical objects in adjacent tiles, while allowing the proper processing of large, extended objects in the overlap areas between two tiles.

The pipeline retrieves all imaging data for the tile of interest and creates coaddends of the calibrated *Euclid* VIS and NISP exposures. As part of this step, the background is subtracted from each exposure and each image is sampled to the native VIS pixel scale of 0.′1 pixel$^{-1}$. Information on flagged pixels and areas around bright stars that are affected by diffraction spikes or blooming is propagated, as is the case for the model PSF in each exposure.

In principle, the information in all bands could be combined to maximise the prospects for object detection, but this would result in complex selection biases for the primary probes. To avoid this, we opt for a staged process, where we detect and subsequently deblend objects in the VIS and NIR bands separately. This ensures a clean selection of the weak lensing source sample, while recording all objects that are visible in the *Euclid* data. For the NIR detections we employ a deep image generated from the combined $Y_E, J_E$, and $H_E$ data. Source detection is done with *SourceXtractor++*<fnref n="17" /> (Bertin et al. 2020), with detection parameters that are optimised for completeness, while keeping the false detection rate below 1%. The *denclue* algorithm (Tramacere et al. 2016) is used for the deblending procedure. The resulting VIS and NIR detections are then combined into a single joint catalogue that is used to perform photometry and to determine a number of other properties. To enable the selection of VIS-detected objects, the combined catalogue contains a flag that indicates whether a source was detected in VIS or in the deep NIR image.

Accurate multi-band photometry is essential for the determination of robust photometric redshifts, but optimising the pipeline for this application might not be ideal for other science

cases. Therefore, to maximise the usefulness of the catalogue, fluxes are measured in the following ways.

- Total object flux within a Kron aperture on the detection image with T-PHOT<fnref n="18" /> (Merlin et al. 2015).

- Isophotal flux measured by summing the flux of the pixels above the detection threshold in the detection image.

- PSF flux measured on the VIS image.

- Aperture flux measured by A-PHOT (Merlin et al. 2019) on images that are PSF-matched to the one with the worst resolution (typically a ground-based image). The circular aperture is set to twice the worst FWHM for each object.

- Template-fitting photometry computed by PSF-convolving the VIS object shape to the different bands and fitting the surface-brightness profiles.

- Single-Sérsic fitting photometry using *SourceXtractor++* (Kümmel et al. 2022) in all available bands.

To increase the scientific value of the catalogues, morphological properties are also determined. A robust separation of point-like and extended objects is achieved by computing the difference between the central surface brightness $\mu_{max}$ and the total brightness. This quantity is translated to a probability of the object being point-like using calibrations based on simulated data (Sect. 6.2). We include the non-parametric estimates for concentration, asymmetry and smoothness (CAS; Tohill et al. 2021), as well as the Gini index (Lotz et al. 2004) for all objects. Moreover, the catalogue includes the best-fit Sérsic models in all bands obtained using *SourceXtractor++* (Kümmel et al. 2022).

## 7.5. Processing of NIR spectroscopy data

The first steps in the processing of the dispersed NIR images are similar to what is done with the imaging data used for photometry (Sect. 7.2): the same detector-level effects need to be accounted for (bad pixel flagging, nonlinearity correction, persistence masking, dark subtraction, and cosmic-ray rejection). Therefore, these steps are based on a common set of processing elements. The different nature of the images does lead to some modifications, for instance in the cosmic ray rejection. The main difference, however, pertains to the instrument models that are adopted for the two observing modes. Specific calibration products are used for each case. The resulting images are then used to extract the spectra.

### 7.5.1. Extraction of spectra

To extract the spectra from the preprocessed NISP spectroscopic images, we need to precisely locate the dispersed image of each object (the so-called ‘spectrogram’) in the spectroscopic frame. Specifically, we are interested in the first order of the grism dispersion, where more than 96% of the object flux is concentrated. To this end, the object catalogues derived from the corresponding direct photometric images are critical.

First, the two-dimensional first-order spectrograms for each object are located on the full NISP spectroscopic frame, based on its coordinates measured from the corresponding direct image. This involves applying the astrometric solution, which translates sky coordinates into detector reference positions of the different dispersion orders. The precise location of the spectrogram is then traced, accounting for any inclination and curvature. Along the

<footnote n="18">¹⁸ The transformation of the PSF from the VIS band with the highest resolution to the NIR and ground-based images, which is necessary for A-PHOT and T-PHOT, is done using convolution kernels as described in Boucaud et al. (2016).</footnote>

<footnote n="17">¹⁷ https://github.com/astrorama/SourceXtractorPlusPlus</footnote>
dispersion direction, the spectrum is re-mapped into wavelength steps, using the wavelength solution computed from reference planetary nebulae emission-line spectra (see Sect. 4.3.1, and Euclid Collaboration: Paterson et al. 2023).

Detector-level pixel-to-pixel variations are corrected from calibrated detector flats, and values for the background over the full focal plane are sampled on detector areas where no spectra are present. These are then averaged per detector, and the values subtracted over the whole frame.

In slitless observations, overlapping spectra from nearby objects represent an important contribution to the noise affecting a given spectrum. Euclid adopts a specific observing sequence to mitigate this, collecting four exposures at varying dispersion directions in each ROS, as described in detail in Sect. 3.5.3. A model of the spurious contribution of each neighbouring source contaminant is built using its spatial extent estimated from the J_E-band photometric image, and its intrinsic SED. If possible, the latter is estimated from uncontaminated portions of the spectrograms extracted from one of the available exposures taken as part of the ROS. If an object happens to be contaminated in all four exposures, a power-law SED interpolation of the available NIR photometric measurements is used instead. These decontamination operations are performed on each identified 2D spectrogram, ultimately resampled on a rectilinear grid along the dispersion and cross-dispersion directions.

At this stage, 1D spectra can then be extracted from the individual exposures in a ROS by properly integrating over the cross-dispersion dimension. This is performed through optimal extraction (Robertson 1986; Horne 1986), using a weighting profile derived from the J_E-band photometric image of the source. This also includes appropriate rotation, matching of both NISP-P and NISP-S PSFs, shear, and resampling, to account, in the case of galaxies, for the inclination of the source with respect to the dispersion direction, so as to minimise effective spectral line-spread function and maximise S/N.

The spectrum extraction is followed by a chromatic relative flux calibration, using bright point sources in the self-calibration field repeatedly observed over the NISP field of view. This process normalises to the same relative flux scale all 1D spectra from different observations, detectors, and location in the focal plane. The absolute flux is then obtained using the overall chromatic sensitivity curves, derived from observations of spectrophotometric standards for each observing setup.

Finally, the flux-calibrated 1D spectra corresponding to each of the four exposures in a ROS<fnref n="19" /> are averaged into a combined 1D spectrum using inverse-variance weighting. During this operation, statistically outlying pixel values (corresponding, e.g. to cosmic ray or contamination residuals) are identified and discarded. This is this 1D spectrum that is then passed to the next step of the spectral analysis, together with the corresponding statistical variance, bit mask – that is per-pixel flagging – and combined effective line-spread function estimate.

Fig. 29 shows the results for a galaxy in the COSMOS field with a known redshift of $z = 1.1770 \pm 0.0005$ (Mainieri et al. 2007). Comparison of the redshift and line flux estimates allows us to quantify the accuracy of the NISP spectroscopic calibrations. The top panel shows the four individual spectrograms extracted for this object over the full red grism (RG)

<figure><img src="image_31.png" /><figcaption>Fig. 29. An example of NISP spectroscopic data for a galaxy in the COSMOS field with $z = 1.1770 \pm 0.0005$ (Mainieri et al. 2007). The top figure shows the four spectrograms, with the $H\alpha$ line clearly visible. The bottom plot shows the corresponding combined and flux-calibrated 1D spectrum (in black) and its associated statistical noise (in orange), while the green line shows the combined continuum and emission line model that fits the data best. The bright $H\alpha$ line is detected with S/N = 14 yielding a redshift of $z = 1.1783 \pm 0.0005$ (vertical blue line), which is in agreement with the previously published value. The flux of the line, $f_{H\alpha} = 2 \times 10^{-15} \text{ erg cm}^{-2} \text{ s}^{-1}$, is approximately ten times higher than the limiting flux for the EWS.</figcaption></figure>

domain, after applying the decontamination procedure necessary to remove the signal from nearby objects, and rectilinear resampling. The bottom panel shows the combined and flux calibrated 1D spectrum for this object, with the $H\alpha$ line clearly detected. We note that this galaxy is not representative for the majority of emission line detections, because its $H\alpha$ flux is approximately ten times brighter than the limiting flux for the EWS.

### 7.5.2. Redshift measurement

The next step is to determine redshifts and measure line fluxes from the extracted spectra, and to provide an estimate of the reliability of the measurements. In the case of EWS observations, the spectra are obtained with the RG, while for the EDFs, the blue grism (BG) spectra are included to enhance the performance of the redshift estimation and to maximise the purity and completeness. Photometric measurements from VIS and NISP could be optionally included and used for the redshift measurement.

The redshift measurement itself is performed through template fitting over a regular grid of redshifts, spaced logarithmically. A first pass uses a coarser logarithmic grid (initial step $\Delta z = 10^{-3}$ logarithmic grid, after which a redshift probability density function (zPDF) is calculated for each model. All those individual zPDFs are combined into one so-called 'first pass zPDF'. A first list of redshift solutions is based on the main peaks (up to 10) of this first zPDF. The redshift measurement is then refined around each of these solutions by refitting over all the

<footnote n="19">¹⁹ These may in fact be more or less than four, depending on the specific position of the object on the detector. In the EWS, the few cases when more than four sub-exposures are available correspond to objects near the borders of the field, which benefit from the overlaps between pointings.</footnote>
models, but with a finer redshift grid; the zPDF is recalculated and the best final solutions (up to five) and corresponding models are then obtained from the strongest peaks. The peaks are ordered following the value of the integral of the zPDF under each peak with a 3 σ window.

The model includes a set of six distinct continuum models representing various star-formation histories from the Bruzual & Charlot (2003) library. Given the resolution and limited wavelength coverage of NISP, we found that it is not necessary to have more detailed continuum parameters, because it only leads to degeneracies between the different templates. We then add nebular emission lines from templates describing various ratios between the different emission lines. These templates have been built from the more than 35 000 galaxies of the VVDS programme (Le Fèvre et al. 2013), which provides an unbiased sample of spectra for a wide range of galaxy types with 0.4 < z < 3, but with magnitude limit of $I_{AB} = 24$. Using these templates is more efficient than fitting all possible emission lines independently, since the number of free parameters becomes too high to provide a reliable solution. Finally the model includes intrinsic reddening based on Calzetti et al. (2000) and intergalactic medium (IGM) absorption as described in Meiksin (2005).

It is also possible to include priors in the calculation of the zPDF, which can be used to favour solutions where a detected emission line is identified as one of the strongest emission lines usually detected in galaxy spectra (e.g. [O II], [O III], $H\alpha$), or just $H\alpha$. Alternatively, one can use the known redshift distribution of $H\alpha$ emitters as a baseline for the zPDF instead of a flat prior. Tests have shown that such an empirical prior gives the best results in terms of purity and completeness (also see Sect. 7.7.1).

In addition to the galaxy model, the pipeline also provides solutions for quasar and stellar models. The quasar model is built in the same way as was done for galaxies, but including a series of double power-law continuum and Lorentzian broad emission lines. The stellar models are built from a set of 36 templates covering all stellar types (Pickles 1998). The selected object class corresponds to the highest statistical evidence over the three categories (galaxy, star, or quasar).

Once the redshift has been calculated, the fluxes of the detected emission lines are measured using both direct integration (DI) and a Gaussian fit (GF). In the DI method, the spectrum is first continuum-subtracted, using the continuum evaluated from a median-iterate filtering smoothing with a variable window. Each line is then integrated, starting from the position of the peak as provided in the previous step, until the flux remains positive; together with the flux, also the S/N, equivalent width (EW), and position of the lines are provided. In the GF method, a multi-Gaussian model plus a constant continuum model is considered, with N Gaussians, depending on the line considered: (i) if a blend with adjacent lines is expected (e.g. we assume three Gaussians for the $H\alpha$+[N II] doublet complex); or (ii) if the lines are separate but close in wavelength, and their ratios can be linked through physics (e.g. two Gaussians for the [O III] and [S II] doublets), or if a line is isolated (e.g. one for $H\beta$). In contrast to the DI method, this allows for us to deconvolve the contribution of different lines in a complex (e.g. $H\alpha$ from [N II]). This model is used to fit the data with a Levenberg-Marquardt algorithm, deriving the flux, S/N, continuum, FWHM, EW, and wavelength for all the lines, assuming the (up to) five redshift solutions (with their uncertainties) obtained in the step before for each galaxy.

Spectra of objects below redshift 0.9 display very few features in NISP spectra, and hence any artefact in the spectra

might be misinterpreted as an emission line. It is therefore essential to identify those spectra that are affected by artefacts: these interlopers could outnumber the targets of the EWS ($0.9 < z < 1.8$) because the redshift distribution peaks around 0.5 for the limiting AB magnitude of $H_e = 24.0$ in the EWS. To discriminate objects with secure redshift measurements from possible interlopers, the pipeline provides a numeric indicator of the 'reliability' of the redshift measurement, which quantifies the quality of the spectrum based on the analysis of the zPDF using a deep-learning algorithm. This algorithm will be trained on EWS-like observations of the EDS, for which a correct redshift at the 99% confidence level is expected. This training will be applied to the set of zPDF of the EWS spectra to quantify the reliability of the redshift measurements.

## 7.6. Measurements for 3×2pt statistics

As discussed in Sect. 2.3 the shape measurements from Euclid are combined with photometric redshifts derived from multi-band photometry. The precise photometric redshifts are required to divide the sample into tomographic bins, while their redshift distributions need to be characterised well for a correct interpretation of the clustering and lensing signals. In Sect. 7.6.1 we describe how the photometric redshifts are determined and we plan to calibrate the corresponding redshift distributions to high accuracy. As described in Sect. 7.6.2, the photometry is also used to classify the objects for further science applications.

The galaxy shape measurements benefit greatly from the sharp diffraction-limited PSF. Nonetheless, our objectives require its size and shape to be determined with unprecedented accuracy. This implies that we have exquisite knowledge of the optical properties of the telescope and understand how the detectors record the incoming photons. This is, however, not sufficient: the PSF is chromatic, which means we need to estimate the appropriate PSF for each galaxy based on its observed SED. This requires a dedicated procedure that is described in Sect. 7.6.3.

The requirements for Euclid were derived by considering the change in the observed unweighted quadrupole moments, $Q_{ij}^{\text{obs}}$, when a galaxy image is convolved by the PSF. For an object with an observed surface-brightness distribution, $I^{\text{obs}}(\mathbf{x})$, and total flux, $F$, its quadrupole moments are defined as

$$Q_{ij} = \frac{1}{F} \int d^2\mathbf{x} \, x_i x_j I^{\text{obs}}(\mathbf{x}). \quad (20)$$

The shape of an object can then be quantified by combining the quadrupole moments into the complex polarisation (or distortion)$^{20}$, $e = e_1 + ie_2$, where

$$e_1 = \frac{Q_{11} - Q_{22}}{Q_{11} + Q_{22}}, \quad \text{and} \quad e_2 = \frac{2 Q_{12}}{Q_{11} + Q_{22}}, \quad (21)$$

while the size of an object is captured by

$$R^2 = Q_{11} + Q_{22}. \quad (22)$$

Although unweighted moments are not practical in the presence of noise and blending, they do provide a convenient framework to quantify the impact of the PSF on shape measurement (Paulin-Henriksson et al. 2009; Massey et al. 2013). In this case, the observed quadrupole moments of a galaxy are given by the

<footnote>$^{20}$ This definition is related to the third flattening (see Sect. 2.2) through $\epsilon = e/(1 + \sqrt{1 - |e^2|})$.</footnote>
sum of the quadrupole moments of the true galaxy image and those of the PSF (Valdes et al. 1983), so that

$$Q_{ij}^{\text{obs}} = Q_{ij}^{\text{gal}} + Q_{ij}^{\text{PSF}}. \qquad (23)$$

Cropper et al. (2013) used this to allocate tolerances for a wide range of instrumental effects, starting from an acceptable level of bias in cosmological parameters. A limitation of this simple ‘flow down’ is that it did not capture spatial variations of sources of bias, resulting in conservative estimates for the impact of residual systematics on the observed lensing signal (Kitching et al. 2019). To explore a more realistic scenario, Euclid Collaboration: Paykari et al. (2020) considered a ‘flow-up’ from perturbations to the defocus of the PSF model to the bias in cosmological parameters. They found that PSF variations at the level of the requirements induced biases of about 4% of the expected statistical uncertainty in dark energy parameters, well within requirements. Further tests, using simulated data, confirmed that the requirements derived by Cropper et al. (2013) are adequate and conservative and thus remain the basis for the calibration of VIS data (Sect. 4.3) and the modelling of the Euclid PSF, which is described in detail in Sect. 7.6.4. The challenging task of accurate shear estimation from the Euclid images is discussed in Sect. 7.6.5.

### 7.6.1. Photometric redshift estimation

As shown in Eq. (14), the amplitude of the lensing signal depends on the redshift distribution of the sources, in addition to the cosmological parameters. Moreover, being able to separate the overall galaxy sample into subsamples that are separated in distance enhances the amount of information that can be extracted from a weak lensing, or 3×2pt analysis. The current baseline configuration uses 13 evenly populated bins in the redshift range 0.2 < z < 2.5.

The mean redshift of these subsamples is of particular importance for cosmic shear, since a bias in the estimated mean redshift directly translates into a bias in cosmological parameter inference. Given the objectives of Euclid, this implies that the mean redshift must be known to $\sigma_{\langle z \rangle} < 0.002(1+z)$ per bin (Ma et al. 2006; Amara & Réfrégier 2007). Importantly, the modelling of the photometric clustering signal in the 3×2pt analysis relies not only on accurate knowledge of the mean redshift, but also of the width of the redshift distributions of the different subsamples (e.g. Tutusaus et al. 2020; Porredon et al. 2022). The need for precise photometric redshifts is therefore two-fold: to place objects in narrow tomographic redshift intervals, minimising the overlap between subsamples; and to ensure that the distribution in redshift of the sources in each subsample is accurately known. These two needs are addressed separately in the photometric redshift pipeline, since the former requires that we optimise the photo-z for precision, while the latter demands very high accuracy.

For individual galaxies, the need to place objects in tight redshift intervals, with $\sigma_z < 0.05(1+z)$ being a key requirement, leads to two considerations: the effective dispersion of the zPDF around its assumed true redshift (‘photo-z scatter’); and the probability that our assumed redshift for an object is far from its true redshift (‘outlier rate’). Further requirements are placed on the choice of the photo-z algorithm by the nature of the EWS. It must be able to run efficiently for the billions of galaxies we will observe and to be able to account for effects that vary on an object-by-object basis, namely the impact of Galactic reddening and the variation of filter transmission curves across the

<figure><img src="image_32.png" /><figcaption>Fig. 30. Photometric redshift performance of the mode of individual probability distributions using NNPZ, taken from Euclid Collaboration: Desprez et al. (2020) who used simulated DES and Euclid NIR data. Regions of photometric redshift space that will be excluded from the weak lensing analyses are shown in grey.</figcaption></figure>

FoV of an instrument (Euclid Collaboration: Paltani et al. 2024). These additional concerns led to the development of the NNPZ (Nearest-Neighbour Photometric Redshifts) algorithm; its performance on a sample of spectroscopic redshifts taken from the COSMOS field is shown in Fig. 30. These results were taken from the Euclid photo-z data challenge of Euclid Collaboration: Desprez et al. (2020), which compared different methods using simulated DES (optical) and Euclid (NIR) data. The photometric redshift quantity that is used to assign tomographic redshift bins is the mode of the NNPZ zPDF, and for the sample shown in Fig. 30 the scatter ($\sigma_z = 0.059$) and outlier fraction ($\eta = 11.1\%$) of this photo-z with respect to the spectroscopic redshifts lie slightly outside of our requirements (< 0.05 and 10%, respectively). This is adequate for the first data release (DR1; see Sect. 7.8), while we expect to achieve a better performance for later releases, once much deeper LSST data become available.

Machine-learning algorithms such as NNPZ rely heavily on how representative the reference training data are of the target data set. Complete and representative spectroscopic samples are impractical to assemble for the EWS, and so instead we will build a sample of galaxies with accurate and precise photometric redshift zPDFs from a carefully selected and curated set of deep reference fields, the EAFs. These reference field photo-zs will be computed using a new custom SED-fitting package called Phosphoros (see Sect. 3.3 in Euclid Collaboration: Desprez et al. 2020, for a brief description). Through tagging the SEDs with properties, such as a mass-to-light ratio or metallicity, Phosphoros is able to additionally produce sampled posterior distributions of physical quantities (e.g. stellar masses) for each reference galaxy. The zPDF and multi-dimensional physical-property posterior information can then be propagated from these reference field galaxies to the EWS galaxies with NNPZ. In this sense, NNPZ operates as an accelerator for the template-fitting approach, and its lack of an explicit training phase allows us to account for effects that vary on an object-by-object basis (Euclid Collaboration: Paltani et al. 2024). The zPDFs we
output at this stage are still expected to be subject to biases, due to our imperfect knowledge of galaxy SEDs and the priors we use, and are therefore used only for tomographic redshift binning in the weak lensing analysis. The galaxy physical properties inferred through this process are not used in the weak gravitational lensing analyses, but will be used for a vast range of non-cosmological applications (see Sect. 10).

Once galaxies have been assigned to their tomographic redshift samples, we must reconstruct their collective distribution in redshift, *n(z)*, for cosmological inference of the weak lensing signal. The strategy to achieve the target accuracy of $\sigma_{\langle z \rangle} < 0.002(1+z)$ per tomographic bin was outlined in Masters et al. (2015) and has evolved only a little since then. Briefly, the 8- or 9-band photometric space of the target galaxy sample is quantised by way of a SOM (Kohonen & Honkela 2007), resulting in a 2-dimensional array of vertices or ‘cells’ (also see Fig. 19). Each cell of the SOM is represented by a vector of values, where each value corresponds to the ratio of flux in a photometric band with respect to the $I_E$ value, and thus each cell represents a possible galaxy SED with free amplitude. Galaxies in the weak lensing sample are each assigned to their closest SOM cell in the flux-ratio space used to construct the SOM. The SOM algorithm preserves locality of the input space, in the sense that similar vectors of flux ratios with respect to $I_E$ will be nearby to one another in the 2-dimensional SOM map, and therefore each cell has a finite size in the input parameter space. The Euclid survey is supported by ground-based observations with different instruments and filters for the northern and southern regions (see Sect. 4.4). To place all galaxies correctly in the SOM we perform a band-standardisation step whereby a per-galaxy colour term is computed and applied to the observed photometry in each band. Corrections for Galactic reddening and bandpass variation are also carried out during this step.

Construction of the $n(z)$ is based on the spectroscopic sample built up through the C3R2 programme (see Sect. 4.5). The sample is designed to cover the diversity of galaxy SEDs in the shear sample, but the measured SED of an object depends on the noise properties of the observations and other survey characteristics. In order to be able to account for varying survey depth, the spectroscopic sample is drawn from the EAFs (Sect. 4.2.2), where the S/N of the data are five times higher than the EWS. We can then draw multiple realisations of the EWS noise properties and apply them to these objects to form a calibration sample that covers both the range in galaxy SEDs and photometric scatter of the shear sample.

With the calibration sample in hand we perform the $n(z)$ estimation for each tomographic redshift subsample separately. Each galaxy in a subsample is assigned to its best-matching SOM cell and the shear weights of objects (see Sect. 7.6.5) are summed within each SOM cell. These sums of weights represent the statistical power of each SOM cell, with an analogous quantity determined by each objects' probability to pass the selection flux cut in the case of photometric clustering. The calibration objects are treated in the same way, with the exception that their ‘shear weights’ are always unity, that is, their zPDFs are measured via NNPZ and their artificially noised flux measurements are dereddened and band-standardised in the same way as for the shear sample. The calibration objects are then assigned to tomographic bins on the basis of the mode of the zPDF, exactly as for the shear sample objects. They are also rendered into the SOM and the number of them per SOM cell, within the given tomographic bin, is computed. The redshift distribution of a tomographic bin is then just a weighted histogram of the spectroscopic redshifts, where the weight of a spectroscopic object is the sum of

shear weights in its cell, divided by the number of spectroscopic objects in that cell. In this way, all trusted redshifts per cell are used, but the relative sampling difference between the shear sample and the selection of objects for which we have spectroscopic redshifts is accounted for.

There are some SOM cells containing galaxies that can be used for the weak lensing analysis, but that lack spectroscopic measurements. Those galaxies are flagged and removed from the analysis because we are unable to represent them in the recovered redshift distribution. With this strategy, the $n(z)$ distribution is subject to sample noise, but Masters et al. (2015) showed that we will be able to meet our target requirement of $\sigma_{\langle z \rangle} < 0.002(1+z)$ on the mean redshift per tomographic bin, provided that each tomographic subsample is represented by a large (> 600) number of cells.

As a further validation of the redshift distributions estimated with the SOM method, we will use clustering redshifts (Newman 2008). This method employs angular cross-correlation measurements of the positions of the sources in the tomographic bins and spectroscopic calibration samples that overlap on the sky. Crucially, this drops the assumption of the calibration sample fully covering the colour and magnitude range of the source sample. As such, different calibration samples that are typically brighter and cover larger areas are used, which makes this approach highly complementary to the SOM method. Measuring the angular cross-correlation amplitude between the tomographic source samples and the spectroscopic calibration samples finely binned in redshift can result in an accurate reconstruction of the tomographic redshift distributions. As shown in Naidoo et al. (2023), an overlap area of only a few hundred square degrees between the Euclid source sample and spectroscopic surveys like BOSS, DESI, and the Euclid NISP-S sample is sufficient to reach the requirement on the accuracy of the mean redshift of the tomographic bins.

### 7.6.2. Classification

The objects detected in Euclid images comprise galaxies, stars, QSOs, globular clusters, Solar System objects, and a diverse array of contaminants and artefacts. In order to avoid introducing selection effects in the weak gravitational lensing analysis, each detected object with a full complement of photometric information available is treated as a possible galaxy at the photometric redshift determination stage. Contaminants and stars are determined during the shape measurement process (see Sect. 7.6.5) by their low shear weights. However, to facilitate the modelling of the PSF we must identify a very pure sample of stars that spans the range of stellar colours. Using only morphological information, such as an object's apparent compactness, risks introducing biases into the PSF model, and so this step is ideally performed using only photometric information.

Our classification uses a set of three pretrained probabilistic random forests (pRFs), each a binary classifier for a type of object, such as star versus not a star. The pRF returns a probability that an object is of that type and thus each object receives three probabilities, one for each of the considered types, namely star, galaxy, and QSO. An object is assigned a given class if its probability to be that class exceeds a predetermined threshold, and a single object can be given multiple classes. Only those objects classed as stars and not any other additional classes are considered suitable for PSF modelling. Class probabilities are retained for use in non-cosmological legacy science processing (see Sect. 10).
### 7.6.3. SED modelling

The convolution of galaxy images by the PSF is the dominant source of bias for weak lensing studies, and an accurate estimate of the PSF is required to obtain unbiased shape measurements. Here, we focus on the challenges that arise from the fact that the *Euclid* telescope is diffraction limited and the VIS passband is very broad (see Cypriano et al. 2010; Eriksen & Hoekstra 2018, for more background). In the analysis we need to take into account that each galaxy is convolved by an effective PSF that depends on its SED. Moreover, to determine the underlying optical PSF model successfully, we must know the SED of each of the stars that is used in this process.

For bright stars ($G \le 16.5$) *Gaia* spectra can be used directly, but for fainter objects we must estimate their SEDs at high accuracy from their broad-band photometry (Eriksen & Hoekstra 2018). To perform this task we re-use NNPZ, the nearest-neighbour method that was developed for the main photometric redshift pipeline, but choose as target quantities the fluxes at different wavelengths. Similar to what is done to determine the redshift probability density functions (PDFs), we construct a reference data set of objects that is in principle representative of the objects for which we need to recover the SED information. Specifically, we use the *Gaia* spectra, assuming that across the *Euclid* footprint we have examples of all relevant stellar types and metallicities. An analysis of the impact of differential Galactic reddening between brighter and fainter stars of the same types (and thus different radial distances) is ongoing. We integrate the *Gaia* spectra through a series of 55 synthetic narrow-band filters of width 10 nm, filling the range $450 \text{ nm} < \lambda < 1000 \text{ nm}$, and recover the weighted mean flux from the 30 closest neighbours. The weight is a pseudo-likelihood computed from the $\chi^2$ distance in flux space between the target object and a reference object, with a free SED amplitude parameter.

The measurement of galaxy SEDs follows a similar procedure, but in this case we lack an equivalent to the *Gaia* spectro-photometric data set. Moreover, existing sets of galaxy SED templates are either incomplete or biased and would thus introduce biases in the PSF construction that would propagate to cosmological parameter estimation (see Eriksen & Hoekstra 2018).

The strategy to create the NNPZ reference sample for galaxy SEDs is summarised as follows. Instead of using spectra, we begin with broad- and intermediate-band photometry from the COSMOS and CDFS fields to provide a coarse sampling of galaxy SEDs through much of the $I_E$ bandpass range. We then apply a combination of Gaussian process (GP) interpolation and template-guided filter colour terms to achieve the finer 55 narrow-band sampling that we also use for the stellar SEDs. To reduce the impact of sample variance we plan to collect additional medium band data across the EAFs.

### 7.6.4. VIS PSF model

As discussed in detail in Cropper et al. (2013), accurate measurement of weak lensing shear imposes stringent requirements on model accuracy for the *Euclid* VIS PSF. In the case of *Euclid*, this means that the PSF model must be known throughout the mission lifetime with a residual temporal-spatial model uncertainty of $\sigma(e) < 2 \times 10^{-4}$ per ellipticity component, while the size needs to be known such that $\sigma(\Delta R_{PSF}^2)/R_{PSF}^2 < 10^{-3}$. Although the *Euclid* VIS system possesses a highly stable, diffraction-limited PSF, meeting these requirements is nontheless challenging because of a number of design choices.

First, to reach the required depth, VIS uses a broad band-pass, but this also results in a strong chromatic dependence for the PSF. Hence, as already discussed in Sect. 7.6.3, the PSF varies between stars and galaxies according to their colour or SED, and within galaxies due to local changes in stellar populations (i.e. colour gradients within galaxies; Semboloni et al. 2013a; Er et al. 2018). Moreover, *Euclid* utilises a dichroic to split the visible and NIR components of the beam, which was designed to produce a hard bandpass edge with minimal out-of-band light (in conjunction with coatings also applied to two of the fold mirrors). This is achieved using a complex multi-layer dielectric coating, which induces wavelength- and polarisation-dependent phase errors on the wavefront and consequent effects in the PSF model (Gaspar Venancio et al. 2016). This is further influenced by chromatic, polarised reflection from the silver optical surfaces and fold mirror coatings, and transmission at the detector. Finally, the pixel size of the *Euclid* VIS CCDs is such that the PSF is undersampled at all wavelengths in the VIS bandpass.

Hence, in addition to capturing the variation across the FoV, the VIS PSF model also needs to allow for the propagation of the model between stars and galaxies with differing SEDs. Together with the need to produce an oversampled model, to address pixel undersampling, this motivates utilising a PSF forward-modelling approach. The Fraunhofer condition links the image-plane PSF to the wavefront error at the exit pupil of the telescope (e.g. Hopkins & Yzuel 1970). Under this condition, one may forward model the PSF on the detector FPA using a combination of chromatic wavefront modelling for the optical contributions and both chromatic and achromatic convolutional kernels for the detector and guiding error contributions (also see Ma et al. 2008). When modelling the PSF, it is essential to analyse image data that have had all linear and nonlinear detector effects corrected, including detector nonlinearity and the brighter-fatter effect (BFE), the latter currently being corrected by a modified version of the algorithm due to Coulton et al. (2018).

An important benefit of forward modelling is that a wide range of effects can be included in a consistent fashion. In the case of the VIS PSF model, these include the following.

- A model for vignetting of the pupil by the secondary optics structure, constructed from industry-supplied CODE V<fnref n="21" /> model inputs and tested against CODE V model outputs.

- Optical path differences due to optical layout, modelled locally on the pupil plane as a weighted sum of Zernike polynomials. Zernike polynomials are a natural basis set for this contribution, since they are defined to be orthogonal on a circular pupil, although the vignetting breaks the orthogonality of the polynomial set in this application. Variation across the FPA is included as a polynomial fit across the FoV for each Zernike polynomial.

- The surface figure error (SFE) contribution to the wavefront error, resulting from imperfections of each reflecting surface, are included by propagating a model for the beam footprint on pre-launch measured SFE maps for each optical surface.

- Telescope jitter, or guiding error, of the telescope during an observation induces a field-dependent, achromatic convolutional effect (Ma et al. 2008). This convolutional kernel is calculated from three-axis measurements of the RPE from the FGS, Wiener filtered to remove noise.

- Detector charge diffusion is modelled as in Niemi et al. (2015), as a chromatic convolutional effect.

<footnote n="21">²¹ CODE V® is a commercial optical design tool.</footnote>
- Chromatic dependencies, including polarisation-dependent telescope throughput, detector quantum efficiency, and source SEDs, are included as spectral weightings in the broadband, chromatic PSF model, which is computed by spline interpolation of monochromatic calculations at multiple wavelengths, followed by integration across the VIS passband and out-of-band optical wavelength region.

- The chromatic response to polarisation from mirror coatings, most notably the layered dielectric coatings of the dichroic, the silver coating on FOM 3, the coatings on FOM 1 and 2 at high angles of incident light, and the hafnium oxide coating on the detector. In particular, the layered dichroic coating induces a strong chromatic and FoV dependence, which must be accounted for in the model. To characterise this dependence and to provide inputs for the PSF modelling, we have commenced an extensive test campaign on the spare model of the dichroic (Baron et al. 2022, 2023; Baron 2024). The polarisation dependence of the VIS optical system is modelled following Chipman et al. (2019, chapters 11, 16). An FoV- and wavelength-dependent set of four 'Jones pupils' is produced from the CODE V model of the optical surfaces and layout. The Jones pupils describe the system's polarised response to orthogonal input polarisation states across the pupil. They are included in the construction of the complex electric field at the exit pupil, and transformed to the monochromatic PSF at a given wavelength for a given set of Stokes parameters, via the construction of an 'amplitude response matrix' and Mueller point-spread matrix. We emphasise that this effect must be included even for unpolarised incident light, which should be considered as a combination of two incoherent orthogonal polarisation states.

- Detector undersampling. The model PSF is produced on oversampled pixels at the FPA, ensuring that the exit pupil is Nyquist sampled. This oversampled image may be down-sampled to the Euclid VIS detector pixel sampling for comparison with data.

The wavefront model for the *Euclid* VIS PSF allows for a fully flexible and modular parameterisation, including chromatic, spatial, and temporal variations. The SFE, dichroic, and detector contributions, with the exception of CTI, are expected to be invariant across the mission lifetime and therefore the models for these are fixed using laboratory measurements. The telescope pointing stability is measured during observations using the FGS, and propagated to the model for any given exposure. However, the optical model is expected to vary throughout the mission, due to thermal variations arising from changing telescope orientation with respect to the Sun and due to variations in spacecraft heat inputs. Distortions to the telescope structure under gravity and perturbations resulting from the launch process mean that laboratory measurements of the optical alignment cannot predict wavefront errors, and instead these must be inferred from in-orbit data.

The PSF model is therefore calibrated across the lifetime of the survey using a hierarchical calibration process consisting of three steps. First, the initial state of the PSF model is inferred from PV (Sect. 4.3.1) and early science operations (ESOP) data across a range of spacecraft attitudes with respect to the Sun vector. This is used to produce a basis set of spatially dependent PSF variations, which can be fit to further observations at other attitudes. Since PSF optical modes are expected to be degenerate with detector modes (e.g. the charge diffusion kernel), this process is conducted iteratively, alongside fits for detector modes, until requirements on validation metrics are met.

Because detector modes are not expected to vary across PV and ESOP calibration data, they are fit in a meta-analysis of the full calibration data set. This basis set is then updated from continuous observations with a monthly cadence, and any perturbations to the basis set modes are fed through to further optimisation. The monthly calibrations are taken on a limited set of defined fields, to simplify the interpretation of PSF variations between observations. This process should capture long-timescale variations. Finally, the resulting calibrated basis set is fit to each science observation, producing short-timescale variation on the timescale of science observation cadence. This enables an accurate PSF model to be produced for shear measurement in each science observation.

In the calibration, the initial telescope optical state is assumed to be maximally unknown, whereas the PSF model contains degenerate effects, both between optical Zernike modes and between detector and guiding error effects. Together, this can make it difficult for initial PSF calibration to robustly determine the underlying optical parameters that are used to infer the PSF across the mission lifetime.

We therefore utilise phase-diverse observations to calibrate the PSF model during the observation of PV and ESOP calibration images (see Section 4.3.2). The phase diversity technique (e.g. Gonsalves 1982; Fienup 1982) adds known aberration to the observations, and for the *Euclid* mission we add defocus by deliberate movement of the M2 mirror along the beam propagation direction. Adding a known defocus serves two purposes: first, it causes optical modes to become distinct from detector modes, due to dispersion of the PSF across a wider scale of pixels; and second, perturbations in optical modes induce differing impacts on the PSF image for each defocus, allowing a joint fit that lifts degeneracies between optical modes. Together, the phase-diverse calibration allows for the model parameter optimisation to robustly find the true best-fit model parameters, by simplifying the posterior that is optimised, through lifting degeneracy and multi-modality. These phase-diverse data are analysed with a model-fitting approach, to account for the wide VIS passband and the image undersampling (cf. Fienup 1999).

Figure 31 shows the pre-calibration VIS PSF model, fit to a PDC test data set taken during the PV phase. This comprised three sets of 10 exposures, taken at nominal in-focus and with two defocus offsets, applied as a movement of the M2 mirror by ±18 µm, corresponding to about 0.8 waves of defocus. Larger defocus amounts are prohibited by the requirement to operate the FGS. Also shown are stacked profiles of stars in that data set, applying pixel-level sigma-clipping to remove contaminating sources. This stacking procedure was carried out on both the stars and the model images fit to those stars. Red contours show isophotes for the model PSF. We see that the features of the PSF, as observed in the stars selected in the data, are well described by the PSF model. There is good agreement on all scales, including in the wings of the PSF, with good alignment of the diffraction spikes caused by the spiders holding the M2 mirror, including lateral translation of those spikes in the defocused images, and the halo caused by SFE on M1. We note that these models were fit to a set of phase-diverse test data with approximately 15 times fewer exposures than the PV and ESOP PDC calibration sets, and the models shown are not yet expected to satisfy the model accuracy requirements discussed above.

### 7.6.5. Shape measurement

Although *Euclid* is designed to provide sharp images of distant galaxies, the observed shapes are biased because of the
<figure><img src="image_33.png" /><figcaption>**Fig. 31.** Comparison between the mean profile of stars and the VIS PSF model during the PV PDC test data. The left panels show stacked images of stars selected from CCD 4-5, while the right panels show the corresponding stacked model PSFs. The top panels show in-focus data and models, while the bottom panels show defocused data and models, obtained with M2 movement of -18 µm. All images are produced as the mean of flux-selected stars, applying 2σ clipping on pixels to remove contaminating objects. Red contours are isophotes of the model PSF, superimposed on both the model and data images. The model is a pre-calibration one, fit to a small PDC test data set.</figcaption></figure>

convolution with the PSF, noise in the images, blending with other galaxies, the presence of cosmic rays, as well as a variety of detector-related sources of bias. These amount to biases that dwarf the lensing signal itself, and thus all need to be carefully accounted for.

Over the years much effort has been spent on developing and testing methods that aim to provide accurate ellipticities for individual galaxies (see e.g. Heymans et al. 2006; Massey et al. 2007; Bridle et al. 2010; Kitching et al. 2012a, 2013, 2012b; Mandelbaum et al. 2015, for the results from several community-wide efforts). The sensitivity of algorithms to certain systematics will differ, and hence it is advantageous to consider several approaches. A well established method is to fit parameterised models to the observed surface brightness. To this end, we have developed LensMC (Euclid Collaboration: Congedo et al. 2024), which is based on forward modelling galaxies with fast Markov chain Monte Carlo (MCMC) sampling and marginalisation over nuisance parameters. It measures the properties of objects from the multiple exposures jointly, while mitigating the bias due to detected neighbours by measuring objects in groups. The robustness has been proven on realistic simulations of EWS images, using galaxies drawn from the EFS as input. The development of LensMC has benefited from the knowledge gained from the ground with lensfit (Miller et al. 2007, 2013; Kitching et al. 2008a), which has been successfully applied to the various data

releases of KiDS (e.g. Fenech Conti et al. 2017; Kannawadi et al. 2019; Li et al. 2023a).

An alternative route uses the observed moments of the surface brightness. For instance, MomentsML (Tewes et al. 2019) is a machine-learning algorithm that trains a neural network to estimate the shape parameters from the moments of the surface-brightness profiles of galaxy images. In principle, this can reach the required accuracy, but it does rely critically on the realism of the simulated training data. This is a general concern for shape measurement (e.g. Hoekstra et al. 2017; Kannawadi et al. 2019), also because biases are already introduced during the object detection stage (e.g. Fenech Conti et al. 2017; Hoekstra et al. 2021). To (partially) circumvent these concerns, MetaCal (Sheldon & Huff 2017; Huff & Mandelbaum 2017) measures the shear response of an object for any shape measurement method by directly distorting its observed images. It has already been extensively applied to the DES lensing analysis (Zuntz et al. 2018; Gatti et al. 2021). Moreover, tests on realistic simulated images have shown its potential for Euclid (Hoekstra 2021; Kannawadi et al. 2021; Hoekstra et al. 2021). Hence, work is ongoing to develop a MetaCal setup for the analysis of Euclid data.

Simulated images are also needed to calibrate the biases introduced by other complications, such as the blending of galaxies at different redshifts (MacCrann et al. 2022;
Li et al. 2023a) or specific detector effects. The latter includes the non-uniformity of the pixel response over the CCD, bleeding due to charge overflow in a pixel, and CTI in the readout process. Other effects that need to be quantified are the various sources of background noise, variation in the star density, cosmic rays, galaxy blending, and contribution of unresolved galaxies to the background (Euclid Collaboration: Martinet et al. 2019). Many of these effects only become relevant because of the precision that Euclid can achieve. Some are very specific to space-based observations where radiation damage leads to CTI. The impact on the VIS images can be modelled following Massey et al. (2014), All these effects must be accurately included in the simulation pipeline described in Sect. 6.2 to reach a successful calibration (e.g. Hoekstra et al. 2015, 2017).

As already discussed in Sect. 2.2, the shear biases can be split into a multiplicative and additive contribution (see Eq. (17)), assuming the applied shear in the simulations is sufficiently small, or the algorithm sufficiently linear. In some cases it can be useful to isolate PSF leakage from the additive bias, adding the term $p_{\text{leak}} \epsilon^{\text{PSF}}$ to Eq. (17), where $p_{\text{leak}}$ quantifies the imperfect correction for the PSF ellipticity $\epsilon^{\text{PSF}}$.

Massey et al. (2013) showed that to reach our scientific goals (also see Sect. 2), the uncertainties in the multiplicative and additive shear bias must remain below $2 \times 10^{-3}$ and $2 \times 10^{-4}$, respectively, for each component of the shear. These numbers represent the total error budget, which then needs to be divided among several contributions (Cropper et al. 2013). Reaching such accuracy on shear measurement requires calibrating on realistic image simulations over a wide enough area so that measurement errors on these biases are negligible compared to the target accuracy. The required number of simulated galaxies, $N_g$, to reach a given uncertainty in the total multiplicative shear bias, $\sigma_m$, is given by (Fenech Conti et al. 2017; Euclid Collaboration: Congedo et al. 2024)

$$N_g = \left( \frac{\sigma_{\epsilon}}{\sigma_m |\gamma|} \right)^2, \qquad (24)$$

where $|\gamma|$ is the modulus of the shear applied to the simulated galaxies,<fnref n="22" /> and $\sigma_{\epsilon}$ is the dispersion of intrinsic ellipticities, taken to be $\sigma_{\epsilon} = 0.26$ in the case of Euclid (Euclid Collaboration: Martinet et al. 2019). However, this number can be lowered by shape (Massey et al. 2007) and pixel (Euclid Collaboration: Martinet et al. 2019) noise cancellation. In the former, pairs of galaxies with 90° rotation are considered, so that the mean intrinsic ellipticity reaches 0 in the absence of noise, and in the latter an extra identical pair with noise of opposite sign is included to lower the impact of shot noise in the calibration process. Together, these cancellations improve the runtime of the calibration simulations by a factor of 7 in the case of Euclid (Jansen et al. 2024).

Since galaxy morphologies depend on redshift, the calibration must be performed for each tomographic bin to avoid any undesired selection effects (e.g. Kannawadi et al. 2019). This process multiplies the number of nuisance parameters related to the calibration that need to be passed to the likelihood by the chosen number of tomographic slices. It also means that the shear and photometric redshift calibration must be performed jointly and that the galaxy morphology and photometry dependence on

redshift must be accurately reproduced in the calibration simulations (Li et al. 2023a). Finally, the chromaticity of the VIS PSF, combined with the SED variation across galaxy bulges and disks, introduces a non-negligible residual bias in shear measurements. This effect is referred to as ‘colour gradient’ bias (Voigt et al. 2012). Semboloni et al. (2013a) showed that it is possible to correct for it on average using estimates for the bias from HST imaging in two narrow filters (F606W and F814W), which has been confirmed to be sufficient for Euclid by Er et al. (2018).

## 7.7. Science-ready data products

The carefully calibrated data products that have been described in the previous sections can be used for scientific exploitation. The enhanced galaxy catalogues contain information about the redshifts and stellar masses of galaxies, as well as detailed morphological information and extensive photometry in all bands. To enable the various cosmological analyses, however, the catalogues are processed further. The cosmological information contained in the catalogues is compressed into the various two-point statistics that have been introduced in Sect. 2. In Sect. 7.7.1 we summarise the steps in the calculation of the clustering statistics of the spectroscopic clustering sample, while Sect. 7.7.2 describes the summary statistics needed for the 3x2pt analysis. Moreover, high-resolution maps of the projected mass distribution (Sect. 7.7.3) and a catalogue of clusters of galaxies with a well defined selection function (Sect. 7.7.4) are provided.

### 7.7.1. Spectroscopic clustering statistics

As discussed in Sect. 2.1, the statistical properties of the galaxy distribution can be quantified in terms of a set of clustering moments, which compress the cosmological and astrophysical information from galaxy surveys. The most fundamental measures of galaxy clustering are two-point statistics, which fully characterise the fluctuations in a Gaussian density field.

The starting point is the spectroscopic galaxy catalogue, described in Sect. 7.5. Since the spectroscopic measurements are made with a slitless spectrograph, it is necessary to apply a selection process to identify the galaxy sample with confident redshift measurements. The EWS will be primarily sensitive to luminous emission line galaxies in the redshift range $0.9 < z < 1.8$ with an $H\alpha$ line flux above the nominal limit of $2 \times 10^{-16} \text{ erg cm}^{-2} \text{ s}^{-1}$. The catalogue will be selected to maximise the galaxy number density, while keeping the fraction of spurious redshift measurements under control. It is important to account for potential systematic effects that can modulate the detection limit and remove genuine galaxies from the sample, or lead to errors in the redshift measurement that add interloper galaxies and reduce the purity. These inaccuracies in the selection have a direct impact on clustering statistics.

The redshift purity and completeness of the sample will be evaluated using the EDFs (Sect. 4.2.1), which enable a detailed assessment of the spectroscopic selection function across the EWS. The selection function will be characterised by a random catalogue of mock, unclustered objects that is constructed to closely trace the mean density of galaxies in the EWS and exhibit the same angular and radial selection effects. Constructing this random catalogue involves a forward-modelling process that relies on the noise level estimated in the NISP exposures and approximate end-to-end simulations of the spectroscopic measurement pipeline.

The resulting selected galaxy catalogue, and its corresponding random catalogue, are provided as data products. Their ‘reconstructed’ versions are also provided; they are obtained by applying a nonlinear transformation to the observed positions

<footnote n="22">²² The simulations may use a different shear for each galaxy to capture the impact of the blending of galaxies at different redshifts (MacCrann et al. 2022; Li et al. 2023a) or a constant shear across the entire scene (Kannawadi et al. 2019; Li et al. 2023a, e.g.).</footnote>
of the galaxies aimed at reducing the effect of the nonlinear evolution of cosmic structures and to enhance signal-to-noise on mildly nonlinear scales. Reconstructions can be performed either using the standard Zeldovich approximation approach (e.g. Padmanabhan et al. 2012) or with an efficient implementation of the cosmological least-action method (Sarpa et al. 2021), both of which provide excellent results at the BAO scale.

Original and reconstructed catalogues are used to compute the clustering statistics. The clustering moments can either be measured in configuration or Fourier space, and while they are theoretically equivalent, in practice their estimation can lead to different types of uncertainties and they can be differently impacted by systematic effects. To ensure accuracy, it is important to measure them in both spaces. In configuration space, we measure the spatial 2PCF of galaxies, which quantifies the excess probability of finding two objects at a given separation in a discrete sample of mass tracers with respect to a random Poisson sample. Estimating this quantity involves counting pairs, which can be computationally intensive. Therefore, we prioritise minimising memory allocation and maximising computational efficiency in our implementation of the estimator, while ensuring accuracy within sub-percent requirements.

Apparent deviations from statistical isotropy in clustering measurements depend on the geometry of the Universe and the growth rate of cosmic structures, breaking the circular symmetry of iso-correlation contours. To quantify these effects accurately, our estimator accounts for deviations by decomposing pair separation vectors into polar or Cartesian coordinate systems, which define the transverse and parallel components to the LoS within a local plane-parallel approximation. To achieve computational efficiency and minimise redundancy, we implement two parallel pair-counting algorithms, namely ‘chain-mesh’ and ‘kd-tree’ techniques, both exact in counting all object pairs within a pre-defined separation range and with varying performance levels in CPU time and memory usage at different scales (Marulli et al. 2016). Our estimator also incorporates the random-splitting technique to enhance computational speed (Keihänen et al. 2019) in the validation tests. This technique involves splitting the random sample, which typically contributes significantly to the computational load, into smaller subsets. The number of subsets can be determined by the user, finding a balance between computational efficiency and accuracy needs. The outputs of our estimator include the first five Legendre multipole moments, see Eq. (9), of the 2PCF, the anisotropic, 2-dimensional 2PCF in polar and Cartesian coordinates, the projected 2PCF, and data and random pair counts. Depending on the option selected by the user, the estimator computes the auto-correlation function of the objects in the catalogue or the cross-correlation of the objects in two different catalogues. Finally, the 2PCFs of the original and of the Zeldovich reconstructed catalogues will be computed.

Various effects render the observed distribution of galaxies non-Gaussian and thus transfer information to higher clustering moments. In order to enable extraction of that information, Euclid will additionally produce measurements of three-point statistics. In particular, we have developed an efficient estimator for the three-point correlation function (3PCF), which employs both the local spherical harmonics approximation technique proposed by Slepian & Eisenstein (2015) and the brute-force triplet-counting approach, which is mainly used for verification purposes. The three-point estimator takes the same input as the two-point estimator, but produces different outputs, namely the connected and reduced 3PCFs for all triangle configurations, as well as the corresponding data and random triplet counts.

The Fourier-space analogue of the two- and three-point correlation functions are the power spectrum and bispectrum, respectively. Their implementation within Euclid follows the methodology of the standard Feldman-Kaiser-Peacock (FKP; Feldman et al. 1994) estimator, which measures correlations of Fourier modes from a weighted fluctuation field. This field is constructed from the difference between the galaxy and random catalogues (see paragraph above), multiplied by weights, which are chosen such that they minimise the variance of the estimator and depend on the number density of the selected galaxy sample. The fluctuation field is then smoothed onto a regular grid that is sufficiently large to encompass the full galaxy catalogue, which enables the use of fast Fourier transforms (FFTs).

In the presence of RSDs, the clustering moments acquire a dependence on the LoS towards a given pair or triplet of galaxies. Measurements of the power spectrum and bispectrum must appropriately account for the variation of the LoS over the extent of the survey footprint, which led Yamamoto et al. (2006) to derive an extension of the FKP estimator (see also Szalay et al. 1998; Yoo & Seljak 2015; Scoccimarro 2015). However, since the Yamamoto estimator presents a computational challenge for data sets as large as those produced by Euclid, in practice one approximates the LoS by that of one of the galaxies in the pair or triplet. With this approximation it is still possible to extract the anisotropic signal using FFTs (Scoccimarro 2015; Bianchi et al. 2015), making the computation efficient, while consequences of this choice of LoS can be addressed at the level of the theoretical models (e.g. Beutler et al. 2019). The implementation for Euclid allows us to measure the first five Legendre multipoles of the power spectrum, Eq. (10), in addition to the monopole, quadrupole, and hexadecapole of the bispectrum.

The clustering two-point statistics moments estimated in Fourier space carry the imprint of the survey selection function (also called the survey window function), which must be accounted for when fitting theoretical models to the measurements. The required correction of the theoretical templates can be computed from the power spectrum of the random catalogue, which is therefore a separate data product of the Euclid pipeline.

### 7.7.2. Photometric 3×2pt statistics

Euclid’s photometric galaxy samples are determined by their detection in the $I_E$ filter. Nominally, this includes all galaxies with $I_E < 24.5$ (10 $\sigma$, extended source) and with valid photometry in all NIR and all available ground-based bands. The resulting samples and their key properties, such as number density and redshift distributions, are further modulated by statistical weights, which typically correlate strongly with S/N in $I_E$. Data processing enables the full galaxy sample to be split into up to 13 tomographic bins (Euclid Collaboration: Pocino et al. 2021) via point estimates for the photometric redshifts. The bin boundaries will be placed such that the statistical constraining power is optimised (see e.g. Zuntz et al. 2021; Euclid Collaboration: Pocino et al. 2021), while also keeping systematic trends and cross-correlation in the redshift distribution calibration to a minimum.

The galaxy samples used for shear measurement and for photometric galaxy clustering will generally differ, due to different statistical weights and the number and placement of tomographic bins, as well as potentially brighter cuts in $I_E$ to reduce spurious clustering patterns caused by spatially varying selection effects. These selection effects are not only caused by systematic variations in the $I_E$ detection efficiency, such as background noise level, foreground stellar density, or Galactic extinction, but also
by trends in all bands entering the photometric redshift estimation, including those from ground-based surveys. This makes it challenging to predict or simulate the full complement of selection effects. By default, we will therefore apply a data-driven approach that learns systematic galaxy density variations from the observed data, closely following methods applied to current imaging surveys (Johnston et al. 2021; Rodríguez-Monroy et al. 2022). Thus, we will construct visibility maps and, equivalently, random catalogues for each tomographic bin in the photometric galaxy clustering sample.

As discussed in Sects. 2.2 and 2.3, we will extract cosmological information from the 3×2pt statistics of the photometric galaxy catalogue. To do so, we measure the two-point correlations of:

1. galaxy ellipticities, tracing cosmic shear;

2. galaxy positions, tracing angular galaxy clustering; and

3. cross-correlations between galaxy positions and ellipticities, tracing galaxy-galaxy lensing.

The 3×2pt statistics are angular statistics, which are measured in projection on the sphere. To obtain information from the cosmological evolution of the galaxy sample, the 3×2pt statistics are computed for all auto- and cross-correlations of the full set of tomographic redshift bins.

We can measure angular statistics on the sphere either in real space, as observed, or in harmonic space, after a spherical harmonic transform of the observations. To maximise the scientific return of *Euclid*, we will measure the 3×2pt statistics in both real and harmonic space. While there is a mathematical relation such as Eq. (12) for each angular correlation function and its power spectrum, the transformation between real and harmonic space requires information on all angular scales, from the full sky to the infinitesimally small. In practice, we therefore cannot exactly transform one measured statistic into the other, even if both ultimately probe the same information, and we will hence measure them separately.

The main real-space measurement consists of the following estimates of the 3×2pt angular correlation functions:

1. the cosmic shear estimator $\xi_{\pm}(\theta)$ is a weighted average of the observed shears (Schneider et al. 2002a);

2. angular clustering is measured separately using estimators $\hat{w}(\theta)$ of Landy & Szalay (1993) and Hamilton (1993), together with catalogues of random positions tracing the survey footprint and systematics;

3. the galaxy-galaxy lensing estimators $\hat{\gamma}_t(\theta)$ and $\hat{\gamma}_\times(\theta)$ correlate galaxy positions of a 'lens' sample with shear estimates from a 'source' sample (see Joachimi et al. 2021; Prat et al. 2022 for recent discussions of these estimators).

Each real-space estimate is computed in linear or logarithmic bins of angular separation $\theta$, and for all combinations of tomographic redshift bins. A secondary set of derived two-point statistics is also measured, consisting of band-power spectra (Schneider et al. 2002a) and Complete Orthogonal Sets of E/B-Integrals (COSEBIs; Schneider et al. 2010), which are different linear combinations of the primary correlation function measurements.

The harmonic-space measurement consists of the following angular power spectra:

1. the cosmic shear angular power spectrum from pairs of galaxy ellipticity maps;

2. the angular clustering power spectrum from pairs of galaxy number density maps, together with visibility maps accounting for the survey footprint and systematics;

3. the clustering-shear cross-correlation power spectrum from combinations of galaxy ellipticity and number density maps.

These spectra are measured from partial-sky data that are observed only within the *Euclid* survey footprint, which imprints the respective angular selection function of each tomographic bin on the measurements (Brown et al. 2005). To account for this effect, every harmonic-space 2-point measurement comes with a so-called mixing matrix, which imprints the same angular selection of the survey on the full-sky angular power spectra obtained from theory.

Error estimates for the measured two-point statistics, in both real and harmonic space, are provided by a delete-one jackknife computation, where the available survey area is divided into a set of smaller regions, and all measurements are repeated while leaving out each region in turn. The resulting sample is then used to estimate the covariance of each set of 3×2pt statistics. Finally, all two-point measurements are repeated for all shape measurement methods (Sect. 7.6.5) that are used in a particular data release.

To illustrate the wealth of data that *Euclid's* 3×2pt statistics will deliver, we measure a synthetic harmonic-space data vector from the EFS. Galaxies from the simulation are selected using an approximate footprint for the northern part of *Euclid's* first data release (DR1), and binned into 13 photometric redshift bins with equal number density. The results are shown in Figs. 32–34.

### 7.7.3. Weak lensing convergence maps

The spatial correlations in the galaxy shapes provide direct information on the projected mass distribution along the LoS: in the weak lensing regime, the observed shear as a function of position can be used to reconstruct the corresponding convergence field. The result can be compared to the distribution of luminous matter, but the main use of these maps is to compute a range of statistics that provide additional cosmological information, beyond the two-point statistics discussed so far. These include the one-point probability distribution function, peak counts, and Minkowski functionals (for an extensive overview of possible estimators, see Euclid Collaboration: Ajani et al. 2023). Hence, convergence or mass maps have become a standard product of weak lensing surveys (e.g Massey et al. 2007; Van Waerbeke et al. 2013; Shan et al. 2014; Jeffrey et al. 2018; Oguri et al. 2018), and *Euclid* is no exception.

Given the substantial survey area of nearly 14 000 deg$^2$, mass mapping$^{23}$ can be performed either in the plane by dividing the *Euclid* survey into small fields, so that the flat-sky approximation remains valid, or directly on the curved sky. The latter, also known as spherical mass mapping, allows for large scales to be probed and is particularly suited for cross-correlations with other observables, such as CMB measurements, the distribution of galaxies, or all-sky cluster catalogues. However, the resolution of the spherical convergence maps and the complexity of the algorithm are limited by the computation time and memory required. As a consequence, planar mass mapping remains important to reconstruct convergence maps with a good resolution and precision, needed in particular to probe the non-Gaussian features of the weak lensing field, for example for higher-order statistics or to study the complex mass distribution in merging clusters (e.g. Clowe et al. 2006; Jee et al. 2012). Since the two methods are complementary, *Euclid* will provide convergence maps over the full survey area based on both approaches, as well as high-resolution maps of smaller patches around the most massive clusters of galaxies. The angular resolution of these maps

<footnote>$^{23}$ We refer to 'mass mapping' when reconstructing the convergence maps from the observed reduced-shear field.</footnote>
<figure><img src="image_34.png" /><figcaption>**Fig. 32.** Measured *E*-mode angular power spectra for cosmic shear from the galaxy ellipticity in the Flagship simulation, after applying the expected survey footprint of the northern part of *Euclid*'s first data release (DR1). Shown for each tomographic redshift bin (numbered panels) are the cosmic shear signal of that bin (black) and the cross-correlations with both lower-numbered (blue) and higher-numbered (orange) bins. The shading of the colour indicates the difference between the two bin numbers. Galaxy ellipticities have intrinsic ellipticity variations ('shape noise'), but no shape measurement error has been added here. Spectra are binned into 32 logarithmic bins. The *y*-axis changes to linear when crossing zero.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>**Fig. 33.** Similar to Fig. 32, but for angular clustering from galaxy positions.</figcaption></figure>
<figure><img src="image_36.png" /><figcaption>**Fig. 34.** Similar to Fig. 32, but for galaxy–galaxy lensing from the positions of galaxies and their ellipticity *E*-mode. Here, cross-correlations in each panel are shown for positions in that bin and foreground or background ellipticities. In harmonic space, the galaxy–galaxy lensing signal is negative; the apparent positive signal at higher redshifts is due to the intrinsic alignment of galaxies.</figcaption></figure>

will depend on the galaxy density and should not exceed 2' for the wide-field and 1' for the small-field maps.

Irrespective of the approach, the reconstruction of the convergence maps from the shear is a difficult task because of shape noise, irregular sampling, complex survey geometry, and the fact that the shear is not a direct observable. Taking these considerations into account, the *Euclid* data releases will include mass reconstructions based on two different planar inversion algorithms. The first method is the standard method proposed by Kaiser & Squires (1993). This Kaiser–Squires (KS) method has several shortcomings, but it is nevertheless commonly used. We therefore include it to allow cross-checks with previous measurements. The second algorithm employs the nonlinear inversion method KS+, described in Pires et al. (2020). This aims at reconstructing the convergence with minimal information loss, while controlling systematic effects. Details of the algorithms and an assessment of their performance are given in Pires et al. (2020). Extensions of these two methods to the curved sky have been implemented (e.g. Kansal 2023), and the results of those algorithms are also included in the *Euclid* data releases. A key application of the results is the study of higher-order statistics, which complement the constraints from the two-point correlation functions (see Sect. 9). For instance, Euclid Collaboration: Ajani et al. (2023) showed that various higher-order weak lensing statistics can tighten constraints on $\Omega_m$ and $\sigma_8$ by a factor of about 2, with prospects for further improvement.

### 7.7.4. Catalogue of clusters of galaxies

Another way to study the peaks in the matter distribution is to identify overdensities of galaxies. Such clusters of galaxies represent an extreme environment, affecting the star-formation

histories and morphological properties of galaxies (e.g. Dressler 1980), and as such can be used to test models of galaxy formation. Moreover, their abundance as a function of mass and redshift is sensitive to the underlying cosmological model (e.g. Allen et al. 2011; Kravtsov & Borgani 2012, also see Sect. 9.2). *Euclid* probes a large cosmological volume, and thanks to the deep NIR imaging it will extend the redshift range for cluster studies considerably, thus complementing large surveys that use the Sunyaev-Zeldovich effect (Planck Collaboration XXVII 2016; Bocquet et al. 2019; Hilton et al. 2021) and X-ray emission (e.g. Liu et al. 2022; Bulbul et al. 2024) to detect and study them.

The investigation of individual clusters, such as the interacting cluster 1E0657–558 (also known as the Bullet Cluster; Clowe et al. 2006), can provide useful insights, but most applications involve large samples that are used in statistical analyses. To correctly interpret the results, understanding how the sample is established is paramount. Moreover, the efficacy of detecting clusters in a multi-band imaging survey depends on a wide number of parameters in the algorithms employed.

Euclid Collaboration: Adam et al. (2019) explored this issue in detail by comparing the performance of various cluster-finding algorithms using synthetic *Euclid* data. The main criteria were their performance in terms of sample purity and completeness. This resulted in the selection of two algorithms that will be used to generate the catalogues of galaxy clusters. We describe their main features below. Importantly, the two codes employ rather different approaches to detect clusters, thus enabling a useful internal cross-check on any potential systematics associated with cluster detection.

The Adaptive Matched Identifier of Clustered Objects (AMICO; Maturi et al. 2005; Bellagamba et al. 2018; Maturi et al. 2019) employs an optimised matched-filter algorithm that can be
trained directly upon the survey data. The baseline cluster model incorporates a NFW radial density profile (Navarro et al. 1997) and a Schechter LF for the cluster members, although different models can be used as well. Clusters are iteratively identified within a 3D (angular position and redshift) significance map. This algorithm has already been used to find clusters in KiDS (Maturi et al. 2019) and the Javalambre-Physics of the Accelerating Universe Astrophysical Survey (miniJ-PAS; González Delgado et al. 2022; Maturi et al. 2023).

The second algorithm, PZWaveL makes minimal assumptions about expected cluster properties. This code uses a difference-of-Gaussians smoothing kernel to detect galaxy overdensities on the physical scale of galaxy clusters. Clusters are identified as statistically significant overdensities within a 3D data cube, with angular pixels with sizes of roughly 12" and bin widths in redshift that are appropriate for the photometric redshift PDFs. Versions of this code have been used for multiple surveys, including the IRAC Shallow Cluster Survey (ISCs; Eisenhardt et al. 2008), and a search in the S-PLUS fields (Werner et al. 2023).

Both codes have been validated and tests on simulated data indicate that they will yield samples with high purity and completeness (e.g. Euclid Collaboration: Adam et al. 2019). We recently used the Flagship simulation (version 2.1.10; Sect. 6.1), processed with Phosphoros, plus AMICO and PZWaveL. The analysis accounted for the uncertainties associated with the photometry, the photometric redshifts, the efficiency of the detection algorithms, and the intrinsic scatter of the cluster astrophysical properties. Based on these findings, we expect that Euclid will detect about $10^6$ clusters within the redshift range $0.1 < z < 2.0$, of which 30% are expected to lie at $z \ge 1$. The final sample selection will be based on the purity level estimated with tailored numerical simulations and the data-driven approach implemented in SinFoniA (Maturi et al. 2019). These results are in agreement with the earlier forecasts presented in Sartoris et al. (2016).

The pipeline for cluster finding will yield a merged cluster catalogue containing detections from both the AMICO and PZWaveL algorithms, which will be matched using both geometrical and membership-based procedures. Information from the detection (e.g. cluster position and significance) will be augmented with additional information from subsequent analyses for a detailed optical characterisation of the sample. This includes refined photometric and spectroscopic redshift estimates, as well as richness and membership estimates based upon the approaches of Castignani & Benoist (2016) and Andreon (2016). For a subset of clusters for which measurements are possible, the catalogue will also include weak lensing masses, cluster density profiles, luminosity functions, and velocity dispersions. Additionally, the size and the sensitivity of the EWS will allow the construction of a shear-selected cluster catalogue (see, e.g. Schneider 1996; Wittman et al. 2001; Schirmer et al. 2007) with more than $10^4$ clusters.

## 7.8. Data releases

The SGS processes Euclid data as soon as they become available. The data acquired on board are downloaded to Earth once a day (Sect. 3.3.1) and received at the mission operation centre (MOC), located at European Science Operations Centre (ESOC) in Darmstadt, Germany. It is then transferred to the Science Operations Centre (SOC) at European Space Astronomy Centre (ESAC) in Spain, where the raw instrument files, together

<figure><img src="image_37.png" /><figcaption>**Fig. 35.** Tentative timeline for public data releases, indicating the three main DRs as well as four smaller quick releases (Q1–Q4). The moment of release is linked to the start of early survey operations, but unforeseen changes to the mission operation may lead to some changes to this nominal schedule.</figcaption></figure>

with the telemetry of the spacecraft are turned into processable flexible image transport system (FITS) files, which are then ingested in the central SGS processing archive. The subsequent processing of the data is distributed among the 9 Science Data Centres (SDCs) which all run the same processing pipeline, but on different parts of the data.

This daily 'on-the-fly' processing uses the latest version of the pipeline and calibration products, and is used to monitor the health of the instrument and the quality of the data product produced by each of the processing steps. Once validated, these processed data are published in the ESA Science Archive. From that moment, these data are available to all members of the Euclid Consortium (EC). In principle, these data can be used for scientific investigation, but special care has to be taken, because the calibrations and processing functions continue to improve. Although such upgrades are integrated into the pipeline in a controlled manner, over time the data become increasing inhomogeneous. This may not pose a problem for many applications, but it can prevent robust high-level tests that are needed to validate the data for further cosmological analysis. To this end, about $500 \text{ deg}^2$ of EWS data are regularly reprocessed to ensure homogeneity and all data products are made available without blinding the summary statistics (see Sect. 8.2.4 about our blinding strategy).

The data processing is also homogenised for each of the three major public DRs. All DRs will occur through the ESA Science Archive, which is already on-line and currently gives access to the Early Release Observations data, presented in accompanying papers. A tentative release schedule, relative to the start of early survey operations, is shown in Fig. 35. The three DRs will contain the complete Euclid data set for roughly the first year of the survey, the first three years of the survey, and the complete survey, respectively. Each DR will comprise all the EWS data products detailed in the preceding sections, similar products obtained on the EDS during the period covered by the DR, as well as the associated calibration products. The SGS software version used for the DR will also be made public at the release time.

As Fig. 35 shows, the DRs are interspersed with so-called quick releases (Q) of smaller volume. At the time of writing, only the contents of Q1 have been planned. This release comprises data for a single visit over the EDFs: $20 \text{ deg}^2$ of the EDF-N, $10 \text{ deg}^2$ of EDF-F, and $23 \text{ deg}^2$ of the EDF-S. We aim to release the imaging and spectroscopic data, as well as catalogues and photometric redshifts.
## 8. Constraining cosmology with the primary probes

A major development in cosmology since the publication of Laureijs et al. (2011) and the selection of *Euclid* has been the release of the results from the analysis of *Planck* data, resulting in cosmological parameter constraints with high precision (Planck Collaboration XVI 2014; Planck Collaboration XIII 2016; Planck Collaboration VI 2020). Moreover, results from other cosmological probes have also continued to improve. As a result, the parameters that describe the concordance ΛCDM model are now well constrained, even though the physical nature of the main ingredients is yet to be explained. Interestingly, in recent years, some inconsistencies between probes have been claimed. For instance, local measurements of the Hubble constant (Riess et al. 2021) differ from the preferred value reported by Planck Collaboration VI (2020). Similarly, weak lensing studies find a lower amplitude for the lumpiness of matter, quantified by $S_8$ (e.g. Asgari et al. 2021). Whatever the origin of these differences, it is clear that the much smaller statistical uncertainties and reduced systematic effects of the *Euclid* results will have a significant impact in resolving the current debate and defining a new standard model of cosmology.

To exploit the impressive statistical power of *Euclid*, it is essential that all sources of systematic errors, whether of instrumental, astrophysical, or theoretical origin, are properly identified and their impact fully assessed. This includes assessing the choice of observables, but also the determination of the likelihood and exploring options for the modelling of various systematic effects. Doing so requires a comprehensive analysis of the performance of the mission, from the pixel-level data to the cosmological inference.

In Sect. 8.1 we present the general methodology that has been adopted to assess the impact of systematic errors and describe some of the tools that have been developed to ensure that cosmological constraints derived from *Euclid* are robust. These efforts are supported by extensive sanity and consistency checks of the data, which in turn drive some of the calibration needs. In Sect. 8.2 we summarise the work done to date to ensure robust parameter estimation, and provide updated predictions for the primary probes. In Sect. 8.3 we explore the prospects of constraining models beyond standard ΛCDM using *Euclid*.

### 8.1. Science performance verification

As discussed in Sect. 2, the science objectives that drove the design of the project are to derive precise constraints on the parameters $w_0$ and $w_a$, which allow us to probe the dynamical nature of dark energy, and to derive constraints on the growth-related $\gamma_g$ parameter. The performance of *Euclid* for the former is quantified by the FoM defined by Eq. (2). The choices for observables and requirements on residual systematic effects all derive from these core objectives.

It is often not clear how experimental design choices are related to changes in the precision and accuracy with which cosmological parameters can be determined. A first attempt was presented in Laureijs et al. (2011) in order to establish the feasibility of the mission. Since then, new insights have led to additional sources of bias that need to be accounted for, while others may now be better understood. Also, as design choices are made, requirements should be updated, as values derived from generic considerations tend to be conservative (for example, see the discussion in Euclid Collaboration: Paykari et al. 2020, for the weak lensing case).

To quantify how the inferred values of cosmological parameters from a set of observables depend on the adopted underlying cosmological model, uncertainties in our understanding of astrophysical sources of bias, or imperfections in the calibration of the instruments, we have developed tools that enable an end-to-end analysis. This allows us to link changes in the mission design to biases in the cosmological parameters of interest, which in turn helps to prioritise their importance, and to derive requirements for specific steps in the analysis pipelines. This is the main objective of the Science Performance Verification (SPV) exercise.

The basis for this SPV are the EFSs, which provide a sample of realistic galaxies in the redshift range of interest, for both the EWS and EDS. As described in Sect. 6.1, there are multiple purposes of such simulations. In the context of SPV, they are first used to define the set of fiducial values of the reference *Euclid* model, resulting in a common input catalogue for realistic synthetic observations. As discussed in Sect. 6.2, such simulated data are used to validate the performance of the elements of the pipeline in terms of purity, redshift precision, and photometric performance, as well as evaluating residuals in shape measurements. Such SPV studies also allow us to take global systematic effects into account. For instance, the star density and Galactic extinction varies across the survey. Similarly, radiation damage or ice build-up introduce large-scale variations that need to be quantified.

The performance of individual pipeline elements can be used to approximate the full end-to-end by a series of catalogue-level operations. Although this may not capture all co-dependencies, it provides a fast way to produce catalogues for large areas, while capturing the various sources of bias in a realistic fashion. As an example, biases in shape measurement are introduced by instrumental effects on the scale of the galaxy image. Hence, it is not necessary to create simulated images for a full *Euclid* survey. Instead, determining the biases as a function of relevant properties for a representative sample is sufficient. These dependencies can then be used to propagate biases at the catalogue level. For instance, Euclid Collaboration: Paykari et al. (2020) used this approach to explore the impact of spatially varying PSF and CTI residuals for the *Euclid* weak lensing measurements.

#### 8.1.1. Reference observables and nuisance parameters

As a preparation for the scientific exploitation of *Euclid*, an evaluation of its performance for the core science goals has been carried out that incorporates our current best knowledge of the mission. Here, we focus on the primary probes (described in Sect. 2) and limit the evaluation using the observables defined in harmonic space.

For the spectroscopic clustering measurements, we used the first three even multipoles of the galaxy clustering spectra $P_ℓ(k)$, defined in Eq. (10). As was done in Euclid Collaboration: Blanchard et al. (2020), we used four redshift bins, with mean redshift values $z \in \{1, 1.2, 1.4, 1.65\}$ and we used the same binning in $k$. Figure 36 shows the predicted signals for the four redshift bins, where the shaded indicate the uncertainties for the EWS based on the analytical covariance matrix (Sect. 8.1.2). Scale cuts in $k$, which need to be applied for the cosmological inference, were chosen to alleviate so-called projection effects (also known as prior volume effects that shift marginalised posterior distributions; Moretti et al. 2023). We used a maximum value of $k_{max} = 0.3 \, h \, \text{Mpc}^{-1}$ in the analysis.
<figure><img src="image_38.png" /><figcaption>**Fig. 36.** Legendre multipoles of the redshift-space power spectrum of galaxy clustering, $P_\ell(k)$, as expected from the spectroscopic survey data within four redshift bins (respectively, $0.9 < z < 1.1$, $1.1 < z < 1.3$, $1.3 < z < 1.5$ and $1.5 < z < 1.8$, where the $P_\ell(k)$ are evaluated at the mean of the redshift intervals). The plots show the monopole ($\ell = 0$, solid line), quadrupole ($\ell = 2$, dashed-dotted line), and hexadecapole ($\ell = 4$, dashed line), together with their error corridors (shaded regions). The latter simply connect the $1-\sigma$ errors from the diagonal values of the analytical covariance matrix, computed for narrow bins of $\Delta k = 0.0017 h \text{ Mpc}^{-1}$. As a result of this fine binning, the shaded areas do not fully reflect the actual constraining power of the measurements.</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>**Fig. 37.** Top: normalised redshift distributions $n(z)$, measured from the EFS, for the 13 equi-populated bins that were used for the 3×2pt analysis for the SPV. Middle: resulting photometric magnification kernels for the 13 redshift bins shown above. Bottom: corresponding shear kernels before (dashed) and after (solid lines) BNT transformation. The latter case gives a better grasp of the tomographic information that can be inferred from WL observations.</figcaption></figure>

The 3×2pt analysis combines the auto- and cross-angular spectra $C_\ell$ for the weak-lensing and photometric galaxy-clustering probes. We use 13 equi-populated redshift bins, which are presented in Fig. 37. These are based on photometric redshifts derived from the EFS assuming LSST-like external photometric data (Sect. 7.6.1). As such, these represent the best-case scenario, because the initial analysis will be based on shallower

ground-based data. In the following, we also assume that the samples used for weak lensing measurements and photometric clustering are the same. For a more exact treatment, an estimation of the joint errors in photometric redshift and multiplicative bias would be needed.

Although magnification does not improve cosmological constraints, we do need to account for it to avoid biased parameter estimates (e.g. Duncan et al. 2022; Mahony et al. 2022). The middle panel of Fig. 37 shows the magnification kernels, which highlight that magnification is particularly important for the highest redshift bins. The bottom panel shows the shear kernels. We also show results after applying the BNT transformation (Bernardeau et al. 2014; Taylor et al. 2018a), which allows us to control the mixing of scales from LoS projections. Figure 38 shows the synthetic angular power spectra for the weak lensing auto-correlations between the redshift bins, while Fig. 39 shows the corresponding photometric clustering signals. We do not show the cross-correlations, but they are included in the SPV analysis. For the binning in $\ell$-space we adopted a set of 32 log-spaced bins, ranging from [10, 5000], for the weak lensing and photometric galaxy clustering measurements. This ensures a manageable size of the data vector, but these choices may be adjusted for the actual data analysis. Cuts in $\ell$ are implemented on the observables in harmonic space after BNT transformation, so that $k > k_{max}$ accounts for less than 20% of each $C(\ell)$. For the 3×2pt analysis, we chose $\ell_{max} = 5000$ for cosmic shear and $\ell_{max} = 3000$ for angular clustering and galaxy-galaxy lensing.

The SPV aims to capture the impact of biases in the data or errors in the modelling itself. In the cosmological inference our ignorance is quantified by nuisance parameters that need to be used consistently, in particular when combining probes. The nuisance parameters can describe astrophysical quantities, such as galaxy biases, or aggregated quantities, such as uncertainties in the purity factor, the mean redshift, or residual multiplicative biases. The nuisance parameters used here are listed and described in Table 5. Furthermore, careful modelling needs to take place at nonlinear scales (see Sect. 8.2.2 for details), to ensure that the maximum information encoded at large $k$-values is used to exploit the full Euclid cosmological constraining power. Consequently, a number of key assumptions have been made, regarding the impact of baryonic feedback processes and the modelling of galaxy bias, magnification bias, interlopers, and
<figure><img src="image_40.png" /><figcaption>**Fig. 38.** Synthetic angular power spectra $C_\ell$ for weak lensing (EE) for the auto-correlation between the 13 photometric redshift bins shown in Fig. 37. The shaded light blue area shows the corresponding uncertainty given by the corresponding analytical covariance matrix, including the super-sample covariance (SSC) term.</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>**Fig. 39.** Similar to Fig. 38, but for the photometric galaxy clustering (gg) for the auto-correlation between the 13 photometric redshift bins shown in Fig. 37.</figcaption></figure>
<table><caption>**Table 5.** Reference values and prior probability distributions for the cosmological and nuisance parameters of the $w_0w_a$ CDM and $\Lambda$CDM + $\gamma_g$ models (adopting a flat geometry).</caption><thead><tr><th>Parameters</th><th></th><th>Fiducial value</th><th>Prior</th></tr></thead><tbody><tr><td colspan="4">Cosmology</td></tr><tr><td>Dimensionless Hubble constant</td><td>h</td><td>0.6737</td><td>U(0.55, 0.91)</td></tr><tr><td>Present-day physical baryon density</td><td>Ω<sub>b</sub>h<sup>2</sup></td><td>0.0227</td><td>N(0.0227, 0.00038)</td></tr><tr><td>Present-day physical cold dark matter density</td><td>Ω<sub>c</sub>h<sup>2</sup></td><td>0.1219</td><td>U(0.01, 0.37)</td></tr><tr><td>Dark energy equation-of-state parameters</td><td>{w<sub>0</sub>, w<sub>a</sub>}</td><td>{-1, 0}</td><td>{U(-3.0, -0.5), U(-3.0, 3.0)}</td></tr><tr><td>Slope primordial curvature power spectrum</td><td>n<sub>s</sub></td><td>0.966</td><td>U(0.87, 1.07)</td></tr><tr><td>Amplitude of the primordial curvature power spectrum</td><td>ln 10<sup>10</sup> A<sub>s</sub></td><td>3.04</td><td>U(1.6, 3.9)</td></tr><tr><td>Growth index</td><td>γ<sub>g</sub></td><td>0.545</td><td>U(0.01, 1.1)</td></tr><tr><td>Baryonic feedback efficiency factor of the HMCode emulator</td><td>log<sub>10</sub>(T<sub>AGN</sub>/K)</td><td>7.75</td><td>N(7.75, 0.17825)</td></tr><tr><td colspan="4">Photometric sample</td></tr><tr><td>Amplitude of intrinsic alignments</td><td>A<sub>IA</sub></td><td>0.16</td><td>U(-2, 2)</td></tr><tr><td>Power-law slope of intrinsic alignment redshift evolution</td><td>η<sub>IA</sub></td><td>1.66</td><td>U(0, 3.0)</td></tr><tr><td>Coefficients of cubic polynomial for clustering bias</td><td>b<sub>gal,i=0...3</sub></td><td>{1.33291, -0.72414, 1.01830, -0.14913}</td><td>U(-3, 3)</td></tr><tr><td>Coefficients of cubic polynomial for magnification bias</td><td>b<sub>mag,i=0...3</sub></td><td>{-1.50685, 1.35034, 0.08321, 0.04279}</td><td>U(-3, 3)</td></tr><tr><td>Per-bin shear multiplicative bias<sup>(‡)</sup></td><td>m<sub>i=1...13</sub></td><td>0.0</td><td>N(0, 0.0005)</td></tr><tr><td>Per-bin mean redshift shift</td><td>Δz<sub>i=1...13</sub></td><td>{-0.025749, 0.022716, -0.026032, 0.012594, 0.019285, 0.008326, 0.038207, 0.002732, 0.034066, 0.049479, 0.066490, 0.000815, 0.049070}</td><td>N[z<sup>fid</sup><sub>i</sub>, 0.002(1 + z<sup>fid</sup><sub>i</sub>)]</td></tr><tr><td colspan="4">Spectroscopic sample</td></tr><tr><td>Per-bin linear bias</td><td>b<sub>1,i=1...4</sub></td><td>{1.412, 1.769, 2.039, 2.496}</td><td>U(1.0, 3.0)</td></tr><tr><td>Per-bin second-order bias</td><td>b<sub>2,i=1...4</sub></td><td>{0.695, 0.870, 1.162, 2.010}</td><td>U(-5.0, 5.0)</td></tr><tr><td>Poissonian shot noise for extra-stochastic parameters</td><td>α<sub>P,i=1...4</sub></td><td>{0.056, 0152, 0.144, 0.309}</td><td>U(-1.0, 2.0)</td></tr><tr><td>Per-bin counter term for Legendre monopole</td><td>c<sub>0,i=1...4</sub></td><td>{11.603, 14.475, 15.667, 26.413}</td><td>Fixed</td></tr><tr><td>Per-bin counter term for Legendre quadrupole</td><td>c<sub>2,i=1...4</sub></td><td>{35.986, 44.914, 43.819, 62.353}</td><td>Fixed</td></tr><tr><td>Per-bin counter term for Legendre hexadecapole</td><td>c<sub>4,i=1...4</sub></td><td>{56.943, 55.443, 44.214, 42.89}</td><td>Fixed</td></tr><tr><td>Per-bin purity factor (assuming Poisson distributed interlopers)</td><td>f<sub>i=1...4</sub></td><td>{0.195, 0.204, 0.306, 0.121}</td><td>N(f<sup>fid</sup><sub>i</sub>, 0.01)</td></tr></tbody></table>

**Notes.** These fiducial values are used to compute the self-generated synthetic data used in Sect. 8.2.3. For the photometric nuisance parameters, we use a polynomial fitting formula for both the galaxy $b_{gal}$ and magnification $b_{mag}$ biases, whose coefficients run from $i = (0, 3)$, whereas we have a constant multiplicative bias $m_i$ and shifts in the bin redshift means $\Delta z_i$ per each of the 13 bins (fiducial values measured from EFS). For the spectroscopic nuisance parameters, we have one per each of the 4 redshift bins. The prior probability distributions are either uniform $U$ or Gaussian $N$. $^{(‡)}$ We assume a constant nuisance parameter for the multiplicative bias in each photometric redshift bin $i = 1...13$, which corresponds to the $m_{0}^{bias}$ parameter appearing in Eq. (17). The other components (i.e. $m_{4}^{bias}$) are assumed to be negligible.

intrinsic alignments, which all will have to be revisited once data are acquired.

We review the model predictions in Sect. 8.2.2, but briefly discuss our choice for the IA model here. For this study case, we adopt the so-called nonlinear linear alignment (NLA) model (Bridle & King 2007) to describe the scale dependence of the signal, while we assume that the amplitude of the signal

scales as

$$f_{\text{IA}}(z) = -A_{\text{IA}} C_{\text{IA}} \frac{\Omega_{\text{m}}}{g_{+}(z)} (1+z)^{\eta_{\text{IA}}} \quad (25)$$

where $g_+(z)$ is the nonlinear growth rate. We assume that the redshift dependence is described by a power law with slope $\eta_{\text{IA}}$.
The amplitude is quantified by the dimensionless parameter $A_{IA}$, which is scaled by the constant $C_{IA} = 5 \times 10^{-14} h^{-2} M_{\odot} M_{\text{pc}}^{-3}$, whose value is determined by the amplitude measured by Brown et al. (2002) in the low-redshift SuperCOSMOS survey (Hambly et al. 2001).

### 8.1.2. Data covariances

The SPV relies on an accurate estimate of the uncertainties. Hence, attention needs to be paid to the computation of covariance matrices. These are also essential to capture the correlations between probes, redshift bins, and angular scales. For the results presented in Sect. 8.2.1, we made a number of assumptions. First of all, we adopted a Gaussian distribution for the data vector, so that its statistical properties can be entirely characterised by the covariance matrix. Although it might be difficult to go beyond this assumption, we plan to assess its validity in future work. Second, we consider the spectroscopic clustering results to be independent of the 3x2pt measurements. This has been shown to be a reasonable approach for a Euclid-like mission (see Taylor & Markovič 2022). Finally, we evaluated the covariance matrix only for the fiducial model. As shown in Carron (2013), this is required in the Gaussian likelihood approximation to avoid the introduction of spurious parameter information. We assumed that the true model is not too far from the fiducial model (see, e.g. Harnois-Déraps et al. 2019 for a discussion of the impact of this assumption on cosmic shear constraints).

The computation of the covariance matrix of the observables is a challenging task. From a theoretical perspective, the computation of its Gaussian part is well understood. The difficulties arise from the fact that the observed modes are not statistically independent, but are coupled. One dominant source of mode coupling is the nonlinear growth of gravitational instabilities. The effect is particularly important on small scales for the observables in the 3x2pt analysis. The mode coupling through super-sample effects is another significant contributor (see e.g. Euclid Collaboration: Sciotti et al. 2024; Beauchamps et al. 2022). For the SPV analysis this term alone is included since we have found that this should give a realistic estimate. Another source of mode coupling is the impact of masks and visibility functions. When finite volumes and masks have to be taken into account, the harmonic components become coupled (Brown et al. 2005). This depends on the details of the survey, as well as the masking procedures. For now, only the size of the footprint has been taken into account in the SPV exercise – via a rescaling of the covariance by the sky fraction $f_{sky} = \Omega_{S}/4\pi$, with $\Omega_{S}$ the solid angle subtended by the survey, in steradians. A more comprehensive treatment is left for future analysis.

The multi-probe Gaussian covariance is given by (Euclid Collaboration: Blanchard et al. 2020):

$$ C_G[C_{ij}^{AB}(\ell), C_{kl}^{CD}(\ell')] = [(2\ell + 1) f_{sky} \Delta\ell]^{-1} \delta_{\ell\ell'}^{K} \\ \times \left\{ [C_{ik}^{AC}(\ell) + N_{ik}^{AC}(\ell)][C_{jl}^{BD}(\ell') + N_{jl}^{BD}(\ell')] \right. \\ \left. + [C_{il}^{AD}(\ell) + N_{il}^{AD}(\ell)][C_{jk}^{BC}(\ell') + N_{jk}^{BC}(\ell')] \right\}. \quad (26) $$

In the above equation, the Kronecker delta $\delta_{\ell\ell'}^{K}$ enforces the aforementioned independence of the different $\ell$ modes in the absence of convolution with the mask. The noise terms $N_{ij}^{AB}(\ell)$ are non-zero only for the auto-correlations between probes and

tomographic bins:

$$ N_{ij}^{AB}(\ell) = \begin{cases} \delta_{ij}^{K} \sigma_{\epsilon}^{2} / \bar{n}_{i}^{S} & A = B = L \\ 0 & A \neq B \\ \delta_{ij}^{K} / \bar{n}_{i}^{L} & A = B = G, \end{cases} \quad (27) $$

where $\sigma_{\epsilon}^{2}$ is the variance of the total intrinsic ellipticity of the sources. Finally, $\bar{n}_{i}^{X}(z)$ are the number densities of sources ($X = S$) and lenses ($X = L$), relevant for cosmic shear and photometric galaxy clustering respectively.

### 8.2. Expected cosmological parameter constraints

Comparing the Euclid measurements to model predictions is not straightforward, owing to the small statistical uncertainties, the need to marginalise over a large number of nuisance parameters (that are needed to quantify residual systematic effects), our limited knowledge of the nonlinear evolution of structure, and the impact of astrophysical processes on the matter distribution. Moreover, an accurate covariance matrix is needed so that all correlations between measurements can be correctly accounted for. Here, we provide an overview of the tools and procedures that have been developed in order to derive cosmological parameter estimates for Euclid. This updates some previously published forecasts presented in Euclid Collaboration: Blanchard et al. (2020).

#### 8.2.1. CLOE: the Cosmology Likelihood for Observables in Euclid

Given the unprecedented precision of the Euclid data, it is important that the comparison with theoretical predictions uses codes that have been tested rigorously. To this end, we have developed the Cosmology Likelihood for Observables in Euclid (CLOE), a highly flexible modular analysis pipeline written in python3. To ensure the fidelity of the results, the development of CLOE has combined the practices of continuous integration and delivery, enforcing automation in its construction, with careful unit testing and deployment of the code against similar pipelines (Euclid Consortium: Martinelli et al., in prep.).

As a baseline, CLOE provides the theoretical predictions for Euclid’s primary cosmological probes<fnref n="24" /> for a given set of cosmological and nuisance parameters. It computes the corresponding likelihood given the measurements (Sect. 7.7) and outputs the posterior probability distributions for the cosmological and nuisance parameters. It relies on the publicly available Boltzmann solvers CAMB (Lewis et al. 2000) and CLASS (Lesgourgues 2011) to compute the theoretical background parameters that are the foundation of the calculations of the primary observables. CLOE computes the predictions for the primary cosmological probes in both harmonic and real space.

A Bayesian approach is used to determine constraints on a given set of cosmological parameters, $\theta$, given the Euclid data vector, $\mathbf{d}$. According to Bayes’ theorem, the key ingredient in the estimation of the posterior distribution of the parameters, $P(\theta|\mathbf{d}, M)$, is the likelihood function $\mathcal{L}(\mathbf{d}|\theta, M)$, which describes the plausibility of a certain parameter value $\theta$, given a model

<footnote n="24">²⁴ The development of CLOE is open to the whole Euclid Consortium, allowing the merging of additional cosmological probes, such as cross-correlations with the CMB or cluster of galaxies.</footnote>
M, after observing a particular outcome. To sample the full posterior distributions of the cosmological parameters of interest, CLOE can be linked to the Bayesian analysis frameworks Cobaya (Torrado & Lewis 2021) and CosmoSIS (Zuntz et al. 2015) as external likelihoods. As a result, CLOE can employ a large number of different MCMC sampling algorithms, such as Metropolis-Hastings, classic nested sampling (e.g. PolyChord; Handley et al. 2015a,b) or advanced nested samplers (e.g. Nautilus; Lange 2023).

Generally, we assume that the likelihood probability distribution $L(d|\theta, M)$ of these measurements $d$ and the underlying physical model $M$, given the Euclid primary observables $t(\theta)$, is Gaussian$^{25}$, with a covariance matrix $C$ that does not depend on cosmology, so that, up to an additive constant,

$$-2 \log L(d|\theta, M) \approx [d - t(\theta)]^T C^{-1} [d - t(\theta)], \quad (28)$$

where $C^{-1}$ is the inverse of the covariance matrix and $t(\theta)$ is the theory vector constructed with the predictions for the Euclid primary observables assuming a cosmological model $M$. CLOE can compute non-Gaussian terms of the likelihood distribution $L(d|\theta, M)$ if the type of covariance matrix $C$ selected for the statistical analysis is numerical, according to Sellentin & Heavens (2016) and Percival et al. (2022).

Specifically, CLOE consists of a series of semi-autonomous python modules that interface with a Bayesian statistical framework tool, such as Cobaya, to read the relevant data vectors and covariance matrices, to compute the theory vectors (see Figs. 36, 38 and 39), and to calculate the likelihood. The modules that contain the relevant recipes that are needed to compute the theoretical predictions interface with another module that includes the modifications arising from nonlinear structure formation (Sect. 8.2.2). In an exercise of Open Science, CLOE participates in the pilot study case of ESA datalabs$^{26}$ (Navarro et al. 2024) as one the selected Euclid software pipelines to directly interface with the Euclid science archive system in the near future. As a demonstration of the capabilities, all the CLOE-related figures in this paper have been computed using ESA datalabs as the reference analysis framework.

### 8.2.2. Nonlinear structure formation

A major strength of CMB measurements is that the signal for a given combination of cosmological parameters can be computed directly, because the primary CMB fluctuations are in the linear regime. This is no longer the case for Euclid, unless the cosmological interpretation is restricted to rather large scales. To exploit the information contained in the smaller scales, nonlinear structure formation and the complexities of galaxy formation need to be taken into account. The challenge is to find an appropriate balance between the desire to minimise the statistical uncertainties, whilst ensuring that the predictions are robust. This involves defining the smallest scales that can reliably be used for a particular observable (e.g. Martinelli et al. 2021). The problem is not limited to the signal itself, but also involves robustly quantifying the covariance between measurements at different scales and the combination of observables that probe common structures. As a minimum, this implies a

coherent description of the primary probes of Euclid, ensuring that theoretical and astrophysical sources of bias are adequately accounted for. These efforts combine analytical and numerical calculations. Below, we summarise the different approaches that have been adopted to model the various summary statistics from Euclid. We do expect further developments, especially in the implementation of emulators (e.g. Eggemeier et al. 2022; Pellejero-Ibanez et al. 2022). As a first step, efforts have focused on analytical prescriptions and emulators that have been extensively validated with simulations and used in some of the most recent data analyses (e.g. Bose et al. 2020; Arnold et al. 2022; Carrilho et al. 2023; Piga et al. 2023).

For the spectroscopic galaxy clustering probe, we need to model the nonlinear galaxy power spectrum in redshift space. We will use a 1-loop perturbation theory model, with counterparts computed using the effective field theory (EFT) of LSS (Ivanov et al. 2020; d'Amico et al. 2020), which allows us to predict the clustering of matter in the mildly nonlinear regime. To this end, we have developed a new code (Moretti et al. 2023) using the FAST-PT algorithm for the loop-integral evaluations (McEwen et al. 2016; Fang et al. 2017). The code has been validated against N-body simulations (Oddo et al. 2021; Tsdrik et al. 2023, also see Fig. 2) and has been used for the re-analysis of BOSS data (Moretti et al. 2023; Carrilho et al. 2023).

For the analytical covariance, the prescription of Wadekar & Scoccimarro (2020) is used, including contributions from nonlinear effects, the window function, super-survey modes, and the integral constraint. As part of the analysis of BOSS DR12, this prescription has been shown to be in excellent agreement with the covariance estimated from over 2000 numerical mocks up to $k = 0.6 h Mpc^{-1}$, leading to negligible differences in cosmological parameters using either approach (Wadekar et al. 2020).

The modelling of the 3×2pt signals presents its own challenges, because of the desire to probe the matter power spectrum on nonlinear scales, and the need to account for baryonic feedback for the shear-shear and galaxy-shear correlations (Semboloni et al. 2011). In this case, two emulators, based on N-body simulations, have been developed and validated. The resulting Euclid Emulator (Euclid Collaboration: Knabenhans et al. 2019) and the bacco emulator (Angulo et al. 2021) can predict the nonlinear power spectrum with an accuracy of about 1% out to $k = 10 h Mpc^{-1}$ in the redshift range $0 < z < 3$, and include predictions for massive neutrinos and $w_0 w_a$ CDM. Ongoing development of these emulators aims to ensure that their accuracy is sufficient for the analysis of the Euclid data. For baryonic feedback effects, we have implemented the BCEmu (Giri & Schneider 2021) and bacco (Aricò et al. 2021) emulators. Other popular prescriptions like halofit (Smith et al. 2003; Takahashi et al. 2012) and HMCode (Mead et al. 2016, 2021) are also available (see Fig. 40). Ongoing simulated data challenges will determine the minimum number of baryonic feedback parameters needed for unbiased parameter inference and the associated scale cuts required, while we also plan to implement nonlinear and non-local bias models for photometric galaxy clustering with prescriptions based on either perturbation theory (Pandey et al. 2020) or hybrid emulators for biased tracers (Zennaro et al. 2023).

As alluded to in Sect. 2.2, we need to account for intrinsic alignments of galaxies (Joachimi et al. 2015; Troxel & Ishak 2015). To capture a wide range of possible alignments, we have implemented the tidal alignment and tidal torque (TATT) model (Blazek et al. 2019), used in the DES cosmic shear analyses (see e.g. Troxel et al. 2018). Progress will come from linking IA

<footnote>$^{25}$ The impact of non-Gaussian terms in the likelihood has been thoroughly studied, and found to depend on the range of scales employed. For our choice of scales for the angular power spectra of photometric probes, non-Gaussian terms have been demonstrated to be negligible for Euclid (Upham et al. 2021; Hall & Taylor 2022).</footnote>

<footnote>$^{26}$ https://datalabs.esa.int</footnote>
<figure><img src="image_42.png" /><figcaption>**Fig. 40.** Ratio of photometric galaxy clustering $C_\ell$ between different nonlinear models and the result for Euclid Emulator 2, for the autocorrelation of the redshift bin centred at $z = 0.83446$. Also shown is the expected Euclid error bar, including the contribution from super-sample covariance.</figcaption></figure>

models to observations, including those made by Euclid, so that the dependence of the IA signal on galaxy properties can be used to reduce the number of nuisance parameters (Fortuna et al. 2021).

Finally, to consistently model the correlations between scales, tomographic bins, and the different probes, we need an accurate covariance matrix for the 3×2pt measurements. Currently, several independent codes have been studied and validated, including the PySSC module (Lacasa & Grain 2019). Thus far, the focus has been on quantifying the impact of super-sample covariance, which is the largest of the expected non-Gaussian contributions (Barreira et al. 2018; Upham et al. 2022; Beauchamps et al. 2022; Euclid Collaboration: Sciotti et al. 2024).

We also want to explore models that include complexity beyond the baseline ΛCDM model (see Sect. 8.3). Although the range of possibilities is vast, we have focused the development on the main science objectives of Euclid. Hence, we have included nonlinear modelling prescriptions for wCDM, $\gamma_g$, and massive-neutrino cosmologies. Ongoing efforts include developing and implementing nonlinear models for a suite of exotic dark energy and modified gravity cases.

### 8.2.3. Parameter estimation

Since the estimates for the performance of Euclid were presented in Laureijs et al. (2011), the fidelity of the predictions has steadily improved. In particular, Euclid Collaboration: Blanchard et al. (2020) presented results from the first collaborative analysis to verify forecasting tools. This study focused on forecasts based on Fisher matrix techniques applied to both primary probes. A key aspect involved the comparison of different numerical implementations. The results showed optimistic and pessimistic scenarios for several cosmological models (flat and non-flat, and different cuts of the nonlinear scales), highlighting the role of the cross-correlations, especially for models beyond a cosmological

constant, potentially increasing the dark energy FoM by at least threefold.

In this section, we update the forecasts from Euclid Collaboration: Blanchard et al. (2020) using the primary probes only (see Sect. 8.1.1 for details about our setup). We limit the discussion to the baseline models of interest: a spatially flat cosmological model, using Eq. (1) to describe the dark energy equation of state ($w_0w_a$CDM); and a model where the parameter $\gamma_g$ is left free ($\Lambda$CDM + $\gamma_g$).

We generated synthetic noiseless data vectors by running CLOE v2.0.2 (Euclid Consortium: Joudaki et al., in prep.) using the fiducial values for the cosmological model parameters presented in Table 5. For this analysis we used the analytical super-sample covariance matrix for the 3×2pt observable (see Euclid Collaboration: Sciotti et al. 2024), and a Gaussian covariance matrix for the spectroscopic probe (GCsp), following the specifications of the SPV exercise for the survey area and other experimental systematics. For the correction of the matter power spectrum at nonlinear scales, we used the latest version of HMcode (Mead et al. 2021) for the photometric probes (allowing the baryonic feedback efficiency factor of the code emulator, $\log_{10}(T_{AGN}/\text{K})$, to be free) and a 1-loop perturbation theory model based on FAST-PT for the spectroscopic measurements, as described in Sect. 8.2.2.

The corresponding theoretical predictions and the calculation of the Euclid likelihood (see Eq. (28)) were computed using CLOE v2.0.2. The sampling of the posterior distributions were obtained using the novel nested sampler Nautilus, interfaced with Cobaya, imposing the priors for all the free cosmological and nuisance parameters presented in Table 5. For forecasting purposes, and to speed up the sampling process, we have used BBN information as a Gaussian prior for the baryon density parameter $\Omega_b h^2$ (Cooke et al. 2018). In this section, we show the results corresponding to three different forecasting cases: 3×2pt, GCsp and 3×2pt + GCsp (see Euclid Consortium: Cañas-Herrera et al., in prep., for a more complete discussion of the forecasts with CLOE). For the GCsp case, we need to sample a total of 23 free parameters, while we need 44 free parameters for the 3×2pt case, and 60 for the joint one. Moreover, we track several derived parameters on the fly, and simultaneously fit 364 (3×2pt) and 12 (GCsp) different spectra. During the analysis, we have fixed the per-bin counter terms (see Table 5) in the 1-loop perturbation theory model used for the GCsp probe, in an effort to mitigate the so-called projection effects. We thus implicitly assume that in future analyses we will be able to impose tighter priors on these model nuisance parameters. To achieve convergence, the parameter space is explored by calculating the Bayesian evidence using 4000 live points<sup>27</sup> in each iteration. To obtain the results presented here, we used approximately 0.2 million CPU hours in total.

Figure 41 shows constraints on the dark energy parameters $w_0$ and $w_a$, as well as the corresponding constraints on the other cosmological parameters, while Fig. 42 shows constraints on the $\gamma_g$ parameter. Regardless of the cosmological models, similar converged distributions are obtained for the nuisance parameters. For all the cases, the fiducial values are recovered for all the sampled parameters, with associated uncertainties that improve by one order of magnitude compared to current surveys (Abbott et al. 2022). The values obtain for the FoM of the dark energy parameters $w_0$ and $w_a$ are consistent with those of

<footnote n="27"><sup>27</sup> In nested sampling, a *live point* is a sample of the likelihood distribution as given by the prior, that is later used to construct the evidence.</footnote>
<figure><img src="image_43.png" /><figcaption>**Fig. 41.** Forecast of the constraints for the $w_0w_a$ CDM cosmological model (adopting a flat geometry) using only the *Euclid* primary probes, as described in Sect. 8.1.1. The sampled parameter space also included the cosmological parameters ($\Omega_bh^2$, $\Omega_ch^2$, $H_0$, $n_s$, $A_s$, $w_0$ and $w_a$) and several nuisance parameters listed in Table 5. The grey dashed lines show the fiducial values of the parameters, that are also listed in Table 5. The posterior distributions were obtained using CLOE v2.0.2 and the sampler Nautilus, with 4000 live points and 16 neural networks. For the photometric probes, we used $\ell_{max} = 5000$ for cosmic shear and $\ell_{max} = 3000$ for photometric angular clustering, and galaxy-galaxy-lensing, while for the spectroscopic probe we used, $k_{max} = 0.3 \ h\ Mpc^{-1}$. We show the 2D-posterior distribution for the parameters $w_0$ and $w_a$ in detail, citing the corresponding FoM obtained for each probe as well as for the combination of both.</figcaption></figure>

Euclid Collaboration: Blanchard et al. (2020), and they are obtained by marginalising over all the sampled parameters.

### 8.2.4. Blinding strategy

Cognitive bias of the scientists undertaking an experiment can lead to priors on data analyses linking together otherwise independent measurements. Blinding strategies are designed to broaden these priors by separating the influence of the scientists' predictions for the measurements from the results themselves. We aim to avoid such biases introduced by the way data are processed, selected, or modelled. This requires some care, because

the tendency to either consciously or unconsciously select, process, or interpret data in such a way as to confirm prior beliefs often leads to biased results (Nickerson 1998). The quest to determine cosmological parameters is no exception (Croft & Dailey 2011) because certain values for cosmological parameters may be preferred based on theoretical grounds.

Many recent cosmological analyses have taken steps to avoid cognitive bias by adopting a so-called ‘blinding’ strategy. This can take many forms, for example, by shifting values or theoretical models in plots (e.g. Wong et al. 2020), modifying the data vectors (e.g. Muir et al. 2020), or adjusting the covariance matrix (Sellentin 2020). The *Euclid* data are required to pass stringent validation tests, while the software that is employed to process
<figure><img src="image_44.png" /><figcaption>**Fig. 42.** Similar to Fig. 41, but for the ΛCDM + γg model (adopting a flat geometry). We show the 1D-posterior distribution for the γg parameter in detail, citing the corresponding 1-sigma uncertainty associated with each probe as well as for the combination of both.</figcaption></figure>

the data is under strict control. This distinguishes it from most
previous cosmological experiments, but there is still the potential
to introduce biases during the scientific analysis.

The blinding strategy we propose for *Euclid* starts with
the extensive use of synthetic data, such as those described
in Sect. 6.2, to build and test pipeline elements. We will then
allow the analysis of the first 500 deg² observed at the start of
the project without any blinding. This represents 1/30 of the
expected final sample. Assuming that variance scales with the
inverse of the volume, the results from the early data should have
an error 5.5 times larger than that from the final data and should
be comparable to current constraints. Thus we can consider that
analysing this sample is equivalent to blinding the signal at the
5.5 σ level. As the pipeline is developed and refined before the

first data release (DR1; Sect. 7.8), development will concentrate
on this sample, and it will be the only sample upon which revi-
sions in the pipeline are retrospectively applied until the DR1
sample is constructed.

As described in Sect. 4.3, the *Euclid* mission has defined a
careful and comprehensive calibration strategy. The observations
taken to facilitate the calibration tend to use the instrument in a
different way to the EWS, so that potential systematic problems
can be identified and quantified. Crucially, the calibration mea-
surements do not strongly depend on cosmology, and hence can
be analysed without additional blinding.

For the core cosmology measurements to come from DR1,
we adopt a strategy similar to that adopted by the DES (Muir
et al. 2020), but without forcing the shift to match between
galaxy clustering and weak lensing measurements. Independent
shifts mean that differences between statistics are blind: while
we cannot then use them to test our analysis methodology, the
impact of *Euclid* on any tensions will be blinded. We choose
to introduce blinding at the level of derived statistics rather than
raw data to facilitate calibration and validation tests and to ensure
we do not impact other analyses of these data.

Specifically, we will shift the two-point data products
(e.g. the correlation function and power spectrum multipole
moments) by the expected difference between two cosmolog-
ical models. This includes shifts in the effects of all of the
key cosmological processes to be measured. The model off-
set applied will be between the best-fit flat ΛCDM model of
Planck Collaboration VI (2020), and one randomly chosen
within a 3 σ interval.

8.3. Beyond ΛCDM models

As shown above, the primary probes of *Euclid* will place tight
constraints on the parameters of the ΛCDM model. However, as
already highlighted in Sect. 2, the data also greatly advance our
ability to explore various extensions to the standard cosmolog-
ical model, potentially shedding light on the dark constituents
of the Universe and the underlying theory of gravity. Here,
we highlight some of these cases, where we note that fur-
ther improvements can be achieved when combining the *Euclid*
results with complementary cosmological probes, such as the
ones summarised in Sect. 9.

The combination of the primary probes over the EWS area
offers the unique opportunity to test gravitational physics on
cosmological scales, placing constraints on several modifica-
tions of GR and models of dark energy (Amendola et al. 2018).
In addition to testing specific models, more agnostic descrip-
tions of modified gravity and dark energy can be explored. For
instance, the phenomenological functions µ<sub>mg</sub>(k, z) and Σ<sub>mg</sub>(k, z)
parameterise, respectively, the relation between the matter den-
sity contrast and the Newtonian and lensing potentials in Fourier
space (e.g. Pogosian et al. 2010), while the EFT of dark energy
(Frusciante & Perenon 2020) provides a framework to explore
deviations from GR consistently. Provided the observables can
be modelled accurately on nonlinear scales, *Euclid* will pro-
vide outstanding constraints on extensions of ΛCDM, espe-
cially when cross-correlations with CMB measurements are also
included (see Sect. 9.4; Casas et al. 2023; Frusciante et al.
2024).

Massive neutrinos suppress the matter power spectrum on
small scales. Hence, a measurement of this subtle signature
allows the sum of neutrino masses, Σm<sub>ν</sub>, to be constrained
using cosmological data. Current CMB and LSS measure-
ments provide stringent upper bounds (Σm<sub>ν</sub> ≤ 0.1 eV, 95% confidence; Planck Collaboration VI 2020; Alam et al. 2021; Palanque-Delabrouille et al. 2020), well below the limits of cur-
rent laboratory experiments studying β decay (upper limit on
effective electron anti-neutrino mass m<sub>ν</sub> < 0.8 eV, 90% confi-
dence, Aker et al. 2022). The absolute neutrino mass sum is
still unknown, but this situation will change thanks to *Euclid*.
In the minimal ΛCDM+Σm<sub>ν</sub> model, *Euclid*'s primary probes
alone can constrain the neutrino mass with a 1 σ error σ(Σm<sub>ν</sub>) =
0.05 eV for a fiducial neutrino mass of 0.06 eV. In combina-
tion with *Planck*, a precision of σ(Σm<sub>ν</sub>) = 0.02 eV can be
reached, implying a 3 σ detection of a non-zero neutrino mass
(Euclid Collaboration: Archidiacono et al. 2024). If the true neu-
trino mass is below 0.08 eV, these constraints will also provide

evidence in favour of the normal neutrino mass ordering<sup>28</sup>. We
note that these estimates are conservative because they assumed
the pessimistic scenario of Euclid Collaboration: Blanchard et al.
(2020). Ongoing efforts to improve the modelling of nonlin-
ear structure formation (Sect. 8.2.2) should ultimately result in
smaller uncertainties.

*Euclid* will also improve the constraints on the effective number of relativistic degrees of freedom $N_{\text{eff}}$, which accounts both for the number of standard model neutrinos ($N_{\text{eff}}^{\text{SM}} = 3.044$, Frouste et al. 2020; Bennett et al. 2021) and for additional light particles, dubbed ‘dark radiation’. The sensitivity of *Euclid*, in combination with current and future CMB surveys, to $N_{\text{eff}}$ will potentially exclude several theoretically well-motivated particles beyond the standard model (Euclid Collaboration: Archidiacono et al. 2024). Therefore, *Euclid* will also shed light on dark matter models predicting a deviation of $N_{\text{eff}}$ from the standard model value, such as models involving interactions between dark matter and dark radiation. Moreover, *Euclid* will inform the models for alternative dark matter scenarios beyond the cold dark matter paradigm, by improving the constraints on warm dark matter, and decaying dark matter (Euclid Collaboration: Lesgourges et al., in prep.).

Current observations are consistent with initial conditions
that correspond to a flat Universe with nearly Gaussian adiabatic
perturbations, whose spectrum is described by a simple power
law (Planck Collaboration I 2020; Alam et al. 2021). Thanks to
the different sensitivity of the primary probes to the expansion
of the Universe and to the growth of structure, *Euclid* provides
an invaluable snapshot of the initial conditions at low redshift.
The measurements will reduce the uncertainty on the spatial cur-
vature Ω<sub>K</sub> an order of magnitude below the current constraints
from galaxy surveys (Alam et al. 2021). A similar improvement
is expected for the uncertainty in the value of the scalar spectral
index n<sub>s</sub> and its running α<sub>s</sub> = dn<sub>s</sub>/d ln k. *Euclid* will also improve
the constraints on features in the primordial power spectrum as
forecast in Ballardini et al. (2024).

*Euclid* will test the statistics of primordial fluctuations beyond the power spectrum. The spectroscopic survey is expected to improve constraints on $f_{\text{NL}}^{\text{local}}$, the local shape of primordial non-Gaussianity, by approximately a factor 8 over current results (Mueller et al. 2022), reaching an uncertainty of about 3.4 when combining power spectrum and bispectrum information and assuming universality for the halo mass function. These uncertainties on initial conditions are comparable to those obtained by Planck (Planck Collaboration X 2020; Planck Collaboration IX 2020), but target a markedly different range in redshift and scale.

9. Additional cosmological probes

*Euclid* is designed with the primary probes in mind, but the data
enable a wide range of additional measurements that can improve
cosmological parameter constraints (Laureijs et al. 2011). For
instance, the cosmological information is not limited to the two-
point statistics that we have focused on so far. In Sect. 9.1
we discuss how higher-order clustering and lensing statistics
can be used to improve cosmological parameter constraints. In
Sect. 7.7.4 we already highlighted the large number of clusters
that *Euclid* will discover. Their use to improve cosmological
parameter constraints is reviewed in Sect. 9.2. As discussed in
Sect. 9.3, the sharp imaging data are ideal for the discovery

<footnote n="28"><sup>28</sup> The minimum mass allowed by neutrino oscillation experiments in normal (inverted) ordering is 0.058 eV (0.100 eV).</footnote>
of strong gravitational lenses, which enable a unique study of the distribution of dark matter on small scales, as well as additional tests of cosmology. The wide area galaxy and matter maps (see Sect. 7.7.3) can be cross-correlated with measurements of the CMB, enabling new probes that are presented in Sect. 9.4. Finally, high-redshift quasars with X-ray data can complement low-redshift cosmological probes in the determination of the cosmological parameters as discussed in Sect. 9.5, while Sect. 9.6 explores the use of passive galaxies as chronometers to provide an independent constraint on the expansion history.

## 9.1. Higher-order statistics

While the two-point statistics would capture all cosmological information in the LSS if it were Gaussian, nonlinear structure formation has introduced non-Gaussian features into the cosmic matter distribution. The full information content can, therefore, only be unlocked with higher-order statistics (HOS). A wide variety of observables that capture the higher-order information have been proposed, which can be roughly grouped into two categories: those that consider N-point correlation functions and Nth-order moments of the density distribution; and those that use topological information of the density distribution. Examples for the first category are higher-order moments (e.g. Gatti et al. 2022; Porth & Smith 2021), higher-order correlation functions (e.g. Heydenreich et al. 2023; Burger et al. 2024) and one-point probability distributions (e.g. Barthelemy et al. 2020; Boyle et al. 2021). The second category includes peak statistics (e.g. Martinet et al. 2018; Harnois-Déraps et al. 2021), Minkowski functionals and persistent homology (e.g. Parroni et al. 2020; Heydenreich et al. 2022), and scattering transforms (e.g. Cheng et al. 2020; Cheng & Ménard 2021). Many estimates can be inferred from WL convergence maps (see Sect. 7.7.3), while some, for example, higher-order correlation functions and aperture mass moments, can be directly measured from shear catalogues (Jarvis et al. 2004; Secco et al. 2022; Porth et al. 2024).

When combined with two-point statistics, HOS enhance cosmological constraints by (partially) resolving parameter degeneracies (e.g. Kayo & Takada 2013; Heydenreich et al. 2023). Euclid Collaboration: Ajani et al. (2023) found that combining each of ten different WL HOS with two-point statistics results in a twofold improvement in constraining $\Omega_m$ and $\sigma_8$ compared to relying solely on two-point statistics. Combining all HOS leads to a factor of about 4.5 improvement.

The power of HOS is illustrated in Fig. 43, which shows expected constraints from a Fisher forecast analysis on $\sigma_8$ and $w_0$ with all other cosmological parameters fixed, using the shear correlation functions $\xi_\pm$, the convergence PDF, or the combination of both. The covariance for the analysis is estimated from the Scinet LiGhtCone Simulations (SLICS) N-body simulations (Harnois-Déraps et al. 2018) and the derivatives of the data vectors are taken either from theoretical predictions (see Boyle et al. 2021 for details on the PDF modelling) or from the DUSTGRAIN-pathfinder simulations (Giocoli et al. 2018), where for the simulations we used the Kaiser-Squires mass reconstruction scheme (see Sect. 7.7.3). The PDF shows a different degeneracy direction between $w_0$ and $\sigma_8$ than the second-order statistics, illustrated by the tilted ellipse in the lower-left corner of Fig. 43. This change in the degeneracy leads to a tightening of the constraints on $\sigma_8$ and $w_0$ when $\xi_\pm$ and the PDF are combined. Consequently, the HOS carry additional cosmological information, which needs to be included to unlock all of Euclid's potential.

<figure><img src="image_45.png" /><figcaption>**Fig. 43.** Constraints on $\sigma_8$ and $w_0$ from a Fisher analysis of $\xi_\pm$ and the convergence PDF, when keeping all other cosmological parameters fixed, normalised by the constraints of second-order statistics alone. We assumed a Euclid-like source redshift distribution to derive the results. The $\xi_+$ and $\xi_-$ values were taken in the range of 1.65 to 201'. The PDF was measured for convergence fields smoothed by a tophat filter of radius 4'.69. Covariances were estimated from the SLICS (Harnois-Déraps et al. 2018), derivatives were either modelled analytically (dashed lines) or estimated from the DUSTGRAIN-pathfinder simulations (Giocoli et al. 2018, solid lines).</figcaption></figure>

HOS can also test for residual systematics and constrain astrophysical effects such as intrinsic alignment (Pyne & Joachimi 2021), baryonic feedback (Semboloni et al. 2013b), or galaxy bias (Huterer et al. 2006). Since the HOS react differently to these effects than two-point statistics, combined analyses allow us to simultaneously constrain cosmological and nuisance parameters without additional data sets.

## 9.2. Clusters of galaxies

Galaxy clusters have long proven to be a valuable cosmological tool (e.g. Bahcall & Fan 1998; Borgani et al. 2001; Haiman et al. 2001; Weller et al. 2002; Reiprich & Böhringer 2002; Vikhlinin et al. 2009; Mantz et al. 2015; Planck Collaboration XXIV 2016; Bocquet et al. 2019; Costanzi et al. 2019; DES Collaboration 2020; Lesci et al. 2022; Chiu et al. 2023; Ghirardini et al. 2024). Arising from the highest density peaks of the initial matter density field, the abundance and spatial distribution of clusters contains information on the growth of structures and expansion history of the Universe (Allen et al. 2011; Kravtsov & Borgani 2012). More specifically, the cluster abundance is a sensitive probe of the parameter $S_8$. The evolution of cluster counts, which effectively measures the growth rate of cosmic structure, constrains dark energy and modified gravity models (e.g. Mantz et al. 2015; Cataneo et al. 2015; Bocquet et al. 2019). The galaxy cluster correlation function, probing the same matter field traced by galaxies, is sensitive to the same cosmological effects and parameters detailed in Sect. 2.1. The relative lower S/N of the cluster clustering measurement, due to the sparser nature of the cluster sample, is mostly compensated by a theoretically predictable halo bias (e.g. Sartoris et al. 2012). The combination of cluster counts and clustering has the potential to deliver
independent, competitive, and complementary cosmological
constraints to those provided by the primary *Euclid* probes, but
depends sensitively on our ability to accurately calibrate their
masses (Sartoris et al. 2016; Fumagalli et al. 2024). Indeed,
while it is possible to predict with percent level accuracy the
abundance of dark matter halos as a function of mass and red-
shift in an arbitrary cosmology (e.g. Tinker et al. 2008; Castro
et al. 2021), halo masses themselves are not directly observable.
In cluster cosmology studies, it is hence crucial to identify and
calibrate observational proxies – such as the number of member
galaxies (richness), X-ray luminosity, or the Sunyaev–Zeldovich
(SZ) signal – against mass. At present, the calibration of these
scaling relations represents the main limiting factor for cluster
cosmology studies at all wavelengths (e.g Planck Collaboration
XXIV 2016; DES Collaboration 2020).

The combination of *Euclid*'s wide sky coverage and high-quality optical data will allow the detection of an order of $10^6$ clusters above ~ $10^{14} M_{\odot}$ out to redshift 2, increasing the number of detected systems by more than an order of magnitude compared to current surveys (Sartoris et al. 2016). In particular, *Euclid* will perform a census of the cluster population above $z \sim 1$ for the first time in the optical-NIR wavelength regime, a critical stage of the Universe's evolution for studying dark energy. The unprecedented large statistics, along with the large volume probed by the survey, will allow *Euclid* to beat down shot noise and sample variance, enabling a statistically significant measurement of the cluster correlation function at $z \gtrsim 0.3$. Along with the exquisite imaging data provided by *Euclid*'s space observations, this will enable a weak lensing calibration of the observable-mass relation out to $z \approx 1$ (Köhlinger et al. 2015), while the sparse spectroscopic data for cluster-member galaxies at $0.9 \lesssim z \lesssim 1.8$ will provide a valuable mass proxy for high-redshift systems (Sartoris et al. 2016).

The large number of clusters and the precision of the lensing
measurements demand stricter control over systematic effects
compared to current cluster surveys. This is especially true
for a photometric cluster survey such as *Euclid*'s, capable of
detecting systems down to group mass scales and becoming
mass-complete above $M > \times 10^{14} M_{\odot}$ (see Sect. 7.7.4). These
systems, which outnumber their more massive descendants by
orders of magnitude, potentially encode valuable cosmological
information, but prove to be difficult to include in abundance
studies (DES Collaboration 2020); the limited resolution that a
photometric cluster-finder algorithm can achieve along the LoS
leads to unavoidable uncertainties and biases in the richness esti-
mate, which become more severe in the low-S/N regime. The
correlation of these systematics with others affecting the lens-
ing measurements, or dynamical mass proxies, further hampers
the characterisation of these systems (Sunayama et al. 2020; Wu
et al. 2022). For *Euclid*, the calibration of such selection effects,
as well as the determination of the threshold for the minimum
cluster richness, will be tackled using a combination of simu-
lated and multi-wavelength data analyses (e.g. Costanzi et al.
2019, 2021; Grandis et al. 2021), along with the inclusion of
the clustering of clusters statistics. The latter, thanks to the mass
dependence of the halo bias, will enable us to break the degener-
acy between cosmological and scaling-relation parameters when
combined with the other cluster observables (e.g. To et al. 2021).
In summary, from the combination of cluster counts, cluster
clustering, and *Euclid*'s mass-proxies, we expect to increase
the precision on the estimation of $S_8$ by an order of magnitude
compared to current galaxy cluster studies based on imaging sur-
veys, such as DES and KiDS, or SZ surveys, such as Atacama
Cosmology Telescope (ACT) and South Pole Telescope (SPT).

<figure><img src="image_46.png" /><figcaption>**Fig. 44.** Simulated *Euclid* observation in the $I_E$ band of the central region of the strong lensing galaxy cluster MACSJ0416.1−2403 ($z = 0.397$, Balestra et al. 2016). The image was obtained with the code Hst2Euclid (Bergamini et al., in prep.), using HST observations taken as part of the Hubble Frontier Fields Survey (Lotz et al. 2017). The image reproduces the depth of the EWS and several giant arcs are clearly visible. The inset shows a zoom into a known galaxy-galaxy strong lensing system, where the lens is a cluster member and the source a background galaxy at redshift $z = 3.222$ (ID14, Vanzella et al. 2017b).</figcaption></figure>

Moreover, *Euclid*'s ability to sample the cluster population at
$z > 1$ will ensure an improvement of a factor of 2 on the DE
FoM, compared to cluster surveys at lower redshift (e.g. Sartoris
et al. 2016; Bocquet et al. 2019).

Furthermore, the inclusion of *Euclid* weak lensing data
is expected to dramatically improve the cosmological con-
straints derived from intra-cluster-medium-selected cluster sam-
ples, such as the ones provided by the extended ROentgen Survey
with and Imaging Telescope Array (eROSITA; Bulbul et al.
2024), or the high-resolution SZ surveys conducted by SPT (e.g.
Bleem et al. 2015; Benson et al. 2014) and ACT (Hilton et al.
2021).

9.3. *Strong gravitational lensing*

Multiple images of a distant source can be produced if light
rays pass sufficiently close to a massive structure along the
LoS. Such favourable configurations are rare, and typically are
not resolved in ground-based observations. Thanks to the sharp
images provided by *Euclid*, vast numbers of strong lenses will
be discovered, and we highlight some of the main applications
here. Strong lensing studies with *Euclid* will cover a wide range
of mass, from galaxies to groups and clusters of galaxies (see
e.g. Fig. 44), while probing the distribution of dark matter on
relatively small scales.

We expect to detect approximately 200 000 galaxy-scale
strong lenses, which thus make up the largest fraction of systems.
These typically feature an $M^*$ lens at a median redshift $z_d = 0.6$
within the EWS. The sources in these systems are expected to
have redshifts ranging from $1 < z_s < 2$ (Collett 2015; Metcalf
et al. 2019; Euclid Collaboration: Leuzzi et al. 2024), resulting
in an average Einstein radius of $\theta_E = 0.5'$, which aligns well
with the capabilities of *Euclid*, but is too small for ground-based surveys (e.g. Petrillo et al. 2019). In the EDS, approximately 3500 lenses are expected, compared to 500 in the equivalent area of the EWS. On average, the unlensed $I_E$ magnitude of the sources is 2 magnitudes fainter than that of the lensed ones, highlighting the significant magnifying effect of lensing (Lotz et al. 2017). This magnification enables the search for ultra-high redshift (lensed) sources up to $z = 8$ in the EDS, down to an unlensed $I_E$ magnitude of 28. This approach is critical for investigating questions about re-ionisation by distant protogalaxies. As *Euclid* targets the Ly α line and benefits from magnifications up to $\mu = 30$ (Mason et al. 2015; Atek et al. 2015; Vanzella et al. 2017a), we anticipate the discovery of 35 to 75 such galaxies in the EDS, while brighter and rarer events may be discovered in the EWS. Moreover, in spite of the lower NISP resolution, in the EWS we expect to find tens of thousands of bright lensed NIR sources in the NISP photometric data that are too faint to be detected in the $I_E$ band (Pearson et al. 2024). Most of these will be sub-millimetre galaxies (at $z \sim 2-4$) or evolved red galaxies.

We expect the identification of up to 2300 lensed quasars, with approximately 16% predicted to be quadruply imaged (Oguri & Marshall 2010). The discovery of lensed supernovae (e.g. Pascale et al. 2025) will complement the lensing of quasars; due to the EDS observations being divided into roughly 40 independent epochs over the six-year mission duration, the detection of several lensed type-Ia supernovae is expected. Both lensed quasars and supernovae will facilitate time-delay cosmography (Treu et al. 2022), allowing for measurements of the Hubble constant (Refsdal 1964; Kelly et al. 2015; Grillo et al. 2018; Wong et al. 2020). This process requires adding the time dimension to the precise *Euclid* data, demonstrating a significant synergy with LSST.

Moreover, microlensing in lensed quasars can be utilised to estimate the fraction of dark matter in compact form in galaxies (Mao & Schneider 1998; Moustakas & Metcalf 2003; Koopmans 2005; Vegetti & Koopmans 2009; Vegetti et al. 2012; Hezaveh et al. 2016; Nierenberg et al. 2017; Gilman et al. 2020; Wagner-Carena et al. 2023; Powell et al. 2023). Analysing the flux ratios between lensed images or their distribution across the sample of lensed quasars provides insights into the quantity and distribution of low-mass dark matter halos in lensing galaxies and along the line of sight (LoS). High spatial resolution follow-up observations, using very-long-baseline interferometry or adaptive optics at large ground-based optical telescopes, will facilitate the detection of such halos through gravitational imaging.

Furthermore, compound lenses, where multiple sources at different redshifts are lensed by the same foreground galaxy, are instrumental in overcoming the mass-sheet degeneracy inherent in lensing. This enables measurements of both the mass slope in lensing galaxies and the distance ratios between lenses and sources. Compound lenses serve as a robust tool for both mass profile measurement and cosmography, with additional sensitivity to dark energy parameters (Gavazzi et al. 2008; Collett & Bacon 2016; Sharma et al. 2023), provided that the multi-lens-plane mass-sheet degeneracy is broken (Schneider 2014).

*Euclid* is set to observe strong lensing features, such as families of multiple images and giant arcs, in thousands of galaxy clusters within the redshift range of $0.2 \le z \le 1.0$ (Boldrin et al. 2012, 2016). These observations will allow for the creation of detailed mass models of clusters' inner regions (e.g. Kneib et al. 1993; Bradač et al. 2005; Diego et al. 2005; Liesenborgs et al. 2006; Coe et al. 2008; Jullo et al. 2007; Zitrin & Broadhurst 2009; Oguri 2010; Zitrin et al. 2013; Lam et al. 2014). The

resulting constraints will test the predictions of ΛCDM and alternative dark matter models, such as self-interacting dark matter, on a cluster scale (Meneghetti et al. 2001; Miralda-Escudé 2002; Meneghetti et al. 2020; Rocha et al. 2013; Peter et al. 2013; Meneghetti et al. 2022, 2023; Granata et al. 2023).

In addition, akin to the aforementioned compound lenses, the simultaneous observation of numerous sources at varying redshifts, all lensed by the same clusters, will provide constraints on cosmological parameters such as $\Omega_m$ and $w$ through the lensing sensitivity to angular-diameter distances (Gilmore & Natarajan 2009; Jullo et al. 2010; Moresco et al. 2022; Caminha et al. 2022; Acebron et al. 2023; Bergamini et al. 2024).

Lastly, *Euclid*'s unique ability to integrate both strong and weak lensing measurements will enable the determination of mass profiles of galaxy clusters from kiloparsec to megaparsec scales (Bartelmann et al. 1996; Bradač et al. 2005; Umetsu et al. 2016). These measurements are essential for accurately determining the total cluster mass and are crucial for constraining the shape and redshift evolution of the cluster concentration-mass relation. Hydrodynamical simulations indicate that the concentration of dark matter halos correlates with the universe's density at the time of their collapse (Navarro et al. 1997; Gao et al. 2008; Ludlow et al. 2013). Therefore, measuring this key relation is vital for validating the ΛCDM cosmological framework (Meneghetti et al. 2011, 2014; Merten et al. 2015).

## 9.4. Cross-correlation with CMB observables

During their journey towards us the CMB photons interact with the large-scale structures of the Universe as they are forming. These structures leave their imprint on the CMB through gravitational lensing (Lewis & Challinor 2006) and via the scattering of CMB photons with electrons having significant thermal and bulk velocities, called the thermal Sunyaev-Zeldovich (tSZ) and kinetic Sunyaev-Zeldovich (kSZ) effects, respectively (Carlstrom et al. 2002). Additionally, the decay of gravitational potentials caused by the accelerated expansion in the late Universe generates new anisotropies in the CMB temperature at large angular scales (integrated Sachs-Wolfe (ISW) effect; Sachs & Wolfe 1967). Maps of the CMB lensing convergence ($\kappa$ hereafter) and of the strength of tSZ (parameterised through the Compton $y$ parameter) can be extracted from high-resolution multi-frequency observations of CMB anisotropies. They are sensitive to the total integrated matter or pressure distribution along the LoS between us and the surface of last scattering, respectively. The *Euclid* survey overlaps on the sky with the major existing CMB data sets, such as *Planck*, SPT-3G and ACT (Planck Collaboration VIII 2020; Benson et al. 2014; Wu et al. 2019; Qu et al. 2024), as well as future ground-based CMB experiments such as the Simons Observatory (SO; Ade et al. 2019) and CMB-Stage 4 (Abazajian et al. 2016) or space-based experiments such as LiteBIRD (Hazumi et al. 2020).

The cross-correlation power spectra between CMB lensing and *Euclid* galaxy clustering ($C_\ell^{kg}$) and weak lensing maps ($C_\ell^{kY}$), together with the CMB lensing auto-correlation $C_\ell^{kk}$, will provide additional observables that are sensitive to cosmological parameters affecting the angular-diameter distances and the growth of the matter perturbations, and as such they can tighten the statistical uncertainties (Sailer et al. 2021). Moreover, they are also free from additive systematic biases (Vallinotto 2012; Schaan et al. 2017). They will thus allow us to break degeneracies and minimise the impact of systematic effects and theoretical uncertainties that might affect the *Euclid* observables when analysed
on their own, or provide new estimators that are less sensitive to systematic effects (Giannantonio et al. 2016; Bianchini & Reichardt 2018). These CMB-*Euclid* cross-correlations will therefore not only add statistical power, but they will also allow us to marginalise over parameters describing systematics with minimal loss of constraining power. Adding these observables on top of the 3×2pt analysis in a combined 6×2pt analysis has become a standard in the field for current surveys (Abbott et al. 2023; Robertson et al. 2021; Marques et al. 2020). Considering also the correlation of galaxy clustering with CMB temperature $C_{\ell}^{kg}$ in a full 7×2pt analysis with multiple tomographic redshift bins, the constraining power of cosmological parameters for the joint analysis of *Euclid* with SO and CMB-S4 data can reduce the statistical uncertainty by a factor of 2–3 and in some cases even more than 10, in particular for generalised cosmological models including curvature or modifications of gravity (Euclid Collaboration: Ilić et al. 2022).

Further improvements are expected if higher-order statistics, involving mixed bispectra are used, correlating one or more CMB lensing fields with *Euclid* probes (Chen et al. 2021; Farren et al. 2023). Recent studies have also shown that cross-correlations between $\kappa$ and biased density tracers (such as galaxies or QSOs) can be used to set competitive constraints on local primordial non-Gaussianity through the scale-dependent galaxy bias, Cross-correlations do this in a more robust way compared to what is achievable with the auto-correlation analysis of the tracers, which is plagued by large-scale survey systematic effects, such as inhomogeneous depth, air mass, or selection effects (Rezaie et al. 2024; Krolewski et al. 2024). *Euclid*'s space-based observations will allow us to carry out these measurements with exquisite precision, as has recently been demonstrated in the context of the Gaia mission (Alonso et al. 2023; Storey-Fisher et al. 2024). The AGNs and QSOs sample detected by *Euclid* (see Sects. 9.5 and 10.4) in particular will probe primordial non-Gaussianity in a redshift range never surveyed before.

*Euclid* probes can also be cross-correlated with SZ maps. Cross-correlation with the tSZ $y$ maps ($C_{\ell}^{y\ell}$, $C_{\ell}^{\gamma\ell}$, $C_{\ell}^{\gamma\gamma}$), on top of being interesting probes in their own right to probe properties of the hot gas in the Universe, are highly sensitive to the physics of baryons. This can be used to constrain (and marginalise) models of feedback or other baryonic effects in 3×2pt analyses of the *Euclid* probes (Pandey et al. 2022; Tröster et al. 2022; Osato et al. 2020; Koukoufilippas et al. 2020; Kou & Bartlett 2023) using a full combination 10×2pt analysis with tSZ and CMB lensing (Fang et al. 2024). Unlike for the tSZ effect, the kSZ signal cannot be separated from CMB temperature maps, since it has the same frequency dependency as the CMB itself and hence can only be seen in combination with an external tracer of the LSS, either in cross-correlation with 2D matter tracers or through velocity-weighted stacking techniques using *Euclid* spectroscopic sample (Hill et al. 2016; Schaan et al. 2021). The cross-correlation between the squared CMB temperature of SO or CMB-S4 maps and *Euclid* galaxy clustering and weak lensing data will enable measurements of the kSZ effect with an overall S/N of about 20, providing statistical constraints on the parameters describing the shape of the gas radial density profile in halos (and thus on the underlying physical mechanisms) at the 10–20% precision level (Bolliet et al. 2023). Stacking techniques will be extended for the first time to higher redshifts, thanks to the capabilities of the NISP instrument. The spectroscopic power of *Euclid* will also allow us to detect the kSZ effect, and more generally the cosmological information encoded in the velocity field, through approaches like the pairwise momenta

(Hand et al. 2012; Planck Collaboration Int. XXXVII 2016), velocity field reconstructions (DeDeo et al. 2005; Planck Collaboration Int. XXXVII 2016; Schaan et al. 2016), or cross-correlation to maps of the so-called ‘angular redshift fluctuations’, which are sensitive to galaxy radial motions (Hernández-Monteagudo et al. 2021; Chaves-Montero et al. 2021). Combined with angular galaxy clustering, the latter are expected to improve constraints on the dark energy equation of state by almost an order of magnitude, compared to angular galaxy clustering alone (Legrand et al. 2021).

## 9.5. Cosmology with high-redshift quasars

Quasars have long been known to obey a nonlinear relation between the rest frame 2500 Å ($L_{UV}$) and the rest frame 2 keV ($L_X$) emission (e.g. Tananbaum et al. 1979; Zamorani et al. 1981), parameterised as $L_X \propto L_{UV}^{\gamma}$, with $\gamma \approx 0.6$. Recently, this relation has been employed to provide an independent measurement of quasar distances, thus turning these objects into *standardisable* candles and extending the distance modulus–redshift relation (the so-called Hubble diagram) of Type Ia supernovae to a redshift range that is still poorly explored ($z > 2$; Risaliti & Lusso 2015, see also Moresco et al. 2022).

The applicability of this technique is based upon two main factors. First, the understanding that most of the observed dispersion in the $L_X$-$L_{UV}$ relation is not intrinsic to the relation itself, but due to observational issues (e.g. X-ray absorption by gas, UV extinction by dust, calibration uncertainties in the X-rays, variability, and selection biases associated with the flux limits of the different samples). Once corrected for and with an optimal selection of *clean* sources (i.e. where the *intrinsic* UV and X-ray quasar emission is measured), the dispersion in the relation becomes rather small, namely $\approx 0.2$ dex (Lusso & Risaliti 2016). Second, the realisation that the slope of the $L_X$-$L_{UV}$ relation does not evolve with redshift up to $z \approx 4$, which is the highest redshift where the source statistics are currently sufficient to verify any possible dependence of the slope with distance. Before *Euclid*, the largest quasar sample that can be used for cosmological analysis is composed of approximately 2400 sources out to $z \approx 7.5$ (Lusso et al. 2020, with about 500 quasars beyond redshift 2). As of today, the precision achieved by the combined use of quasars and Type Ia supernovae is on the order of 25% on $\Omega_m$ and 20% on $w_0$, assuming an evolving dark energy equation-of-state parameter as in Eq. (1) – see Table 1 in Bargiacchi et al. (2022).

We expect a quasar sample of around $2 \times 10^6$ sources in the EWS out to $z \approx 5$ detected in all four *Euclid* bands ($5\sigma$ detection following the colour-colour AGN selection described in Sect. 10.4), with X-ray emission above the eROSITA limiting depth (at a 2–10 keV flux limit higher than $5.2 \times 10^{-18} \text{ W m}^{-2}$; see Selwood et al., in prep. for details)$^{29}$. We forecast that roughly 20% of *Euclid* quasars will fulfil the selection criteria described in Lusso et al. (2020), that is, unobscured at both UV and X-ray energies, radio quiet, with negligible contamination from the host galaxy, with several thousands of objects at $z > 2$ surviving the selection cuts. Approximately, 10% of this sample will have spectroscopic redshifts from the NISP instrument (see also Euclid Collaboration: Lusso et al. 2024). The *Euclid* quasar sample, complemented with eROSITA data and additional X-ray

<footnote>$^{29}$ eROSITA sources detected in the 0.2–2.3 keV energy range with a detection threshold $> 6$, corresponding to a point source flux limit of $6.5 \times 10^{-18} \text{ W m}^{-2}$, assuming a power law with photon index $\Gamma = 1.9$ (Brunner et al. 2022).</footnote>