© The Authors 2025

1. Introduction

Galaxies in the Universe tend to cluster as they evolve within the large-scale structure, itself growing under the influence of gravity and universal expansion. Their observed positions are possibly the best tracers of the overall three-dimensional (3D) matter spatial distribution, given their high and unrivalled observed number density, among all observable objects on cosmological scales. This makes galaxy clustering and its features crucial for studying the expansion and structure growth in the late Universe (e.g. Amendola et al. 2018). In particular, the baryon acoustic oscillations (BAO) feature can be used as a standard ruler, making galaxy clustering sensitive to the geometry and energy content of the Universe (e.g. Alam et al. 2021). Moreover, peculiar velocities contribute to observed galaxy redshifts and distort the genuine spatial distribution along the line-of-sight direction. On large scales, these velocities map the coherent flows towards over-densities induced by the growth of structure, whose strength is dictated by the laws of gravity. This apparent feature is of paramount importance in testing standard gravity and the cosmological model (Kaiser 1987; Guzzo et al. 2008).

Galaxy clustering can be quantified in a statistical manner by determining the n-point statistics of the spatial distribution of galaxies. This is achieved in configuration space by estimating the two-point correlation function (2PCF), as well as higher order n-point correlation functions, which enable the exploration of non-Gaussian features of the galaxy spatial distribution. The 3D 2PCF, alongside the power spectrum, is the most commonly used tool for analysing galaxy spectroscopic survey data and inferring cosmological model parameters. For this purpose, it is crucial to estimate the contributions of the 2PCF along the parallel and transverse directions to the line of sight, as the apparent spatial distribution of galaxies is anisotropic.

The Euclid mission is a medium-sized European Space Agency (ESA) space mission to unravel the dark sector of the Universe (Euclid Collaboration: Mellier et al. 2025). It is primarily devoted to the investigation of the nature of dark energy and the dark matter distribution. The Euclid satellite was launched on 1 July 2023 and will survey about 14 000 deg² of the extragalactic sky, performing one of the largest galaxy surveys ever made. It will probe the last 10 billion years of the universal expansion history via its main cosmological probes: weak gravitational lensing and galaxy clustering. In particular, Euclid’s Near-Infrared Spectrometer and Photometer (NISP, Euclid Collaboration: Jahnke et al. 2025) will measure the spectroscopic redshifts of about 20–30 million Hα-emitting galaxies using grism spectroscopy. The resulting 3D galaxy map will serve to characterise the clustering of matter and, in turn, determine the cosmological parameters with an unprecedented precision. In this regard, the estimation of the 2PCF of spectroscopic galaxies will play a crucial role. Nonetheless, the estimation of the 2PCF from such a large dataset presents a number of challenges, particularly in terms of the efficiency of the estimation. To address those challenges, we have designed a robust, amply tested, and efficient software that can be used within the Euclid Science Ground Segment (SGS). This pipeline element, which provides one of the critical end-products for Euclid scientific exploitation, is the processing function 2PCF-GC.

This paper presents the method and 2PCF-GC processing function software devised to estimate the 3D 2PCF within the Euclid SGS pipeline. It describes the ensemble of tests that have been performed to ensure the highest accuracy and efficiency in the estimation. Furthermore, it provides a discussion on how the 2PCF-GC processing function addresses Euclid computational challenges and provides forecasts on real survey 2PCF estimation. In the following, we simply refer to the 2PCF-GC processing function and software as 2PCF-GC. This paper is part of a series that describes all Euclid galaxy clustering processing functions, which also includes those for three-point correlation function (Veropalumbo et al., in prep.), as well as the power spectrum and bispectrum estimations (Salvalaggio et al., in prep.).

The paper is structured as follows. Section 2 presents the 2PCF definition and estimators. Section 3 discusses the optimisation of the software. Section 4 describes the validation and performance tests as well as associated results. Section 5 shows predictions of 2PCF measurements in the Euclid Wide Survey and Sect. 6 provides our conclusions.

2. Two-point correlation function

2.1. Definition and estimators

The 2PCF is defined as the probability of finding two objects at volume elements, dV₁ and dV₂, separated by a vector, r, with respect to a Poisson random distribution. This excess probability is defined as (Peebles 1980):

$$\begin{array}{c}\hfill \mathrm{d}P=n^2[1+\xi (\mathit{r})]\mathrm{d}{V}_1\mathrm{d}{V}_2,\end{array}$$(1)

where n is the object mean number density; therefore, ξ(r) measures the clustering in excess or in deficit compared to a random Poisson point distribution in space. In the previous equation, statistical homogeneity is implicitly assumed in that the 2PCF is only a function of the separation vector and not on the positions of the two objects.

In 2PCF-GC, we used the minimum-variance estimator from Landy & Szalay (1993), which allows the estimation of the two-point auto-correlation function for a set of spatially distributed objects, based on pair counts as a function of separation. As opposed to other proposed estimators, it minimises the estimator variance and mitigates discreteness and edge effects. This estimator is defined as

$$\begin{array}{c}\hfill \xi (\mathit{r})=\frac{\mathrm{DD}(\mathit{r})-2\mathrm{DR}(\mathit{r})+\mathrm{RR}(\mathit{r})}{\mathrm{RR}(\mathit{r})},\end{array}$$(2)

where DD, DR, and RR correspond to data-data, data-random, and random-random normalised (distinct) pair counts, respectively. It requires a catalogue of unclustered objects (random catalogue hereafter) that randomly samples the effective volume of the data catalogue, for which we want to estimate the 2PCF. Schematically, the Landy & Szalay (1993) estimator can be derived by taking the two-point correlation of the catalogue density contrast δ defined from data and random catalogues, such that δ = (D − R)/R, where D and R represent the data and associated random catalogue counts, respectively. The raw pair counts need to be normalised to the total number of pairs. With N_D and N_R as the number of objects in the data and random catalogues, respectively, the normalised counts are obtained as

$$\begin{array}{cc}\hfill \mathrm{DD}(\mathit{r})& =\frac{2\widehat{\mathrm{DD}}(\mathit{r})}{N_{\mathrm{D}}(N_{\mathrm{D}}-1)},\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill \mathrm{DR}(\mathit{r})& =\frac{\widehat{\mathrm{DR}}(\mathit{r})}{N_{\mathrm{D}}{N}_{\mathrm{R}}},\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill \mathrm{RR}(\mathit{r})& =\frac{2\widehat{\mathrm{RR}}(\mathit{r})}{N_{\mathrm{R}}(N_{\mathrm{R}}-1)},\hfill \end{array}$$(5)

where $\widehat{\mathrm{DD}}$, $\widehat{\mathrm{DR}}$, $\widehat{\mathrm{RR}}$ correspond to raw counts. These factors result from the fact that the number of distinct pairs that can be drawn from N objects is N(N − 1)/2 and the number of cross pairs that can be drawn from two sets of N₁ and N₂ objects is N₁N₂. In the presence of object weights, w, as in the case when we are using Feldman et al. (1994) weights, each pair in the counts contribute with the multiplication of individual object weights and the normalised counts become

$$\begin{array}{cc}\hfill \mathrm{DD}(\mathit{r})& =\frac{2\widehat{\mathrm{DD}}(\mathit{r})}{{\left({\displaystyle {\sum }_{i=1}^{N_{\mathrm{D}}}}{w}_i\right)}^2-{\displaystyle {\sum }_{i=1}^{N_{\mathrm{D}}}}{w}_i^2},\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill \mathrm{DR}(\mathit{r})& =\frac{\widehat{\mathrm{DR}}(\mathit{r})}{\left({\displaystyle {\sum }_{i=1}^{N_{\mathrm{D}}}}{w}_i\right)\left({\displaystyle {\sum }_{i=1}^{N_{\mathrm{R}}}}{w}_i\right)},\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill \mathrm{RR}(\mathit{r})& =\frac{2\widehat{\mathrm{RR}}(\mathit{r})}{{\left({\displaystyle {\sum }_{i=1}^{N_{\mathrm{R}}}}{w}_i\right)}^2-{\displaystyle {\sum }_{i=1}^{N_{\mathrm{R}}}}{w}_i^2},\hfill \end{array}$$(8)

where the index, i, goes over the objects in the catalogues.

We could also want to measure the cross-correlation function between two sets of objects in an overlapping volume. For this, we can devise a cross-correlation estimator in a similar fashion as for the auto-correlation case. By defining the density contrasts for the two populations, δ₁ = (D₁ − R₁)/R₁ and δ₂ = (D₂ − R₂)/R₂, where now D₁ (D₂) and R₁ (R₂) stand for the data and random catalogue counts of the population 1 (2), the two-point cross-correlation function estimator is expressed as

$$\begin{array}{c}\hfill {\xi }_{12}(\mathit{r})=\frac{{\mathrm{D}}_1{\mathrm{D}}_2(\mathit{r})-{\mathrm{D}}_1{\mathrm{R}}_2(\mathit{r})-{\mathrm{R}}_1{\mathrm{D}}_2(\mathit{r})+{\mathrm{R}}_1{\mathrm{R}}_2(\mathit{r})}{{\mathrm{R}}_1{\mathrm{R}}_2(\mathit{r})},\end{array}$$(9)

where D₁D₂, D₁R₂, R₁D₂, and R₁R₂ are the data 1-data 2, data 1-random 2, random 1-data 2, and random 1-random 2 normalised pair counts, respectively.

For some purposes, it can also be useful to define a modified density contrast that uses two different random catalogues. Such a modification allows the subtraction of known correlations from the data catalogue, for instance, for the purpose of building BAO reconstruction estimators (Padmanabhan et al. 2012) or correcting for observational systematic errors (e.g. Paviot et al. 2022). If we define an auxiliary random catalogue, S, and modified density contrast, δ^m = (D − S)/R, the estimator for the auto-correlation function associated with that density contrast becomes

$$\begin{array}{c}\hfill {\xi }^{\mathrm{m}}(\mathit{r})=\frac{\mathrm{DD}(\mathit{r})-2\mathrm{DS}(\mathit{r})+\mathrm{SS}(\mathit{r})}{\mathrm{RR}(\mathit{r})}·\end{array}$$(10)

Similarly, the cross-correlation function estimator associated with the modified density contrast is expressed as

$$\begin{array}{c}\hfill {\xi }_{12}^{\mathrm{m}}(\mathit{r})=\frac{{\mathrm{D}}_1{\mathrm{D}}_2(\mathit{r})-{\mathrm{D}}_1{\mathrm{S}}_2(\mathit{r})-{\mathrm{S}}_1{\mathrm{D}}_2(\mathit{r})+{\mathrm{S}}_1{\mathrm{S}}_2(\mathit{r})}{{\mathrm{R}}_1{\mathrm{R}}_2(\mathit{r})}.\end{array}$$(11)

The pair counts and estimated 2PCF with these estimators involve a binning scheme in the pair separation. In 2PCF-GC, the latter can be set to either linear or logarithmic base 10. For instance, in the case of the angle-averaged¹ 2PCF with linear binning in r, the notation ξ(r) will refer to the 2PCF for separations between r − Δr/2 and r + Δr/2, where Δr is the constant linear bin size. In the case of logarithmic binning, it will refer to the 2PCF for logarithms of the separation between log₁₀(r/u)−Δlog₁₀(r/u)/2 and log₁₀(r/u)+Δlog₁₀(r/u)/2, where Δlog₁₀(r/u) is the constant logarithmic bin size and u is the fiducial unit of length.

Finally, it is important to emphasise that all previous estimators are exact and unbiased in the limit of infinite random catalogues, for Poisson and binomial point processes and generally for a point process with structures much smaller than the size of the sample (Kerscher 1999). From the computational perspective, those estimators reduce to counting pairs as a function of separation from different catalogues. It is clear from the definitions that the main technical challenge in estimating the 2PCF is the ability of counting all pairs as a function of separation from large catalogues. Naively, this process would scale as N², where N is the number of objects in the catalogue. This can become intractable in the case of very large catalogues, typically with numbers of object above 10⁶, but efficient algorithms can reduce the amplitude of the scaling.

2.2. Scale dependence

The 2PCF is a function of the separation vector r and under the assumption of an isotropic universe, the 2PCF may only depend on the norm of the separation vector r = |r|. Nonetheless, correlations are affected by apparent physical anisotropies, such as redshift-space distortions, which we want to quantify. For this, the 3D separation vector can be decomposed into transverse and parallel to the line-of-sight directions. We define the transverse and parallel separation vectors, ⃗r_⊥ and ⃗r_∥, such that

$$\begin{array}{c}\hfill \mathit{r}=r_{\perp }+r_{\Vert }\end{array}$$(12)

and

$$\begin{array}{cc}\hfill r_{\Vert }& \equiv r_{\Vert }{u}_{\Vert }=(\mathit{r}·u_{\Vert })u_{\Vert },\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill r_{\perp }& \equiv r_{\perp }{u}_{\perp }=r-(\mathit{r}·u_{\Vert })u_{\Vert },\hfill \end{array}$$(14)

where r = ⃗x₂ − ⃗x₁, ⃗x₁ and ⃗x₂ are the 3D positions of the two objects, x is the pair line-of-sight vector (defined below in Eq. (15)), ⃗u_∥ = x/|x|, and ⃗u_⊥ is a unit vector perpendicular to ⃗u_∥ and belonging to the plane determined by the non-collinear vectors r and x. With these definitions, we can define the anisotropic 2PCF ξ(r_⊥, r_∥).

In the decomposition presented in Eq. (12), we made the implicit assumption that all galaxy lines of sight are parallel and we defined a mean pair line-of-sight direction, x. This choice is motivated in the distant-observer limit, when r_⊥ ≪ r_∥. In full generality, one should consider the dependence on the triangular configuration comprising the observer and the two objects of the pair. This would require the use of three scalars to define the pair configuration. The systematic error in the correlation function estimated within the distant-observer approximation, usually referred to as wide-angle effects is, however, only relevant at large transverse separations and is generally neglected. Usual prescriptions² for the pair line-of-sight definition are the mid-point, bisector, or end-point definitions, the mid-point being the most commonly used. These correspond to using

$$\begin{array}{c}\hfill \mathit{x}=\{\begin{array}{cc}{x}_1/2+x_2/2\hfill & \text{for mid-point},\hfill \\ x_1\hfill & \text{for end-point},\hfill \\ x_1/|x_1|+x_2/|x_2|\hfill & \text{for bisector}.\hfill \end{array}\end{array}$$(15)

The geometry and different definitions are made explicit in Fig. 1. Given the symmetries present in the mid-point and bisector definitions, they are less sensitive to wide-angle distortions than the end-point one. Those three definitions use the so-called local plane-parallel approximation, while using a single constant line-of-sight direction for all pairs in a catalogue would lead to the so-called global plane-parallel approximation (Samushia et al. 2015). The standard approach used in 2PCF-GC is to use the mid-point definition but the other definitions can be chosen for specific cases.

Fig. 1. Geometry and separation definition for end-point (left), mid-point (middle), and bisector (right) pair line-of-sight conventions. The observer is at the point ‘O’ and the two other points represent the objects of the pair. In those figures, all vectors and segments are contained in the plane determined by the vectors ⃗x and ⃗r.

Instead of directly using a decomposition of the separation vector around the line-of-sight direction, it can be useful to express the separation vector in the associated polar basis with coordinates

$$\begin{array}{cc}\hfill r& =|\mathit{r}|,\hfill \end{array}$$(16)

$$\begin{array}{cc}\hfill \mu & =r_{\Vert }/|\mathit{r}|,\hfill \end{array}$$(17)

where μ is the cosine angle between the separation and line-of-sight directions (see Fig. 1). From this we can further define the 2PCF multipole moments of the order, ℓ,

$$\begin{array}{c}\hfill {\xi }_{\ell }(r)=\frac{2\ell +1}{2}{\displaystyle {\int }_{-1}^1}\mathrm{d}\mu \xi (r,\mu )L_{\ell }(\mu ),\end{array}$$(18)

where L_ℓ(μ) is the Legendre polynomial of order ℓ. In the auto-correlation function case, there is no physically expected odd multipole signal. However, in the case of the cross-correlation between different populations, a relativistic signal can arise in the odd multipole correlation functions (e.g. Breton et al. 2019). 2PCF-GC has been built such that this type of signal can be extracted and provides in output all even and odd moments up to ℓ = 4.

Finally, an estimate of real-space clustering, namely, without the effect of redshift-space distortions, can be obtained by integrating ξ(r_⊥, r_∥) over r_∥. This leads to the projected 2PCF

$$\begin{array}{c}\hfill w_{\perp }(r_{\perp })={\displaystyle {\int }_{-\infty }^{\infty }}\mathrm{d}{r}_{\Vert }\xi (r_{\perp },r_{\Vert }),\end{array}$$(19)

where, in practice, the integral is definite and finite integration limits have to be defined to mitigate the effect of large uncorrelated pairs in the integral. In 2PCF-GC, both integrals in Eqs. (18) and (19) are evaluated as Riemann sums over the linearly binned anisotropic 2PCF.

2.3. Data spatial partitioning

The implementation of the 3D 2PCF estimation in Euclid uses specific data partitionings to enable an efficient estimation from the huge dataset that Euclid will produce. An overall spectroscopic sample of about 20–30 million galaxies with redshifts is expected in the completed Euclid Wide Survey (Euclid Collaboration: Mellier et al. 2025). To address this task, three efficient pair-counting algorithms have been developed, based on different data spatial partitionings: linked-list, k-d tree, and octree. They exploit the observation that, at fixed requested scale range, not all possible distinct pairs have to be computed and stored. Those data spatial partitionings allowed us to explore efficiently all pairs where the separation falls within the requested scale range and prune the exploration of irrelevant data. The purpose of developing different pair-counting methods is to have different ways of assessing the same quantity and identify the fastest and more reliable method. After the development and optimisation of the methods, we found that in the end all three are very efficient, as discussed in the next section. It is worth emphasising that those methods are exact and no approximation is involved at the pair-counting level.

2.3.1. Linked-list algorithm

The linked-list algorithm for range searching, also sometimes referred to as the chained-mesh algorithm, implements a 3D regular pixellation scheme. Our implementation builds on the work of Marulli et al. (2016) and a similar algorithm has been used in other implementations (e.g. Alonso 2012; Donoso 2019; Sinha & Garrison 2020). The data bounding volume was divided into a regular Cartesian mesh and the indexes of objects residing in each cell are stored in a list: the linked list. In practice, the elements (or nodes) of the list are not explicitly linked, instead they are stored in a vector of indexes that map to their 3D positions.

The pair counting was performed through two nested loops: a first loop goes through all cells containing objects, while the second explores the neighbouring cells whose maximum distance is below the requested maximum separation. In the inner loop, all pairs of objects are counted as in the naive nested-loop algorithm. This strategy allows the pruning of irrelevant pair counts that are outside of the required separation range. The efficiency of the algorithm depends on the cell size and mean number of object per cell. Optimally, one would like to have a cell volume corresponding to a multiple of the search sphere and the most appropriate mean number of object per cell. The latter has to be large enough to avoid having to explore too many cell-cell pairs but small enough to avoid having a too large cell volume, and in turn the pruning to be inefficient. After several trials, an optimal choice of setting the number of cells such that there are about 100 objects per cell on average is used in our implementation. Practically, this is done by first taking the maximum separation as cell size and estimating the averaged number of objects per cell. Then, the cell size is divided by the integer value that allows us to reach approximately 100 objects per cell on average.

2.3.2. k-d tree algorithm

The k-d tree range-search algorithm (Bentley 1975) is based on partitioning the set of object positions into a binary tree. Starting from the smallest axis-aligned bounding volume encompassing all data points, namely, the tree root, tree nodes are obtained by recursively dividing into two equipopulated subsets until reaching the leaves of the tree. Each subsequent partitioning forms two children nodes that contain subsets of the parent node set and in which objects are close in space. Each binary split increases by unity the depth of the tree. In the classical k-d tree, the node splitting is performed along the dimension with largest spread, at the median object position. The leaves correspond to the highest-depth nodes, when the number of objects reaches a minimum value. In our implementation, we used a minimum value of 100 objects and the sliding-midpoint method for the node splitting rule, which is more adapted for clustered objects (Maneewongvatana & Mount 1999). This choice has proven to be the most efficient.

A crucial aspect of k-d tree is that the bounding volume coordinates of each node are kept in the data structure. This information is then used to search through the tree. For the purpose of counting pairs, we use the dual-tree approach that is a generalisation of the single-tree range-search algorithm for pairs (Moore et al. 2001; Zhao 2023). The search is performed by spanning two trees simultaneously and testing each node pair recursively. Depending of the minimum and maximum distance between the tested nodes, the search can be stopped or passed through to the children nodes. This process stops when reaching the leaves. With this method, one can efficiently prune pairs, when separations are outside of the required separation range. The pair counting is effectively performed at the level of the leaves, using a nested loop going over all pairs, whose complexity goes as 𝒪(n²), where n is the number of objects in the leaves.

2.3.3. Octree algorithm

In the octree algorithm (Meagher 1980), the partitioning of the data volume is organised in a tree structure in which each internal node has exactly eight children, as opposed to two children for the k-d tree. The octree is built from the smallest cubic bounding volume encompassing all the data points. This root node is then recursively subdivided in octants, obtained by dividing the side in each of the three dimensions in two equal parts. Each subsequent partitioning increases the depth of the tree and this process stops when reaching the leave nodes. The leaves can be defined by imposing the maximum depth of the tree, or a minimum value for the number objects in the leaves. Similarly as for the k-d tree, and after several trials, we set the latter minimum value to 100 objects.

Our implementation used a hashed octree structure (Warren & Salmon 1993), where tree nodes are stored in a hash table and each node is identified by its Morton binary code (Morton 1966). This allows for the optimisation of the data structure storage and memory access, while also leading to the most optimal performance in tree traversing. Because of the regularity of the octree spatial structure, the node bounding volume coordinates can efficiently be deduced from the Morton code by using bitwise operations. This information is then used to search through the tree. The pair counting is performed similarly as for the k-d tree, by spanning two octrees simultaneously and testing each node pair recursively. Depending of the minimum and maximum distance between the tested nodes, the search can be stopped or passed through to the children nodes, and the process stops when reaching the leaves. The pair counting is performed at the level of leaves, in a similar manner as for the k-d tree.

2.4. Software architecture

The 2PCF-GC processing function is a processing element of the Euclid SGS, which carries out the entire data processing up to cosmological parameter extraction. There are ten Organisational Units (OU) within the SGS, each one having the responsibility to define, design, and validate a specific analysis of the SGS workflow. 2PCF-GC belongs to OU-Level 3 that is in charge of producing the highest-level scientific data products. The data processing within the Euclid SGS is performed in a distributed system across the Science Data Centers (SDC) from Finland, France, Germany, Italy, Netherlands, Spain, Switzerland, United Kingdom, and the United States. While most of the SDC have a High Throughput Computing (HTC) design, the underlying infrastructure can vary across SDC. For high-performance computation, as particularly required by this processing function, the overall infrastructure allows for parallelisation.

The development of the 2PCF-GC code was performed in C++11 within the framework of the Euclid SGS common tools and guidelines to ensure homogeneity in terms of development, storage, and computing, independently of the location (Frailis et al. 2019). This development followed the common Euclid Development ENvironment (EDEN), which establishes the set of libraries and associated versions to be used by any of the Euclid software and prevents inconsistencies or changes in the functionality of different libraries between development and production. 2PCF-GC was integrated in the COllaborative DEvelopment ENvironment (CODEEN), a continuous integration and delivery (CI/CD) platform that automates the building, unit testing, and distribution of all the scientific software in the SGS. The source code is stored in a Version Control System (Gitlab) and can be run through a CI/CD pipeline to be finally deployed on a distributed file system available on all SDCs. This system design allows SGS operations to be performed smoothly and efficiently across all SDC, providing the extra advantage of increased computing power and storage capacity. 2PCF-GC has been developed to run within the SGS infrastructure, but can also be run in standalone mode, which has demonstrated to be crucial for testing and validation activities.

2.4.1. Process overview

The overall 2PCF-GC processing follows four main steps as illustrated in Fig. 2. In the first step, it reads the inputs: a configuration file, data and random catalogues, and pre-computed pair counts if this option is selected in the configuration file. In pipeline mode, those catalogues are provided by the preceding processing function in the SGS pipeline chain, the SEL-ID processing function. The latter extracts a catalogue from the Euclid Wide Survey using certain selection criteria and provides it with the associated random catalogue to 2PCF-GC. The input galaxy and random catalogues are then read and decomposed into an internal spatial representation: linked-list, k-d tree, or octree, depending on the chosen pair-counting method. After building the internal data structure, the counting algorithm roams over it to identify the pairs as a function of separation. Weighted pair counts in each separation bin are stored in arrays. 2PCF-GC performs the necessary pair counts in series depending on the requested estimator. Alternatively, they can be read from input files. Those pair counts are finally combined to obtain the 2PCF estimate. At the end of the process the 2PCF and individual pair counts products are prepared and delivered in the form of FITS files (Pence et al. 2010).

Fig. 2. 2PCF-GC process overview. The different boxes illustrate the different steps involved in the processing.

2.4.2. Inputs and outputs

The inputs and outputs of 2PCF-GC are defined in the Euclid Common Data Model, which defines the format of the input and output data and metadata, and ensures the stability of interfaces between pipelines and the Euclid Archive System (EAS). The latter is the database containing all data products and metadata processed for the Euclid mission (Williams et al. 2019). The input products of 2PCF-GC are a configuration file, a set of data and random catalogues, and possible pre-computed pair counts from a previous run. The input catalogues contain the 3D spatial information of objects as well as optional statistical weights. The latter can be used for instance to up-weight or down-weight objects in pairs to account for variations in the spatial sampling of certain objects. The input celestial coordinate system can be either the equatorial Cartesian or spherical system, where in the spherical case, the radial coordinate is either a redshift or a comoving radial distance. In the case that the redshift is provided, 2PCF-GC first converts the redshift to a comoving radial distance from a provided fiducial cosmological model. After that, Cartesian coordinates are used in the pair counting. The output products are tabulated 2PCF measurements and associated pair counts. Apart from the input catalogues, through the parameter list, the configuration file allows for the specification of:

• the correlation function estimator: auto-correlation, cross-correlation, modified auto-correlation, modified cross-correlation;

• the correlation function type: angle-averaged ξ(r), anisotropic ξ(r_⊥, r_∥) and w_⊥(r_⊥), anisotropic ξ(r, μ) and ξ_ℓ;

• the pair-counting method: linked-list, k-d tree, or octree;

• the type of binning (linear or logarithmic) and definition of bins;

• the pair line-sight definition for ξ_ℓ;

• the upper limit of integration along r_∥ for the projected correlation function;

• an option to enable reusing pre-computed pair counts;

• the number of splits of the random catalogue when using the random split option.

3. Optimisation

The main challenge in estimating the 2PCF from Euclid data is to perform pair counting from massive galaxy and random catalogues in a reasonable amount of time. Euclid scientific accuracy requirements impose a 10% accuracy relative to statistical uncertainty. This translates at estimation level into imposing a number of random catalogue objects of at least 50 times that of the galaxy catalogue, since the estimator variance depends on the number of objects in the data and random catalogues (Landy & Szalay 1993; Keihänen et al. 2019). The choice of using random catalogues with at least 50 times more objects than in the data catalogue is a formal requirement of the Euclid mission endorsed by ESA, allowing us to reach only a percent-level contribution to the 2PCF standard error from the finiteness of the random catalogue. A significant effort has been put on optimising 2PCF-GC for speed and memory consumption. This was performed by making use of optimal data spatial partitioning and pair-counting method as previously discussed, but also by implementing parallelisation and a specific treatment of random-random pair counts based on a the random split technique. We describe those two aspects in the following.

3.1. Parallelisation

The parallelisation of pair counting can be achieved straightforwardly for the three considered algorithms. For the linked-list algorithm, this is done by splitting the loop over mesh cells in chunks and performing the computation for each chunk in parallel. In the case of the k-d tree and octree algorithms, instead of starting from the tree root, the pair counting starts in parallel from all internal tree nodes at a given depth (greater than zero). In all cases, the parallel instances computes partial counts that are then summed up in order to obtain the final counts. Those parallelisation strategies have been implemented using the shared-memory Open Multi-Processing (OpenMP) application programming interface. The scaling and performance of these implementations are presented in Sect. 4.

3.2. Treatment of random-random pair counts

Random-random pair counts dominate the overall computation time of the 2PCF as it involves the largest number of pairs. A major gain in runtime can be obtained by using the so-called random split technique in the computation of those pairs. This technique, described in Keihänen et al. (2019), relies on first splitting the random catalogue in N_S sub-catalogues, calculating all DR_i and RR_i pair counts, and finally summing up sub-catalogue pair counts to obtain the final random-random and data-random pair counts. Therefore the estimator for the 2PCF (auto-correlation case) becomes

$$\begin{array}{c}\hfill \xi (\mathit{r})=\frac{\mathrm{DD}(\mathit{r})-2{\mathrm{DR}}^{\prime }(\mathit{r})+{\mathrm{RR}}^{\prime }(\mathit{r})}{{\mathrm{RR}}^{\prime }(\mathit{r})},\end{array}$$(20)

where

$$\begin{array}{cc}\hfill {\mathrm{DR}}^{\prime }(\mathit{r})& =\frac{1}{N_{\mathrm{S}}}{\displaystyle \sum _{i=1}^{i=N_{\mathrm{S}}}}{\mathrm{DR}}_i(\mathit{r}),\hfill \end{array}$$(21)

$$\begin{array}{cc}\hfill {\mathrm{RR}}^{\prime }(\mathit{r})& =\frac{1}{N_{\mathrm{S}}}{\displaystyle \sum _{i=1}^{i=N_{\mathrm{S}}}}{\mathrm{RR}}_i(\mathit{r}).\hfill \end{array}$$(22)

With this strategy the maximum number of random-random pairs to be computed is smaller by a factor of $N_{\mathrm{S}}^{-1}$ compared to the case without random split, leading to a significant gain in computational time. An optimal N_S is such that each random sub-catalogue is of the same size as the data catalogue. Keihänen et al. (2019) studied the bias and variance of the 2PCF estimator with this treatment and found that, for a random catalogue 50 times larger than the galaxy catalogue, this technique reduces the computation time by a factor of more than ten without affecting estimator variance or bias. In the following, we adopt as baseline random catalogues that have 50 times the number of galaxies in the data catalogue, and when using the random split technique, an optimal value of N_S = 50. We show in the next sections that this effectively allows us to reach the accuracy requirement of Euclid.

4. Validation and performance

The testing and validation of 2PCF-GC were performed at the different stages of development of the code. 2PCF-GC has successfully passed six maturity level gates, each one involving a series of validation tests including 2PCF calculations on mock data. Significant efforts have been put in making the code meet high coding standards following the quality requirement defined by SGS. We present in this section the results from the most significant validation tests, with a focus on accuracy, runtime, and memory tests.

4.1. Benchmark catalogues

To perform benchmarks of 2PCF-GC, we made use of different sets of mock data catalogues. Each one was selected to perform a specific test and has different characteristics. We review in the following those characteristics.

• Euclid Large Mocks: The ELM suite has been designed for studying observational systematic errors on galaxy clustering. Those are realisations of Hα-emitting galaxies in a lightcone of radius 30 deg on the sky. They are based on the Pinocchio approximate method for efficiently generating halo catalogues and lightcones (Monaco et al. 2013). Pinocchio halos are populated with emission-line galaxies using a halo occupation distribution model calibrated on the flagship galaxy mock presented in the next section. As for the flagship galaxy mock, the simulated galaxy catalogue corresponds to Hα emitters with flux above 2 × 10⁻¹⁶ erg s⁻¹ cm⁻² and distributed over the redshift range 0.9 < z < 1.8. To perform validation tests, we only considered galaxies at 0.9 < z < 1.1 in the first mock, leading to an effective volume of V = 1.5 h⁻³ Gpc³. Details of the ELM mocks are provided in Monaco et al. (in prep.).

• CoxMock suite: The CoxMock suite has been designed to produce mocks for which the true 2PCF is known, and which can be produced massively in a reasonable amount of time. To do this, we used an isotropic line-point Cox process where lines of a given length are randomly placed in a periodic cubical volume and points are randomly scattered on those lines (Stoyan et al. 1995). We generated 10 000 catalogue realisations of a cube of side 1.74 h⁻¹ Gpc containing 10⁶ objects. For each realisation, 100 000 lines of length 800 h⁻¹ Mpc and 1 000 000 points were drawn. More precisely, each line was drawn at random position and orientation by randomly selecting one of the endpoints inside the simulation box. Each point was then randomly placed on one line by first randomly choosing a line and second uniformly distributing it on the chosen line. If lines or points fell outside of the cube boundaries, they were mirrored by the opposite face, effectively applying periodic boundary conditions. The expected 2PCF of points in these mocks resembles a cosmological ξ(r) as it is described by a damped power-law with an index of −2 (see Sect. 4.3).

4.2. Runtimes and scaling with number of objects

We tested the runtime of 2PCF-GC with the different pair-counting algorithms and as a function of the number of objects in the catalogues. For this, we extracted four sub-catalogues with N_o = [10⁴, 10⁵, 10⁶, 5 × 10⁶] randomly chosen objects from the ELM mock. We then built the associated random catalogues with 50 times more objects than in the data. The redshift distribution of random objects was drawn from the estimated mean distribution over 1000 mocks. Since the volume of the mock is the same, the sub-catalogues probe different number densities of objects. In this test, we considered the estimation of the multipole correlation functions in 40 linear bins in r spanning the interval [0, 200] h⁻¹ Mpc and 200 in μ. The runtimes are obtained using a computer cluster node equipped with 32 physical central processing units (CPU) Intel(R) Xeon(R) Silver 4216 at 2.10 GHz, and by performing 5 runs with 4, 8, 16, or 32 parallel threads, except for the longest runs with N_o > 10⁶ where only one run with 32 threads was performed.

Figure 3 presents the times in CPU-hour that we obtained when estimating the multipole correlation functions with the different algorithms, as a function of the number of objects. The total runtimes, denoted by 2PCF in the figure, are divided in different parts corresponding to the time spent on data structure construction, DD and RR pair counting. While the data structure construction time scales linearly with the number of objects, the DD, DR, and RR pair-counting times scale quadratically for a constant volume, as expected. Overall, the data structure construction has a subdominant contribution of one up to seven orders of magnitude smaller than pair counting. The DD, DR, and RR pair-counting times are similar for the tree-based algorithms, while the linked-list algorithm shows a steeper slope with the number of objects, that is, smaller times for low N_o and slightly longer times above 10⁶ objects. The full calculation of the 2PCF is dominated by the RR counts and the overall runtimes reflect this, with the tree-based algorithms and particularly the octree performing best. Furthermore, the random split option allows us to drastically reduce the runtimes by up to a factor of 10 for a data catalogue of 5 × 10⁶ objects, as shown in Fig. 3. In that case, the 2PCF multipoles are computed in about 7 hours on 32 threads, while without this lasts 4 days. The scaling with the number of objects is slightly different for the three pair-counting algorithms. Since the individual random split sub-catalogues have approximately the same size as the data catalogue, the times are dominated by smaller-size catalogues for which the linked-list algorithm proves to be slightly faster. This is true up to approximately N_o = 5 × 10⁵, and at larger N_o the runtimes are similar for the three algorithms, with only the octree being marginally faster for 5 × 10⁶ objects. For the full 2PCF calculation, the linked-list algorithm runtimes scale very accurately as N_o², while for tree-based algorithms as $N_{\mathrm{o}}^{1.96}$ in the regime of N_o > 10⁵.

Fig. 3. Run times for the calculation of the multipole correlation of galaxies obtained from the ELM mock. The times are expressed in CPU-hour and as a function of the data or random catalogue size. The various symbols represent the time spent on the data structure construction, DD calculation, DR calculation, RR calculation, and the overall 2PCF runtime with and without random split option. The DR calculation times are provided as a function of the number of objects in the data catalogue and assuming a fifty times larger random catalogue. The different curves represent the runtimes obtained with the linked-list (solid), k-d tree (dashed), and octree (dotted) algorithms. The blue vertical band shows the range of expected number densities in the spectroscopic sample at redshifts within 0.9 < z < 1.8. The abscissa refers to the number of object in the data catalogue except for RR calculation where it refers to that in the random catalogue.

The volume of the ELM mock that we used for runtime tests corresponds approximately to the volume probed spectroscopically by the Data Release 1 of Euclid in the redshift interval 0.9 < z < 1.1. We show with the vertical blue band in Fig. 3 the expected number density of galaxies over 0.9 < z < 1.8, which gives an idea of the expected runtime to measure the multipole correlation functions in the first redshift interval: about 30 min on 32 threads. It is worth noting that the runtimes also scale in an independent manner with volume and maximum scale. Up to approximately 5 millions of Hα-emitting galaxies 0.9 < z < 1.1 are expected in the spectroscopic catalogue by the end of the Euclid Wide Survey, and according to our tests, the calculation would last a bit less than a day on 32 threads using the random split option. This represents a reasonable time and meets Euclid mission requirement, given that more than 32 threads will be usable in the SDC for the actual Euclid data processing. In particular, some SDC can provide up to 128 usable parallel threads.

We compare in Fig. 4 the runtimes for the correlation function multipole moments measured by 2PCF-GC with those obtained with the publicly available Corrfunc code v2.5.3 (Sinha & Garrison 2020). In this comparison, we use the same binning configuration and input catalogues. In addition to OpenMP, Corrfunc makes use of Single Instruction Multiple Data (SIMD) parallelism that permits vectorising some operations at the CPU level. This can allow a speed-up in the pair counting from about 50% to up to a factor of 2–3 depending on the CPU used (Sinha & Garrison 2020; Zhao 2023). Because computing infrastructures can vary across Euclid SDC and the constraint on portability of the different PF, no attempt was made to include SIMD parallelism in 2PCF-GC. Therefore, for a fair comparison we only used Corrfunc without SIMD vectorisation enabled. As shown in Fig. 4, in this configuration 2PCF-GC is slightly faster than Corrfunc by up to 30%, depending on the number of objects.

Fig. 4. Comparison of the runtimes for the calculation of the multipole correlation of galaxies in the ELM mocks using 2PCF-GC and Corrfunc publicly available code.

4.3. Accuracy tests

Within the validation of 2PCF-GC, we conducted a series of tests on the accuracy of the estimated 2PCF. We tested both the absolute accuracy and relative accuracy to estimates obtained from an external software. Those tests are presented in the following.

4.3.1. Absolute accuracy

We made a comprehensive analysis that uses a series of mocks for which the underlying 2PCF is perfectly known, namely the CoxMock suite of mocks. We generated 10 000 realisations allowing us to reach extremely precise summary statistics on the 2PCF. Pons-Bordería et al. (1999) give the following expression for the 2PCF of the line-point Cox process when r ≤ L:

$$\begin{array}{c}\hfill {\xi }_{\mathrm{cox}}(r)=\frac{1}{2\pi L{n}_{\mathrm{L}}}(\frac{1}{r^2}-\frac{1}{\mathit{rL}}),\end{array}$$(23)

where L is the line length and n_L is the line number density. Due to the finite number of lines and the fact that each point is placed on a randomly chosen line, for a given point, the probability that another point falls on another line is given by 1 − 1/N_L. This leads to an additional term $-N_{\mathrm{L}}^{-1}$ to the above expression for ξ_cox and in turn, the reference bin-averaged CoxMock isotropic correlation function becomes

$$\begin{array}{c}\hfill {\xi }_{\mathrm{ref}}(r_{\mathrm{low}},r_{\mathrm{up}})=\frac{3}{4\pi L^2{n}_{\mathrm{L}}}\frac{(r_{\mathrm{up}}-r_{\mathrm{low}})(2L-r_{\mathrm{up}}-r_{\mathrm{low}})}{r_{\mathrm{up}}^3-r_{\mathrm{low}}^3}-\frac{1}{N_{\mathrm{L}}},\end{array}$$(24)

where r_low and r_up are the lower and upper limits of the bin in r.

We measured the angle-averaged and multipole correlation functions in the mock realisations using 200 and 40 linear bins in r spanning the range [0, 200] h⁻¹ Mpc and including the random split option. In those mocks there is no anisotropic clustering and we recover vanishing quadrupole and hexadecapole moments with residual random variation around zero, as expected. To quantify the accuracy on the estimation, we calculate the mean bias relative to the standard deviation Δξ/σ and mean relative bias Δξ/ξ_ref, respectively defined as

$$\begin{array}{cc}\hfill \frac{\mathrm{\Delta }\xi (r)}{\sigma (r)}& =\frac{⟨\xi (r)⟩-{\xi }_{\mathrm{ref}}(r)}{\sigma (r)},\hfill \end{array}$$(25)

$$\begin{array}{cc}\hfill \frac{\mathrm{\Delta }\xi (r)}{{\xi }_{\mathrm{ref}}(r)}& =\frac{⟨\xi (r)⟩-{\xi }_{\mathrm{ref}}(r)}{{\xi }_{\mathrm{ref}}(r)},\hfill \end{array}$$(26)

where ⟨ξ(r)⟩ refers to the mean ξ(r) over the realisations, and σ corresponds to the reference statistical uncertainty in Euclid. The detailed calculation of the latter is given in Appendix A. The mean over the 10 000 realisations suppresses field stochasticity and sample variance. Both quantities are shown in Fig. 5. We find that the mean bias in the estimated correlation is extremely low, always below 10% of the statistical uncertainties for all considered scales. There are residual sample variance fluctuations, but overall we do not find any trend of bias. On closer inspection, we find that, except on scales below 10 h⁻¹ Mpc, the mean bias is always below about 5%. Similarly, the mean relative bias is always below 1% and can reach 0.1 per cent at the smallest scales. Finally, we can also estimate the mean bias relative to the uncertainty in the CoxMocks, where σ(r) is replaced by σ_cox(r) in Eq. (25). This is presented in Fig. 6, where we can see that intrinsically the relative accuracy on the 2PCF is essentially always within 3% for a volume of 5.268 h⁻³ Gpc³.

Fig. 5. Accuracy of the estimated real-space correlation function from the CoxMocks. While the top panel presents the mean bias relative to the standard deviation Δξ/σ, the bottom panel exhibits the absolute value of the mean relative bias |Δξ/ξref|. These estimates are obtained by averaging over 10 000 mock realisations. In the top panel, the horizontal coloured band represents ±10% of σ(r).

Fig. 6. Same as in the top panel of Fig. 5, but relative to the uncertainty in the CoxMocks. This estimate is obtained by averaging over 10 000 mock realisations. The horizontal coloured band represents ±3% of σcox(r).

4.3.2. Relative accuracy

We further compared the standard and modified correlation function multipole moments measured by 2PCF-GC with those obtained with Corrfunc. For this test we used the Flagship Galaxy Mock v1 galaxies at 0.9 < z < 1.1 (see Sect. 5) and considered 40 linear bins in r spanning the range [0, 200] h⁻¹ Mpc. We used the same catalogues as input of 2PCF-GC and Corrfunc. The difference in the monopole, quadrupole, and hexadecapole relative to the expected statistical error expected in Euclid, are presented in Fig. 7. For each multipole moment, the relative accuracy between the two estimates is of the order of the machine precision, both for the standard and modified 2PCF. In the modified 2PCF estimation, the reconstructed data and auxiliary random catalogues were obtained using the recSym BAO reconstruction algorithm, as detailed in Sarpa et al. (in prep.). The latter reconstruction algorithm uses Lagrangian Perturbation Theory to predict bulk-flow displacements and substracts them from the random point and galaxy positions (Padmanabhan et al. 2012; Burden et al. 2015).

Fig. 7. Difference in the monopole (top panel), quadupole (central panel), and hexadecapole (bottom panel) moments of the standard (ξℓ) and modified (ξℓm) correlation function between the 2PCF-GC and Corrfunc estimates, relative to the standard deviation. These measurements are obtained from a single mock realisation.

Finally, we tested the impact of using the random split option on the accuracy of the estimated 2PCF multipole moments. For this, we compared the measurements obtained in the same mock catalogue with the same setup as before while using or not this option. The difference in the monopole, quadrupole, and hexadecapole relative to the expected statistical error is shown in Fig. 8.

Fig. 8. Impact of the random split option on the accuracy of the measured 2PCF. The curve shows the difference in the monopole (top panel), quadupole (central panel), and hexadecapole (bottom panel) moments of the correlation function, while using or not the random split option with NS = 50. The shaded areas represent the regions encompassing respectively 10% and 20% of the expected statistical uncertainty. These measurements are obtained from a single mock realisation.

The results show that the random split option with N_S = 50 introduces random errors below 10% of the statistical uncertainties for all considered scales and no systematic error. It is worth emphasising that those uncertainties can be reduced by reducing the number of random splits, at the expense of longer runtimes.

Overall, these tests validates 2PCF-GC accuracy and the choice of using 50 times more random points than objects in the data catalogue and the usage of the random split option. We would like to emphasise that, even if not presented here, we also performed other accuracy and consistency tests on the cross-correlation function, and in all cases, we found similar, very good accuracies. This includes for instance consistency tests of the cross-correlation function where the two data catalogues are taken as two random sub-samples of a data catalogue. 2PCF-GC has been used to estimate cross-correlation functions in Risso et al. (in prep.).

4.4. Memory usage

The random-access memory (RAM) usage of 2PCF-GC is dominated by the storage of the input data and random object positions and weights. Furthermore, the memory required to process the data can increase depending on the amount of information comprised by the linked-list or tree structures. Although object properties are not directly copied into the spatial partitioning structures, additional metadata associated with the latter must be stored in memory. The linked-list structure, which comprises an ordered array of integers that convey information about the objects present in each cell, formally requires the least memory. The tree structures generally require the storage of more information. Each tree node contains several pieces of information including the coordinates of the bounding volume associated with the node, the indices of the objects contained in the node, and pointers to child nodes. This increases the amount of information to be stored in memory. However, the use of hash tables allows the pointer information to be removed, and in the case of the octree, the Morton code associated with each node allows the bounding volume coordinates to be computed directly, thus saving further memory.

In Fig. 9 we present estimates of the memory usage of 2PCF-GC from the ELM mock, both accounting for the RAM resident set size and swap space, for the three methods and as a function of the number of objects in the data catalogue. The memory usage increases with the number of objects and the storage of the input catalogues dominates the overall budget. The linked-list method tends to use less memory, but thanks to the implementation of hash tables for both k-d tree and octree, the memory usage is only marginally larger for the latter methods. Overall, for a data catalogue of five million objects and 50 times more in the random catalogue, about 19 GB of memory is required. In the case where the random split option is used, only a fraction of the random catalogue needs to be stored in memory at a time and thus the necessary memory for storing the random object information is significantly lower. In that case only about 0.9 GB of memory is required for a data catalogue containing five million objects.

Fig. 9. Memory usage as function of the number of objects in the data catalogue, for the considered data partitionings and while using or not the random split option. The different lines shows the results for the different pair-counting methods, while the different symbols indicate whether or not the random split option is used.

5. Expectations for Euclid Wide Survey

We present in this section some forecasts on 2PCF estimation and showcase some interesting features of 2PCF-GC using a realistic full-sky mock catalogue of the Euclid Wide Survey. We made use of the Flagship Galaxy Mock (FGM) v2.1 (Carretero et al. 2017; Tallada et al. 2020), which has been constructed from the Flagship 2 N-body simulation by populating dark matter halos with emission-line galaxies as described in Euclid Collaboration: Castander et al. (2025). A sub-sample of the galaxies with an Hα flux greater than 2 × 10⁻¹⁶ erg s⁻¹ cm⁻² and a redshift in the range 0.9 < z < 1.8 has been extracted, which represents the targeted galaxies for the spectroscopic sample of the Euclid Wide Survey. The mock consist of a light-cone covering an octant of the sky, and for the purpose of performing realistic forecasts, we replicated the octant in order to have full-sky catalogue. We note that this procedure does not introduce spurious clustering features on the scales smaller than 200 h⁻¹ Mpc, which are of interest for our analysis. This mock represents current most realistic expectation for the intrinsic clustering of Hα-emitting galaxies in the Euclid Wide Survey available within the Euclid collaboration.

5.1. Galaxy two-point correlation function across survey timeline

To illustrate the performance of 2PCF-GC and provide forecasts on the expected measured 2PCF at different stages of Euclid observations, we estimated the correlation function multipole moments in the FGM after 1, 3, and 6 years of observations, as defined in the Euclid Reference Survey³ (Euclid Collaboration: Scaramella et al. 2022). The angular extent of year 1, year 3, and year 6 observations as defined in the Euclid Reference Survey is presented in Fig. 10. We defined four redshift intervals: 0.9 < z < 1.1, 1.1 < z < 1.3, 1.3 < z < 1.5, 1.5 < z < 1.8, as used for Euclid cosmological forecasts (e.g. Euclid Collaboration: Blanchard et al. 2020).

Fig. 10. Sky coverage of Euclid Wide Field after one, three, and six years of observations, as defined in the Euclid Reference Survey (Euclid Collaboration: Scaramella et al. 2022).

The correlation function monopole, quadrupole, and hexadecapole are presented in Figs. 11, 12, and 13, respectively. Even if these measurements are extracted from a specific mock realisation of the observable universe, with its own sample variance, it gives a sense of the improvement on the estimation of the correlation function and high precision recovered after 6 years of observation. In particular, one can see from the monopole that the BAO peak in the lowest redshift interval is not very pronounced with the first year coverage, while at the end of the survey, the signal is much more significant. On scales beyond 100 h⁻¹ Mpc, we can see significant variations of the amplitude between year 1 and year 6, which can be attributed to sample variance, the fact that the over-density distribution associated with galaxies differs significantly. Similar trends are seen in the quadrupole and hexadecapole. Particularly for the latter, which is the most uncertain, one can see that the end-of-survey amplitude tends to vanish on large scales, as expected.

Fig. 11. Monopole correlation function estimated from the FGM mock for galaxies with Hα flux above 2 × 10−16 erg s−1 cm−2 at different epochs of observations. The different colours show the monopole in the redshift intervals: 0.9 < z < 1.1, 1.1 < z < 1.3, 1.3 < z < 1.5, 1.5 < z < 1.8.

Fig. 12. Same as Fig. 11, but for the quadrupole correlation function.

Fig. 13. Same as Fig. 11, but for the hexadecapole correlation function.

5.2. Impact of pair line-of-sight definition

In order to showcase the capabilities of 2PCF-GC, we estimated the multipole correlation functions in the FGM using the mid-point and end-point line-of-sight definition. Even if the estimator assumes the local plane-parallel approximation, the choice of the line-of-sight definition has different sensitivity to wide-angle effects. From its maximally-symmetric properties, the mid-point definition minimises the latter effects, while the end-point definition is more affected (e.g. Reimberg et al. 2016; Beutler et al. 2019). This can be seen in Fig. 14, which shows the absolute difference between using the mid-point and end-point definitions in the estimator, for the monopole and quadrupole correlation functions in the year 6 dataset. In particular, it is instructive to see that the absolute difference has no strong redshift dependence in the monopole. This difference remains within the 1σ statistical error expected in the completed Euclid Wide Survey at r < 100 h⁻¹ Mpc, shown with the dotted line in the figure. In the quadrupole instead, the effect can be more significant. In the two highest redshift intervals, the difference goes beyond the expected 2σ statistical error on the quadrupole. Overall, these results gives some insights on the typical wide-angle effects expected on large scales in the final Euclid Wide Survey sample when using the end-point line-of-sight definition.

Fig. 14. Impact of the choice of pair line-of-sight definition on the estimated monopole (upper panel) and quadrupole (lower panel) correlation functions in FGM mock. The curves show the absolute value of the difference between using mid-point and end-point definitions in the estimator for the monopole and quadrupole correlation functions. The different colours show this quantity for the redshift intervals: 0.9 < z < 1.1, 1.1 < z < 1.3, 1.3 < z < 1.5, 1.5 < z < 1.8. The dotted (dashed) line shows the expected 1σ (2σ) statistical error in the completed Euclid Wide Survey in the interval 1.1 < z < 1.3.

5.3. Integral constraint

The normalisation of the pair counts to the observed number of galaxies in the LS estimator imposes a constraint on the integral of the observed over-density field in the survey, which is that the latter vanishes. However, because of the finite size of the sample, the integral over the observed over-density field does not necessarily vanish and this leads to biasing negatively the correlation function on large scales (Peebles 1980). This effect, commonly referred to as the integral constraint, is particular important for small surveys. In the Euclid Wide Survey, which will be among the largest spectroscopic survey of galaxies, the expected effect is very small. Nonetheless, given the unprecedented statistical precision in the correlation function measurements expected in Euclid, it is important to assess the level and scales at which the integral constraint can impact the measurements.

The integral constraint introduces a bias on the observed correlation function and, at first order, the model for the latter is (e.g. de Mattia & Ruhlmann-Kleider 2019):

$$\begin{array}{c}\hfill \xi (r)-2⟨\delta ϵ⟩+⟨{ϵ}^2⟩\simeq \xi (r)-⟨{ϵ}^2⟩,\end{array}$$(27)

where ξ(⃗r) is the theoretical correlation function,

$$\begin{array}{c}\hfill ϵ={\displaystyle {\int }_{\mathrm{V}}}{\mathrm{d}}^{3}x\delta (x)W(x)/{\displaystyle {\int }_{\mathrm{V}}}{\mathrm{d}}^{3}xW(x),\end{array}$$(28)

δ(⃗x) is the over-density field, W(⃗x) is the survey window function, and ∫_Vd³x denotes an integral over the survey volume. The leading integral-constraint term in the correlation function is (de Mattia & Ruhlmann-Kleider 2019):

$$\begin{array}{c}\hfill ⟨{ϵ}^2⟩=\frac{{\displaystyle {\int }_{\mathrm{V}}}{\mathrm{d}}^3\mathrm{\Delta }\xi (\mathrm{\Delta })\mathcal{W}(\mathrm{\Delta })}{{\displaystyle {\int }_{\mathrm{V}}}{\mathrm{d}}^3\mathrm{\Delta }\mathcal{W}(\mathrm{\Delta })}=\frac{{\displaystyle {\int }_0^{\infty }}\mathrm{d}\mathrm{\Delta }{\mathrm{\Delta }}^2{\displaystyle {\sum }_{\ell =0}^{\infty }}\frac{{\xi }_{\ell }(\mathrm{\Delta })}{2\ell +1}{\mathcal{W}}_{\ell }(\mathrm{\Delta })}{{\displaystyle {\int }_0^{\infty }}\mathrm{d}\mathrm{\Delta }{\mathrm{\Delta }}^2{\mathcal{W}}_0(\mathrm{\Delta })},\end{array}$$(29)

where the survey window function correlation function and associated multipole moments are

$$\begin{array}{cc}\hfill \mathcal{W}(\mathrm{\Delta })& ={\displaystyle {\int }_{\mathrm{V}}}{\mathrm{d}}^{3}xW(x)W(x+\mathrm{\Delta }),\hfill \end{array}$$(30)

$$\begin{array}{cc}\hfill {\mathcal{W}}_{\ell }(\mathrm{\Delta })& =\frac{2\ell +1}{2}{\displaystyle {\int }_{-1}^1}\mathrm{d}{\mu }_{\mathrm{\Delta }}\mathcal{W}(\mathrm{\Delta })L_{\ell }({\mu }_{\mathrm{\Delta }}),\hfill \end{array}$$(31)

respectively. In the previous equations, ⃗Δ denotes a separation vector. By expanding Eq. (27) in multipole moments, one can see that the constant integral constraint term will only survive in the monopole correlation function. Moreover, by realising that the survey window function correlation function is directly related to random-random counts by RR(⃗Δ)∝Δ²𝒲(⃗Δ) (Wilson et al. 2017), Eq. (29) can be rewritten as

$$\begin{array}{cc}\hfill ⟨{ϵ}^2⟩& =\frac{{\displaystyle {\int }_0^{\infty }}\mathrm{d}\mathrm{\Delta }{\displaystyle {\sum }_{\ell =0}^{\infty }}{(2\ell +1)}^{-1}{\mathrm{RR}}_{\ell }(\mathrm{\Delta }){\xi }_{\ell }(\mathrm{\Delta })}{{\displaystyle \int \mathrm{d}}\mathrm{\Delta }{\mathrm{RR}}_0(\mathrm{\Delta })}\hfill \end{array}$$(32)

$$\begin{array}{cc}& =\frac{{\displaystyle {\sum }_{i=1}^{N_{\mathrm{b}}}}{\displaystyle {\sum }_{\ell =0}^{\infty }}{(2\ell +1)}^{-1}{\mathrm{RR}}_{\ell }({\mathrm{\Delta }}_i){\xi }_{\ell }({\mathrm{\Delta }}_i)}{{\displaystyle {\sum }_{i=1}^{N_{\mathrm{b}}}}{\mathrm{RR}}_0({\mathrm{\Delta }}_i)},\hfill \end{array}$$(33)

where in Eq. (33) the continuous integral is approximated by a discrete sum over bins⁴Δ_i, and N_b is the number of bins. In the previous equation, RR_ℓ is defined as

$$\begin{array}{c}\hfill {\mathrm{RR}}_{\ell }(r)=\frac{2\ell +1}{2}{\displaystyle {\int }_{-1}^1}\mathrm{d}\mu \mathrm{RR}(s,\mu )L_{\ell }(\mu ).\end{array}$$(34)

Therefore, by using Eq. (33) we can directly estimate the impact of the integral constraint from the random-random pair counts and a prior knowledge on the true underlying multipole correlation functions.

We compute the integral constraint correction expected in Euclid from Eq. (33) by using the RR counts estimated by 2PCF-GC and performing the discrete sum up to Δ = 1000 h⁻¹ Mpc. The correlation function multipoles in the integral are taken from the measured multipole (symbols) or by using an analytical model for the redshift-space multipole correlation functions. The latter was obtained by fitting the FGM measurements with most accurate full-shape models for the correlation function multipole moments, as described in Kärcher et al. (in prep.). The main advantages of using a model ξ_ℓ is to have a prediction up to the largest scales and that is not itself affected by the integral constraint. Overall, we found values of ⟨ϵ²⟩ varying within (3–8) × 10⁻⁵, depending on the epoch of observations and the redshift interval.

We present in Fig. 15 the leading integral-constraint term ⟨ϵ²⟩ estimated from the FGM, for the different redshift intervals and observing periods considered. We can see that the impact of the integral constraint depends on the size of the surveyed volume and the considered redshift interval. The effect is the most important at year 1 for all redshift intervals, as expected. At year 3 and then year 6, the effect diminishes. The predictions obtained from using the analytical ξ_ℓ in the integral constraint calculation for year 6 are shown in Fig. 15 with the solid curves. Those predictions are perfectly in line with those obtained from taking the measured ξ_ℓ, although on the largest scales, they provide more robust predictions as they are not affected by sample variance. In essentially all cases, the effect remains within the 1σ statistical uncertainty expected at year 6, as shown with the dotted line in Fig. 15. In future Euclid cosmological analyses, this effect can be straightforwardly accounted for at likelihood level, when confronting measurements to model predictions, or simply neglected if the analysis is limited to scales below about r = 150 h⁻¹ Mpc.

Fig. 15. Impact of the integral constraint on the monopole correlation function in the Euclid Wide Survey. The curves with the different different symbols show the leading integral-constraint term ⟨ϵ2⟩, after 1, 3, and 6 years of observation respectively. The solid curves show the same quantity for year 6 observations, but when the integral constraint effect is derived from a model correlation function. The different colours show this quantity for the redshift intervals: 0.9 < z < 1.1, 1.1 < z < 1.3, 1.3 < z < 1.5, 1.5 < z < 1.8. The dotted curves delineate the expected 1σ statistical uncertainty on the monopole correlation function in the completed Euclid Wide Survey in the interval 1.1 < z < 1.3.

As a final remark, we note that we did not consider here the impact of the sample purity and completeness, as well as possible radial integral constraint effects (de Mattia & Ruhlmann-Kleider 2019), which are beyond the scope of this paper and are the subject of forthcoming Euclid Collaboration papers.

6. Conclusions

This paper describes the development and validation of the processing function of the Science Ground Segment for estimating the 3D 2PCF, as part of the Euclid mission. The methodology employed in the software is both robust and effective. The data handling process has been optimised through innovative spatial partitioning strategies, which drastically reduce computation time while maintaining high accuracy. These optimisations are crucial for meeting the Euclid mission’s stringent scientific requirement, which demand speed and precision. This involves advanced algorithms such as k-d tree, octree, and linked-list, employed for pair counting and that ensure an efficient, yet exact, estimation of the 2PCF. The other key feature of the software is its optimisation and parallelisation capabilities. By leveraging the OpenMP application programming interface, the software achieves significant improvements in computational efficiency. These methods are indispensable for managing the extensive and complex dataset anticipated from the Euclid survey, which involves billions of galaxy pairs.

The software has undergone extensive validation using a variety of mock galaxy catalogues, which were designed to simulate the conditions and challenges of Euclid survey data. Rigorous tests have demonstrated that the software can handle the anticipated data volume and complexity, ensuring that the Euclid mission’s scientific goals can be met. The validation process has confirmed the software’s high accuracy and reliability in estimating the 2PCF, leading to an accuracy below that required by the Euclid mission. Moreover, the software’s integration into the Euclid SGS pipeline exemplifies its readiness for large-scale deployment and its pivotal role in the mission’s data processing framework for the production of highest-level scientific products for cosmology. This integration ensures that the vast amount of data collected by the Euclid mission can be processed efficiently.

We have presented forecasts of 2PCF-GC performance in measuring the 2PCF at different stages of the Euclid Wide Survey. These forecasts are based on a realistic mock dataset and illustrated the robustness and efficiency in estimating the 2PCF. This mock dataset has also allowed us to investigate the impact of different pair line-of-sight definitions on the estimation of the multipole correlation functions, highlighting differences in sensitivity to wide-angle effects over the scales that will be reached by Euclid after six years of observations. We found that deviations in the monopole and quadrupole correlation functions are less than 5% below 100 h⁻¹ Mpc, but can exceed 25% at larger scales. We have also estimated the integral constraint bias affecting Euclid Wide Survey measurements. While this is expected to be minimal in the Euclid Wide Survey, due to the unprecedented precision required in the Euclid measurements, the assessment of the magnitude and extent of the integral constraint effect remains critical. Using both measured and analytically modelled multipole correlation functions, we estimated that the integral constraint effect is of the order of 5 × 10⁻⁵, depending on the redshift interval and observing epoch considered. Our results indicate that the effect is most pronounced in the first year of the survey, reaching 10% at the BAO scale and up to 30% at 200 h⁻¹ Mpc, but it diminishes in subsequent years. Nonetheless, future works will need to consider additional factors in the estimation of the 2PCF, such as sample purity, completeness, and potential radial integral constraint effects, which are beyond the scope of this paper. In particular, a first assessment of the impact of expected redshift errors on the 2PCF indicate that the results shown in this paper should be considered optimistic (Risso et al., in prep.).

Acknowledgments

The authors would like to thank Jean-Charles Lambert for his help in using and extracting reliable computational measures from the LAM computer cluster. This work was supported by the ASI/INAF agreement n. 2018-23-HH.0 “Scientific activity for Euclid mission, Phase D”, the MIUR, PRIN 2017 research grant ‘From Darklight to DM: understanding the galaxy/matter connection to measure the Universe’ and the INFN project “InDark”. FM acknowledges the financial contribution from the grant PRIN-MUR 2022 20227RNLY3 ‘The concordance cosmological model: stress-tests with galaxy clusters’ supported by Next Generation EU and from the grant ASI n. 2024-10-HH.0 ‘Attività scientifiche per la missione Euclid – fase E’. This work has made use of CosmoHub. CosmoHub is developed and maintained by PIC, IFAE, CIEMAT, in collaboration with ICE-CSIC. It is partially financed by the European Union NextGenerationEU(PRTR-C17.I1) and by Generalitat de Catalunya, as well as by the grant EQC2021-007479-P funded by MCIN/AEI/10.13039/501100011033 and by the “European Union NextGenerationEU/PRTR”. The Euclid Consortium acknowledges the European Space Agency and a number of agencies and institutes that have supported the development of Euclid, in particular the Agenzia Spaziale Italiana, the Austrian Forschungsförderungsgesellschaft funded through BMK, the Belgian Science Policy, the Canadian Euclid Consortium, the Deutsches Zentrum für Luft- und Raumfahrt, the DTU Space and the Niels Bohr Institute in Denmark, the French Centre National d’Etudes Spatiales, the Fundação para a Ciência e a Tecnologia, the Hungarian Academy of Sciences, the Ministerio de Ciencia, Innovación y Universidades, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Research Council of Finland, the Romanian Space Agency, the State Secretariat for Education, Research, and Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom Space Agency. A complete and detailed list is available on the Euclid web site (https://www.euclid-ec.org).

References

Appendix A: Reference statistical uncertainty on the multipole two-point correlation functions in Euclid

The expected statistical uncertainty on the multipole 2PCF is based on the expected Euclid survey characteristics, number density, and clustering of observed Hα-emitting galaxies. Based on the assumptions detailed in Table 1–3 of Euclid Collaboration: Blanchard et al. (2020), 4 redshift intervals have been considered for which the observed volume, galaxy mean number density, and galaxy linear bias are given in Table A.1. With this information and a theoretical prediction for the matter power spectrum, we can calculate the expected Gaussian covariance matrix for the galaxy correlation function monopole, quadrupole, and hexadecapole using the method of Grieb et al. (2016). For the purpose of testing the accuracy on the 2PCF, we considered the case of the interval at 1.1 < z < 1.3 as representative of Euclid. The standard deviations σ₀(r), σ₂(r), and σ₄(r) were then obtained by taking the square root of the diagonal elements of the covariance matrices of ξ₀(r), ξ₂(r), and ξ₄(r), respectively. We assumed σ(r) = σ₀(r), where σ(r) is associated with the real-space 2PCF. Those reference uncertainties are shown in Fig. A.1.

Fig. A.1. Reference statistical 1σ uncertainty on the monopole, quadrupole, and hexadecapole 2PCF for the completed Euclid Wide Survey.

All Tables

All Figures

	Fig. 1. Geometry and separation definition for end-point (left), mid-point (middle), and bisector (right) pair line-of-sight conventions. The observer is at the point ‘O’ and the two other points represent the objects of the pair. In those figures, all vectors and segments are contained in the plane determined by the vectors ⃗x and ⃗r.
In the text

	Fig. 2. 2PCF-GC process overview. The different boxes illustrate the different steps involved in the processing.
In the text

Fig. 3. Run times for the calculation of the multipole correlation of galaxies obtained from the ELM mock. The times are expressed in CPU-hour and as a function of the data or random catalogue size. The various symbols represent the time spent on the data structure construction, DD calculation, DR calculation, RR calculation, and the overall 2PCF runtime with and without random split option. The DR calculation times are provided as a function of the number of objects in the data catalogue and assuming a fifty times larger random catalogue. The different curves represent the runtimes obtained with the linked-list (solid), k-d tree (dashed), and octree (dotted) algorithms. The blue vertical band shows the range of expected number densities in the spectroscopic sample at redshifts within 0.9 < z < 1.8. The abscissa refers to the number of object in the data catalogue except for RR calculation where it refers to that in the random catalogue.

In the text

	Fig. 4. Comparison of the runtimes for the calculation of the multipole correlation of galaxies in the ELM mocks using 2PCF-GC and Corrfunc publicly available code.
In the text

	Fig. 5. Accuracy of the estimated real-space correlation function from the CoxMocks. While the top panel presents the mean bias relative to the standard deviation Δξ/σ, the bottom panel exhibits the absolute value of the mean relative bias |Δξ/ξref|. These estimates are obtained by averaging over 10 000 mock realisations. In the top panel, the horizontal coloured band represents ±10% of σ(r).
In the text

	Fig. 6. Same as in the top panel of Fig. 5, but relative to the uncertainty in the CoxMocks. This estimate is obtained by averaging over 10 000 mock realisations. The horizontal coloured band represents ±3% of σcox(r).
In the text

	Fig. 7. Difference in the monopole (top panel), quadupole (central panel), and hexadecapole (bottom panel) moments of the standard (ξℓ) and modified (ξℓm) correlation function between the 2PCF-GC and Corrfunc estimates, relative to the standard deviation. These measurements are obtained from a single mock realisation.
In the text

Fig. 8. Impact of the random split option on the accuracy of the measured 2PCF. The curve shows the difference in the monopole (top panel), quadupole (central panel), and hexadecapole (bottom panel) moments of the correlation function, while using or not the random split option with NS = 50. The shaded areas represent the regions encompassing respectively 10% and 20% of the expected statistical uncertainty. These measurements are obtained from a single mock realisation.

In the text

	Fig. 9. Memory usage as function of the number of objects in the data catalogue, for the considered data partitionings and while using or not the random split option. The different lines shows the results for the different pair-counting methods, while the different symbols indicate whether or not the random split option is used.
In the text

	Fig. 10. Sky coverage of Euclid Wide Field after one, three, and six years of observations, as defined in the Euclid Reference Survey (Euclid Collaboration: Scaramella et al. 2022).
In the text

	Fig. 11. Monopole correlation function estimated from the FGM mock for galaxies with Hα flux above 2 × 10−16 erg s−1 cm−2 at different epochs of observations. The different colours show the monopole in the redshift intervals: 0.9 < z < 1.1, 1.1 < z < 1.3, 1.3 < z < 1.5, 1.5 < z < 1.8.
In the text

	Fig. 12. Same as Fig. 11, but for the quadrupole correlation function.
In the text

	Fig. 13. Same as Fig. 11, but for the hexadecapole correlation function.
In the text

Fig. 14. Impact of the choice of pair line-of-sight definition on the estimated monopole (upper panel) and quadrupole (lower panel) correlation functions in FGM mock. The curves show the absolute value of the difference between using mid-point and end-point definitions in the estimator for the monopole and quadrupole correlation functions. The different colours show this quantity for the redshift intervals: 0.9 < z < 1.1, 1.1 < z < 1.3, 1.3 < z < 1.5, 1.5 < z < 1.8. The dotted (dashed) line shows the expected 1σ (2σ) statistical error in the completed Euclid Wide Survey in the interval 1.1 < z < 1.3.

In the text

Fig. 15. Impact of the integral constraint on the monopole correlation function in the Euclid Wide Survey. The curves with the different different symbols show the leading integral-constraint term ⟨ϵ2⟩, after 1, 3, and 6 years of observation respectively. The solid curves show the same quantity for year 6 observations, but when the integral constraint effect is derived from a model correlation function. The different colours show this quantity for the redshift intervals: 0.9 < z < 1.1, 1.1 < z < 1.3, 1.3 < z < 1.5, 1.5 < z < 1.8. The dotted curves delineate the expected 1σ statistical uncertainty on the monopole correlation function in the completed Euclid Wide Survey in the interval 1.1 < z < 1.3.

In the text

	Fig. A.1. Reference statistical 1σ uncertainty on the monopole, quadrupole, and hexadecapole 2PCF for the completed Euclid Wide Survey.
In the text