EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 707 (March 2026) A&A, 707 (2026) A228 Full HTML
Open Access
Issue
A&A
Volume 707, March 2026
Article Number A228
Number of page(s) 21
Section Cosmology (including clusters of galaxies)
DOI https://doi.org/10.1051/0004-6361/202556177
Published online 17 March 2026

© The Authors 2026

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. This email address is being protected from spambots. You need JavaScript enabled to view it. to support open access publication.

1. Introduction

The Euclid mission (Laureijs et al. 2011) is, along with the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration 2016), Vera Rubin Observatory (Ivezic et al. 2009), and Nancy Grace Roman Telescope initiatives (Dore et al. 2019), a stage IV (Albrecht et al. 2006) observational campaign expected to significantly advance our understanding of the Universe by mapping cosmological perturbations on an unprecedented volume. Euclid will combine two major probes of the large-scale structure: galaxy clustering and weak lensing. The scientific potential of this combination, along with a description of the instrumentation that will enable it, is described in the recent overview paper Euclid Collaboration: Mellier et al. (2025). The spectroscopic galaxy sample, in particular, will cover the redshift range 0.9 ≤ z ≤ 1.8, collecting redshift measurements for millions of Hα-emitting galaxies across a total area of 14 000 deg2.

Galaxy clustering exploits the statistical properties of the fluctuations in the galaxy distribution at large scales, measuring and analysing, in the standard approach, its correlation functions, starting with the two-point correlation function (2PCF) or its Fourier-space counterpart, the power spectrum. A major contribution to cosmological studies over the past two decades came from using the baryonic acoustic oscillations (BAOs) present in the 2PCF as a standard ruler to reconstruct the background expansion (Seo & Eisenstein 2003; Eisenstein et al. 2005; Percival et al. 2007; Wang et al. 2017; Zhao et al. 2017; Adame et al. 2025). Additional constraints come from the anisotropy induced on two-point correlators by redshift-space distortions (Peacock et al. 2001; Guzzo et al. 2008; Beutler et al. 2017; Grieb et al. 2017; Pezzotta et al. 2017; Hou et al. 2018) and, more generally, from the full-shape analysis of their redshift-space multipoles, with the aim of extracting all available information from the clustering measurements (Sánchez et al. 2013, 2017; D’Amico et al. 2020; Ivanov et al. 2020; Tröster et al. 2020; Adame et al. 2025).

In recent years, the inclusion in galaxy clustering analyses of higher-order correlation functions – such as the three-point correlation function (3PCF) and, in Fourier space, the bispectrum – has become increasingly standard. These statistics quantify the non-Gaussian properties of the galaxy distribution as a random field. Different sources of non-Gaussianity are directly related to non-linearities in the evolution of matter perturbations (Fry 1984), galaxy bias (Fry & Gaztañaga 1993; Fry 1994; Frieman & Gaztanaga 1994), and redshift-space distortions (Hivon et al. 1995; Scoccimarro et al. 1999), or possibly to a primordial non-Gaussian component due to inflationary physics (Verde et al. 2000; Scoccimarro 2000; Scoccimarro et al. 2004). The joint analysis of the galaxy power spectrum and bispectrum has been performed over the last few years in several studies involving Baryon Oscillations Spectroscopic Survey (BOSS) datasets (Dawson et al. 2013) datasets, leading to improvements in constraints on cosmological parameters in the context of the standard model and its extensions (see e.g. Gil-Marín et al. 2017; D’Amico et al. 2020; Philcox & Ivanov 2022) and better constraints on initial conditions (D’Amico et al. 2025; Cabass et al. 2022).

The full exploitation of the information provided by the 3PCF in configuration space, despite its clear advantages over the bispectrum in terms of accounting for survey geometry effects, has long been hampered by the computational cost of its estimation in large datasets. Early studies focused on a subset of all potentially measurable configurations, focusing on the galaxy properties and non-linear bias (Jing & Börner 2004; Gaztañaga et al. 2005; McBride et al. 2011; Marín 2011; Marín et al. 2013; Moresco et al. 2017) and the detection of acoustic features (Gaztañaga et al. 2009; Moresco et al. 2021). The introduction of an estimator based on a spherical harmonic expansion (Slepian & Eisenstein 2015, 2018) has substantially reduced the computational burden, enabling a full exploitation of the information potentially encoded in the 3PCF (Slepian et al. 2017, 2018).

Yet, in a likelihood analysis, an additional disadvantage of the 3PCF compared to the bispectrum is still present at the level of the evaluation of the theoretical model. Predictions for configuration-space statistics are typically obtained from a first evaluation of their Fourier-space counterpart. For the 2PCF, this mapping is efficiently handled using the FFTLog algorithm (Hamilton 2000) in one dimension. In the case of the 3PCF, while some works applied the 1D fast Fourier transform (FFT) to leading-order (LO) perturbative predictions (Slepian & Eisenstein 2017; Sugiyama et al. 2021; Veropalumbo et al. 2022), more general methods instead adopted a 2D FFT (Fang et al. 2020) to transform any bispectrum model, when expanded into spherical harmonics, into its configuration-space equivalent (Umeh 2021; Guidi et al. 2023; Farina et al. 2026; Pugno et al. 2025).

Cosmological perturbation theory (PT) plays a crucial role in our understanding of galaxy clustering at large scales (see Bernardeau et al. 2002, for a classical review). Significant effort has been made in the past decade to identify all potential effects that could lead to biased estimates of cosmological parameters if ignored. This includes the impact of small-scale matter fluctuations on large-scale clustering; this motivated the development of the effective field theory of large-scale structure (EFTofLSS; Baumann et al. 2012; Carrasco et al. 2012), non-local (Chan et al. 2012; Baldauf et al. 2012), re-normalised (McDonald 2006; McDonald & Roy 2009; Assassi et al. 2014; Eggemeier et al. 2019), higher derivative bias (Bardeen et al. 1986; Fujita et al. 2020), and infrared re-summation (Blas et al. 2016b). Here we set aside, for the sake of brevity, redshift-space distortion effects, which are extremely important in actual data analysis but not relevant for this work. The state-of-the-art model, which is a perturbative prediction at the next-to-leading order (NLO) for the power spectrum and at the LO for the bispectrum, is at the basis of current full-shape analyses of the galaxy power spectrum, bispectrum, and 2PCF. A full-shape analysis involving the 3PCF (on equal footing) is still lacking due to the computational burden of its numerical evaluation.

This work is part of a series of papers in preparation for the Euclid mission. The series provides an extensive validation of analytical and, in some cases, numerical models for galaxy clustering statistics. Paper I focused on the real-space power spectrum (Euclid Collaboration: Pezzotta et al. 2024).

Here we instead explore the range of validity of perturbative models in describing real-space measurements of the 2PCF and 3PCF in synthetic catalogues that reproduce the density and clustering amplitude expected for the Euclid spectroscopic sample. In addition, we apply such models to a full-shape joint analysis of the two statistics. To do so, we developed an emulator for both 2PCF and 3PCF predictions as a function of three cosmological parameters, the amplitude of scalar perturbations (As), the cold dark matter (CDM) density parameter (ωcdm), and the reduced Hubble parameter (h). This constitutes an important step towards a complete, configuration-space analysis of galaxy clustering datasets in general and towards the goals of the Euclid mission more specifically.

This paper is organised as follows. In Sects. 2 and 3 we introduce the basic definitions of Fourier-space and configuration-space correlators along with a brief overview of their theoretical model in PT. Section 4 presents the synthetic galaxy catalogues and the simulation on which they are based, the 2PCF and 3PCF measurements, and the theoretical model for the relative covariance matrices. In Sect. 5 we describe the fitting procedure and the performance metrics employed to assess model accuracy and precision. Our results are presented in Sects. 6 and 7: we first fix the cosmological parameters to assess the goodness of fit (GoF) of the model and then explore the full-shape analysis of 2PCF and 3PCF in order to constrain the same parameters. Lastly, in Sect. 8 we present our conclusions. In Appendix B we include a more detailed description of the 2PCF and 3PCF emulator and its validation.

2. Definitions and conventions

For convenience, we use the following notation to refer to the integration over the infinite volume of a loop variable, q,

∫ q : = ∫ d 3 q ( 2 π ) 3 , Mathematical equation: $$ \begin{aligned} \int _\mathbf{q } := \int \frac{\mathrm{{d^3 q}}}{(2\pi )^3} , \end{aligned} $$(1)

and we also adopt the following convention for the Fourier transform and its inverse

δ ( k ) : = ( 2 π ) 3 ∫ q δ ( q ) e i k · q , Mathematical equation: $$ \begin{aligned} \delta (\mathbf k )&:= (2\pi )^3 \int _\mathbf{q } \,\delta (\mathbf q ) \, \mathrm{e} ^{\mathrm{i}\mathbf{k \cdot \mathbf{q }}}, \end{aligned} $$(2)

δ ( x ) : = ∫ k δ ( k ) e − i k · x . Mathematical equation: $$ \begin{aligned} \delta (\mathbf x )&:= \int _\mathbf{k } \,\delta (\mathbf k ) \, \mathrm{e} ^{- \mathrm{i}\mathbf k \cdot \mathbf{x }}. \end{aligned} $$(3)

Under the assumption of statistical homogeneity and isotropy, the power spectrum P(k) of a generic density contrast δ(k) is therefore defined as

⟨ δ ( k 1 ) δ ( k 2 ) ⟩ : = ( 2 π ) 3 δ D ( k 1 + k 2 ) P ( k ) , Mathematical equation: $$ \begin{aligned} \langle \delta (\mathbf k _1)\,\delta (\mathbf k _2) \rangle := \, (2 \pi )^3 \delta _{\mathrm{D} }(\mathbf k _1 \!+\! \mathbf k _2)\,P(k), \end{aligned} $$(4)

where k := |k1|. The corresponding two-point function ξ(r) := ⟨δ(x)δ(x + r)⟩, where r := |r|, can be obtained as its Fourier transform, is

ξ ( r ) = ∫ k P ( k ) e i k · r = ∫ d k 2 π 2 k 2 j 0 ( k r ) P ( k ) , Mathematical equation: $$ \begin{aligned} \xi (r) = \int _\mathbf{k } P(k)\ \mathrm{e} ^{\mathrm{i} \mathbf k \cdot \mathbf r } = \int \frac{{\mathrm{d} } k}{2 \pi ^2} k^2 \, j_0(kr) \, P(k) , \end{aligned} $$(5)

where jn represent the spherical Bessel functions of n-th order. Similarly, the bispectrum B(k1, k2, k3), i.e. the 3PCF of the Fourier-space density contrast δ(k), is defined as

⟨ δ ( k 1 ) δ ( k 2 ) δ ( k 3 ) ⟩ : = ( 2 π ) 3 δ D ( k 1 + k 2 + k 3 ) B ( k 1 , k 2 , k 3 ) . Mathematical equation: $$ \begin{aligned} \langle \delta (\mathbf k _1) \delta (\mathbf k _2) \delta (\mathbf k _3)\rangle := (2 \pi )^3 \, \delta _\mathrm{D} (\mathbf k _1 \!+ \! \mathbf k _2 \!+ \! \mathbf k _3)\,B(k_1, k_2, k_3). \end{aligned} $$(6)

Here the Dirac delta distribution forces the three wave vectors {k1, k2, k3} to form a closed triangle, so that k3 can be written as a function of k1, k2 and the cosine of the angle between k1 and k2, denoted as μ12 := k1 ⋅ k2/k1k2. We can therefore consider an expansion with such an angle dependence in Legendre polynomials,

B ( k 1 , k 2 , k 3 ) = ∑ ℓ B ℓ ( k 1 , k 2 ) L ℓ ( μ 12 ) , Mathematical equation: $$ \begin{aligned} B(k_1,k_2,k_3)=\sum _\ell B_\ell (k_1,k_2)\,\mathcal{L} _{\ell }(\mu _{12}), \end{aligned} $$(7)

with coefficients defined as

B ℓ ( k 1 , k 2 ) = 2 ℓ + 1 2 ∫ − 1 1 d μ 12 B ( k 1 , k 2 , k 3 ) L ℓ ( μ 12 ) . Mathematical equation: $$ \begin{aligned} B_\ell (k_1, k_2) = \frac{2 \ell + 1}{2} \int _{-1}^{1} \mathrm{d} \mu _{12} \, B (k_1, k_2, k_3) \,\mathcal{L} _{\ell }(\mu _{12}) . \end{aligned} $$(8)

The 3PCF, in configuration space, is defined, adopting the notation rij := ri − rj, as

ζ ( r 12 , r 13 , r 23 ) : = ⟨ δ ( r 1 ) δ ( r 2 ) δ ( r 3 ) ⟩ , Mathematical equation: $$ \begin{aligned} \zeta (r_{12}, r_{13}, r_{23}) := \langle \delta (\mathbf r _1)\,\delta (\mathbf r _2)\, \delta (\mathbf r _3)\rangle , \end{aligned} $$(9)

and, in the same way, the dependence on the angle between r12 and r13 can be expanded in Legendre polynomials as

ζ ( r 12 , r 13 , r 23 ) = ∑ ℓ ζ ℓ ( r 12 , r 13 ) L ℓ ( μ 12 , 13 ) , Mathematical equation: $$ \begin{aligned} \zeta (r_{12}, r_{13}, r_{23}) = \sum _{\ell }\,\zeta _\ell (r_{12}, r_{13} ) \mathcal{L} _\ell (\mu _{12,13}) , \end{aligned} $$(10)

where now we define the cosine μ12, 13 := r12 ⋅ r13/(r12r13) and the coefficients are given by

ζ ℓ ( r 12 , r 13 ) : = 2 ℓ + 1 2 ∫ − 1 1 d μ 12 , 13 ζ ( r 12 , r 13 , r 23 ) L ℓ ( μ 12 , 13 ) . Mathematical equation: $$ \begin{aligned} \zeta _{\ell }(r_{12}, r_{13}) := \frac{2 \ell + 1}{2} \!\int _{-1}^{1} \mathrm{d} \mu _{12,13}\, \zeta (r_{12}, r_{13}, r_{23}) \,\mathcal{L} _{\ell }(\mu _{12,13}). \end{aligned} $$(11)

The relation between the 3PCF and the bispectrum can be expressed in terms of their respective multipoles Bℓ and ζℓ as

ζ ℓ ( r 12 , r 13 ) = ( − 1 ) ℓ ∫ d k 1 d k 2 ( 2 π ) 6 B ℓ ( k 1 , k 2 ) j ℓ ( k 1 r 12 ) j ℓ ( k 2 r 13 ) . Mathematical equation: $$ \begin{aligned} \zeta _{\ell }(r_{12}, r_{13}) = (-1)^{\ell }\int \frac{\mathrm{d} k_1 \, \mathrm{d} k_2}{(2 \pi )^6} B_{\ell }(k_1, k_2)\,j_\ell (k_1r_{12}) \,j_\ell (k_2r_{13}). \end{aligned} $$(12)

As we will see, the expansion in multipoles will be crucial to defining an efficient estimator for the 3PCF and for a fast evaluation of its theoretical prediction.

3. Theoretical models

This section provides a brief introduction to the real-space modelling of the 2PCF and 3PCF in cosmological PT (Bernardeau et al. 2002). We considered specifically a one-loop expression for the 2PCF accounting for a general bias expansion (Desjacques et al. 2018), corresponding to the EFTofLSS expression described in a companion paper (Euclid Collaboration: Kärcher et al., in prep.), to which we refer the reader for further details. The prediction for the 3PCF is instead at tree level in PT.

3.1. Eulerian perturbation theory and galaxy bias

In Eulerian PT, the equations describing the evolution of the matter density and velocity field, which are continuity and Euler equations, are solved perturbatively assuming the density contrast δ(x) to be small at large scales. The non-linear solution is given by the sum of the solution to the linearised equations plus higher-order corrections.

In addition, the relation between the galaxy density contrast δg(x) and the matter field δ(x) is also given in terms of a perturbative expansion. The terms relevant for the one-loop power spectrum and 2PCF models are given by

δ g ( x ) = & b 1 δ ( x ) + b 2 2 δ 2 ( x ) + b ∇ 2 ∇ 2 δ ( x ) + b G 2 G 2 ( Φ v | x ) + b Γ 3 Γ 3 ( x ) , Mathematical equation: $$ \begin{aligned} \delta _{\rm g}(\mathbf x ) = \&b_1\,\delta (\mathbf x ) +\frac{b_2}{2}\delta ^{\,2}(\mathbf x ) + b_{\nabla ^2} \nabla ^2 \delta (\mathbf x ) \nonumber \\&+ b_{\mathcal{G} _2}\mathcal{G} _2\left(\Phi _\mathrm{v} \,|\,\mathbf x \right) + b_{\Gamma _3}\Gamma _3(\mathbf x ), \end{aligned} $$(13)

where 𝒢2 and Γ3 are non-local operators, defined as

G 2 ( Φ | x ) : = [ ∂ i ∂ j Φ ( x ) ] 2 − [ ∂ 2 Φ ( x ) ] 2 , Mathematical equation: $$ \begin{aligned}&\mathcal{G} _2 (\Phi |\mathbf x ) := \big [\partial _i \partial _j \Phi (\mathbf x ) \big ]^2 - [\partial ^2 \Phi (\mathbf x )]^2, \end{aligned} $$(14)

Γ 3 ( x ) : = G 2 ( Φ | x ) − G 2 ( Φ v | x ) , Mathematical equation: $$ \begin{aligned}&\Gamma _3(\mathbf x ) := \mathcal{G} _2(\Phi |\mathbf x ) - \mathcal{G} _2(\Phi _v |\mathbf x ), \end{aligned} $$(15)

and Φ(x) and Φv(v) represent the gravitational and velocity potential. The bias relation in Eq. (13) consists of all operators built from second derivatives of the gravitational and velocity potential. It includes the linear (b1) and quadratic (b2), local bias contributions (Kaiser 1984; Szalay 1988; Coles et al. 1993; Fry & Gaztañaga 1993), omitting the cubic local operator that leads to a correction that can be absorbed in a re-normalised linear bias parameter (Szalay 1988; McDonald 2006; McDonald & Roy 2009; Assassi et al. 2014; Eggemeier et al. 2019). It also includes non-local contributions (b𝒢2, bΓ3) induced by non-linear evolution (Catelan et al. 1998; Chan et al. 2012; Baldauf et al. 2012) and higher derivative correction to linear bias (b∇2; Bardeen et al. 1986; Desjacques 2008; Fujita et al. 2020). For a detailed review, see Desjacques et al. (2018).

The most general and conservative assumption in fitting the galaxy correlation models to their measurements consist in considering all bias parameters as free, independent parameters to be determined by the data along with the cosmological parameters. On the other hand, we do know that the bias parameters are correlated, and large portions of the parameter space are in fact non-physical. We can take advantage of specific relations among the bias parameters to reduce the number of nuisance parameters. We considered two such relations. The first, given by

b G 2 ex − set = 0.524 − 0.547 b 1 + 0.046 b 1 2 , Mathematical equation: $$ \begin{aligned} b_{\mathcal{G} _2}^\mathrm{{ex-set}} = 0.524 - 0.547\ b_1 +0.046 \ b_1^2 , \end{aligned} $$(16)

is a quadratic fit obtained by Eggemeier et al. (2020) to the excursion-set prediction of Sheth et al. (2013). The second relation, given instead by

b Γ 3 coev = − 1 6 ( b 1 − 1 ) − 5 2 b G 2 + b Γ 3 L , Mathematical equation: $$ \begin{aligned} b_{\Gamma _3}^\mathrm{{coev}} = -\frac{1}{6}(b_1 -1) - \frac{5}{2} b_{\mathcal{G} _2} + b_{{\Gamma _3}}^{\mathcal{L} }, \end{aligned} $$(17)

is derived in Eggemeier et al. (2019) assuming the evolution of conserved galaxy number density (hence, ‘coevolution’) after formation, with the subscript ℒ indicates the corresponding Lagrangian quantities at the formation moment.

The resulting, perturbative expression for the galaxy power spectrum in real space, Pg(k), omitting stochastic contributions irrelevant for the 2PCF prediction, then reads

P g ( k ) = P g tree ( k ) + P g 1 − loop ( k ) + P g ctr ( k ) . Mathematical equation: $$ \begin{aligned} P_{\rm g}(k)=P_{\rm g}^{\,\mathrm {tree}}(k) + P_{\rm g}^{\,\mathrm {1{-}loop}}(k) + P_{\rm g}^{\,\mathrm {ctr}}(k) . \end{aligned} $$(18)

The linear, tree level, leading term is simply proportional to the linear matter power spectrum,

P g tree ( k ) = b 1 2 P L ( k ) , Mathematical equation: $$ \begin{aligned} P_{\rm g}^{\,\mathrm {tree}}(k)=b_1^{\,2} \, P_{\rm L}(k), \end{aligned} $$(19)

while the one-loop correction in standard PT is

P g 1 − loop ( k ) = P g , 22 ( k ) + P g , 13 ( k ) Mathematical equation: $$ \begin{aligned} P_{\rm g}^{\,\mathrm {1{-}loop}}(k) =&\ P_{\mathrm{g}, 22}(k)+P_{\mathrm{g}, 13}(k) \end{aligned} $$(20)

= 2 ∫ q K 2 2 ( q , k − q ) P L ( | k − q | ) P L ( q ) + 6 b 1 P L ( k ) ∫ q K 3 ( q , − q , k ) P L ( q ) , Mathematical equation: $$ \begin{aligned} =&\ 2\int _\mathbf{q }K_2^{\,2}\left( \mathbf q ,\mathbf k -\mathbf q \right) \, P_{\rm L} \left( \left| \mathbf k -\mathbf q \right| \right) \, P_{\rm L}(q) \nonumber \\&+ 6 \, b_1 \, P_{\rm L}(k) \int _\mathbf{q } K_3\left( \mathbf q ,-\mathbf q ,\mathbf k \right) \, P_{\rm L}(q) , \end{aligned} $$(21)

where the kernels describing non-linear matter density evolution and non-linear bias are

K 2 ( k 1 , k 2 ) = b 1 F 2 ( k 1 , k 2 ) + 1 2 b 2 + b G 2 S ( k 1 , k 2 ) , Mathematical equation: $$ \begin{aligned} K_2(\mathbf k _1, \mathbf k _2) =&\, b_1 F_2(\mathbf k _1, \mathbf k _2) + \frac{1}{2}\,b_2 + b_{\mathcal{G} _2} S(\mathbf k _1, \mathbf k _2) , \end{aligned} $$(22)

K 3 ( k 1 , k 2 , k 3 ) = b 1 F 3 ( k 1 , k 2 , k 3 ) + b 2 F 2 ( k 1 , k 2 ) + 2 b G 2 S ( k 1 , k 12 ) F 2 ( k 2 , k 3 ) + 2 b Γ 3 S ( k 1 , k 12 ) [ F 2 ( k 2 , k 3 ) − G 2 ( k 2 , k 3 ) ] , Mathematical equation: $$ \begin{aligned} K_3(\mathbf k _1, \mathbf k _2, \mathbf k _3) =&\, b_1 F_3(\mathbf k _1, \mathbf k _2, \mathbf k _3) + b_2 F_2(\mathbf k _1, \mathbf k _2) \nonumber \\&+ 2 b_{\mathcal{G} _2}\,S(\mathbf k _1, \mathbf k _{12})\,F_2(\mathbf k _2, \mathbf k _3) \nonumber \\&+ 2b_{\Gamma _3}\,S(\mathbf k _1, \mathbf k _{12})\,\big [F_2(\mathbf k _2, \mathbf k _3) - G_2(\mathbf k _2, \mathbf k _3) \big ] , \end{aligned} $$(23)

with

S ( k 1 , k 2 ) = ( k 1 · k 2 ) 2 k 1 2 k 2 2 − 1 , Mathematical equation: $$ \begin{aligned} S(\mathbf k _1, \mathbf k _{2})&= \frac{(\mathbf k _1 \cdot \mathbf k _2)^2}{k^2_1k^2_2} - 1, \end{aligned} $$(24)

F 2 ( k 1 , k 2 ) = 5 7 + 1 2 k 1 · k 2 k 1 k 2 ( k 2 k 1 + k 1 k 2 ) + 2 7 ( k 1 · k 2 ) 2 k 1 2 k 2 2 , Mathematical equation: $$ \begin{aligned} F_2(\mathbf k _1, \mathbf k _{2})&= \frac{5}{7} + \frac{1}{2} \frac{{\mathbf k _1 \cdot \mathbf k _2}}{k_1k_2}\ \Bigg (\frac{k_2}{k_1} + \frac{k_1}{k_2}\Bigg ) + \frac{2}{7}\frac{(\mathbf k _1 \cdot \mathbf k _2)^2}{k^2_1k^2_2}, \end{aligned} $$(25)

G 2 ( k 1 , k 2 ) = 3 7 + 1 2 k 1 · k 2 k 1 k 2 ( k 2 k 1 + k 1 k 2 ) + 4 7 ( k 1 · k 2 ) 2 k 1 2 k 2 2 . Mathematical equation: $$ \begin{aligned} G_2(\mathbf k _1, \mathbf k _{2})&= \frac{3}{7} + \frac{1}{2} \frac{\mathbf k _1 \cdot \mathbf k _2}{k_1k_2} \ \Bigg (\frac{k_2}{k_1} + \frac{k_1}{k_2}\Bigg ) + \frac{4}{7} \frac{(\mathbf k _1 \cdot \mathbf k _2)^2}{k^2_1k^2_2}. \end{aligned} $$(26)

The last term in the power spectrum model accounts for the fully degenerate contributions from the effective field theory (EFT) matter counter-terms depending on the effective speed of sound, cs (Baumann et al. 2012; Carrasco et al. 2012) and the higher derivative correction to linear bias,

P g ctr ( k ) = − 2 b 1 ( b 1 c s 2 + b ∇ 2 δ ) k 2 P L ( k ) Mathematical equation: $$ \begin{aligned} P_{\rm g}^{\,\mathrm {ctr}}(k)&=-2 \, b_1\left( b_1\,c_{\rm s}^{\,2}+b_{\nabla ^{\,2}\delta } \right) \, k^2 \, P_{\rm L}(k) \end{aligned} $$(27)

: = − 2 c 0 k 2 P L ( k ) , Mathematical equation: $$ \begin{aligned}&:= -2 \, c_0 \, k^2 \, P_{\rm L}(k), \end{aligned} $$(28)

with the parameter c0 := b1cs2 + b∇ 2δ. Finally, the expression for the galaxy bispectrum at tree-level in PT is simply

B g tree ( k 1 , k 2 , k 3 ) = 2 b 1 2 K 2 ( k 1 , k 2 ) P L ( k 1 ) P L ( k 2 ) + cyc . , Mathematical equation: $$ \begin{aligned} B^\mathrm{{tree}}_{\rm g}(k_1, k_2, k_3) = 2b_1^2 \, K_2(\mathbf k _1, \mathbf k _2)\,P_{\rm L}(k_1)\,P_{\rm L}(k_2) + \mathrm{cyc.} , \end{aligned} $$(29)

where, again, we neglect shot-noise contributions.

3.2. Infrared re-summation

The expressions described above do not account for the additional smearing of the baryonic features due to non-linear evolution (Eisenstein et al. 2007; Smith et al. 2007; Crocce & Scoccimarro 2008; Matsubara 2008; Desjacques et al. 2010; Baldauf et al. 2015; Senatore & Zaldarriaga 2015). We modeled this effect in the power spectrum following the approach of Blas et al. (2016a) in terms of a wiggle-no wiggle splitting of the linear power spectrum,

P L ( k ) = P nw ( k ) + P w ( k ) , Mathematical equation: $$ \begin{aligned} P_{\rm L}(k)=P_{\rm {nw}}(k)+P_{\rm w}(k), \end{aligned} $$(30)

obtained adopting the specific procedure of Vlah et al. (2016) (see Appendix C in Euclid Collaboration: Pezzotta et al. 2024, for a detailed description).

The linear galaxy power spectrum in real space is then replaced by a LO contribution,

P g IR , LO ( k ) = b 1 2 [ P nw ( k ) + e − k 2 Σ 2 P w ( k ) ] , Mathematical equation: $$ \begin{aligned} P_{\rm g}^\mathrm{{IR,\, LO}}(k)=b_1^2\left[P_{\rm {nw}}(k)+ \mathrm{e} ^{-k^2\Sigma ^2}P_{\rm w}(k)\right], \end{aligned} $$(31)

with the constant Σ2 representing the anisotropic variance of the relative displacement field (Eisenstein et al. 2007),

Σ 2 = 1 6 π 2 ∫ 0 k S d q P nw ( q ) [ 1 − j 0 ( q k osc ) + 2 j 2 ( q k osc ) ] , Mathematical equation: $$ \begin{aligned} \Sigma ^2 =&\frac{1}{6\pi ^2}\int _0^{k_S} {\mathrm{d} } q\,P_{\rm {nw}}(q)\,\left[1-j_0\left(\frac{q}{k_{\rm {osc}}}\right)+2j_2\left(\frac{q}{k_{\rm {osc}}}\right)\right], \end{aligned} $$(32)

where kosc = 1/ℓosc is the wavenumber associated with the BAOs scale ℓosc = 110 h−1 Mpc, while kS = 0.2 h Mpc−1 is an upper limit of integration whose specific choice does not lead to significant differences on the final result. To give an idea, for the redshift values considered in this work, z ∈ {0.9, 1.2, 1.5, 1.8}, we obtain Σ ∈ {8.72, 6.82, 5.22, 4.34}.

The one-loop prediction is now included in the NLO correction as

P g IR , LO + NLO ( k ) = b 1 2 P g , nw ( k ) + ( 1 + k 2 Σ 2 ) e − k 2 Σ 2 b 1 2 P w ( k ) + P g IR , 1 − loop ( k ) , Mathematical equation: $$ \begin{aligned} P_{\rm g}^\mathrm{{IR,\, LO+NLO}}(k) =&\ b_1^2\,P_{\rm {g},\,\mathrm {nw}}(k) +\left(1+k^2\Sigma ^2\right)\mathrm{e} ^{-k^2\Sigma ^2} \,b_1^2P_{\rm w}(k) \nonumber \\&+ P_{\rm g}^\mathrm{{IR,\, 1-loop}}(k) , \end{aligned} $$(33)

where the last term denotes the galaxy power spectrum one-loop correction evaluated in terms of the LO matter power spectrum, which is, schematically,

P g , nw IR , 1 − loop ( k ) : = P g 1 − loop [ P nw ( k ) + e − k 2 Σ 2 P w ( k ) ] . Mathematical equation: $$ \begin{aligned} P_{\rm {g},\, \mathrm{{nw}}}^\mathrm{{IR,\,1-loop}} (k) := P_{\rm g}^\mathrm{{1-loop}} \, [P_{\rm {nw}} (k) +\mathrm{e} ^{-k^2\Sigma ^2}P_{\rm {w}}(k)]. \end{aligned} $$(34)

The tree-level bispectrum expression is replaced instead by the LO prediction (Ivanov & Sibiryakov 2018)

B g tree , LO ( k 1 , k 2 , k 3 ) = 2 b 1 2 K 2 ( k 1 , k 2 ) × [ P nw ( k 1 ) P nw ( k 2 ) + e − k 1 2 Σ 2 P w ( k 1 ) P nw ( k 2 ) + e − k 2 2 Σ 2 P w ( k 2 ) P nw ( k 1 ) ] + 2 perm . Mathematical equation: $$ \begin{aligned} B_{\rm g}^\mathrm{{tree, LO}}(k_1, k_2, k_3) =&\ 2\, b_1^2\,K_2(\mathbf k _1, \mathbf k _2) \nonumber \\&\times \left[P_{\rm {nw}}(k_1)\,P_{\rm {nw}}(k_2) \right. +\mathrm{e} ^{-k_1^2\Sigma ^2}\,P_{\rm {w}}(k_1)\,P_{\rm {nw}}(k_2) \nonumber \\&\left.+\mathrm{e} ^{-k_2^2\Sigma ^2}\,P_{\rm w}(k_2)P_{\rm {nw}}(k_1)\right] + \mathrm{{2\ perm}}. \end{aligned} $$(35)

In what follows we retain, for clarity, the notation of Sect. 3.1 for all contributions to each correlation function, while we assume as implicit the implementation of infrared re-summation as described here.

3.3. Evaluation of 2PCF and 3PCF

The 2PCF was computed from the power spectrum of Eq. (18) according to the Hankel transform of Eq. (5), implemented using the FFTLog approach of Hamilton (2000). The decomposition in Legendre polynomials of the 3PCF defined in Eq. (10) allows us to apply an analogous procedure to this statistic. We can in fact consider the 2D extension of the FFTLog algorithm proposed by Fang et al. (2020) for an efficient evaluation of the integral in Eq. (12) over the Bℓ(k1, k2) functions, (see Umeh 2021; Guidi et al. 2023; Farina et al. 2026; Pugno et al. 2025, for further details).

However, the implementation of the 2DFFTLog approach still constitutes a computational burden when considered in the context of Monte Carlo sampling of a likelihood function over a large parameter space. For this reason, we developed an emulator for the 3PCF for a full set of cosmological and nuisance parameters. This allowed us to efficiently extend the full-shape analysis of the 2PCF to include the next, higher-order statistic. In fact, our emulator provides all contributions to the 2PCF and 3PCF where a combination of bias parameters can be factorised and describes the cosmological dependence of each term on the scalar amplitude parameter (As), the matter density (ωcdm), and the Hubble parameter (h). A detailed description of the emulator, along with validation tests, can be found in Appendix B.

4. Data

We validated our theoretical model against a set of synthetic galaxy catalogues obtained from the Euclid Flagship I N-body simulation. The catalogues adopt a halo occupation distribution (HOD) prescription to describe the population of Hα galaxies expected for the Euclid spectroscopic sample. We provide below a concise description of the galaxy catalogues, referring the reader to Euclid Collaboration: Pezzotta et al. (2024) for further details.

4.1. Euclid simulation

The Flagship I simulation uses the PKDGRAV3 code (Potter et al. 2017) to track the evolution of two trillion dark matter particles in a box of size L = 3780 h−1 Mpc, with a mass resolution, mp ∼ 2.398 × 109h−1 M⊙. It assumes a flat, ΛCDM cosmology with the fiducial parameters given in Table 11.

Table 1.

Fiducial parameters of the flat ΛCDM model.

We considered four comoving snapshots at redshift z = 0.9, 1.2, 1.5, and 1.8. In each snapshot, a friends-of-friends halo catalogue with a minimum mass corresponding to ten particles, is constructed. Galaxies are then assigned to these halos using a HOD prescription derived from the Flagship I light cone catalogue, designed to reproduce the Model 3 distribution of Pozzetti et al. (2016). Specifically, the HOD parameters are are chosen to match the expected selection function of the Euclid spectroscopic sample with a Hα flux limit of fHα > 2 × 10−16 erg cm−2 s−1, as outlined by Euclid Collaboration: Scaramella et al. (2022). Table 2 shows for each snapshot the total number of galaxies, the number density and a value for the linear bias obtained from measurements of the galaxy and matter power spectra. For a more detailed description of the mock galaxy catalogue we refer the reader, again, to Euclid Collaboration: Pezzotta et al. (2024).

Table 2.

Key properties of the synthetic galaxy catalogues.

We fitted our model to measurements from the full simulation volume, which is approximately three to six times larger than the effective volume of the Euclid redshift bins assumed for the forecasts analysis of Euclid Collaboration: Blanchard et al. (2020). The large simulation volume enables a stringent test of the model, which ensures that any residual theory systematic error is well within the expected precision of future measurements.

The catalogue we used does not take into account observational systematic effects such as target incompleteness, sample purity, and the impact of the angular footprint or radial selection function. A comprehensive exploration of such effects is left to other works (see Euclid Collaboration: Monaco et al., (in prep.); Euclid Collaboration: Risso et al. 2026).

4.2. Clustering measurements

We estimated the 2PCF using the natural estimator (Peebles 1973):

ξ ̂ ( r ) : = N DD ( r ) N RR ( r ) − 1 , Mathematical equation: $$ \begin{aligned} \hat{\xi }(r) \, := \, \frac{N_{\rm {DD}}(r)}{N_{\rm {RR}}(r)} - 1, \end{aligned} $$(36)

where NDD(r) and NRR(r) represent the pair counts, as a function of separation r, in the data and a random distribution with, in this ideal case, constant density. This is equivalent to the usual Landy–Szalay estimator (Landy & Szalay 1993) in the case of a cubic box with periodic boundary conditions (they should both lead to the same variance). We took advantage of this property to also estimate analytically the number of pairs from the random catalogue as

N RR ( r ) = 4 π 3 n ¯ 2 V [ ( r + Δ r 2 ) 3 − ( r − Δ r 2 ) 3 ] , Mathematical equation: $$ \begin{aligned} N_{\rm {RR}}(r) \, = \, \frac{4\pi }{3}\bar{n}^2 V \left[ \left(r + \frac{\Delta r}{2}\right)^3 - \left(r - \frac{\Delta r}{2}\right)^3 \right], \end{aligned} $$(37)

where we assume a separation bin of size Δr centred on r, n ¯ Mathematical equation: $ \bar{n} $ being the mean galaxy number density and V the volume of the box.

We measured the 2PCF in the four snapshots for separations ranging from rmin = 0 h−1 Mpc to rmax = 150 h−1 Mpc. We used linearly spaced bins with a width of Δr = 5 h−1 Mpc. The results are displayed in the top panels of Fig. 1, where each column represents one of the four redshift snapshots. In the analysis presented in Sects. 6 and 7, we restricted the 2PCF data vector to separations r ≤ rmax = 140 h−1 Mpc.

Thumbnail: Fig. 1. Refer to the following caption and surrounding text. Fig. 1.

Top panels: 2PCF measurements for four comoving snapshots from the Flagship I simulation of the Model 3 HOD. Bottom panels: Same but for the 3PCF measurements. The triangular configurations are ordered by increasing values of r12, r13, and r23 under the condition r12 ≤ r13 ≤ r23, as depicted in Fig. A.1. The circles are colour-coded according to the η, while the vertical grey lines mark a group of triangles that share the same value of the smaller side (r12).

The Landy–Szalay estimator for the 2PCF can be easily extended to the analogous estimator for the 3PCF, constructed in terms of triplet counts (see e.g. Szapudi & Szalay 1998). However, since such an estimator scales with the number of galaxies Ng as 𝒪(Ng3), its computational cost quickly becomes prohibitive when it comes to very large datasets such as those we considered here, where Ng is around 108. To resolve this issue, Slepian & Eisenstein (2015) introduced an estimator based on the spherical harmonics decomposition (SHD), which improves significantly the computational efficiency. In this approach, the dependence of the 3PCF on the angle between the triangle sides r12 and r13 is expanded in terms of Legendre polynomials according to Eq. (10). As a result, the full estimator ζ ̂ ( r 12 , r 13 , r 23 ) Mathematical equation: $ \hat{\zeta}(r_{12}, r_{13}, r_{23}) $ is written as a function of an estimator ζ ̂ ℓ ( r 12 , r 13 ) Mathematical equation: $ \hat{\zeta}_{\ell} (r_{12}, r_{13}) $ of the expansion coefficients as

ζ ̂ ( r 12 , r 13 , r 23 ) = ∑ ℓ = 0 ℓ max ζ ̂ ℓ ( r 12 , r 13 ) L ( μ 12 , 13 ) , Mathematical equation: $$ \begin{aligned} \hat{\zeta }(r_{12}, r_{13}, r_{23}) = \sum _{\ell =0}^{\ell _{\rm {max}}} \hat{\zeta }_{\rm \ell } (r_{12}, r_{13}) \, \mathcal{L} (\mu _{12,13}) , \end{aligned} $$(38)

thereby reducing the computational complexity to 𝒪(Ng2). The coefficients are estimated as

ζ ̂ ℓ ( r 12 , r 13 ) = N DDD , ℓ − 3 N DDR , ℓ + 3 N DRR , ℓ − N RRR , ℓ N RRR , 0 , Mathematical equation: $$ \begin{aligned} \hat{\zeta }_\ell (r_{12}, r_{13}) \, = \, \frac{N_{\rm {DDD},\,\ell }-3N_{\rm {DDR},\,\ell }+3N_{\rm {DRR},\,\ell }-N_{\rm {RRR},\,\ell }}{N_{\rm {RRR},\, 0}}, \end{aligned} $$(39)

where NDDD,  ℓ, NDDR,  ℓ, NDRR,  ℓ, and NDRR,  ℓ are, in turn, the coefficients of the Legendre expansion of the triplets counts from the galaxy catalogue, from the random and the two mixed terms, directly evaluated in the approach of Slepian & Eisenstein (2015). Again, since we are working with a simulation box with periodic boundary conditions, it is possible to simplify the evaluation of the coefficients ζ ̂ ℓ ( r 12 , r 13 ) Mathematical equation: $ \hat{\zeta}_{\ell} (r_{12}, r_{13}) $ and derive an analytical expression for the monopole of the random triplet counts NRRR,  0.

The SHD estimator is known to be less efficient for isosceles or nearly isosceles triangle configurations (i.e. with r12 ≃ r13). In these cases, a large value of ℓmax is required to recover results consistent with those of the standard, direct-counting, estimator. Indeed, for isosceles triangles each coefficient of the expansion in Eq. (10) is the result of an integration over the bispectrum down to small, non-linear scales, beyond the range of validity of the perturbative model. On the other hand, using a fixed ℓmax ensures internal consistency between measurements and theory. To mitigate this issue, we introduced the quantity (see e.g. Veropalumbo et al. 2022)

η : = r 13 − r 12 Δ r , Mathematical equation: $$ \begin{aligned} \eta := \frac{r_{13} - r_{12}}{\Delta r}, \end{aligned} $$(40)

with Δr being the separation bin size, to parametrise the proximity of a given configuration to the isosceles case. Fixing a lower limit (ηmin) amounts to excluding isosceles, or nearly isosceles configurations.

In the bottom panels of Fig. 1 we present the 3PCF measurements for the four snapshots. The triangle index is defined by increasing value of r12, r13, and r23 under the condition r12 ≤ r13 ≤ r23, as depicted in Fig. A.1. The measurements include all configurations with side lengths from rmin = 0 h−1 Mpc to rmax = 150 h−1 Mpc in bins of width Δr = 5 h−1 Mpc. We expanded the 3PCF up to ℓmax = 10, which strikes an ideal balance between computational efficiency and information content. A more detailed description of the numerical implementation and validation of the SHD estimator for the 3PCF will be given in Euclid Collaboration: Veropalumbo et al. (in prep.), while a similar presentation of the 2PCF estimator is given in Euclid Collaboration: de la Torre et al. (2025).

4.3. Covariance

We made use of a covariance matrix obtained from an analytical model for both two- and three-point functions. We worked in the Gaussian approximation, so all covariance depends simply on the galaxy power spectrum, and we further assumed a linear model for this spectrum.

For the 2PCF covariance we have (Grieb et al. 2016; Lippich et al. 2019)

C ξ ( r , r ′ ) = 1 π 2 V ∫ 0 ∞ d k k 2 P tot 2 ( k ) j ¯ 0 ( k r ) j ¯ 0 ( k r ′ ) , Mathematical equation: $$ \begin{aligned} \mathcal{C} _{\xi }(r, r\prime ) = \frac{1}{\pi ^2 V} \, \int _0^{\infty } {\mathrm{d} } k \, k^2 \, P^2_{\rm {tot}}(k) \, \bar{j}_{0}(kr) \, \bar{j}_{0}(kr\prime ), \end{aligned} $$(41)

where V is the box volume, j ¯ 0 Mathematical equation: $ \bar{j}_0 $ are the bin-averaged, zeroth-order Bessel functions (see Eq. A19 in Grieb et al. 2016), while

P tot ( k ) : = b 1 2 P lin ( k ) + 1 n ¯ g Mathematical equation: $$ \begin{aligned} P_{\rm {tot}}(k) := b_1^2 \, P_\mathrm{lin} (k) + \frac{1}{\bar{n}_g}\, \end{aligned} $$(42)

is the linear galaxy power spectrum including a Poisson shot-noise contribution, n ¯ g Mathematical equation: $ \bar{n}_g $ being the galaxy number density. For each redshift snapshot, we computed the covariance matrix using the bias and density values listed in Table 2. As noted by Smith (2009), a more detailed treatment of discrete tracers introduces additional terms in the covariance – particularly relevant at low number density – even under Gaussian assumptions. These include contributions from combinatoric sampling effects, which are not included in our model; this is sufficient given the number density of the simulated dataset.

For the 3PCF covariance, we followed the expression proposed in Slepian & Eisenstein (2015), rigorously validated in Veropalumbo et al. (2022). This again takes advantage of the SHD writing the covariance between two measurements ζ ̂ ( r 12 , r 13 , r 23 ) Mathematical equation: $ \hat{\zeta}(r_{12},r_{13},r_{23}) $ and ζ ̂ ( r 12 ′ , r 13 ′ , r 23 ′ ) Mathematical equation: $ \hat{\zeta}(r_{12}{\prime},r_{13}{\prime},r_{23}{\prime}) $ as

C ζ ( r 12 , r 13 , r 23 ; r ′ 12 , r ′ 13 , r ′ 23 ) = = ∑ ℓ , ℓ ′ ℓ max C ζ , ℓ ℓ ′ ( r 12 , r 13 ; r ′ 12 , r ′ 13 ) L ℓ ( μ 12 , 13 ) L ℓ ′ ( μ 12 , 13 ) , Mathematical equation: $$ \begin{aligned} \begin{aligned} \mathcal{C} _{{\zeta }}(&r_{12}, r_{13}, r_{23};r\prime _{12}, r\prime _{13}, r\prime _{23}) = \\&= \sum _{\ell , \ell \prime }^{\ell _{max}} \mathcal{C} _{{\zeta }, \ell \ell {\prime }}(r_{12}, r_{13}; r\prime _{12}, r\prime _{13}) \, \mathcal{L} _{\ell }(\mu _{12,13}) \, \mathcal{L} _{\ell \prime }(\mu _{12,13}), \end{aligned} \end{aligned} $$(43)

where 𝒞ζ, ℓℓ′(r12, r13; r′12, r′13) represents the covariance between the coefficients ζ ̂ ℓ ( r 12 , r 13 ) Mathematical equation: $ \hat{\zeta}_\ell(r_{12},r_{13}) $ and ζ ̂ ℓ ′ ( r 12 r 12 , r 13 r 12 ) Mathematical equation: $ \hat{\zeta}_{\ell{\prime}}(r_{12}r_{12},r_{13}r_{12}) $. This, in turn, is given by

C ζ , ℓ ℓ ′ ( r i , r j ; r ′ i , r ′ j ) = = ( 2 ℓ + 1 ) ( 2 ℓ ′ + 1 ) ( − 1 ) ℓ + ℓ ′ ( 2 π ) 6 V ∫ d p 1 d p 2 d p 3 δ D ( p 123 ) × P tot ( p 1 ) P tot ( p 2 ) P tot ( p 3 ) j ¯ ℓ ( p 1 r i ) j ¯ ℓ ( p 2 r j ) L ℓ ( μ 12 ) × [ j ¯ ℓ ′ ( p 1 r i ′ ) j ¯ ℓ ′ ( p 2 r j ′ ) L ℓ ′ ( μ 12 ) + 5 perm . ] . Mathematical equation: $$ \begin{aligned} \mathcal{C} _{\zeta ,\ell \ell \prime }&(r_i, r_j;r\prime _i, r\prime _j )= \\ \nonumber =&\frac{(2\ell +1)\, (2\ell \prime +1)\, (-1)^{\ell +\ell \prime }}{(2\pi )^6\,V}\!\int \!{\mathrm{d} } p_1 \, {\mathrm{d} } p_2 \, {\mathrm{d} } p_3 \, \delta _{\rm {D}}(\mathbf p _{123}) \\ \nonumber&\times P_{\rm {tot}}(p_1) \, P_{\rm {tot}}(p_2) \, P_{\rm {tot}}(p_3) \bar{j}_\ell (p_1r_i) \, \bar{j}_\ell (p_2r_j) \, {\mathcal{L} }_\ell (\mu _{12}) \\ \nonumber&\times \left[\bar{j}_{\ell \prime }(p_1r_i\prime ) \, \bar{j}_{\ell \prime }(p_2r_j\prime ) \, {\mathcal{L} }_{\ell \prime }(\mu _{12})+5\ \mathrm{{perm}.}\right]. \end{aligned} $$(44)

We refer the reader to Slepian & Eisenstein (2015) for a detailed description of the practical evaluation of this expression.

In the Gaussian approximation, the cross-covariance between the two- and three-point functions vanishes, as the only contributions arise from the bispectrum of the density field and the connected 5-point function (see e.g. Sefusatti et al. 2006). While these terms are expected to introduce non-negligible correlations – potentially at the ∼20% level for amplitude-related parameters – we neglected them in this work for consistency with the Gaussian covariance framework. This choice is justified by the scope of our analysis, which focuses on validating the modelling pipeline rather than delivering optimal cosmological constraints.

5. Likelihood analysis and performance metrics

We define in the section the likelihood function we adopted to fit the parameters of our model to the simulation measurements. We also introduce the statistical tools we used to determine the range of validity of the model in terms of the GoF, the level of theoretical uncertainties introduced by the model, and the ability to constrain cosmological parameters.

5.1. Likelihood function

We adopted a Gaussian likelihood function for both the 2PCF and the 3PCF measurements. We thus have

− 2 ln L ξ ( θ ) = χ ξ 2 ( θ ) , Mathematical equation: $$ \begin{aligned} -2\ln \mathcal{L} _\xi (\mathbf {\theta} ) = \chi _\xi ^2 (\mathbf {\theta} ), \end{aligned} $$(45)

where vector θ denotes the model parameters while the χ2 for the 2PCF is given by

χ ξ 2 ( θ ) = ∑ i , j = 1 N r [ ξ i ( θ ) − ξ ̂ i ] C ξ , i j − 1 [ ξ j ( θ ) − ξ ̂ j ] , Mathematical equation: $$ \begin{aligned} \chi ^2_{\rm \xi }(\mathbf {\theta} ) = \sum _{i,j=1}^{N_{r}} \big [ \xi _i(\mathbf {\theta} ) - \hat{\xi }_i \big ] \ C^{-1}_{\xi , ij} \ \big [ \xi _j(\mathbf {\theta} ) - \hat{\xi }_j \big ], \end{aligned} $$(46)

with ξ(θ) and ξ ̂ Mathematical equation: $ \hat{\xi} $ representing respectively the 2PCF theoretical prediction and its measurements and with the sum running over the separation bins ri and rj, Nr being the total number of bins, corresponding to the size of the data vector. Similarly, for the 3PCF we have

− 2 ln L ζ ( θ ) = χ ζ 2 ( θ ) , Mathematical equation: $$ \begin{aligned} -2\ln \mathcal{L} _\zeta (\mathbf {\theta} ) = \chi _\zeta ^2 (\mathbf {\theta} ), \end{aligned} $$(47)

with

χ ζ 2 ( θ ) = ∑ i , j = 1 N t [ ζ i ( θ ) − ζ ̂ i ] C ζ , i j − 1 [ ζ j ( θ ) − ζ ̂ j ] , Mathematical equation: $$ \begin{aligned} \chi ^2_{\zeta }(\mathbf {\theta} ) = \sum _{i,j=1}^{N_{\rm t}} \big [ \zeta _i(\mathbf {\theta} ) - \hat{\zeta }_i \big ] \ C^{-1}_{\zeta , ij} \ \big [ \zeta _j(\mathbf {\theta} ) - \hat{\zeta }_j \big ], \end{aligned} $$(48)

with the sum now extending over all triangular configuration up to their total number (Nt).

The likelihood of joint measurements of the 2PCF and 3PCF is simply given by

L ξ + ζ ( θ ) = L ξ ( θ ) L ζ ( θ ) , Mathematical equation: $$ \begin{aligned} \mathcal{L} _{\xi +\zeta }(\mathbf {\theta} ) = \mathcal{L} _{\xi }(\mathbf {\theta} ) \, \mathcal{L} _{\zeta }(\mathbf {\theta} ), \end{aligned} $$(49)

and, correspondingly, the joint χ2 will be

χ ξ + ζ 2 ( θ ) = χ ξ 2 ( θ ) + χ ζ 2 ( θ ) Mathematical equation: $$ \begin{aligned} \chi ^2_{\rm \xi +\mathrm {\zeta} }(\mathbf {\theta} ) = \chi ^2_{\xi }(\mathbf {\theta} ) + \chi ^2_{\zeta }(\mathbf {\theta} ) \, \end{aligned} $$(50)

since, as mentioned above, we neglected the cross-covariance between the two statistics.

We sampled the posterior using the emcee code2 (Foreman-Mackey et al. 2013), implementing an affine-invariant ensemble sampler, which is particularly well suited to high-dimensional parameter spaces. We employed 50 walkers in each run and ensured that the chains are terminated only after exceeding 100 integrated autocorrelation times, thereby ensuring full convergence.

The model parameters, along with uniform priors assumed in all likelihood evaluations, are outlined in Table 3. The choice of large intervals is meant to ensure uninformative priors.

Table 3.

Cosmological and nuisance parameters of the 2PCF and 3PCF models.

5.2. Performance metrics

Our main goal is the assessment of the range of validity of the model, determining the minimum scale, rmin, down to which the model can still accurately describe the measurements without introducing significant systematic uncertainties to the recovered cosmological parameters. As in the companion paper (Euclid Collaboration: Pezzotta et al. 2024), this task is carried out using three performance metrics: a GoF, a figure of bias (FoB), and figure of merit (FoM) as defined below.

5.2.1. Goodness of fit

For an assessment on the GoF, we simply used the standard χ2 test. The χ2 definitions for the 2PCF, 3PCF, and joint at a given point θ in the parameter space are given in Eqs. (46), (48), and (50). The posterior-averaged value, ⟨χ2⟩post, is then compared to the thresholds from the χ2 distribution at the 68% and 95% confidence levels, determined in terms of the appropriate number of degrees of freedom, computed as the difference between the size of the data vector and the number of model parameters.

5.2.2. Figure of bias

We quantified the uncertainties in recovering a subset of fiducial, cosmological parameters (θfid) due to model systematic uncertainties in terms of an FoB defined as

FoB ( θ ) = [ ( ⟨ θ ⟩ post − θ fid ) T S − 1 ( θ ) ( ⟨ θ ⟩ post − θ fid ) ] 1 / 2 , Mathematical equation: $$ \begin{aligned} \mathrm{{FoB}}(\mathbf {\theta} ) = \left[ \left( \langle \mathbf {\theta} \rangle _{\rm {post}} - \mathbf {\theta} _{\rm {fid}} \right)^\mathrm{T} S^{-1}(\mathbf {\theta} ) \left( \langle \mathbf {\theta} \rangle _{\rm {post}} - \mathbf {\theta} _{\rm {fid}} \right) \right]^{1/2}, \end{aligned} $$(51)

where ⟨θ⟩post represents the parameter averaged over the posterior distribution, while S(θ) is the parameters’ covariance matrix, evaluated from the same posterior distribution.

When θ consists of a single parameter, the FoB measures the deviation of the posterior mean from the fiducial value in units of the parameter’s standard deviation, with FoB values of one and two corresponding to the 68% and 95% confidence levels, respectively. However, when dealing with multiple parameters, these thresholds must be recalculated by integrating a multivariate normal distribution over the appropriate number of dimensions. For instance, in the case of three parameters relevant for the results of Sect. 7, the 68% and 95% confidence levels correspond to FoB values of 1.87 and 2.80, respectively.

5.2.3. Figure of merit

Finally, we defined an FoM (Albrecht et al. 2006; Wang 2008) to quantify the constraining power of our model over a specific range of scales and under given bias assumptions:

FoM ( θ ) = 1 det [ S ( θ ) ] , Mathematical equation: $$ \begin{aligned} \mathrm{{FoM}}(\mathbf {\theta} ) = \frac{1}{\sqrt{\det [S(\mathbf {\theta} )]}} , \end{aligned} $$(52)

where det[S(θ)] is the determinant of the covariance matrix of the parameters of interest. This represents the hyper-volume enclosed by the hyper-surface defined by the covariance matrix S(θ). Therefore, a larger FoM indicates tighter constraints on the model parameters.

6. Galaxy bias

In this section we assess the range of validity of our bias model for both the 2PCF and 3PCF. Our focus will be then on the GoF test, keeping fixed the cosmological parameters at their fiducial values. We considered both the case where all bias parameters are free and the case where the bias relations mentioned above is assumed to reduce the parameter space. We also needed to make sure that no significant running of the bias parameters was present, i.e. there was no significant variation in the recovered parameters, as a function of the minimal scale defining the data vectors.

6.1. Maximal model

We started with the maximal model (see e.g. Oddo et al. 2021), with all bias parameters free, with the priors given in Table 3. In Fig. 2 we show the χ2 averaged over the posterior distribution for the four redshift snapshots and for the 2PCF, 3PCF, and joint 2PCF+3PCF data vector.

Thumbnail: Fig. 2. Refer to the following caption and surrounding text. Fig. 2.

GoF test for the bias 2PCF and 3PCF models with fixed cosmological parameters as a function of the minimal scale, r min 2 PCF Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{2PCF}} $ and r min 3 PCF Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{3PCF}} $, respectively, assumed for the analysis. The posterior-averaged χ2 is compared to the 68% and 95% confidence intervals for the χ2 distribution, denoted by the grey shaded areas, given the corresponding Nd.o.f. (fixing the case ηmin = 1 for the 3PCF confidence levels; see the main text for more details). For the 3PCF, we considered three different values for the ηmin parameter, Eq. (40), progressively excluding near-isosceles configurations. In the joint analysis, the minimal scale for the 2PCF was kept fixed at r min 2 PCF = 25 h − 1 Mpc Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{2PCF}}=25\, h^{-1} \, \mathrm{{Mpc}} $, with the results shown as a function of r min 3 PCF Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{3PCF}} $ alone. The four columns span the four simulation snapshots at redshifts (z) of 0.9, 1.2, 1.5, and 1.8.

We notice in the first place that while the one-loop 2PCF model provides a good fit for all scales considered, the tree-level 3PCF prediction fails for separations below 20 h−1 Mpc at z = 0.9, with a better performance at higher redshift. Under the choice of ηmin = 1, Eq. (40), including configurations closer to the isosceles case, the validity of the model is further restricted to rmin > 20 h−1 Mpc at z = 0.9 and rmin > 10 h−1 Mpc for z = 1.2 and above.

In Fig. 3 we show the marginalised 2D constraints on the nuisance parameters (bias and counter-term) from the 2PCF analysis against the joint 2PCF plus 3PCF analysis for several values of r min 3 PCF Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,3PCF}} $. This shows in the first place the peculiar and significant constraining power of higher-order statistics on linear and non-linear bias. This is particularly evident in the reduction of the posterior uncertainty on the non-linear parameter b2 and non-local parameters b𝒢2 and bΓ3. It further shows that down to r min 3 PCF = 30 h − 1 Mpc Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,3PCF}}=30\, h^{-1} \, \mathrm{{Mpc}} $ no running of the bias parameter is evident, with all analyses at different rmin consistently reducing the posterior uncertainty.

Thumbnail: Fig. 3. Refer to the following caption and surrounding text. Fig. 3.

Left panel: Marginalised 2D constraints at confidence intervals of 68.3% and 99.5%, on the nuisance parameters (bias and counter-term) from the 2PCF analysis (grey shaded areas) against the joint 2PCF and 3PCF analysis for several values of r min 3 PCF Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,3PCF}} $ at z = 1.5 (dashed blue lines, dashed orange lines, and green shaded areas). All cases assume η = 3; see Eq. (40). The vertical grey band is an estimate of the linear bias b1 obtained from a comparison of the galaxy power spectrum with measurements of the matter power spectrum in Euclid Collaboration: Pezzotta et al. (2024). Right panel: Marginalised 2D constraints at confidence intervals of 68.3% and 99.5%, on the nuisance parameters (bias and counter-term) from the joint 2PCF-3PCF analysis at z = 1.5 for r min 3 PCF = 30 h − 1 Mpc Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,3PCF}}=30\, h^{-1} \, \mathrm{{Mpc}} $ (red contours) compared to the same set-up but with the assumption of the b𝒢2(b1) bias relation alone (red contours) and the combination of the relations b𝒢2(b1) and bΓ3(b1) (violet contours).

The plots include, as a vertical grey band, an estimate of the value of the linear bias b1 obtained from a comparison of the galaxy power spectrum with measurements of the matter power spectrum in Euclid Collaboration: Pezzotta et al. (2024). We find that the inclusion of the 3PCF consistently yields posterior distributions for b1 that are visually closer to the reference value across all the minimum scale cuts explored. While the 2PCF-only results remain statistically consistent, they exhibit broader posteriors, reflecting a stronger degeneracy with other model parameters. We revisit this aspect in the context of cosmological parameter constraints in the following section.

6.2. Bias relations

The right panel of Fig. 3 shows the same results for the joint analysis with rmin = 30 h−1 Mpc and ηmin = 3 (green contours), here compared to the same set-up but with the additional assumptions of two bias relations. We considered, in particular, the b𝒢2(b1) relation alone, Eq. (16), (red contours) and its combination with the relation for bΓ3 of Eq. (17). We display only the results for the z = 1.5 snapshots, noticing that for the measurements at other redshifts we obtain qualitatively similar results.

We notice that in both cases, the bias relations provide a reduction of the parameter space by one and two parameters, respectively, without affecting the b1, b2 constraints. On this specific contour in fact we do not observe any induced systematic shift as we do not notice any appreciable reduction in the uncertainty. This is noticeable instead on the c0 parameter, describing the combined effect of the EFT counter-term and higher derivative bias. These results provide some evidence that the bias relations can be used to speed up the likelihood evaluation without compromising the main results. We return to this point when discussing cosmological parameters constraints.

We close the section by noticing that the adoption of one or two bias relations does not lead to any worsening of the GoF, in terms of the χ2 test. We avoid a dedicated figure since the difference with the maximal model case would be barely visible.

7. Cosmological parameters from the joint, full-shape analysis

In this section we present the results for the full-shape, joint analysis of 2PCF and 3PCF aiming at constraining three cosmological parameters: the amplitude of scalar fluctuations (As), the CDM density (ωcdm), and the Hubble constant (h). This minimal set captures the main physical effects to which large-scale structure observables are directly sensitive: As determines the amplitude of fluctuations and is probed by the power spectrum and bispectrum; ωcdm governs the growth history and the shape of the clustering signal; and h fixes the physical scales through the matter–radiation equality scale and the sound horizon, affecting the shape of correlation functions in real space. Limiting our analysis to these parameters highlights the constraining power of the full-shape information.

To the best of our knowledge, this is the first attempt at estimating the improvement due to adding the 3PCF to the full-shape analysis of the 2PCF, although still limited to real space. Until now, 3PCF analyses have been restricted to template fitting, i.e. using fixed-shape models to extract cosmological information from the BAO scale position or from the anisotropy induced by redshift-space distortions, rather than performing full-shape parameter inference, due to the considerable computational demands of a full cosmology-dependent 3PCF model.

A key ingredient, to achieve this goal, is the emulator described in Appendix B. Constructed using a PyTorch3 framework, the emulator has undergone extensive testing and has been fully integrated into the same sampler. The emulator extends to the 2PCF prediction, further enhancing the overall efficiency of our analysis.

Figure 4 shows marginalised constraints on cosmological parameters As, ωcdm, h, and on the linear bias b1 derived from the joint 2PCF and 3PCF analysis as a function of the minimal scale rmin that characterises the 3PCF data vector, with the minimal separation in the 2PCF measurements fixed at r min 2 PCF = 25 h − 1 Mpc Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,2PCF}}=25\, h^{-1} \, \mathrm{{Mpc}} $. The four columns span the different redshifts. As in Fig. 2, varying shades of these colours correspond to different choices for the parameter ηmin. The dashed black lines denote the fiducial values assumed for the simulations, together with the grey line that shows the same fiducial values for b1 considered in Fig. 3. Overall, the results show good agreement between the emulated predictions and the expected values. For large values of rmin, the estimates for As and b1 tend to deviate further from the expected values, largely due to projection effects caused by a large degeneracy between these parameters. On the other hand, moving to the small-scale regime, a possible failure of the model to recover the fiducial values is only evident for ωcdm at the lowest redshift of z = 0.9. Below rmin = 40 h−1 Mpc, the estimates for As and b1(z) exhibit a slight dependence on the choice of ηmin, with the case of ηmin = 1 leading to more biased values for these parameters, particularly at low redshift.

Thumbnail: Fig. 4. Refer to the following caption and surrounding text. Fig. 4.

Marginalised constraints on the cosmological parameters and linear bias from a joint 2PCF and 3PCF analysis. The posterior mean along with its uncertainties is displayed as a function of the minimal scale ( r min 3 PCF Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,3PCF}} $) for different choices of ηmin and fixing r min 2 PCF = 25 h − 1 Mpc Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,2PCF}}=25\, h^{-1} \, \mathrm{{Mpc}} $ and rmax = 140 h−1 Mpc for both statistics. The dashed black line refers to the fiducial values assumed in the simulation, while the dashed grey line refers to the fiducial b1 adopted in Euclid Collaboration: Pezzotta et al. (2024), also shown in Fig. 3.

In Fig. 5 we present the performance metrics corresponding to the results of Fig. 4, all as a function of the minimal scales r min 3 PCF Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,3PCF}} $, with r min 2 PCF Mathematical equation: $ r^{\mathrm{2PCF}}_{\mathrm{min}} $ fixed at 25 h−1 Mpc, and for the same three different values of ηmin. The top panels show the posterior-averaged ⟨χ2⟩post. The grey bands denote the 68% and 99.7% confidence levels for the χ2 distribution. We find a good fit for r min 3 PCF ≳ 20 h − 1 Mpc Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,3PCF}}\gtrsim 20\, h^{-1} \, \mathrm{{Mpc}} $ except for the lowest redshift where r min 3 PCF ≳ 30 h − 1 Mpc Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,3PCF}}\gtrsim 30\, h^{-1} \, \mathrm{{Mpc}} $ ensures a safer result. This is naturally consistent with the results at fixed cosmology. The middle panels show the FoB, defined in terms of the parameters As, ωcdm, and h. The uncertainties are shown again in two shades of grey, corresponding to 68% and 95% confidence levels. The model provides overall unbiased results except for the first separation bin at rmin = 10 h−1 Mpc, for all redshifts and for all values of ηmin. Finally, the bottom panels of Fig. 5 show the FoM, defined again in terms of the three cosmological parameters. The FoM increases significantly across the whole separation range (notice the log-scale on the y-axis) and presents also larger values in the lowest redshift snapshots. In addition, the choice of a low value for ηmin can also make a large difference, since this increase the size of the data vector and relative S/N.

Thumbnail: Fig. 5. Refer to the following caption and surrounding text. Fig. 5.

Performance metrics for the full-shape, joint analysis of the 2PCF and the 3PCF. All quantities are shown as a function of the minimal scale ( r min 3 PCF Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,3PCF}} $), with r min 2 PCF Mathematical equation: $ r^{\mathrm{2PCF}}_{\mathrm{min}} $ fixed at 25 h−1 Mpc, and for the three different values of ηmin in different colour shades. From top to bottom we show the posterior-averaged ⟨χ2⟩post for the GoF, the FoB, and the FoM, these last two both defined in terms of the three cosmological parameters As, ωcdm, and h. The grey bands denote the 68% and 99.7% confidence levels for the χ2 distribution and the FoB.

The left plot of Fig. 6 shows the constraints on the cosmological parameters at z = 0.9 plus the linear bias obtained from the combined statistics for different choices of r min 3 PCF Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,3PCF}} $ and ηmin, compared to the 2PCF-only constraints where we include an additional, informative Gaussian prior on b1 for illustrative purposes defined by a 10% error on the best-fit value obtained from the joint analysis. This is required by the strong degeneracy between As and b1 in the 2PCF likelihood. This degeneracy is still present in the joint analysis, although significantly reduced without imposing any prior on b1. We obtain a discrepancy with the fiducial value of As at confidence interval of 68.3% for the configuration with r min 3 PCF = 30 h − 1 Mpc Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,3PCF}}=30\, h^{-1} \, \mathrm{{Mpc}} $ and ηmin = 2, while the discrepancy is smaller for the other parameters, despite the additional degeneracy between h and ωcdm. Interestingly, no tension is visible on the h and ωcdm constraints.

Thumbnail: Fig. 6. Refer to the following caption and surrounding text. Fig. 6.

Left panel: Marginalised 2D and 1D posterior distributions at z = 0.9 from the joint 2PCF and 3PCF analysis of the cosmological parameters h, ωcdm, and As and the linear bias parameter b1 for various 3PCF scale cuts ( r min 3 PCF Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,3PCF}} $) and the ηmin parameter defined in Eq. (40), fixing r min 2 PCF = 25 h − 1 Mpc Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,2PCF}} = 25 \, h^{-1} \, \mathrm{{Mpc}} $. Grey contours correspond to the 2PCF-only analysis assuming a Gaussian prior on b1 given by the 10% uncertainties around the best-fit value obtained from the joint analysis. Each panel shows the results for a different snapshot. Vertical and horizontal lines denote the fiducial values of the cosmological parameters. Right panel: Marginalised 1D and 2D posterior distributions of the same subset of model parameters for the combined analysis of 2PCF and 3PCF at z = 0.9 fixing r min 2 PCF = 25 h − 1 Mpc Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,2PCF}} = 25 \, h^{-1} \, \mathrm{{Mpc}} $. The reference case of the maximal model with r min 3 PCF = 30 h − 1 Mpc Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,3PCF}} = 30 \, h^{-1} \, \mathrm{{Mpc}} $ and ηmin = 3 is compared to the results obtained adopting the bias relations of Eq. (16), in red, and its combination with Eq. (17), in brown.

The same results are found with no relevant difference for all snapshots, as described in Fig. C.1. We remind the reader that the posterior distributions are not statistically independent since each realisation is a different snapshot of the same evolving dark matter perturbations, sharing the same initial seeds. We observe that at lower redshifts the constraints are tighter due to the larger density of the galaxy distribution and that, in general, both h and ωcdm exhibit accurate estimates of the fiducial values. The largest discrepancy, at the limit of the confidence interval of 99.5% contour, is obtained for As at z = 1.2. We notice that now also the linear bias parameter is correctly recovered at all redshifts.

In the right plot of Fig. 6, we compare the maximal model with all free parameters, assuming r min 3 PCF = 30 h − 1 Mpc Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,3PCF}} = 30 \, h^{-1} \, \mathrm{{Mpc}} $ and ηmin = 3, to the results obtained assuming the bias relation of Eq. (16) for b G 2 Mathematical equation: $ b_{\mathcal{G}_2} $, and its combination with the relation for bΓ3 of Eq. (17), showing the same set of parameters as in the right plot. Overall, we notice how the adoption of the bias relations enhances the precision of the cosmological parameter constraints. Notably, the 2D posterior distribution of As and b1, affected by a strong degeneracy, benefits the most from this additional information, although, again, we notice a slightly biased estimate for As. As presented in Fig. C.2, which shows the full set of parameters, the bΓ3 relation greatly reduces the uncertainty on the higher derivatives and/or counter-term parameter c0, as in the bias-only case.

8. Conclusions

In this work we provide a validation of PT models for the 2PCF and 3PCF in real space across a redshift range and against a synthetic galaxy catalogue representative of the Euclid spectroscopic galaxy survey (see e.g. Euclid Collaboration: Mellier et al. 2025).

We took advantage of galaxy clustering measurements obtained from four redshift snapshots of the Flagship I simulation populated with the HOD following the Model 3 prescription of Pozzetti et al. (2016). In particular, we determined HOD parameters by selecting galaxies with a Hα flux limit of fHα = 2 × 10−16 erg cm−2 s−1. The 3PCF measurements in particular, presented in Sect. 4, are based on the SHD introduced by Slepian & Eisenstein (2015); see also Euclid Collaboration: Veropalumbo et al. (in prep). This estimator allowed us to measure all 3PCF configurations while keeping a manageable computational cost.

The galaxy 2PCF is modelled at the NLO within the framework of EFTofLSS, involving up to five nuisance parameters: b1, b2, b𝒢2, bΓ3, and the parameter c0, which describes an EFT counter-term and higher-derivative bias. The 3PCF, on the other hand, was modelled at the LO in PT, with the tree-level expression depending only on the bias parameters b1, b2, and b𝒢2.

We first assessed the GoF provided by the model via a χ2 test and explored its range of validity as a function of the minimal scale included in the data vector and of the additional parameter ηmin Eq. (40), progressively excluding nearly isosceles triangles from the analysis (Sect. 6). The predictions for these configurations are, in fact, affected by systematic uncertainties when the Legendre expansion of Eq. (10) is limited to a relatively small number of multipoles, as in any practical application. We find that the tree-level model for the 3PCF fails below separations of 20 h−1 Mpc when ηmin ≥ 2, while a larger minimum scale should be considered for when ηmin = 1. Commonly adopted relations between bias parameters, such as those for b𝒢2(b1) and bΓ3(b1), can help reduce the parameter space without affecting the determination of the other bias parameters, except for the c0 parameter.

The main results of the paper, presented in Sect. 7, are from a joint, full-shape analysis of the 2PCF and 3PCF, carried out with the aim of constraining three cosmological parameters: As, ωcdm, and h. Cosmological analyses involving the 3PCF have so far been based on template fitting (see e.g. Slepian et al. 2017; Slepian & Eisenstein 2018; Veropalumbo et al. 2021; Sugiyama et al. 2021; Farina et al. 2026), as a direct evaluation of the 3PCF model is computationally too expensive. To overcome this hurdle, we developed an emulator able to provide fast and accurate full predictions for both 2PCF and 3PCF, thereby enabling a complete joint analysis of the two statistics.

In this case, in addition to the χ2 test, we adopted an FoB and an FoM to quantify the model’s systematic uncertainties along with its ability to constrain cosmological parameters. Focusing only on the results for the joint analysis, since only this configuration allows us to properly constrain all three parameters without informative priors, we find that a conservative choice for the scale cuts is given by r min 3 PCF = 40 h − 1 Mpc Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,3PCF}} = 40 \, h^{-1} \, \mathrm{{Mpc}} $ and ηmin = 3. A more aggressive configuration with r min 3 PCF = 30 h − 1 Mpc Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,3PCF}} = 30 \, h^{-1} \, \mathrm{{Mpc}} $ and ηmin = 2 is still acceptable, although a tension at the confidence interval of 68.3% with the fiducial value As begins to show. This finding aligns with conclusions drawn from previous methodological studies on the galaxy bias model (see e.g. Veropalumbo et al. 2022, for further details). On the other hand, selecting ηmin = 1 led to more unstable results, particularly at lower redshifts.

Adopting bias relations provides only a limited advantage in a bias-only analysis. However, when cosmological parameters are included, it can lead to tighter constraints, most notably on As, particularly in the case of Eq. (16).

In conclusion, we confirm that, as expected, the combination of 2PCF and 3PCF in real space provides greater constraining power with respect to the 2PCF alone, effectively breaking the degeneracy between As and b1(z). We also expect the improvement to be less pronounced in redshift space, where the 2PCF multipoles typically show reduced degeneracies among cosmological and bias parameters. This work represents a first step towards a full-shape joint analysis in redshift space that, out of necessity, will require tools such as the emulator developed here. In this perspective, this study paves the way for the application of emulated predictions in comparisons to real data, setting the stage for a complete analysis of galaxy clustering data in configuration space that will be on the same footing as power spectrum and bispectrum joint analyses.

This work is part of a series of preparation papers aimed at validating theoretical models for the analysis of galaxy clustering statistics from the Euclid spectroscopic sample. A first assessment of the real-space model for the galaxy power spectrum has been presented in Euclid Collaboration: Pezzotta et al. (2024). The cases of the redshift space power spectrum and 2PCF will be presented respectively in Euclid Collaboration: Camacho et al. (in prep.) and Euclid Collaboration: Kärcher et al. (in prep.). In addition, the redshift-space, joint power spectrum and bispectrum, and 2PCF-3PCF cases will appear in Euclid Collaboration: Pardede et al. (in prep.) and Euclid Collaboration: Pugno et al. (in prep.). Finally, another study will explore the BAO signal in the redshift-space 3PCF.

Acknowledgments

We thank Krister Nagainis, Adi Nusser and Andrea Enia for useful discussions. MG and MM acknowledge support from the research grant ‘Optimizing the extraction of cosmological information from Large Scale Structure analysis in view of the next large spectroscopic surveys’ from MIUR, PRIN 2022 (grant 2022NY2ZRS 001). MM acknowledges the financial contribution from the grant ASI n. 2024-10-HH.0 “Attività scientifiche per la missione Euclid – fase E. 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 (www.euclid-ec.org).

References

  1. Adame, A. G., Aguilar, J., Ahlen, S., et al. 2025, JCAP, 02, 021 [CrossRef] [Google Scholar]
  2. Albrecht, A., Bernstein, G., Cahn, R., et al. 2006, arXiv e-prints [arXiv:astro-ph/0609591] [Google Scholar]
  3. Assassi, V., Baumann, D., Green, D., & Zaldarriaga, M. 2014, JCAP, 08, 056 [CrossRef] [Google Scholar]
  4. Baldauf, T., Seljak, U., Desjacques, V., & McDonald, P. 2012, Phys. Rev. D, 86, 083540 [NASA ADS] [CrossRef] [Google Scholar]
  5. Baldauf, T., Mirbabayi, M., Simonović, M., & Zaldarriaga, M. 2015, Phys. Rev. D, 92, 043514 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bardeen, J. M., Bond, J. R., Kaiser, N., & Szalay, A. S. 1986, ApJ, 304, 15 [Google Scholar]
  7. Baumann, D., Nicolis, A., Senatore, L., & Zaldarriaga, M. 2012, JCAP, 07, 051 [CrossRef] [Google Scholar]
  8. Bernardeau, F., Colombi, S., Gaztañaga, E., & Scoccimarro, R. 2002, Phys. Rep., 367, 1 [Google Scholar]
  9. Beutler, F., Seo, H.-J., Saito, S., et al. 2017, MNRAS, 466, 2242 [Google Scholar]
  10. Blas, D., Garny, M., Ivanov, M. M., & Sibiryakov, S. 2016a, JCAP, 07, 052 [Google Scholar]
  11. Blas, D., Garny, M., Ivanov, M. M., & Sibiryakov, S. 2016b, JCAP, 07, 028 [CrossRef] [Google Scholar]
  12. Cabass, G., Ivanov, M. M., Philcox, O. H. E., Simonović, M., & Zaldarriaga, M. 2022, Phys. Rev. Lett., 129, 021301 [Google Scholar]
  13. Carrasco, J. J. M., Hertzberg, M. P., & Senatore, L. 2012, J. High Energy Phys., 9, 82 [Google Scholar]
  14. Catelan, P., Lucchin, F., Matarrese, S., & Porciani, C. 1998, MNRAS, 297, 692 [CrossRef] [Google Scholar]
  15. Chan, K. C., Scoccimarro, R., & Sheth, R. K. 2012, Phys. Rev. D, 85, 083509 [CrossRef] [Google Scholar]
  16. Coles, P., Moscardini, L., Lucchin, F., Matarrese, S., & Messina, A. 1993, MNRAS, 264, 749 [Google Scholar]
  17. Crocce, M., & Scoccimarro, R. 2008, Phys. Rev. D, 77, 023533 [NASA ADS] [CrossRef] [Google Scholar]
  18. D’Amico, G., Gleyzes, J., Kokron, N., et al. 2020, JCAP, 05, 005 [Google Scholar]
  19. D’Amico, G., Lewandowski, M., Senatore, L., & Zhang, P. 2025, Phys. Rev. D, 111, 063514 [Google Scholar]
  20. Dawson, K. S., Schlegel, D. J., Ahn, C. P., et al. 2013, AJ, 145, 10 [Google Scholar]
  21. DESI Collaboration (Aghamousa, A., et al.) 2016, arXiv e-prints [arXiv:1611.00036] [Google Scholar]
  22. Desjacques, V. 2008, Phys. Rev. D, 78, 103503 [CrossRef] [Google Scholar]
  23. Desjacques, V., Crocce, M., Scoccimarro, R., & Sheth, R. K. 2010, Phys. Rev. D, 82, 103529 [CrossRef] [Google Scholar]
  24. Desjacques, V., Jeong, D., & Schmidt, F. 2018, Phys. Rep., 733, 1 [Google Scholar]
  25. Dore, O., Hirata, C., Wang, Y., et al. 2019, BAAS, 51, 341 [Google Scholar]
  26. Eggemeier, A., Scoccimarro, R., & Smith, R. E. 2019, Phys. Rev. D, 99, 123514 [NASA ADS] [CrossRef] [Google Scholar]
  27. Eggemeier, A., Scoccimarro, R., Crocce, M., Pezzotta, A., & Sánchez, A. G. 2020, Phys. Rev. D, 102, 103530 [NASA ADS] [CrossRef] [Google Scholar]
  28. Eisenstein, D. J., Zehavi, I., Hogg, D. W., et al. 2005, ApJ, 633, 560 [Google Scholar]
  29. Eisenstein, D. J., Seo, H.-J., & White, M. 2007, ApJ, 664, 660 [Google Scholar]
  30. Euclid Collaboration (Blanchard, A., et al.) 2020, A&A, 642, A191 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Euclid Collaboration (Scaramella, R., et al.) 2022, A&A, 662, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Euclid Collaboration (Pezzotta, A., et al.) 2024, A&A, 687, A216 [Google Scholar]
  33. Euclid Collaboration (de la Torre, S., et al.) 2025, A&A, 700, A78 [Google Scholar]
  34. Euclid Collaboration (Mellier, Y., et al.) 2025, A&A, 697, A1 [Google Scholar]
  35. Euclid Collaboration (Risso, I., et al.) 2026, A&A, 707, A233 [Google Scholar]
  36. Fang, X., Eifler, T., & Krause, E. 2020, MNRAS, 497, 2699 [NASA ADS] [CrossRef] [Google Scholar]
  37. Farina, A., Veropalumbo, A., Branchini, E., & Guidi, M. 2026, J. Cosmol. Astropart. Phys., 2026, 028 [Google Scholar]
  38. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
  39. Frieman, J. A., & Gaztanaga, E. 1994, ApJ, 425, 392 [Google Scholar]
  40. Fry, J. N. 1984, ApJ, 279, 499 [NASA ADS] [CrossRef] [Google Scholar]
  41. Fry, J. N. 1994, Phys. Rev. Lett., 73, 215 [Google Scholar]
  42. Fry, J. N., & Gaztañaga, E. 1993, ApJ, 413, 447 [NASA ADS] [CrossRef] [Google Scholar]
  43. Fujita, T., Mauerhofer, V., Senatore, L., Vlah, Z., & Angulo, R. 2020, JCAP, 01, 009 [CrossRef] [Google Scholar]
  44. Gaztañaga, E., Norberg, P., Baugh, C. M., & Croton, D. J. 2005, MNRAS, 364, 620 [Google Scholar]
  45. Gaztañaga, E., Cabré, A., Castander, F., Crocce, M., & Fosalba, P. 2009, MNRAS, 399, 801 [Google Scholar]
  46. Gil-Marín, H., Percival, W. J., Verde, L., et al. 2017, MNRAS, 465, 1757 [Google Scholar]
  47. Grieb, J. N., Sánchez, A. G., Salazar-Albornoz, S., & Dalla Vecchia, C. 2016, MNRAS, 457, 1577 [NASA ADS] [CrossRef] [Google Scholar]
  48. Grieb, J. N., Sánchez, A. G., Salazar-Albornoz, S., et al. 2017, MNRAS, 467, 2085 [NASA ADS] [Google Scholar]
  49. Guidi, M., Veropalumbo, A., Branchini, E., Eggemeier, A., & Carbone, C. 2023, JCAP, 08, 066 [Google Scholar]
  50. Guzzo, L., Pierleoni, M., Meneux, B., et al. 2008, Nature, 451, 541 [Google Scholar]
  51. Hamilton, A. J. S. 2000, MNRAS, 312, 257 [NASA ADS] [CrossRef] [Google Scholar]
  52. Hivon, E., Bouchet, F. R., Colombi, S., & Juszkiewicz, R. 1995, A&A, 298, 643 [Google Scholar]
  53. Hou, J., Sánchez, A. G., Scoccimarro, R., et al. 2018, MNRAS, 480, 2521 [NASA ADS] [CrossRef] [Google Scholar]
  54. Ivanov, M. M., & Sibiryakov, S. 2018, JCAP, 07, 053 [Google Scholar]
  55. Ivanov, M. M., Simonović, M., & Zaldarriaga, M. 2020, JCAP, 05, 042 [CrossRef] [Google Scholar]
  56. Ivezic, Z., Tyson, J. A., Axelrod, T., et al. 2009, Am. Astron. Soc. Meeting Abstr., 213, 460.03 [Google Scholar]
  57. Jing, Y. P., & Börner, G. 2004, ApJ, 607, 140 [Google Scholar]
  58. Kaiser, N. 1984, ApJ, 284, L9 [NASA ADS] [CrossRef] [Google Scholar]
  59. Landy, S. D., & Szalay, A. S. 1993, ApJ, 412, 64 [Google Scholar]
  60. Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints [arXiv:1110.3193] [Google Scholar]
  61. Lazeyras, T., Wagner, C., Baldauf, T., & Schmidt, F. 2016, JCAP, 02, 018 [Google Scholar]
  62. Lippich, M., Sánchez, A. G., Colavincenzo, M., et al. 2019, MNRAS, 482, 1786 [NASA ADS] [CrossRef] [Google Scholar]
  63. Marín, F. 2011, ApJ, 737, 97 [Google Scholar]
  64. Marín, F. A., Blake, C., Poole, G. B., et al. 2013, MNRAS, 432, 2654 [Google Scholar]
  65. Matsubara, T. 2008, Phys. Rev. D, 77, 063530 [NASA ADS] [CrossRef] [Google Scholar]
  66. McBride, C. K., Connolly, A. J., Gardner, J. P., et al. 2011, ApJ, 739, 85 [Google Scholar]
  67. McDonald, P. 2006, Phys. Rev. D, 74, 103512 [CrossRef] [Google Scholar]
  68. McDonald, P., & Roy, A. 2009, JCAP, 08, 020 [CrossRef] [Google Scholar]
  69. Moresco, M., Marulli, F., Moscardini, L., et al. 2017, A&A, 604, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  70. Moresco, M., Veropalumbo, A., Marulli, F., Moscardini, L., & Cimatti, A. 2021, ApJ, 919, 144 [NASA ADS] [CrossRef] [Google Scholar]
  71. Oddo, A., Rizzo, F., Sefusatti, E., Porciani, C., & Monaco, P. 2021, JCAP, 11, 038 [CrossRef] [Google Scholar]
  72. Peacock, J. A., Cole, S., Norberg, P., et al. 2001, Nature, 410, 169 [Google Scholar]
  73. Peebles, P. J. E. 1973, ApJ, 185, 413 [Google Scholar]
  74. Percival, W. J., Cole, S., Eisenstein, D. J., et al. 2007, MNRAS, 381, 1053 [Google Scholar]
  75. Pezzotta, A., de la Torre, S., Bel, J., et al. 2017, A&A, 604, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Philcox, O. H. E., & Ivanov, M. M. 2022, Phys. Rev. D, 105, 043517 [CrossRef] [Google Scholar]
  77. Potter, D., Stadel, J., & Teyssier, R. 2017, Comput. Astrophys. Cosmol., 4, 2 [NASA ADS] [CrossRef] [Google Scholar]
  78. Pozzetti, L., Hirata, C. M., Geach, J. E., et al. 2016, A&A, 590, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  79. Pugno, A., Eggemeier, A., Porciani, C., & Kuruvilla, J. 2025, JCAP, 01, 075 [Google Scholar]
  80. Sánchez, A. G., Kazin, E. A., Beutler, F., et al. 2013, MNRAS, 433, 1202 [Google Scholar]
  81. Sánchez, A. G., Scoccimarro, R., Crocce, M., et al. 2017, MNRAS, 464, 1640 [Google Scholar]
  82. Scoccimarro, R. 2000, ApJ, 542, 1 [Google Scholar]
  83. Scoccimarro, R., Couchman, H. M. P., & Frieman, J. A. 1999, ApJ, 517, 531 [NASA ADS] [CrossRef] [Google Scholar]
  84. Scoccimarro, R., Sefusatti, E., & Zaldarriaga, M. 2004, Phys. Rev. D, 69, 103513 [Google Scholar]
  85. Sefusatti, E., Crocce, M., Pueblas, S., & Scoccimarro, R. 2006, Phys. Rev. D, 74, 023522 [NASA ADS] [CrossRef] [Google Scholar]
  86. Senatore, L., & Zaldarriaga, M. 2015, JCAP, 02, 013 [NASA ADS] [CrossRef] [Google Scholar]
  87. Seo, H.-J., & Eisenstein, D. J. 2003, ApJ, 598, 720 [Google Scholar]
  88. Sheth, R. K., Chan, K. C., & Scoccimarro, R. 2013, Phys. Rev. D, 87, 083002 [NASA ADS] [CrossRef] [Google Scholar]
  89. Slepian, Z., & Eisenstein, D. J. 2015, MNRAS, 454, 4142 [NASA ADS] [CrossRef] [Google Scholar]
  90. Slepian, Z., & Eisenstein, D. J. 2017, MNRAS, 469, 2059 [Google Scholar]
  91. Slepian, Z., & Eisenstein, D. J. 2018, MNRAS, 478, 1468 [Google Scholar]
  92. Slepian, Z., Eisenstein, D. J., Brownstein, J. R., et al. 2017, MNRAS, 469, 1738 [Google Scholar]
  93. Slepian, Z., Eisenstein, D. J., Blazek, J. A., et al. 2018, MNRAS, 474, 2109 [Google Scholar]
  94. Smith, R. E. 2009, MNRAS, 400, 851 [NASA ADS] [CrossRef] [Google Scholar]
  95. Smith, R. E., Scoccimarro, R., & Sheth, R. K. 2007, Phys. Rev. D, 75, 063512 [NASA ADS] [CrossRef] [Google Scholar]
  96. Sugiyama, N. S., Saito, S., Beutler, F., & Seo, H.-J. 2021, MNRAS, 501, 2862 [Google Scholar]
  97. Szalay, A. S. 1988, ApJ, 333, 21 [NASA ADS] [CrossRef] [Google Scholar]
  98. Szapudi, I., & Szalay, A. S. 1998, ApJ, 494, L41 [NASA ADS] [CrossRef] [Google Scholar]
  99. Tröster, T., Sánchez, A. G., Asgari, M., et al. 2020, A&A, 633, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  100. Umeh, O. 2021, JCAP, 05, 035 [Google Scholar]
  101. Verde, L., Wang, L., Heavens, A. F., & Kamionkowski, M. 2000, MNRAS, 313, 141 [Google Scholar]
  102. Veropalumbo, A., Sáez Casares, I., Branchini, E., et al. 2021, MNRAS, 507, 1184 [Google Scholar]
  103. Veropalumbo, A., Binetti, A., Branchini, E., et al. 2022, JCAP, 09, 033 [Google Scholar]
  104. Vlah, Z., Seljak, U., Yat Chu, M., & Feng, Y. 2016, JCAP, 03, 057 [CrossRef] [Google Scholar]
  105. Wang, Y. 2008, Phys. Rev. D, 77, 123525 [NASA ADS] [CrossRef] [Google Scholar]
  106. Wang, Y., Zhao, G.-B., Chuang, C.-H., et al. 2017, MNRAS, 469, 3762 [Google Scholar]
  107. Zhao, G.-B., Wang, Y., Saito, S., et al. 2017, MNRAS, 466, 762 [Google Scholar]
1

For the slight difference with respect to Potter et al. (2017), we refer the reader to the explanation in Euclid Collaboration: Pezzotta et al. (2024).

Appendix A: Dataset

In Fig. A.1 we illustrate how the triangle sides, rij, vary as a function of the triangle index, with each side shown with dashed lines of different colours.

In Fig. A.2 we show the cumulative signal-to-noise ratio (S/N) for the two statistics as a function of the minimal separation rmin computed using the expressions above for the covariance and the ξ ̂ Mathematical equation: $ \hat{\xi} $ and ζ ̂ Mathematical equation: $ \hat{\zeta} $ measurements as

( S N ) ξ 2 : = ∑ r , r ′ ≥ r min ξ ̂ ( r ) C ξ − 1 ( r , r ′ ) ξ ̂ ( r ′ ) Mathematical equation: $$ \begin{aligned} \left(\frac{\mathrm S}{\mathrm N}\right)_\xi ^2 :=\sum _{r,r^{\prime }\ge r_{\rm min}}\hat{\xi }(r)\,C_\xi ^{-1}(r,r^{\prime })\,\hat{\xi }(r^{\prime })\, \end{aligned} $$(A.1)

and

( S N ) ζ 2 : = ∑ r ij , r ij ′ ≥ r min ζ ̂ ( t ) C ζ − 1 ( t , t ′ ) ζ ̂ ( t ′ ) , Mathematical equation: $$ \begin{aligned} \left(\frac{\mathrm S}{\mathrm N}\right)_\zeta ^2 :=\sum _{r_{\rm ij},r_{\rm ij}^{\prime }\ge r_{\rm min}}\hat{\zeta }(t)\,C_\zeta ^{-1}(t,t^{\prime })\,\hat{\zeta }(t^{\prime }), \end{aligned} $$(A.2)

with t and t′ corresponding respectively to the triplets {r12,r13,r23} and {r′12,r′13,r′23}. The S/N for the 2PCF is, as might be expected, higher than the 3PCF, but only slightly, and less so at lower redshifts. The two indeed become comparable at mildly non-linear scales. We notice as well that the total signal for the subset of 3PCF configurations defined by ηmin = 2 can (see Eq. 40) be larger by a factor of a few with respect to the subset defined by ηmin = 2, and, again, it increases at low z.

Thumbnail: Fig. A.1. Refer to the following caption and surrounding text. Fig. A.1.

Triangle sides (rij) as a function of the triangle index, following the condition r12 ≤ r13 ≤ r23.

Thumbnail: Fig. A.2. Refer to the following caption and surrounding text. Fig. A.2.

Cumulative S/N for the 2PCF and 3PCF measurements computed with the predicted Gaussian covariance according to Eqs. (A.1) and (A.2), as a function of the minimal scale included, rmin.

Appendix B: 2PCF and 3PCF emulators

As detailed in Sect. 5, to efficiently sample the likelihood function as across the three cosmological parameters, we developed an emulator to obtain model predictions rapidly and accurately, in the first place, for the 3PCF, but also extending it to include the 2PCF. Indeed, the computational cost of the predictions for the coefficients ζℓ via Eq. (12) as well as the loop integrals coupled to the Hankel transform of Eq. (5) for the 2PCF makes it impractical to explore a broad range of cosmological parameters directly.

Table B.1.

Neural network architecture of the emulator.

To address this issue, we constructed an emulator implemented in PyTorch4. The neural network, whose architecture is described in Table B.1, was trained using input samples drawn from a Latin hypercube sampling of the 3D, cosmological parameter space. This spans the ranges 109As ∈ [1.045, 3.135], h ∈ [0.63, 0.73], and ωcdm ∈ [0.09, 0.15]. The full models of both 2PCF and 3PCF can be written as sums of contributions each factorising the dependence on bias and nuisance parameters. This means that we can keep the analytical dependence on these parameters emulating the cosmology-dependence of each individual contribution. The emulator’s predictions were subsequently validated against a test set drawn from an independent sample to ensure their accuracy and reliability.

The metric chosen for assessing the emulator’s performance compares the systematic error Δξ on each separation bin for which r ≥ 7.5 h−1 Mpc for the 2PCF and Δζ on each triangular configuration for the 3PCF to their statistical error estimated from the covariance computed in Sect. 4 and with the binning choice adopted there. The statistical error is therefore relative to an observed volume of about 54 h−3 Gpc3, larger than any redshift bin expected for the Euclid spectroscopic sample (Euclid Collaboration: Blanchard et al. 2020). In particular, we required that the ratios Δξ/σξ and Δζ/σζ stayed, on average, below the 10% level. The distribution of this metric depends on both the separation (r) and the triangle configuration under consideration, the bias contributions, as well as on the sampled values of the cosmological parameters. This is illustrated by the normalised histograms shown in Fig. B.1. In this test, involving the full model predictions, we assumed a unitary linear bias b1 and c0, the bias relations presented in Eqs. (16) and (17), and the fitting function for b2(b1) obtained from numerical simulations in Lazeyras et al. (2016). All distributions, for both 2PCF and 3PCF at all four redshifts, are characterised by a median value for Δξ/σξ or Δζ/σζ well below 10%, decreasing slightly with redshift.

Thumbnail: Fig. B.1. Refer to the following caption and surrounding text. Fig. B.1.

Normalised histogram of the ratio between the difference between emulated and exact modelling 2PCF predictions and the corresponding uncertainty, Δξ/σξ (top four panels). Different colours refer to different redshift snapshot, while the dashed black line represents the median of the distribution. The same quantity for the 3PCF, Δζ/σζ, is shown in the bottom four panels, for the different redshifts.

This assessment ensures that the emulation process for both the 2PCF and the 3PCF is sufficiently accurate for the analysis presented in this work, but, in practice, also for a direct application to Euclid data or measurements from Stage IV surveys. We estimated a possible effect of systematic errors in the emulated prediction by adding in quadrature the mean value of Δξ and Δζ to the diagonal of the two covariance matrices finding no appreciable difference.

Finally, we notice that the emulator reduces the computational time required for the 3PCF evaluation from approximately 5 minutes per cosmological model to the order of 10−3 seconds.

Appendix C: Cosmological constraints

As a complement to the left panel of Fig. 6, in the in the top right, top left and bottom left panels of Fig. C.1 we present here the corresponding results at higher redshifts, namely z = 1.2, 1.5, and 1.8. We find that the same qualitative conclusions drawn at z = 0.9 continue to hold across redshift, with a general improvement in the precision of parameter inference at lower redshifts. A mild bias in the inferred values of As and b1 is observed, though still well within statistical compatibility with the expected values. We interpret this as a residual effect of projection, and note that all snapshots share the same seed. Moreover, we combined all redshift snapshots to derive constraints on the cosmological parameters As, ωcdm, h, as shown in bottom right panel Fig. C.1. In this combined analysis, we treated the cosmological parameters as shared across all redshifts while allowing the bias parameters to vary independently for each redshift bin, reflecting the expected evolution of galaxy bias with redshift. The combined constraints show significant improvement in precision compared to individual redshift analyses, demonstrating the enhanced constraining power achievable by combining multiple redshifts snapshots. We note that these snapshots are correlated as they share the same simulation seed, which should be taken into account when interpreting the statistical significance of the improvement.

Finally, in Fig C.2 we show the full extension of the right panel of Fig. 6, including all model parameters and, in particular, the non-linear galaxy bias terms, b2, b𝒢2, and bΓ3, and the EFT parameters, c0. This figure also illustrates how the impact of consistency relations on the bias parameter constraints becomes increasingly important as cosmological parameters are marginalised over, with the strongest effects visible for the degeneracy between As and b1. This highlights the critical role of consistency relations in enabling accurate cosmological inference once the cosmological parameter space is opened.

Thumbnail: Fig. C.1. Refer to the following caption and surrounding text. Fig. C.1.

Top right, top left, and bottom left: Joint 2PCF+3PCF cosmological constraints and the linear bias b1 for individual redshift snapshots at z = 1.2, z = 1.5, and z = 1.8, following the format of Fig. 6. Bottom right: Combined cosmological constraints from all redshift bins (z = 0.9, 1.2, 1.5, and 1.8).

Thumbnail: Fig. C.2. Refer to the following caption and surrounding text. Fig. C.2.

Same as the right panel of Fig. 6 but showing the full set of model parameters.

All Tables

Table 1.

Fiducial parameters of the flat ΛCDM model.

In the text
Table 2.

Key properties of the synthetic galaxy catalogues.

In the text
Table 3.

Cosmological and nuisance parameters of the 2PCF and 3PCF models.

In the text
Table B.1.

Neural network architecture of the emulator.

In the text

All Figures

Thumbnail: Fig. 1. Refer to the following caption and surrounding text. Fig. 1.

Top panels: 2PCF measurements for four comoving snapshots from the Flagship I simulation of the Model 3 HOD. Bottom panels: Same but for the 3PCF measurements. The triangular configurations are ordered by increasing values of r12, r13, and r23 under the condition r12 ≤ r13 ≤ r23, as depicted in Fig. A.1. The circles are colour-coded according to the η, while the vertical grey lines mark a group of triangles that share the same value of the smaller side (r12).
In the text
Thumbnail: Fig. 2. Refer to the following caption and surrounding text. Fig. 2.

GoF test for the bias 2PCF and 3PCF models with fixed cosmological parameters as a function of the minimal scale, r min 2 PCF Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{2PCF}} $ and r min 3 PCF Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{3PCF}} $, respectively, assumed for the analysis. The posterior-averaged χ2 is compared to the 68% and 95% confidence intervals for the χ2 distribution, denoted by the grey shaded areas, given the corresponding Nd.o.f. (fixing the case ηmin = 1 for the 3PCF confidence levels; see the main text for more details). For the 3PCF, we considered three different values for the ηmin parameter, Eq. (40), progressively excluding near-isosceles configurations. In the joint analysis, the minimal scale for the 2PCF was kept fixed at r min 2 PCF = 25 h − 1 Mpc Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{2PCF}}=25\, h^{-1} \, \mathrm{{Mpc}} $, with the results shown as a function of r min 3 PCF Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{3PCF}} $ alone. The four columns span the four simulation snapshots at redshifts (z) of 0.9, 1.2, 1.5, and 1.8.
In the text
Thumbnail: Fig. 3. Refer to the following caption and surrounding text. Fig. 3.

Left panel: Marginalised 2D constraints at confidence intervals of 68.3% and 99.5%, on the nuisance parameters (bias and counter-term) from the 2PCF analysis (grey shaded areas) against the joint 2PCF and 3PCF analysis for several values of r min 3 PCF Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,3PCF}} $ at z = 1.5 (dashed blue lines, dashed orange lines, and green shaded areas). All cases assume η = 3; see Eq. (40). The vertical grey band is an estimate of the linear bias b1 obtained from a comparison of the galaxy power spectrum with measurements of the matter power spectrum in Euclid Collaboration: Pezzotta et al. (2024). Right panel: Marginalised 2D constraints at confidence intervals of 68.3% and 99.5%, on the nuisance parameters (bias and counter-term) from the joint 2PCF-3PCF analysis at z = 1.5 for r min 3 PCF = 30 h − 1 Mpc Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,3PCF}}=30\, h^{-1} \, \mathrm{{Mpc}} $ (red contours) compared to the same set-up but with the assumption of the b𝒢2(b1) bias relation alone (red contours) and the combination of the relations b𝒢2(b1) and bΓ3(b1) (violet contours).
In the text
Thumbnail: Fig. 4. Refer to the following caption and surrounding text. Fig. 4.

Marginalised constraints on the cosmological parameters and linear bias from a joint 2PCF and 3PCF analysis. The posterior mean along with its uncertainties is displayed as a function of the minimal scale ( r min 3 PCF Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,3PCF}} $) for different choices of ηmin and fixing r min 2 PCF = 25 h − 1 Mpc Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,2PCF}}=25\, h^{-1} \, \mathrm{{Mpc}} $ and rmax = 140 h−1 Mpc for both statistics. The dashed black line refers to the fiducial values assumed in the simulation, while the dashed grey line refers to the fiducial b1 adopted in Euclid Collaboration: Pezzotta et al. (2024), also shown in Fig. 3.
In the text
Thumbnail: Fig. 5. Refer to the following caption and surrounding text. Fig. 5.

Performance metrics for the full-shape, joint analysis of the 2PCF and the 3PCF. All quantities are shown as a function of the minimal scale ( r min 3 PCF Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,3PCF}} $), with r min 2 PCF Mathematical equation: $ r^{\mathrm{2PCF}}_{\mathrm{min}} $ fixed at 25 h−1 Mpc, and for the three different values of ηmin in different colour shades. From top to bottom we show the posterior-averaged ⟨χ2⟩post for the GoF, the FoB, and the FoM, these last two both defined in terms of the three cosmological parameters As, ωcdm, and h. The grey bands denote the 68% and 99.7% confidence levels for the χ2 distribution and the FoB.
In the text
Thumbnail: Fig. 6. Refer to the following caption and surrounding text. Fig. 6.

Left panel: Marginalised 2D and 1D posterior distributions at z = 0.9 from the joint 2PCF and 3PCF analysis of the cosmological parameters h, ωcdm, and As and the linear bias parameter b1 for various 3PCF scale cuts ( r min 3 PCF Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,3PCF}} $) and the ηmin parameter defined in Eq. (40), fixing r min 2 PCF = 25 h − 1 Mpc Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,2PCF}} = 25 \, h^{-1} \, \mathrm{{Mpc}} $. Grey contours correspond to the 2PCF-only analysis assuming a Gaussian prior on b1 given by the 10% uncertainties around the best-fit value obtained from the joint analysis. Each panel shows the results for a different snapshot. Vertical and horizontal lines denote the fiducial values of the cosmological parameters. Right panel: Marginalised 1D and 2D posterior distributions of the same subset of model parameters for the combined analysis of 2PCF and 3PCF at z = 0.9 fixing r min 2 PCF = 25 h − 1 Mpc Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,2PCF}} = 25 \, h^{-1} \, \mathrm{{Mpc}} $. The reference case of the maximal model with r min 3 PCF = 30 h − 1 Mpc Mathematical equation: $ r_\mathrm{{min}}^\mathrm{{\,3PCF}} = 30 \, h^{-1} \, \mathrm{{Mpc}} $ and ηmin = 3 is compared to the results obtained adopting the bias relations of Eq. (16), in red, and its combination with Eq. (17), in brown.
In the text
Thumbnail: Fig. A.1. Refer to the following caption and surrounding text. Fig. A.1.

Triangle sides (rij) as a function of the triangle index, following the condition r12 ≤ r13 ≤ r23.
In the text
Thumbnail: Fig. A.2. Refer to the following caption and surrounding text. Fig. A.2.

Cumulative S/N for the 2PCF and 3PCF measurements computed with the predicted Gaussian covariance according to Eqs. (A.1) and (A.2), as a function of the minimal scale included, rmin.
In the text
Thumbnail: Fig. B.1. Refer to the following caption and surrounding text. Fig. B.1.

Normalised histogram of the ratio between the difference between emulated and exact modelling 2PCF predictions and the corresponding uncertainty, Δξ/σξ (top four panels). Different colours refer to different redshift snapshot, while the dashed black line represents the median of the distribution. The same quantity for the 3PCF, Δζ/σζ, is shown in the bottom four panels, for the different redshifts.
In the text
Thumbnail: Fig. C.1. Refer to the following caption and surrounding text. Fig. C.1.

Top right, top left, and bottom left: Joint 2PCF+3PCF cosmological constraints and the linear bias b1 for individual redshift snapshots at z = 1.2, z = 1.5, and z = 1.8, following the format of Fig. 6. Bottom right: Combined cosmological constraints from all redshift bins (z = 0.9, 1.2, 1.5, and 1.8).
In the text
Thumbnail: Fig. C.2. Refer to the following caption and surrounding text. Fig. C.2.

Same as the right panel of Fig. 6 but showing the full set of model parameters.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.