spatially resolved colour information from HST is needed to quantify shape-
measurement biases introduced by colour gradients within a galaxy; these
biases are inevitable for a diffraction-limited PSF in Euclid’s wide
$I_{\scriptscriptstyle\rm E}$-band (Voigt et al. 2012; Semboloni et al. 2013a;
Er et al. 2018). Extensive spectroscopic redshift surveys of the EAFs make
them very suitable for the calibration of photometric redshifts. The EAFs are
observed with the ROS (Sect. 4.1), providing full VIS and NISP data sets. In
total, 17% of the time during routine operations is used on the EDFs, EAFs,
and for self-calibration (Fig. 17). These data form a scientific cornerstone
of the Euclid mission with substantial legacy value (Sect. 10) owing to their
great depth and large number of revisits.

A central pillar of Euclid’s calibration scheme is the set of monthly visits
of the self-calibration field (Sect. 4.2.3). These last about 23 hours and
obtain a large number of images and spectra to recalibrate Euclid’s
spectrophotometric response with high accuracy, S/N, and spatial resolution.
In this way we counter any adverse effects from space weathering and
outgassing (Euclid Collaboration: Schirmer et al. 2023), and from
imperfections in the preceding calibration chain, thus enabling a consistent
flux calibration over the full 6-year mission duration. The self-calibration
observations also provide data for regular updates of 20 further calibration
products.

The NISP hardware calibration plan foresees monthly lamp flats and
nonlinearity calibrations, and a lower cadence of reciprocity failure and
wavelength-calibration checks. The VIS calibration plan includes biases,
flats, darks, trap-pumping, and charge-injection lines on a daily basis, as
these are part of the ROS. On a monthly scale we will calibrate the VIS
nonlinearity chain, check for radiation damage using trap-pumping, and obtain
numerous flat-fields for a high-S/N characterisation of the brighter-fatter
effect, and the conversion gain or photon-transfer curve (PTC). The VIS PSF
model is recalibrated every 2–6 weeks, depending on when the survey moves to a
new large survey patch (Euclid Collaboration: Scaramella et al. 2022). In
total, 7% of the time is used for dedicated hardware calibrations that are not
part of the ROS (Fig. 17).

###  4.4 Complementary ground-based photometry

Euclid relies on optical ground-based imaging that complements the VIS and
NISP imaging for photometric redshift estimation (Abdalla et al. 2008) and to
assign the correct SED-weighted PSF to each galaxy in the lensing analysis
(Eriksen & Hoekstra 2018). To this end, a large coordinated campaign of
ground-based observations with different observatories will provide the
necessary multi-band photometry to matching depths across the EWS and EDS
areas. Here, we summarise the characteristics of these data.

The DES provides a good starting point for the southern sky. This survey,
completed in 2019, covers about 3750 deg2 of the EWS, with achieved depths of
$g=24.5$, $r=24.1$, $i=23.6$, and $z=23.4$ ($10\sigma$ for a point source in a
$2^{\prime\prime}$ diameter aperture). These depths meet the requirements in
the $gri$ bands, while it comes close in $z$. These data are sufficient for
the cosmological parameter estimates based on the first data release (see
Sect. 7.8). To enable more precise measurements for the final analyses, the
DES data will be superseded by deeper derived data products from the LSST (Guy
et al. 2022), which will overlap with 7534 deg2 of the EWS.

The northern sky, however, lacked an equivalent data set: KiDS (de Jong et al.
2015) and the HSC survey (Aihara et al. 2018, 2022) do not cover sufficient
area, and largely target regions of the northern sky that are closer to the
ecliptic, whilst the DESI Legacy Imaging Survey (Dey et al. 2019) is one
magnitude too shallow. To address the need for additional complementary
imaging data a new collaboration was set up in 2017: the Ultraviolet Near
Infrared Optical Northern Survey (UNIONS), a wide field $ugriz$ survey of the
northern extragalactic sky that is a ‘collaboration of collaborations’. The
Canada-France Imaging Survey collaboration provides $u$\- and $r$-band imaging
using the Canada-France-Hawaii Telescope (CFHT), $i$-band and part of the
$z$-band data are obtained using the Panchromatic Survey Telescope and Rapid
Response System (Pan-STARRS; Chambers et al. 2016), while the Wide Imaging
with Subaru-Hyper Suprime-Cam Euclid Sky (WISHES) team, a collaboration of
Japanese scientists, acquires $z$-band imaging with the HSC (Miyazaki et al.
2018). In addition, HSC $g$-band data are collected through PI time, via a
Canadian Gemini-Subaru exchange and time from the Institute for Astronomy,
University of Hawaii: the Waterloo-Hawaii-IfA $g$-band Survey (WHIGS). UNIONS
is becoming the definitive broadband optical survey of the northern sky and
once fully acquired and combined, these will provide the required $ugriz$
coverage of the northern part of the EWS over 5711 deg2, joining with the LSST
at a declination of $+15$ degrees over the north Galactic Cap. For reference,
the achieved depths by UNIONS are: $u=23.6$, $g=24.5$, $r=24.1$, $i=23.7$,
$z=23.4$ ($10\sigma$ for a point source in a $2^{\prime\prime}$ diameter
aperture). As it remains a challenge to ensure that at every Euclid data
release the footprints covered from space and from the ground overlap as much
as possible, the first year of the survey will be prioritised towards the
southern sky, while the data collection in the north by UNIONS continues and
balances out the north and the south for the second Euclid data release.

Figure 19: The galaxy multicolour-space to $i=25$ AB mag, encoded in a 2D map
with a 150$\times$75 binning using the self-organising map algorithm (Masters
et al. 2015). On the left is the distribution of spectroscopic coverage of the
map prior to the C3R2 effort. The white regions are those parts of galaxy-
colour space lacking high-confidence spectroscopic redshifts for calibration.
On the right is the current map, after incorporating the $>$5800 C3R2 faint
galaxy spectra. The map coverage has increased from about 51% to $>$90%, with
many colour cells calibrated with multiple galaxies. Spectra to calibrate the
remaining empty cells may be obtained as next-generation spectroscopic
facilities come online, or they can be addressed with clustering redshift
approaches (e.g., Newman 2008). We note that the remaining empty regions
correspond to lower-density (less occupied) parts of the galaxy-colour space.

The EDFs (Sect. 4.2.1) and EAFs have their own challenges, with limiting-
magnitude requirements of around 26 AB mag over 56 square degrees. The Cosmic
Dawn Survey (Euclid Collaboration: McPartland et al., in prep.) is a
complementary, UV–IR multiwavelength survey for the EDFs and EAFs that aims to
optimise the legacy science returns from these fields, with a primary focus on
the high-redshift Universe. The DAWN survey combines dedicated and archival
observations from CFHT MegaCam, Subaru HSC, LSST, Spitzer IRAC and other
ancillary data available in the EDFs and EAFs, to depths matching the Euclid
observations. Photometric catalogues of the imaging data were consistently
produced using The Farmer (Weaver et al. 2022). The Spitzer observations and
data reduction are described in Euclid Collaboration: Moneti et al. (2022). A
description of the ground-based observations from the Hawaii 20 deg2 (H20)
Survey and the first catalogue data release can be found in Euclid
Collaboration: Zalesky et al. (in prep.).

###  4.5 Complementary ground-based spectroscopy

As discussed in Sect. 2.2, the correct interpretation of the photometric
clustering and lensing signals depends critically on accurate estimates of
their redshift distributions. This, in turn, relies on large samples of robust
spectroscopic redshifts. To calibrate the redshift distributions of the lenses
and sources for the weak lensing and photometric clustering measurements, we
have collected a substantial amount of complementary deep spectroscopy.

Our baseline approach is outlined by Masters et al. (2015) and aims to
directly calibrate the relation between galaxy colours and the redshift.
Specifically, Masters et al. (2015) quantified the expected distribution of
galaxy colours using the self-organising map (SOM) algorithm (Kohonen &
Honkela 2007), showing that sizeable regions of the colour space were lacking
high-confidence spectroscopic redshifts (see the left panel in Fig. 19). They
also demonstrated that a targeted campaign of spectroscopic follow-up could
obtain the redshifts needed for calibration.

This analysis motivated the Complete Calibration of the Colour-Redshift
Relation (C3R2) programme, a coordinated effort between Keck (Masters et al.
2017, 2019; Stanford et al. 2021), the Very Large Telescope (Euclid
Collaboration: Guglielmo et al. 2020), and the Large Binocular Telescope
(Euclid Collaboration: Saglia et al. 2022), to measure redshifts for faint
galaxies down to $i=25$ AB mag spanning the full galaxy colour space. The Keck
programme, using DEIMOS, MOSFIRE, and LRIS, resulted in $>$5100 deep spectra
with secure redshifts, while the other observations added over 600 more. The
improvement in colour-space coverage as a result of the C3R2 efforts is shown
in Fig. 19.

We emphasise that the C3R2 surveys built upon extensive existing spectroscopy
to calibrate the colour-redshift relation. Other large-scale deep
spectroscopic surveys that contributed substantially to the external redshift
calibration sample for Euclid include DEIMOS 10k (Hasinger et al. 2018),
DEEP2&3 (Cooper et al. 2012), MOSDEF (Kriek et al. 2015), UDSz (Bradshaw et
al. 2013), VANDELS (McLure et al. 2018), VIPERS (Garilli et al. 2014;
Scodeggio et al. 2018), VUDS (Tasca et al. 2017), VVDS (Le Fèvre et al. 2013),
and zCOSMOS (Lilly et al. 2007).

All deep spectra were collected in a database for careful source-by-source
redshift validation, given the need for high purity in the spectroscopic
calibration sample (Euclid Collaboration: Ilbert et al. 2021). We are
currently testing the performance for the planned tomographic redshift binning
for Euclid, using the spectroscopic calibration sample together with the
Flagship mock galaxy simulation (Sect. 6.1).

##  5 Early results from commissioning and PV

Comprehensive results from the commissioning and PV phases will be published
once the ongoing data analyses have concluded. The estimated amount of raw
data produced by Euclid during the PV phase between August and November 2023
is 21 TB, compared to 12.1 TB of raw data created by HST up until November
2023.171717File names in the HST archive that end in _raw.fits (Matthew
Burger, Space Telescope Science Institute, priv. comm.) Therefore, in this
paper we present some preliminary results, showing that Euclid’s performance
is sufficient for its core scientific goals. More details can be found in the
instrument-specific papers for VIS (Euclid Collaboration: Cropper et al. 2024)
and NISP (Euclid Collaboration: Jahnke et al. 2024).

###  5.1 Photometric throughput and molecular contamination

For initial estimates of the total system throughput we used the Gaia (Gaia
Collaboration: Vallenari et al. 2023), 2MASS (Cutri et al. 2003), and VHS
(VISTA Hemispheric Survey; McMahon et al. 2021) magnitudes of field sources.
Considerable colour terms exist in the transformations from the four wide
Euclid bands to these external bands (Fig. 12). Our transformations are based
on the known passbands and stellar SEDs from our simulations (Sect. 6); in the
case of NISP they are given in Euclid Collaboration: Schirmer et al. (2022).

We find the measured NISP and VIS zero points (ZPs) for sources with a
frequency-flat SED to be considerably better, by about 0.4 mag, than required.
Accurate ZPs and updated Euclid photometric systems will be based on
observations of stable white-dwarf spectrophotometric standards that we
already established with HST, jointly for Euclid and Roman (proposal ID
16702).

Thin layers of water ice formed on optical surfaces due to outgassing,
resulting in throughput modulations from interference and scattering.
Contamination is expected and typical for spacecraft, and can be countered
with thermal decontamination (for details see Euclid Collaboration: Schirmer
et al. 2023).

Throughput monitoring has shown characteristic variations that have been
linked to several nanometres of ice on the optics; selective heating of one of
the folding mirrors in March 2024 completely restored the transmission to
immediate post-launch levels. Because the instrument cavity has little venting
area to the outer space, and because outgassing is a continuous process,
further decontamination activities are expected over the operation period.

Figure 20: VIS image quality. The figure shows a stacked data PSF near the
centre of the VIS FPA, from an observation of the self-calibration field,
averaging over source SEDs. The FWHM is approximately $None$ in this data set.
The effect of trefoil (Sect. 5.2) is evident in this log-scale representation.
Figure 21: Point-like cosmic-ray density in VIS during a low M-class solar
flare. The cosmic rays are caused by X-rays impinging onto the detectors after
penetrating Euclid’s sunshield in some gaps between the solar cells, causing
characteristic geometric patterns. During rare but bright X-class flares, up
to 25% of the VIS detector area must be masked. The location of the pattern
and its shape depends strongly on the LoS of VIS towards the Sun through the
sunshield, and thus on the spacecraft’s attitude.

###  5.2 VIS point-spread function

For the weak lensing measurements we need a detailed and chromatic PSF model
that we introduce in Sect. 7.6.4. Here we present an initial and coarse
evaluation of the PSF. With a plate scale of $None$ pixel-1 the VIS optical
PSF is considerably undersampled (Fig. 20). To avoid aliasing in resampled
images, the sampling frequency would need to be at least $2.0$–$2.5$ pixels
per FWHM (Bernstein 2002). The requirements on the VIS PSF are an ellipticity
less than 0.13, and an FWHM smaller than $None$ at 800 nm. The typical
ellipticity and FWHM we measure in VIS science exposures are 0.04 and $None$,
respectively. This includes jittering effects from the FGS for representative
guide-star densities and a background of cosmic-ray hits. The jittering is
accounted for in the PSF model using the time-series of the guiding
corrections during every exposure. From an opto-mechanical perspective alone,
Euclid’s PSF has great stability owing to its SiC components and a thermally
stable environment at L2.

In pre-launch testing, the polished M1 was found to have a small residual
amount of astigmatism, which would result in a strong dependence of the PSF
ellipticity on the telescope’s focus. To reduce this astigmatism, a mechanical
correction was applied to M1 to compensate the astigmatism prior to launch. It
is thought that this resulted in a small, but noticeable, amount of trefoil in
the PSF, visible as an approximately triangular shape in the PSF contours at
$None$–$None$ from the PSF core (Fig. 20). Like PSF ellipticity, the trefoil
is an inherent part of the PSF model (Sect. 7.6.4) and thus accounted for in
the shape measurements. The trefoil is described – to leading order – as a
spin-3 contribution to the PSF. We therefore do not expect that it has a
significant impact on the multiplicative and additive shear biases in Eq.
(17), as those are described by spin-0 and spin-4 components.

Figure 22: Impact of the spacecraft orientation on the VIS background. Left:
At ${\rm AA}=0$ considerable stray light levels are present that exceed the
zodiacal background by more than one order of magnitude. Right: For ${\rm
AA}<-$$$ the stray light is reduced to a few percent of the zodiacal
background. It still needs to be modelled for some calibrations and low-
surface-brightness science.

###  5.3 X-ray contamination from solar flares

Euclid’s instruments have radiation shields, and are additionally protected by
the spacecraft itself. This also includes protection from X-rays produced
during solar flares, which are absorbed by the silicon in Euclid’s sunshield.
However, gaps in the sunshield’s solar cells do let X-rays pass, some of which
intersect with the VIS FPA at an angle-of-incidence of about $60^{\circ}$.
These are then detected as excess cosmic rays in the VIS images (Fig. 21).
NISP is not affected by this.

Contrary to protons that cause displacement damage in the detector’s atomic
lattice, X-rays harmlessly create electron-hole pairs in the CCD’s depletion
region, like optical photons. The increased density of cosmic rays renders a
fraction of the image unusable for scientific analysis. Using the X-ray
sensors (XRS; Hanser & Sellers 1996) onboard the Geostationary Operational
Environmental Satellites (GOES-16 and GOES-18), we find that considerable data
loss occurs once the solar X-ray flux approaches about $1\times 10^{-5}$ W m-2
in the GOES long band (1.5–12.4 keV).

Typical flares last about 10–60 minutes, so that one or several subsequent VIS
images can be affected. As long as particularly active sunspot regions are
visible on the Sun, up to 10% of the VIS images can be substantially affected
by flares, judging from the GOES solar activity recorded between March and
November 2023. During these periods, we expect to lose data from up to 4% of
all VIS pixels while passing through the solar maximum in 2024–2025. Any area
lost could possibly be recovered during periods of currently unallocated time
(Sect. 4.1.2).

Figure 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
${\rm\lx@glossaries@gls@link{acronym}{AA}{\leavevmode AA}}=$$$. 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.

###  5.4 Optical stray light in VIS

During commissioning it was found that VIS can be affected by considerable
amounts of stray light that exceed the zodiacal background by a factor of 10
or more (Fig. 22). The root cause has not been unambiguously identified, but
is thought to be a thruster nozzle that is illuminated by the Sun. The light
enters the instrument cavity following a triple scattering process, from the
nozzle to the backside of the Sun shield, through a hypothesised opening in
the thermal multi-layer insulation, to a mounting leg of the VIS shutter, and
from there to the VIS focal plane. Even closed-shutter VIS observations are
affected, and several parasitic light paths exist. NISP is not affected, as it
is enclosed in black multi-layer thermal insulation (MLI), and parasitic light
entering through the dichroic is blocked by baffles.

Rotating the spacecraft around its $Z$-axis (Fig. 5) imposing $\rm AA<$$$
moves the nozzle into the shadow, effectively reducing the stray light to
levels of a few percent of the zodiacal background (Fig. 22). Rotating up to
$\rm AA=$$$ is safe following a post-launch evaluation. The reduced AA range
now available for the survey is $[$$,$$]$ (previously $[$$,$$]$ ) with a
margin of $None$ for orbit uncertainty. The stray light is then negligible for
Euclid’s core science. Figure 23 shows how the latest survey configuration
adapts the spacecraft orientation to the stray light constraint. Low-surface
brightness science and some calibrations still require the construction of
stray light models from the large number of survey fields and calibration
data.

Figure 24: Reach in ecliptic longitude around transit for the leading side of
the survey, for ${\rm\lx@glossaries@gls@link{acronym}{SAA}{\leavevmode
SAA}}\in[$$,$$]$. 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 $\rm{\lx@glossaries@gls@link{acronym}{AA}{\leavevmode AA}}\in[$$,$$]$ with
much reduced visibility (grey) that would not permit the completion of the
survey. By allowing the fields to rotate by up to $None$ 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 $None$ around
the origin. Figure 25: EWS coverage and colour-coded yearly progress in an
all-sky Mollweide projection. The blue borders enclose the $16\,000\,\deg^{2}$
RoI that contains the $13\,416\,\deg^{2}$ observed sky of the EWS. Small dark
regions within the EWS are masks for stars brighter than $4$ AB mag.

###  5.5 Latest survey strategy

The reduction of the AA range has the double effect of thinning and skewing
the window of visibility, which no longer contains the transit meridian and
has a zero reach at some latitudes (Fig. 24). This inhibits the scheduling of
Wide Survey patches of any substantial size. Furthermore, half the time it
would be impossible to schedule fields with $|\beta|\leq$$$.

Hence, we allowed the FOV to rotate by up to ${\pm}\,$$None$ with respect to
the tessellation tiles, recovering a large fraction of the tiles’ original
visibility, albeit with a different shape. Visibility was regained at all
latitudes. A post-launch check showed that we could also safely increase the
maximum SAA from $None$ to $None$. The new enlarged window of visibility,
shown for the leading side in Fig. 24, has two distinct parts in the northern
and southern ecliptic hemispheres due to the now asymmetric AA range. The
visibilities for the leading and trailing sides are east-west and north-south
mirrored. In the leading side, the larger reach in longitude favours
observations in the southern hemisphere, while the trailing side favours the
northern hemisphere.

This relaxed tessellation constraint implies that neighbouring FoVs are no
longer necessarily aligned. To avoid gaps between FoVs we tessellate the sky
with smaller tiles, increasing the mean overlap between fields, which reduces
the covered sky area. In the latest survey computation, Fig. 25, the EWS
covers 13 416 deg2 of which 171 deg2 are lost due to bright stars. The
effective sky area is 13 245 deg2, representing a decrease of 5% as compared
to the target of 14 000 deg2 in 6 years (Sect. 4.1). The implications of this
reduction in survey area for the core science objectives, as well as possible
mitigation strategies are being explored.

##  6 Simulated data

The calibration observations discussed in Sect. 4.3 provide important
information about the performance of the telescope and the data that are
collected. To prepare the pipeline and to interpret the results, simulated
data are needed as well. This includes large realistic input universes that
can capture survey characteristics, such as the Euclid Flagship Simulation
(EFS; Sect. 6.1) as well as sophisticated pixel-level instrument simulators
(Sect. 6.2). The latter are essential for exploring the sensitivity of the
measurements to instrumental effects and test our ability to correct these.

###  6.1 The Euclid Flagship Simulation

The optimal exploitation of the Euclid data demands the development of large-
volume and high-mass resolution numerical simulations that reproduce the
large-scale galaxy distribution that the mission will observe with high
fidelity. Not only do these help to assess the performance with a realism that
cannot be achieved otherwise, but such simulations are also an essential tool
for the development of the data processing and science analysis pipelines. The
major advance that the Euclid data will bring implies the need for a dedicated
effort. To this end, we developed the EFS, which is described in detail in
Euclid Collaboration: Castander et al. (2024). Here, we summarise its main
characteristics.

The EFS features a simulation box of 3600 $h^{-1}$Mpc on a side with
$16\,000^{3}$ particles, leading to a mass resolution of $m_{\rm
p}=10^{9}\,h^{-1}M_{\odot}$. This 4 trillion particle simulation is the
largest $N$-body simulation performed to date and matches the basic science
requirements of the mission, because it allows us to include the faintest
galaxies that Euclid will observe, while sampling a cosmological volume
comparable to what the satellite will survey. The simulation was performed
using PKDGRAV3 (Potter & Stadel 2016) on the Piz Daint supercomputer at the
Swiss National Supercomputer Center (CSCS). The input cosmology differs
slightly181818The EFS uses the following values for the density parameters:
$\Omega_{\rm m}=0.319$; $\Omega_{\rm b}=0.049$; and
$\Omega_{\Lambda}=0.681-\Omega_{\rm rad}-\Omega_{\nu}$, with a radiation
density $\Omega_{\rm rad}=0.00005509$, and a contribution from massive
neutrinos $\Omega_{\nu}=0.00140343$. Additional parameters are: the equation
of state of dark-energy $w=-1.0$; the reduced Hubble constant $h=0.67$; the
scalar spectral index of the initial fluctuations $n_{\rm s}=0.96$; and the
scalar power spectrum amplitude $A_{\rm s}=2.1\times 10^{-9}$ (corresponding
to $\sigma_{8}=0.813$) at $k=0.05\,{\rm Mpc}^{-1}$. from the one listed in
Table 1, but this has no material impact on the applications.

The initial conditions were realised at $z=99$ with first-order Lagrangian
perturbation theory (1LPT) displacements from a uniform particle grid. The
transfer functions for the density field and the velocity field were generated
at this initial redshift by CLASS (Lesgourgues 2011) and CONCEPT (Dakin et al.
2022). As the usual scaling of the linear power spectrum at $z=0$ to the
initial redshift of the simulation (known as back-scaling), to generate the
initial conditions, was not used, all linear contributions from radiation,
massive neutrinos, and metric perturbations (in the $N$-body gauge, see Fidler
et al. 2015) were included via a lookup table and applied as a small
corrective PM (particle-mesh) force at each timestep. This ensures a match to
the linear evolution of the matter density field at all redshifts when
including these additional linear terms. The main data product was produced on
the fly during the simulation and is a continuous full-sky particle light cone
out to $z=3$, where each particle was output exactly when the shrinking light
surface sweeps by it. This resulting ball of particles contains 31 trillion
particle positions and peculiar velocities (700 TB of data). The 3D particle
lightcone data were used to identify roughly 150 billion dark-matter halos
using Rockstar (Behroozi et al. 2013), and to create all-sky dark-matter 2D
maps in 200 tomographic redshift shells between $z=0$ and $z=99$, with a
HEALPix (Górski et al. 2005) tessellation resolution $N_{\rm side}=8192$,
corresponding to $None$ per pixel.

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. (2024), 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).

Figure 26: The image on the left shows the lensing convergence for sources
with $z_{\rm s}=1$ for a simulated patch of sky covering 50 deg2. 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 $\pm 3\sigma$, where
$\sigma$ 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.

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 (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 using observed
relations, taking the H$\alpha$ 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\,\mathrm{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. 2024); 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: Vallenari et al. 2023) 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_{\scriptscriptstyle\rm 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 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,
Imagem 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 aXeSIM191919http://axe-
info.stsci.edu/axesim/ (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 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 2-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

Figure 27: VIS view of a $$$\times$$$ 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.

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: Vallenari et al. 2023) 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
$$$\times$$$ 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_{\scriptscriptstyle\rm E}$, $J_{\scriptscriptstyle\rm E}$, and
$H_{\scriptscriptstyle\rm 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
pre-processing, 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 catalogues, application of flat-fielding and photometric
calibration based on the pre-computed 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).

Figure 28: False-colour NISP image of a $$$\times$$$ area of Euclid’s self-
calibration field (Sect. 4.2.3). Filters $Y_{\scriptscriptstyle\rm E}$,
$J_{\scriptscriptstyle\rm E}$, and $H_{\scriptscriptstyle\rm 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.

###  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, e.g., 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 $None$ to $None$),
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. 2016; Waters et al. 2016). Similarly the $g$\-
and $z$-band SEFs from the WHIGS and WISHES surveys 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).

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: Prusti et al. 2016) photometry
and spectroscopy that are stable across the sky with a systematic uncertainty
of $\sim$2 mmag (Gaia Collaboration: Vallenari et al. 2023). 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 SED s across the Euclid
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
$\sim$2 mmag systematic uncertainties in Gaia photometry. Therefore, with the
availability of the Gaia DR3 (Gaia Collaboration: Vallenari et al. 2023) $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ṁmag 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
coadds 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 $None$ 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_{\scriptscriptstyle\rm E}$, $J_{\scriptscriptstyle\rm E}$, and
$H_{\scriptscriptstyle\rm E}$ data. Source detection is done with
SourceXtractor++202020https://github.com/astrorama/SourceXtractorPlusPlus
(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-PHOT212121The 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). (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_{\rm 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 the
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 pre-processed 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 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_{\scriptscriptstyle\rm 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_{\scriptscriptstyle\rm E}$-band
photometric image of the source. This also includes appropriate rotation,
matching of both NISP-P and NISP-S PSF s, 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 222222These 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. 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) 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.

Figure 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 $\text{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_{{\rm H}\alpha}=2\times
10^{-15}$ erg cm-2 s-1, is approximately ten times higher 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 ($z$PDF) is calculated for each
model. All those individual $z$PDFs are combined into one so-called ‘first
pass $z$PDF’. A first list of redshift solutions is based on the main peaks
(up to 10) of this first $z$PDF. The redshift measurement is then refined
around each of these solutions by refitting over all the models, but with a
finer redshift grid; the $z$PDF 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 $z$PDF
under each peak with a $3\,\sigma$ 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_{\rm 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 $z$PDF, 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 $z$PDF 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_{\scriptscriptstyle\rm 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 $z$PDF 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 $z$PDF 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^{\rm obs}_{ij}$, when a galaxy
image is convolved by the PSF. For an object with an observed surface-
brightness distribution, $I^{\rm obs}(\mathbf{x})$, and total flux, $F$, its
quadrupole moments are defined as

$Q_{ij}=\frac{1}{F}\int\,{\rm d}^{2}\mathbf{x}\,\;x_{i}\,x_{j}\,I^{\rm
obs}(\mathbf{x})\;.$ (20)

The shape of an object can then be quantified by combining the quadrupole
moments into the complex polarisation (or distortion)232323This definition is
related to the third flattening (see Sect. 2.2) through
$\epsilon=e/(1+\sqrt{1-|e^{2}|})$. , $e=e_{1}+{\rm i}e_{2}$, where

$e_{1}=\frac{Q_{11}-Q_{22}}{Q_{11}+Q_{22}}\;,\quad{\rm
and}\;\;e_{2}=\frac{2\,Q_{12}}{Q_{11}+Q_{22}}\;,$ (21)

while the size of an object is captured by

$R^{2}=Q_{11}+Q_{22}.$ (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 sum
of the quadrupole moments of the true galaxy image and those of the PSF
(Valdes et al. 1983), so that

$Q_{ij}^{\rm obs}=Q_{ij}^{\rm gal}+Q_{ij}^{\rm PSF}.$ (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.

Figure 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.

For individual galaxies, the need to place objects in tight redshift
intervals, with $\sigma_{z}<0.05\,(1+z)$ being a key requirement leading to
two considerations: the effective dispersion of the $z$PDF 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 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 $z$PDF, 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 (DR) (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 $z$PDFs from a carefully selected and curated set of deep
reference fields, the EAFs. These reference field photo-$z$s 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 $z$PDF 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 $z$PDFs 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_{\scriptscriptstyle\rm 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_{\scriptscriptstyle\rm
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 ${\rm
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, i.e., their
$z$PDFs 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 $z$PDF, 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 pre-trained 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 pre-determined 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\,{\la}\,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\,{\rm nm}\,{<}\,\lambda\,{<}\,1000\,{\rm 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_{\scriptscriptstyle\rm 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_{\rm PSF}^{2})/R_{\rm
PSF}^{2}<10^{-3}$. Although the Euclid VIS system possesses a highly stable,
diffraction-limited PSF, meeting these requirements is nonetheless challenging
because of a number of design choices.

First, to reach the required depth, VIS uses a broad bandpass, 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 V242424CODE V® is a commercial optical
design tool. 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.

  * •

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). 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 downsampled 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.

Figure 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\,$\mathrm{\SIUnitSymbolMicro m}$$. All images are produced as the mean of
flux-selected stars, applying 2 $\sigma$ 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.

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 $\pm 18\,\mu$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 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, b; 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; Kitching et al. 2008a; Miller et al. 2013), 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_{\rm leak}\,\epsilon^{\rm PSF}$ to Eq. 17,
where $p_{\rm leak}$ quantifies the imperfect correction for the PSF
ellipticity $\epsilon^{\rm 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_{\rm 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_{\rm g}=\left(\frac{\sigma_{\epsilon}}{\sigma_{m}\,|\gamma|}\right)^{2}\;,$
(24)

where $|\gamma|$ is the modulus of the shear applied to the simulated
galaxies,252525The 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.,) 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^{\circ}$ 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 2-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
3×2pt 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 2-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\text{\times}{10}^{-16}$\,\mathrm{erg}\,\mathrm{cm}^{-2}\,\mathrm{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 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 providing 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 predefined 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 2PCF s 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 3-point statistics. In particular, we have developed an efficient estimator
for the 3-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 3-point estimator takes the same input as the
2-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 2- and 3-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 2-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

Figure 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.  Figure 33: Similar to Fig. 32,
but for angular clustering from galaxy positions.  Figure 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.

Euclids photometric galaxy samples are determined by their detection in the
$I_{\scriptscriptstyle\rm E}$ filter. Nominally, this includes all galaxies
with $I_{\scriptscriptstyle\rm 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_{\scriptscriptstyle\rm 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_{\scriptscriptstyle\rm 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_{\scriptscriptstyle\rm E}$ detection
efficiency, e.g., due to 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 2-point correlations of:

  1. 1.

galaxy ellipticities, tracing cosmic shear;

  2. 2.

galaxy positions, tracing angular galaxy clustering; and

  3. 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. 1.

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

  2. 2.

angular clustering is measured separately using estimators $\hat{w}(\theta)$