# A joint 2- and 3-point clustering analysis of the VIPERS PDR2 catalogue at
$z\sim 1$: breaking the degeneracy of cosmological parameters

Alfonso Veropalumbo1,2, Iñigo Sáez Casares3, Enzo Branchini1,2,4, Benjamin R.
Granett5,6, Luigi Guzzo5,6, Federico Marulli7,8,9, Michele Moresco7,8, Lauro
Moscardini7,8,9, Andrea Pezzotta10,11,12, Sylvain de la Torre13

1Dipartimento di Matematica e Fisica, Università degli studi Roma Tre, Via
della Vasca Navale, 84, 00146 Roma, Italy
2INFN - Sezione di Roma Tre, via della Vasca Navale 84, I-00146 Roma, Italy
3 École Normale Supérieure Paris-Saclay, 4 avenue des Sciences, 91190 Gif-sur-
Yvette, France
4 INAF - Osservatorio Astronomico di Roma, via Frascati 33, I-00040 Monte
Porzio Catone (RM), Italy
5 INAF - Osservatorio Astronomico di Brera, Via Brera 28, 20122 Milano, via E.
Bianchi 46, 23807 Merate, Italy
6 Università degli Studi di Milano, via G. Celoria 16, 20133 Milano, Italy
7 Dipartimento di Fisica e Astronomia “Augusto Righi” - Alma Mater Studiorum
Università di Bologna, via Piero Gobetti 93/2,
I-40129 Bologna, Italy
8 INAF - Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, via
Piero Gobetti 93/3, I-40129 Bologna, Italy
9 INFN - Sezione di Bologna, viale Berti Pichat 6/2, I-40127 Bologna, Italy
10 Max-Planck-Institut für extraterrestrische Physik, Postfach 1312,
Giessenbachstr., 85741 Garching, Germany
11 Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans,
s/n, 08193 Barcelona, Spain
12 Institut d‘Estudis Espacials de Catalunya (IEEC), 08034 Barcelona, Spain
13 Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We measure the galaxy 2- and 3-point correlation functions at $z=[0.5,0.7]$
and $z=[0.7,0.9]$, from the final data release of the VIPERS survey (PDR2). We
model the two statistics including a nonlinear 1-loop model for the 2-point
function and a tree-level model for the 3-point function, and perform a joint
likelihood analysis. The entire process and nonlinear corrections are tested
and validated through the use of the 153 highly realistic VIPERS mock
catalogues, showing that they are robust down to scales as small as 10
$h^{-1}\,\mathrm{Mpc}$. The mocks are also adopted to compute the covariance
matrix that we use for the joint 2- and 3-point analysis. Despite the limited
statistics of the two (volume-limited) sub-samples analysed, we demonstrate
that such a combination successfully breaks the degeneracy existing at 2-point
level between clustering amplitude $\sigma_{8}$, linear bias $b_{1}$ and the
linear growth rate of fluctuations $f$. For the latter, in particular, we
measure $f(z=0.61)=0.64^{+0.55}_{-0.37}$ and $f(z=0.8)=1.0\pm 1.0$, while the
amplitude of clustering is found to be $\sigma_{8}(z=0.61)=0.50\pm 0.12$ and
$\sigma_{8}(z=0.8)=0.39^{+0.11}_{-0.13}$. These values are in excellent
agreement with the extrapolation of a Planck cosmology.

###### keywords:

large-scale structure of Universe – galaxies: statistics – cosmology:
observations

††pubyear: 2021††pagerange: A joint 2- and 3-point clustering analysis of the
VIPERS PDR2 catalogue at $z\sim 1$: breaking the degeneracy of cosmological
parameters–C

##  1 Introduction

Galaxy clustering has emerged over the past three decades as one main
observational pillar supporting the current standard model of cosmology.
Cosmological constraints were so far derived mostly from two-point statistics
of the galaxy distribution, either in configuration or in Fourier space (see
e.g. Alam et al., 2021, for a comprehensive recent review). Two-point
statistics will remain a central probe also for the next generation of
redshift surveys, which have just started or are about to start. This includes
in particular the measurement of specific features, as baryonic acoustic
oscillations (BAO) (e.g., Cole et al., 2005; Eisenstein et al., 2005) and
redshift space distortions (RSD) (e.g., Peacock et al., 2001; Guzzo et al.,
2008), which are part of the standard “dark energy" probes in projects like
DESI (DESI Collaboration et al., 2016a, b), the Euclid (Laureijs et al., 2011)
and Roman (Akeson et al., 2019) space telescopes or, in the radio, the Square
Kilometre Array (SKA, Maartens et al., 2015).

For a Gaussian random field, this would be enough to fully characterize the
field statistically. However, significant information exists in the galaxy
distribution, beyond the two-point correlation function (2PCF hereafter) or
the power spectrum. This is locked in the n-point correlation functions, or
higher-moments, of the field. The simplest of these, the 3-point correlation
function (3PCF hereafter), is sensitive to non-Gaussian features in the
primordial density perturbations, non-linear effects in their evolution and
galaxy biasing, i.e. the relationship between galaxy tracers and the
underlying distribution of matter. Since these effects are characterised by a
different higher-order signal (see e.g. Fry, 1994; Matarrese et al., 1997),
the measurements of the 3PCF (or equivalently in Fourier space, the
bispectrum) has the potential to separate them and constrain their relative
strength and evolution. In addition, BAO and RSD also leave an imprint on
higher-order statistics that can be exploited to remove parameter degeneracy
and improve cosmological constraints.

The advantage of combining 2- and 3-point statistics to improve cosmological
constraints has been realised in recent years (Sefusatti et al., 2006;
Yankelevich & Porciani, 2019). These works showed that the combination of the
3D power spectrum and bispectrum can significantly reduce statistical errors
in the determination of cosmological parameters, while breaking degeneracies
among these. The benefit is remarkable when using clustering data alone. These
considerations have inspired the recent combined 2- and 3-point analyses on
SDSS data carried out by Gil-Marín et al. (2015); Gil-Marín et al. (2017).

In configuration space, the measurement of the 3PCF is computationally more
demanding. Examples in the literature include estimates from the 2dFGRS (Jing
& Börner, 2004), various releases of the SDSS main sample (Kayo et al., 2004;
Nichol et al., 2006; Kulkarni et al., 2007; McBride et al., 2011a, b; Marín,
2011; Guo et al., 2014), the BOSS CMASS (Guo et al., 2015) and the first
Public Data Release (PDR1) of the VIPERS survey (Moresco et al., 2017). Only
recently, however, thanks to a breakthrough in the efficiency of available
estimators, it has become possible to push measurements of the 3PCF to scales
comparable to those of the BAO peak (Slepian & Eisenstein, 2015; Slepian et
al., 2017a, b), also using galaxy clusters as tracers (Moresco et al., 2020).
Finally, combined 2PCF and 3PCF analyses in configuration space were performed
by Marín et al. (2013) using the WiggleZ spectroscopic galaxy survey.

In this paper, we perform for the first time a joint 2- and 3-point
correlation analysis of the final PDR2 catalogue of the VIPERS survey
(Scodeggio et al., 2018; Guzzo et al., 2014; Garilli et al., 2014), applying
the most advanced estimators and extracting joint cosmological constraints out
to $z=0.9$. One main goal of VIPERS was to constrain the growth rate of
structure $f$ out to $z\sim 1$, a result that was successfully achieved
through a series of complementary estimates (de la Torre et al., 2013, 2017;
Pezzotta et al., 2017; Mohammad et al., 2018). All these analyses were
developed in configuration space, for which the survey geometry and footprint
can be more easily handled, if compared to the complex window function
convolution one has to deal with in Fourier space (Rota et al., 2017). The
main drawback when working in configuration space is the large covariance of
2PCF and 3PCF data and their errors. However, thanks to the availability of
efficient clustering estimators and of the large ensemble of realistic mock
catalogues made available by the VIPERS collaboration (de la Torre et al.,
2017), we can solve the problem by numerically estimating the covariance
matrix of the data.

We thus build upon previous VIPERS analyses to perform a joint 2PCF and 3PCF
measurement, which allows us to break the degeneracy which affects some key
cosmological parameters when these are estimated from the 2PCF alone.
Specifically, we obtain separate estimates of the linear growth rate of
structure, $f$, the clustering amplitude, $\sigma_{8}$, and some of the
parameters entering the galaxy biasing relation. This work, in particular,
expands upon the 3PCF analysis performed by Moresco et al. (2017) on the first
data release of VIPERS (PDR1), by 1) using the final, larger PDR2 release; 2)
exploring all triangle configurations, rather than just a subset; 3) combining
2- and 3-point statistics to obtain joint cosmological constraints. It is also
complementary to Cappi et al. (2015) and Di Porto et al. (2016), where the
non-linearity and evolution of galaxy bias in VIPERS was first studied.

The layout of the paper is as follows. In Sect. 2 we briefly describe the
VIPERS PDR2 catalogue and the mock data used to estimate errors. In Sect. 3 we
describe the 2PCF and 3PCF estimators and perform validation tests using the
mock samples. Covariance matrices and their estimates are discussed in Sect.
4. In Sect. 5 we present the 2PCF and 3PCF models used for the likelihood
analysis of Sect. 6. The results of the 2PCF, the 3PCF and their joint
analyses are presented in Sect. 7. Finally, we discuss the results and draw
our conclusions in Sect. 8.

Throughout the work, unless otherwise specified, we assume a flat
$\Lambda$-cold dark matter ($\Lambda$CDM) cosmological model characterized by
the following parameters:
$\\{\Omega_{M},\Omega_{b},n_{s}\\}=\left(0.3,0.045,0.96\right)$. The Hubble
constant is defined as
$H_{0}=100\,h\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1}$.

##  2 Datasets

###  2.1 VIPERS data

Table 1: Definition of the sub-samples of the VIPERS PDR2 catalogue that we have considered in this work. P1 and P2 are magnitude-limited samples that match those analysed in Pezzotta et al. (2017) and are used here for validation and consistency tests. G1 and G2 are instead volume-limited samples, upon which the joint clustering analysis specific of this work is based. Name |  $z-$range | Magnitude Cut | Number of objects
---|---|---|---
G1 | $0.5\leq z<0.7$ | $M_{B}<-19.3+(0.7-z)$ | 23352
G2 | $0.7\leq z<0.9$ | $M_{B}<-20.3+(0.9-z)$ | 13046
P1 | $0.5\leq z<0.7$ | $i_{\rm AB}\leq 22.5$ | 30764
P2 | $0.7\leq z<1.2$ | $i_{\rm AB}\leq 22.5$ | 35734

The VIMOS Public Extragalactic Redshift Survey (VIPERS) has been completed as
one of the ESO Large Programmes. It was designed to build a spectroscopic
sample of about 100,000 galaxies, aiming at an optimal combination of depth
(reaching beyond $z\simeq 1$), volume and sampling density (Guzzo et al.,
2014). This is obtained covering $\sim$ 24 deg2 over the W1 and W4 fields of
the Canada-France-Hawaii Telescope Legacy Survey, which provides accurate
photometry in five bands. This area was tiled with a mosaic of 288 pointings
with the VIMOS spectrograph at the ESO VLT, measuring moderate-resolution
spectra ($R\simeq 220$) for galaxies brighter than $i_{\rm AB}=22.5$. A colour
pre-selection in the $(r-i)$ vs $(u-g)$ plane was applied prior to the
spectroscopic observations, efficiently and accurately excluding objects with
$z<0.5$ and boosting the spectroscopic sampling to nearly $50\%$. The root
mean square (rms) redshift measurement error of these data is
$\sigma_{z}=5\times 10^{-4}(1+z)$, corresponding to 167 km s-1. For
consistency with previous clustering analyses of the VIPERS data, we only
consider objects with a redshift confirmation rate larger than 96.1%,
corresponding to quality flags from 2 to 9. More details on the survey design
and the final data release can be found in Guzzo et al. (2014) and Scodeggio
et al. (2018), respectively.

For the work presented here, we extract from the PDR2 catalogue four different
sub-samples, whose characteristics are summarised in Table 1. The samples P1
and P2 are magnitude-limited samples corresponding to two redshift bins,
$z=[0.5,0.7]$ and $[0.7,1.2]$, following the selection criteria used in
Pezzotta et al. (2017). They will be used to validate our analysis against
previous measurements from the same data. This selection maximises the number
of available tracers, at the price of a redshift-dependent selection function
and a mean density and bias that vary with redshift, especially in the outer
redshift bin, which reaches $z=1.2$ with a rather sparse bright population.

The samples named G1 and G2 are those used for the specific new analyses of
this paper. They are two volume-limited, non-overlapping bins at $z=[0.5,0.7]$
and $[0.7,0.9]$, coinciding with the samples named L1 and L3 in the SHAM
analysis of PDR2 by Granett et al. (2019), (see also Davidzon et al., 2016).
The absolute magnitude thresholds guarantee completeness above 90%, while
evolution (both in luminosity and bias) is minimised by the limited redshift
size of the bins. Galaxies within G1 and G2 are brighter on average than those
in P1 and P2. As a consequence, the measured clustering amplitude is expected
to be larger due to the higher galaxy bias (Marulli et al., 2013), hence
partly compensating for the larger Poisson noise.

###  2.2 Weights

The completeness of the spectroscopic sample is quantified in the PDR2
catalogue by a direction-dependent target sampling rate (TSR) and a
spectroscopic success rate (SSR) (Scodeggio et al., 2018). To account for
these, we have used the same weighting scheme as in Pezzotta et al. (2017). To
each object, a statistical weight $w_{i}$ is assigned as

$w_{i}=w_{i}^{TSR}\cdot w_{i}^{SSR}\,\,\,.$ (1)

Here, $w_{i}^{TSR}$ is defined as the ratio of the local surface densities of
target and parent galaxies (i.e. before and after applying the target
selection) within an aperture of $60\times 100$ arcsec2. $w_{i}^{SSR}$ is
instead the local fraction of observed spectra with reliable redshift
measurement with respect to the target sample. We also correct for the small-
scale bias introduced by slit “collisions" by up-weighting each galaxy-galaxy
pair at a given angular separation according to de la Torre et al. (2013). On
the scales of our analysis, this last correction impacts the estimated 2PCF
quadrupole moment, but not the monopole. We refer the interested reader to
Sect. 4 of Pezzotta et al. (2017) for further details on the weighting scheme.

###  2.3 Mock VIPERS data

To validate our analysis and quantify statistical errors, we use the publicly
available VIPERS mock catalogues111The VIPERS mock catalogues are publicly
available at http://www.vipers.inaf.it/rel-pdr2.html#mocks . These were
generated from lightcones extracted from the Big MultiDark N-body simulation
(Klypin et al., 2016). The cosmology of that simulations is a flat
$\Lambda$CDM characterised by the set of parameters
$(\Omega_{M},\Omega_{\Lambda},\Omega_{b},h,n_{s},\sigma_{8})=(0.307,0.693,0.0482,0.678,0.960,0.823)$,
which are slightly different from those used in this work. Dark matter halos
were identified in the parent simulation and populated with synthetic galaxies
down to the faint magnitude limit of the survey, as detailed in de la Torre et
al. (2013, 2017). This resulted in 153 independent mock VIPERS W1+W4
catalogues, on which the same footprint and selection function of the real
survey were applied. These mocks were designed to match the luminosity
function, number density and redshift distribution of VIPERS galaxies.
Previous analyses have shown that they also reproduce the VIPERS 2PCF (Granett
et al., 2019). We will see in Sect. 3.2 that, within the errors, they also
match the 3PCF of VIPERS galaxies.

For each mock catalogue, four subsamples corresponding to the G1, G2, P1 and
P2 selections of Table 1 were created.

##  3 Clustering Measurements

Figure 1: Monopole (blue circles) and quadrupole (orange squares) moments of
the 2PCF measured from the G1 (left) and G2 (right) VIPERS sub-samples. The
corresponding measurements from the 153 VIPERS mock catalogues are also
plotted (thin solid lines), together with their average (thick solid line).
The error bars on the data points correspond to the rms dispersion of the
mocks. The consistency between the mocks and the real data highlights the very
good fidelity in clustering properties of the mock galaxies.

The estimators used in this work belong to the class of unbiased, minimum
variance N-point estimators proposed by Szapudi & Szalay (1998). Their general
form is

$\xi^{N}=\frac{(D-R)^{N}}{R^{N}},$ (2)

where $D$ is the data catalogue, $R$ is the so-called “random" sample, and $N$
is the order of the correlation statistics. The estimate relies on counting
and binning all N-tuples $D^{q}R^{p}$ formed by $q$ data and $p$ random
objects, with $q+p=N$.

The "random" catalogue of objects is a synthetic sample with the same geometry
and selection function as the real survey but with no spatial clustering. We
assume that the selection function, i.e. the probability to observe a galaxy
at the spatial position $(\alpha,\delta,z)$, can be factorised as

$P(\alpha,\delta,z)=f(\alpha,\delta)N(z)\,,$ (3)

where $f(\alpha,\delta)$ accounts for the angular footprint of the survey and
$N(z)$ is the redshift distribution of the sources in the catalogue.

All measurements presented here have been performed using the estimators
described below and publicly available in the library
CosmoBolognaLib222https://gitlab.com/federicomarulli/CosmoBolognaLib (Marulli
et al., 2016).

###  3.1 2PCF estimator

Figure 2: Top panels: estimates of the three-point correlation function for
the G1 sample using both the triplet counting method (black dots) and the SHD
method (solid lines), for different $l_{max}$ as indicated by the legend (top
panels). The left and right panels correspond to isosceles
($r_{12}=r_{13}=10\,h^{-1}\,\mathrm{Mpc}$) and non-isosceles
($r_{12}=10\,h^{-1}\,\mathrm{Mpc},r_{13}=30\,h^{-1}\,\mathrm{Mpc}$)
configurations, respectively. The bottom panels show the difference between
the two estimators, for the different choices of $l_{max}$ (solid lines). The
dashed black curves correspond to $\pm 10\%$ of the statistical scatter of the
mocks, highlighting the relative importance of the systematic errors
introduced by the choice of $l_{max}$.

To estimate the anisotropic 2PCF we use the Landy & Szalay (1993) estimator;

$\xi(s,\mu)\,=\frac{DD(s,\mu)-2DR(s,\mu)-RR(s,\mu)}{RR(s,\mu)},$ (4)

where $DD$, $DR$ and $RR$ are the data-data, data-random and random-random
pairs of objects, respectively. The separation vector of two objects has
modulus $s$ and forms a cosine angle $\mu=\cos(\theta)$ with the line of sight
to the pair, measured at the midpoint of the separation vector itself.

We bin the pair counts using a constant logarithmic bin in $s$,
$\Delta\log(s)=0.1$, and a linear one for $\mu$, $\Delta\mu=0.05$, as in
Pezzotta et al. (2017).

For each $s$ bin we compute the 2PCF multipoles

$\xi_{l}(s)=\frac{2l+1}{2}\int_{-1}^{1}\mathrm{d}\mu\mathcal{P}_{l}(\mu)\xi(s,\mu)\,,$
(5)

where $\mathcal{P}_{\ell}(\mu)$ are the Legendre polynomials. We only consider
the monopole $\ell=0$ and quadrupole $\ell=2$, since odd multipoles are zero
by design, and the measured hexadecapole turns out to be consistent with zero
within the (large) errors.

To validate our estimator, we measure the 2PCF multipoles for the subsamples
of Table 1 and compare the results with those of previous VIPERS studies. Fig.
1 shows the 2PCF monopole and quadrupole moments of the VIPERS G1 and G2
samples (blue circles and the orange squares, respectively). The solid thin
lines correspond to measurements of each mock catalogue, with the thicker line
and errorbars on data points showing the average and rms of the mocks
respectively. This plot shows how the mock VIPERS samples reproduce faithfully
the clustering of the real data, suggesting that they are adequate to estimate
their errors and covariance. We have repeated the test using the P1 and P2
samples and verified that the results match those of Pezzotta et al. (2017,
see Fig. 12)

The maximum separation considered in our analysis is set by the signal-to-
noise of the 2PCF. The VIPERS footprint is highly elongated along the
(equatorial) longitude, which significantly reduces the number of distant
pairs in the transverse direction. As a result, the estimated quadrupole
moment becomes noisy beyond $\sim 40\,h^{-1}\,\mathrm{Mpc}$. We therefore set
$s_{max}=40\,h^{-1}\,\mathrm{Mpc}$ in our 2PCF analysis. In Sect. 7.1 we will
also test the sensitivity of the results to the smallest separation scale
considered, showing that the results are robust down to
$s_{min}=15\,h^{-1}\,\mathrm{Mpc}$, which we then assume as the default value.
In summary, our baseline range for the 2PCF analysis is
$[15,40]\,h^{-1}\,\mathrm{Mpc}$.

###  3.2 3PCF estimator

Figure 3: Monopole moments (black circles) of the three-point correlation
function for the G1 and G2 samples (top and bottom panels, respectively),
measured for some selected sets of triangles of different side lengths, as
indicated on the x-axis. The corresponding measurements from the 153 VIPERS
mock catalogues are also shown by the thin grey lines, together with their
averages (solid black line). Error bars correspond to the rms scatter among
the mocks.

To measure the 3PCF of the VIPERS first data release catalogue (PDR1), Moresco
et al. (2017) used the estimator of Szapudi & Szalay (1998) because of its
ability to account for complicated survey geometry (Kayo et al., 2004). This
estimator, however, is computationally demanding since it relies on brute-
force galaxy triplet counting, whose computing time scales as $O(N^{3})$. This
becomes unbearable if applied to 153 mock catalogues with the size of VIPERS.

A major breakthrough in this respect has been the introduction by Slepian &
Eisenstein (2015), of a more efficient estimator based on a spherical
harmonics decomposition (hereafter SHD). For this algorithm, the computational
cost scales as $O(N^{2})$, making it successfully applicable to large galaxy
samples such as the SDSS-DR12 CMASS catalogue (Slepian et al., 2017b). The
method relies on the 3PCF Legendre polynomials expansion proposed by Szapudi
(2004)

$\zeta(s_{12},s_{13},\mu)\,=\sum_{l=0}^{l_{max}}\zeta_{l}(s_{12},s_{13})\mathcal{P}_{l}(\mu)\,,$
(6)

where the 3PCF $\zeta(s_{12},s_{13},\mu)$ is parametrized by two triangle
sides $s_{12}$, $s_{13}$ and their cosine angle
$\mu=\cos(\hat{s}_{12}\cdot\hat{s}_{13})$. This expansion offers two
advantages. The first one is that the multipole moments
$\zeta_{l}(s_{12},s_{13})$ can be efficiently estimated by locally expanding
the density field at distances $s_{12}$ and $s_{13}$ in spherical harmonics
and then by cross-correlating the expansion coefficients using the spherical
harmonics addition theorem. The second advantage is that Eq. (6) typically
requires a limited number of multipoles to converge. These properties
dramatically reduce the computational cost, allowing us to measure the 3PCF in
all VIPERS mocks and consider all triangle configurations. We refer the reader
to Slepian & Eisenstein (2015) for a detailed description of the SHD
estimator.

The main disadvantages of this estimator are that, for a given configuration,
all triangles are mixed together; this is particularly relevant in all cases
where the third side spans a large range of scales. For example, in the case
of isosceles configurations with $s_{12}=s_{13}$, the third side $s_{23}$
varies from 0 to $2\times s_{12}$, i.e. well into the highly nonlinear regime,
where theoretical predictions cannot be trusted. For this reason, here we will
only consider triangles with all side lengths above a minimum value
$s_{min}=15\,h^{-1}\,\mathrm{Mpc}$ and up to a maximum length
$s_{max}=40\,h^{-1}\,\mathrm{Mpc}$. This matches the scale used in the 2PCF
analysis. In Sect. 7.2, we shall test the robustness of our results to such a
choice.

####  3.2.1 Convergence of the SHD method

We first assess the sensitivity of the SHD method to the choice of $l_{max}$
in Eq. 6. We thus compare the 3PCF of the G1 catalogue measured with the SHD
method to that, supposedly exact, estimated with the same brute-force triplet
counting technique used by Moresco et al. (2017). We show in Fig. 2 the
results for two triangle configurations: isosceles triangles with
$s_{12}=s_{13}=10\,h^{-1}\,\mathrm{Mpc}$ (left panel) and non-isosceles
triangles with $s_{12}=10\,h^{-1}\,\mathrm{Mpc}$ and
$s_{13}=30\,h^{-1}\,\mathrm{Mpc}$ (right panel). In both cases, the same $\mu$
binning is used.

In the top panels, we compare the reference 3PCF measured with the triplet
counting technique (black dots), with the one measured with the SHD method
using $l_{max}=10,\,20$ and $30$ (lines with different colours). These
particular choice of configurations cover both the cases where highly non-
linear scales are included (isosceles, left panel) or not (non-isosceles,
right panel). In the bottom panels, we show the difference of the 3PCFs
estimated with the two methods. As an indication of the relative value of the
potential systematic error with respect to statistical errors, the black
dashed lines indicate 10% of the rms scatter among the mock catalogues. For
non-isosceles configurations (right panel), the SHD matches well the expected
result, with differences of the order (or below) 1% of the random errors,
almost independently of the choice of $l_{max}$. This implies that a limited
number of multipoles $\ell$ is sufficient for Eq. (6) to converge to the
correct result, typically $l<10$.

On the contrary, the results of the SHD method for the isosceles
configurations (left panel) depend on the choice of $l_{max}$ and $\mu$. For
$l_{max}=30$, the results are generally satisfactory (the difference between
SHD and triplet counting is of the order of a few percent of the statistical
error), except for $\mu=1$, i.e. when $s_{23}\rightarrow 0$. This is not
surprising, since elongated isosceles triangle shapes are difficult to
reproduce in harmonics space.

Based on these results, we adopt a conservative approach in which: 1) we set
$l_{max}=30$ in the 3PCF estimate, and 2) we exclude the case $\mu=1$ for
isosceles configurations. We also choose to parameterise the 3PCF in terms of
the three side lengths $\zeta(s_{12},s_{13},s_{23})$ instead of two sides and
the cosine angle, $\zeta(s_{12},s_{13},\mu)$. Since the SHD method involves
binning in the radial direction, we perform an additional step, described in
Appendix A, to use the same binning along the three sides $s_{12},s_{13}$ and
$s_{23}$.

Finally, we sort the triangles in increasing size $s_{12}\leq s_{13}\leq
s_{23}$, to avoid repetitions. This choice also facilitates triangle
selection: for example, we can exclude scales smaller than a chosen value by
setting a threshold for one triangle side only.

####  3.2.2 Comparison with the 3PCF of mock samples

The second test is analogous to the one performed to validate the 2PCF
estimator: we measure the 3PCF in the G1, G2, P1 and P2 VIPERS samples and
compare the results with the same quantity measured in the mocks.

In Fig. 3 we show the results for the G1 and G2 samples (top and bottom
panels, respectively). We detect a non-zero 3PCF signal up to the largest
scales considered in our analysis, in agreement with Moresco et al. (2017).
The VIPERS 3PCF agrees well with that of the mock samples, demonstrating that
the VIPERS mocks are indeed adequate to estimate the errors for the 3PCF
analysis, too.

##  4 Covariant errors

Figure 4: The global correlation matrix used in our joint likelihood analysis
of the 2PCF and 3PCF, estimated from the 153 VIPERS mock catalogues using Eq.
7. We show here the one for the G1 sample, as an illustrative example. The
matrix is composed by two sub-blocks, corresponding to the 2PCF monopole and
quadrupole ($C_{\xi,\xi}$, bottom left), and the 3PCF monopole
($C_{\zeta,\zeta}$, top right), respectively, together with their cross-
covariance ($C_{\xi,\zeta}$). The sub-arrays are used for the separate 2PCF
and 3PCF analyses, while the full $45\times 45$ matrix is used in the joint
likelihood analysis.

To estimate the errors, we will not adopt an analytic, Gaussian model (like
in, e.g., Slepian & Eisenstein 2017). Instead, supported by the results of the
previous section, we directly estimate errors and their covariance matrix
$C_{i,j}$ from the mock VIPERS catalogues, as

$C_{i,j}=\frac{1}{N_{m}-1}\sum_{k=0}^{N_{m}}\left(d_{i}^{k}-\overline{d}_{i}\right)\left(d_{j}^{k}-\overline{d}_{j}\right)\,,$
(7)

where $d_{i}^{k}$ indicates the data vector in the $i-{th}$ bin and in the
$k-{th}$ mock catalog. The size of the data vector and of the covariance
matrix depend on the specific clustering analysis being performed. In the 2PCF
analysis, performed in the range $15\leq s_{min}^{2pt}\leq
40\,h^{-1}\,\mathrm{Mpc}$, the size of the data vector is 10 and the
corresponding $10\times 10$ covariance matrix is shown in the lower left part
of Fig. 4. For the 3PCF analysis, performed over the same range of scales, the
size of the data vector is 35 and the covariance matrix corresponds to the
upper-right part of the same figure. In the joint 2PCF and 3PCF analysis, the
data vector has 45 elements and the covariance matrix, which also contains the
2PCF-3PCF cross-terms, corresponds to the full matrix of Fig. 4.

Both precision and accuracy of our covariance matrix depend on the number of
mock catalogues used for its numerical estimation. Since this number is
limited, a correction for the expected systematic errors should be considered
(see, e.g. discussion in Hartlap et al., 2007; Percival et al., 2014). To this
end, we follow Sellentin & Heavens (2016) in using a modified likelihood
function, which we discuss in detail in Sect. 6.

##  5 Modelling the galaxy 2- and 3- point correlation functions

In this section we describe our model for the anisotropic galaxy 2PCF,

$\displaystyle\xi_{g}(s,\mu)\,$ $\displaystyle=$
$\displaystyle\left<\delta_{g}(\mathbf{x}))\delta_{g}(\mathbf{x}+\mathbf{s})\right>,$
(8)

and 3PCF,

$\displaystyle\zeta_{g}(s_{12},s_{13},s_{23})$ $\displaystyle=$
$\displaystyle\left<\delta_{g}(\mathbf{x})\delta_{g}(\mathbf{x}+\mathbf{s_{12}})\delta_{g}(\mathbf{x}+\mathbf{s_{13}})\right>\,.$
(9)

In these expressions, $\delta_{g}(\mathbf{x})$ is the galaxy density contrast
at the position $\mathbf{x}$, while $<...>$ indicates ensemble average.

The model relies on the approach of Scoccimarro et al. (1999); Scoccimarro
(2004); McDonald & Roy (2009); Saito et al. (2014), and has already been used
to perform 2PCF, 3PCF and joint 2PCF+3PCF analyses (Gil-Marín et al., 2015; de
la Torre et al., 2017; Slepian & Eisenstein, 2017).

###  5.1 Galaxy bias

The first ingredient of the model is the galaxy biasing relation of McDonald &
Roy (2009)

$\displaystyle\delta_{g}=b_{1}\delta+\frac{b_{2}}{2}\left[\delta^{2}-<\delta^{2}>\right]+\frac{b_{s^{2}}}{2}\left[s^{2}-<s^{2}>\right]+O\left(\delta^{3}\right).\,\,$
(10)

In this expression, the galaxy density contrast $\delta_{g}$ is a function of
the matter non-linear over-density $\delta$ and of the non-local tidal field
$s^{2}$, which accounts for dependence of $\delta_{g}$ on the local potential
$\phi$ (see e.g. Sect. II of McDonald & Roy, 2009). The terms $<\delta^{2}>$
and $<s^{2}>$ guarantee that $<\delta_{g}>=0$. This bias model is then
characterised by the four parameters $b_{1}$, $b_{2}$, $b_{s^{2}},b_{3nl}$.

###  5.2 Distortion effects

Using the observed spectroscopic redshifts as a distance proxy introduces
spurious distortions in the clustering properties that need to be accounted
for. Three types of distortions need to be considered: those introduced by the
departures from the Hubble flow, those induced by assuming an incorrect
cosmological model in computing distances from redshifts and those due to the
errors on the measured redshifts.

####  5.2.1 Peculiar motions

Peculiar velocities introduce a Doppler shift that modifies the mapping
between spatial positions, $\mathbf{x}$, and the measured positions,
$\mathbf{s}$:

$\mathbf{s}=\mathbf{x}+\frac{v_{\parallel}}{aH(a)}\mathbf{e}_{\parallel}\,$
(11)

where $a$ is the expansion factor, $v_{\parallel}$ is the radial component of
the peculiar velocity vector, $\mathbf{e}_{\parallel}$ is the unit radial
vector, and $\mathbf{s}$ is the redshift-space vector position of the object.

This mapping introduces a relation between the measured (redshift-space)
galaxy over-density, $\delta_{g}(\mathbf{s})$, and the true, real-space one,
$\delta_{g}(\mathbf{x})$

$\delta_{g}(\mathbf{s})=\left[1+\delta_{g}(\mathbf{x})\right]\left|\frac{\mathrm{d}^{3}s}{\mathrm{d}^{3}x}\right|^{-1}-1\,,$
(12)

where $\left|\frac{\mathrm{d}^{3}s}{\mathrm{d}^{3}x}\right|^{-1}$ is the
Jacobian of the map in Eq. (11).

In this work, we make the plane-parallel hypothesis, i.e. we assume that the
relative separations between galaxy pairs are much smaller than the distance
to the observer. This assumption is fully justified, since the maximum scale
considered in our analyses ($40\,h^{-1}\,\mathrm{Mpc}$) is much smaller than
any distance to a VIPERS galaxy, which lies at $z\geq 0.5$. Under this
approximation, Eq. (12) becomes

$\delta_{g}(\mathbf{s})=\frac{\delta_{g}(\mathbf{x})+f\partial_{\parallel}\mathbf{u}}{1-f\partial_{\parallel}\mathbf{u}}\,,$
(13)

where the partial derivative is taken along the radial direction, $\mathbf{u}$
is the peculiar velocity, and $f$ is the linear growth rate of density
fluctuations,

$f=\frac{\mathrm{d}\log D}{\mathrm{d}\log a}\approx\Omega_{M}^{0.545}(z)\,,$
(14)

where $D(z)$ is the linear growth factor and the second relation is a good
approximation in the flat $\Lambda$CDM model assumed here (see e.g. Wang &
Steinhardt, 1998; Huterer & Linder, 2007).

####  5.2.2 Geometrical distortions

Choosing an incorrect fiducial cosmology to estimate distances from redshifts
generates another kind of anisotropy, the well-known Alcock-Paczynski effect
(Alcock & Paczynski, 1979). The detection of this effect on intrinsically
isotropic clustering features as the BAO peak in the 2PCF has been used to
trace the expansion history of the Universe (Eisenstein et al., 2005; Kazin et
al., 2013). Here, we will treat this effect as a nuisance to be marginalised
over, since our analysis is limited to scales well below the BAO peak. As
shown in de la Torre et al. (2017), the impact of this marginalisation in the
error budget of the VIPERS clustering analysis is negligible. Therefore, we
will set our fiducial cosmology equal to the flat $\Lambda$CDM model specified
in Sect. 1 and safely ignore the Alcock-Paczynski effect.

####  5.2.3 Redshift measurement errors

Redshift measurement errors introduce a noise in the redshift-to-distance
relation analogous to that of incoherent motions (Marulli et al., 2012),
effectively erasing information on scales below

$\sigma_{\pi}=\frac{c\sigma_{z}}{H(z)}\,h^{-1}\,\mathrm{Mpc}\,,$ (15)

where $\sigma_{z}$ is the $rms$ error typical of the adopted instrumental
setup. For VIPERS, the errors are very well described by a Gaussian
distribution with $\sigma_{z}=5\times 10^{-4}(1+z)$, corresponding to a length
scale $\sigma_{\pi}\sim 1.5\,h^{-1}\,\mathrm{Mpc}$ (Scodeggio et al., 2018;
Sereno et al., 2015).

###  5.3 2PCF model

To model the galaxy 2PCF we start from the matter power spectrum $P(k)$

$\left<\delta(\mathbf{k})\delta(\mathbf{k^{\prime}})\right>=2\pi^{3}\delta_{D}(\mathbf{k}-\mathbf{k^{\prime}})P(|\mathbf{k}|)\,.$
(16)

Using Eq. (13) to account for redshift space distortions, we obtain an
expression of the redshift-space power spectrum

$\displaystyle P^{s}(k,\nu)$ $\displaystyle=$ $\displaystyle\int
d^{3}\mathbf{r}e^{-i\mathbf{k}\cdot\mathbf{r}}\left<e^{-ifk\nu\Delta
u_{\parallel}}\right.$
$\displaystyle\left.\left[\delta(\mathbf{x})+f\partial_{\parallel}u_{\parallel}(\mathbf{x})\right]\left[\delta(\mathbf{x}^{\prime})+f\partial_{\parallel}u_{\parallel}(\mathbf{x}^{\prime})\right]\right>\,,$

where $\nu=k_{\parallel}/k$, $u_{\parallel}=-v_{\parallel}/(faH(a))$,
$v_{\parallel}$ and $k_{\parallel}$ are the radial components of the peculiar
velocity and wavenumber vectors, respectively, $\delta$ is the matter over-
density,
$u_{\parallel}=u_{\parallel}(\mathbf{x})-u_{\parallel}(\mathbf{x}^{\prime})$
and $\mathbf{r}=\mathbf{x}-\mathbf{x}^{\prime}$. The term in square brackets
accounts for the effect of coherent motions that increase the clustering
amplitude, whereas the exponential factor encodes the “Fingers of God" effect
of incoherent motions. For a more detailed description, see Taruya et al.
(2010).

Equation (5.3) is exact but of impractical use. Rather, using the
approximations introduced by Scoccimarro (2004) and Taruya et al. (2010) for
the case of biased mass tracers, we obtain the simpler expression

$P_{g}^{s}(k,\nu)=D(\sigma_{12},\sigma_{z})\left[P_{\rm{gg}}(k)+2\nu^{2}fP_{\rm{g\theta}}(k)+\nu^{4}f^{2}P_{\rm{\theta\theta}}(k)\right],$
(18)

where

$\displaystyle P_{gg}(k)=$ $\displaystyle
b_{1}^{2}P_{\delta\delta}(k)+2b_{2}b_{1}P_{b2,\delta}(k)+2b_{s^{2}}b_{1}P_{bs2,\delta}(k)$
(19)
$\displaystyle+b_{2}^{2}P_{b22}(k)+2b_{2}b_{s^{2}}P_{b2s2}(k)+b_{s^{2}}^{2}P_{bs22}(k)$
$\displaystyle+2b_{1}b_{3\rm nl}\sigma_{3}^{2}(k)P_{\rm lin}(k),$
$\displaystyle P_{g\theta}(k)=$ $\displaystyle
b_{1}P_{\delta\theta}(k)+b_{2}P_{b2,\theta}(k)+b_{s^{2}}P_{bs2,\theta}(k)$
(20) $\displaystyle+b_{3\rm nl}\sigma_{3}^{2}(k)P_{\rm lin}(k)\,.$

In these equations, $P_{\rm{\delta\delta}}$, $P_{\rm{\delta\theta}}$ and
$P_{\rm{\theta\theta}}$ are, respectively, the nonlinear matter density-
density, density-(velocity divergence), and (velocity divergence)-(velocity
divergence) power spectra. $P_{\rm lin}$ is the matter linear power spectrum,
and $P_{b2,\delta}$, $P_{bs2,\delta}$, $P_{b2,\theta}$, $P_{bs2,\theta}$,
$P_{b22}$, $P_{b2s2}$, $P_{bs22}$, $\sigma_{3}^{2}$ are 1-loop bias integrals;
their expressions can be found in e.g. Gil-Marín et al. (2015); de la Torre et
al. (2017). The linear and nonlinear power spectra of density and velocity
divergence are specified at the effective redshift of the sample. To model
$P_{\rm lin}$ we use the CAMB Boltzmann solver code (Lewis et al., 2000),
whereas for $P_{\delta\delta}(k),P_{\delta\theta},P_{\theta\theta}$ we use the
SPT 1-loop model of Scoccimarro (2004, Eqs. 63-65). In this model, the degree
of non-linearity in $P_{\rm{\delta\delta}}$, $P_{\rm{\delta\theta}}$ and
$P_{\rm{\theta\theta}}$ is quantified by the same parameter that measures the
clustering amplitude, i.e. $\sigma_{8}(z)$.

In Eq. (18) we have introduced the term $D(\sigma_{12},\sigma_{z})$ to model
the combined damping effect of incoherent motions and redshift errors. Its
explicit expression is

$D(\sigma_{12},\sigma_{z})=\left(1+k^{2}\nu^{2}\sigma^{2}_{12}\right)^{-1}\cdot\exp\left(-k^{2}\nu^{2}\sigma^{2}_{z}\right)\,,$
(21)

The first term, modelled as a Lorentzian damping, accounts for the effect of
random motions within dark matter halos, quantified by the pairwise velocity
dispersion $\sigma_{12}$. The second term accounts for the redshift
measurement errors and is modelled as a Gaussian function. In our analysis, we
fix $\sigma_{z}$ to the VIPERS estimated value and leave $\sigma_{12}$ as a
free factor we can fit for and then treat it as a nuisance parameter.

Since the 2PCF model fits well the corresponding measurements from the VIPERS
mocks on all scales considered, we decided to ignore the model correction
terms $C_{A}(k,\nu,f,b_{1})$ and $C_{B}(k,\nu,f,b_{1})$ of Taruya et al.
(2010), that are instead considered in other analyses (e.g. Gil-Marín et al.,
2015).

To test the robustness of our model, we also considered an alternative 2PCF
model, in which we used (1) the HALOFIT semi-analytical prescription
calibrated by Takahashi et al. (2012) to model the matter power spectra and,
(2) the fitting functions of Bel et al. (2019) to model $P_{\delta,\theta}(k)$
and $P_{\theta,\theta}(k)$. We find that the alternative model provide results
very similar to the original one, the differences being much smaller than the
statistical errors.

We then use Eq. (18) to model the multipole moments of the anisotropic power
spectrum

$P^{s}_{\ell}(k)=\frac{2\ell+1}{2}\int_{-1}^{1}P_{g}^{s}(k,\nu)\mathcal{P}_{l}(\nu)\,\mathrm{d}\nu$
(22)

and, from these, the multipole moments of the anisotropic 2PCF

$\xi^{s}_{\ell}(s)=i^{\ell}\int\frac{k^{2}}{2\pi^{2}}P^{s}_{\ell}(k)j_{\ell}(ks)\,\mathrm{d}k,$
(23)

where $j_{\ell}$ indicates the spherical Bessel functions. Having neglected
the $C_{A}(k,\nu,f,b_{1})$ terms, in our model the parameters $b_{1}$ and $f$
are fully degenerate with $\sigma_{8}$. The two bias parameters $b_{1}$ and
$b_{2}$ are also highly degenerate.

###  5.4 3PCF model

The 3PCF model we adopt is the same one used by Slepian & Eisenstein (2017)
and Slepian et al. (2017a) to detect the BAO feature in the 3PCF of SDSS DR12
galaxies and to measure their biasing relation in redshift space.

Our 3PCF model is derived from galaxy bispectrum
$B(|\mathbf{k_{1}}|,|\mathbf{k_{2}}|,|\mathbf{k_{3}}|)$, defined as

$\left<\delta_{g}(\mathbf{k_{1}})\delta_{g}(\mathbf{k_{2}})\delta_{g}(\mathbf{k_{3}})\right>=\delta_{D}(\mathbf{k_{1}}+\mathbf{k_{2}}+\mathbf{k_{3}})B(|\mathbf{k_{1}}|,|\mathbf{k_{2}}|,|\mathbf{k_{3}}|).$
(24)

In particular, we assume the redshift-space tree-level galaxy bispectrum
derived in Scoccimarro et al. (1999). From this, we obtain the 3PCF model by
Fourier transform, following the prescriptions of Slepian et al. (2017b).

The 3PCF model can be expressed as

$\displaystyle\zeta(s_{12},s_{13},s_{23})=$
$\displaystyle\sum_{l=0}^{l=4}A_{l}(b_{1},\gamma,\gamma^{\prime},\beta)\cdot
f_{l}\left(s_{12},s_{13},s_{23}\right)+$ (25) $\displaystyle
B(b_{1},\beta)\cdot\sum_{l=0}^{l_{max}}k_{l}\left(s_{12},s_{13},s_{23}\right)\,,$

where $\gamma\equiv{b_{2}}/{b_{1}}$ and $\gamma^{\prime}\equiv
b_{s^{2}}/b_{1}$ are combinations of the bias parameters. The first term of
Eq. (25) accounts for contributions separable in $k_{1},k_{2}$ in the
Scoccimarro et al. (1999) bispectrum model. Its Fourier transform can be
computed analytically as the product of two quantities. The first one,
$A_{l}$, depends on the bias parameters, $b_{1},b_{2}$, $b_{s^{2}}$, and on
the linear distortion parameter $\beta=f/b_{1}$. The second one, $f_{l}$,
depends on the linear power spectrum, $P_{\rm lin}(k)$. The explicit
expression of the $A_{l}$ and $f_{l}$ terms is given in Appendix B.

The second term in Eq. (25) accounts for the non separable terms (i.e. the
terms in which $k_{3}$ appears explicitly). It is also made-up of two terms.
The multiplicative factor $B=b_{1}^{3}(7\beta^{2}+2\beta^{3})$ depends on
$b_{1}$ and $\beta$ while the terms $k_{l}\left(s_{12},s_{13},s_{23}\right)$
in the sum depend on $P_{\rm lin}(k)$. Their explicit expressions are given in
Equations 17, 18 and 19 of Slepian et al. (2017b). The contribution of non-
separable terms is small compared to that of the separable one and therefore
we will ignore them in the 3PCF model, similarly to what we did with the
$C_{A}$ and $C_{B}$ terms in the 2PCF model.

Scoccimarro et al. (1999); Gil-Marín et al. (2015) have modeled the effect of
incoherent motions with a damping term in the galaxy bispectrum

$D^{B}(\alpha\sigma_{12})=\left\\{1+\alpha^{2}\left[\left(k_{1}\nu_{1}\right)^{2}+\left(k_{2}\nu_{2}\right)^{2}+\left(k_{3}\nu_{3}\right)^{2}\right]^{2}\frac{\sigma_{12}^{2}}{2}\right\\}^{-2}.$
(26)

The magnitude of the damping effect is determined by the pairwise velocity
dispersion, $\sigma_{12}$ and modulated, in the model, by a multiplicative
factor, $\alpha$ (with $\alpha=0$ describing a purely coherent flow). In this
work, however, we decided to drop this damping term since it has a negligible
impact on scales larger than $10\,h^{-1}\,\mathrm{Mpc}$ as we show in a
dedicated test in Sect. 7.2, where we assess the sensitivity of our results to
the choice of $s_{min}$. This choice has the further advantage of avoiding
computationally demanding 2-dimensional Fourier transforms.

###  5.5 Integral constraint

The so-called integral constraint effect originates from the assumption that
the mean density estimated from the data coincides with the true one. This
assumption biases the estimate of the correlation function on scales
comparable to the size of the survey (Hui & Gaztañaga, 1999), where it cannot
be neglected, in particular at BAO scales (Slepian & Eisenstein, 2017).
However, the maximum scales considered in our analysis (40
$h^{-1}\,\mathrm{Mpc}$) are significantly smaller than the size of the survey.
This guarantees that the integral constraint is negligible, as shown by Cappi
et al. (2015), and will be ignored in our analysis.

###  5.6 Summary of model parameters

Our 2PCF and 3PCF models are specified by the set of free parameters listed in
Table 2. These include: the clustering rms amplitude, $\sigma_{8}$, the linear
redshift distortion parameter $\beta\equiv f/b_{1}$, the pairwise velocity
dispersion, $\sigma_{12}$, and the two galaxy bias parameters, $b_{1}$ and
$\gamma\equiv b_{2}/b_{1}$.

The other model parameters that are kept fixed to their fiducial values are:
the shape parameters of $P_{lin}(k)$ which are set to the Planck values
(Planck Collaboration et al., 2020); the VIPERS rms redshift error,
$\sigma_{z}=5\times 10^{-4}(1+z)$; the bias parameters $b_{s^{2}}$ and
$b_{3nl}$, which are derived from local Lagrangian theory, i.e.

$\displaystyle b_{s^{2}}$ $\displaystyle=-\frac{4}{7}\left(b_{1}-1\right)\,,$
(27) $\displaystyle b_{3nl}$
$\displaystyle=\frac{32}{325}\left(b_{1}-1\right)\,.$ (28)

This choice is motivated by the fact that their measured values would be too
noisy, as we have verified by running a dedicated likelihood analysis, aimed
at testing the local Lagrangian bias hypothesis. Finally, we set $\alpha=0$; a
choice that has a negligible impact on the results, as we show in Sect. 7.2.

##  6 Parameter Inference

To estimate the free parameters of the model, described as a vector
$\boldsymbol{\theta}$, and their uncertainties, we compute the posterior
probability, $P\left(\boldsymbol{\mu}(\boldsymbol{\theta})|\mathbf{d}\right)$,
where $\boldsymbol{\mu}(\boldsymbol{\theta})$ is the model prediction vector
and $\mathbf{d}$ is the data vector.

We perform three different analyses: one in which the data vector is
represented by the VIPERS 2PCF monopole and quadrupole measured in the two
samples G1 and G2 in the range $[15,40]\,h^{-1}\,\mathrm{Mpc}$; one where we
consider the measured 3PCF monopole of the same galaxies in the same range of
scales; and one considering both statistics. In each case, we use the
corresponding numerical covariance matrix $C$ as described in Sect. 4.

To evaluate the posterior probability, we multiply the likelihood,
$\mathcal{L}(\mathbf{d}|\boldsymbol{\mu},C)$, by the prior probability,
$P\left(\boldsymbol{\theta}\right)$. For each free parameter, the prior is
modelled as a step function with zero value outside the ranges specified in
Table 2. To compute the likelihood
$\mathcal{L}(\mathbf{d}|\boldsymbol{\mu},C)$, we use the covariance matrices
described in Sect. 4. The latter have been numerically estimated from a large
(153) but finite number of mock VIPERS catalogues, which introduces a
systematic error. To correct for this, we adopt the approach of Sellentin &
Heavens (2016) and use the modified likelihood function

$\mathcal{L}\left(\mathbf{d}|\boldsymbol{\mu},C,N_{m}\right)=\bar{c}_{\mathrm{p}}\left(\det{C}\right)^{-1/2}\left[1+\frac{\chi^{2}(\boldsymbol{\mu},C)}{N_{m}-1}\right]^{-\frac{N_{m}}{2}}\,,$
(29)

where

$\chi^{2}(\boldsymbol{\mu},C)=\left(\mathbf{d}-\boldsymbol{\mu}\right)^{T}C^{-1}\left(\mathbf{d}-\boldsymbol{\mu}\right)\,$
(30)

$\bar{c}_{\mathrm{p}}=\frac{\Gamma\left(\frac{N_{m}}{2}\right)}{\pi\left(N_{m}-1\right)^{p/2}\Gamma\left(\frac{N_{m}-\ell}{2}\right)}\,,$
(31)

with $N_{m}$ being the number of mock catalogues, and $\Gamma$ the Gamma
function.

To sample the posterior
$P\left(\mathbf{\mu}(\mathbf{\theta})|\mathbf{d}\right)$ we use a Monte Carlo
Markov Chain (MCMC) approach that generates a set of chains, i.e. collections
of points in the parameter space
$\\{\mathbf{\theta}_{1},\mathbf{\theta}_{2},...,\mathbf{\theta}_{N}\\}$. Some
of the parameters used in our analysis are expected to be (partly) degenerate.
In these cases, we do not consider the individual parameters, but a suitable
combination, as e.g. for $b_{1}\sigma_{8}$ and $f\sigma_{8}$ in the case of
the 2PCF analysis alone.

Table 2: Priors for the free parameters used in our likelihood analyses. The pairwise velocity dispersion, $\sigma_{12}$ is expressed in $h^{-1}\,\mathrm{Mpc}$. Analysis Parameters | $b_{1}$ | $\gamma$ | $\beta$ | $\sigma_{8}(z)$ | $\sigma_{12}$
---|---|---|---|---|---
2PCF | [0.1, 5] | [-5., 5] | [0.1, 5] | [0, 5] | [0, 10]
3PCF | [0.1, 5] | [-5, 5] | [0.1, 5] | [0.1, 5] | –
2PCF+3PCF | [0.1, 5] | [-5, 5] | [0.1, 5] | [0, 5] | [0, 10]
|
---|---
Figure 5:  Marginalized 2D posterior probability levels for the parameters $b_{1}\sigma_{8}$, $f\sigma_{8}$ and $\sigma_{12}$, from the 2PCF analysis of the P1 and P2 “control" samples (left/right respectively), as a function of the minimum scale included in the analysis (see legend). The units of $\sigma_{12}$ are $h^{-1}\,\mathrm{Mpc}$ and the darker/lighter shades for each set correspond to 68% and 95% confidence levels, respectively. 1D distribution functions on the diagonal give the 1D marginalized probability distributions of each parameter, with the 68% interval indicated by the blue-shaded area for the $s_{min}=15\,h^{-1}\,\mathrm{Mpc}$ reference case only. The horizontal and vertical dashed lines give the best-fit values obtained by Pezzotta et al. (2017) on similar sub-samples of the VIPERS PDR2 catalogue.  |
---|---
Figure 6: Marginalized 2D and 1D posterior distributions for the parameters
$b_{1}\sigma_{8}$, $\beta$ and $\gamma$, obtained from the 3PCF analysis of
the G1 and G2 VIPERS samples (left and right panels, respectively). The plots
show the results for the two choices of $s_{min}$ as indicated by the legend.
As in previous figures, darker/lighter shades correspond to 68 % and 95 %
confidence levels, respectively. The green dashed lines report the
marginalized 1D distributions obtained for the same parameters from the 2PCF
analysis, using $s_{min}=15\,h^{-1}\,\mathrm{Mpc}$.

##  7 Results

In this section, we present the results of three analyses based, respectively,
on the VIPERS 2PCF only, 3PCF only, and on the combination of the two
statistics.

The first analysis is very similar to those of de la Torre et al. (2013,
2017); Pezzotta et al. (2017); Mohammad et al. (2018) and we use it to
validate our likelihood pipeline. The main goal of the second analysis, which
builds upon the work of Moresco et al. (2017), is to measure the bias
parameters $b_{1}$ and $\gamma\equiv b_{2}/b_{1}$. Finally, the joint 2PCF and
3PCF analysis is aimed at breaking parameter degeneracies to obtain
independent constraints on the galaxy bias, the evolution of cosmic structures
and the clustering amplitude.

For all such analyses, the G1 and G2 VIPERS sub-samples are used (see Table
1). For the 2PCF study only, we also consider the P1 and P2 samples to allow
us a direct comparison to Pezzotta et al. (2017).

###  7.1 2PCF

The results of the 2PCF analysis are shown in Fig. 5. In the panels we plot
the 2D and 1D posterior probability distributions for the combinations of
degenerate parameters $f\sigma_{8}$, $b_{1}\sigma_{8}$ and for $\sigma_{12}$,
obtained from the analysis of the VIPERS P1 ($\bar{z}=0.61$) and P2
($\bar{z}=0.87$) samples and for different $s_{min}$ values indicated in the
plots. We considered the 2PCF monopole and quadrupole moments and used the
$10\times 10$ covariance sub-matrix $C_{\xi,\xi}$ of Fig. 4.

Our analysis is very similar, though not identical, to that of Pezzotta et al.
(2017). To minimize the differences, we compare our results obtained by using
an SPT 1-loop 2PCF model in redshift space with their “TNS" model (Taruya et
al., 2010). For the same reason, we push the analysis down to
$s_{min}=5\,h^{-1}\,\mathrm{Mpc}$, but then show the results obtained with
larger values up to $s_{min}=15\,h^{-1}\,\mathrm{Mpc}$ (our reference case).

The 2D and 1D posterior distributions shown in Fig. 4 agree well with the
results of Pezzotta et al. (2017) except for the pairwise velocity dispersion,
$\sigma_{12}$, which is significantly ($2\sigma$) smaller. This is not
surprising since, as shown in Fig. 18 of Pezzotta et al. (2017), the
Scoccimarro (2004) 2PCF model tends to underestimate the value of
$\sigma_{12}$.

The 2D contours in Fig. 5 show a mild anti-correlation between $f\sigma_{8}$
and $b_{1}\sigma_{8}$, which is expected since larger redshift distortions
that boost up the 2PCF monopole can be compensated by reducing the clustering
amplitude $b_{1}\sigma_{8}$. In addition, we find a positive correlation
between $f\sigma_{8}$ and $\sigma_{12}$, which is explained by the fact that
damping effects can be compensated by increasing the clustering amplitude.

To assess the sensitivity of our results to some of the model parameters, we
have performed two robustness tests. In the first one we have repeated the
likelihood analysis with different values of
$s_{min}=5,\,10\,\mbox{and}\,15\,h^{-1}\,\mathrm{Mpc}$. The results are shown
in the same Fig. 5. As expected, increasing the value of $s_{min}$ amplifies
the errors, especially for the nonlinear $\sigma_{12}$ parameter. However, no
systematic error is introduced as all results agree with those of Pezzotta et
al. (2017) and with each other.

As a second test we used the alternative 2PCF model introduced in Sect. 5.3 in
which we 1) use HALOFIT (Takahashi et al., 2012) to model $P_{\delta\delta}$
and 2) adopt the fitting formula of Bel et al. (2019) to model
$P_{\delta\theta}$ and $P_{\theta\theta}$. The results are very similar to
those obtained with the reference 2PFC model.

Finally, we have repeated the likelihood analysis on G1 and G2 samples that we
use for the 3PCF analysis too. The results, which are qualitatively similar to
those obtained with P1 and P2, are summarized in Table 3 where we list the
best-fit parameters and their $1\sigma$ uncertainties.

We confirm that the 2PCF analysis successfully constrains the parameter
combinations $f\sigma_{8}$, $b_{1}\sigma_{8}$ and $\beta\equiv f/b_{1}$ but
leave $\gamma$ unconstrained. This is not surprising since this parameter is
sensitive to the small scale clustering that we ignore, having set
$s_{min}=15\,h^{-1}\,\mathrm{Mpc}$. Indeed, when the analysis is extended down
to $s_{min}=5\,h^{-1}\,\mathrm{Mpc}$, we do measure a $\gamma$ value
significantly different from zero, yet mildly inconsistent with the results of
previous VIPERS analyses (Cappi et al., 2015; Di Porto et al., 2016). This
value is also in tension with the one yielded by the 3PCF analysis presented
in the next section. It probably indicates the minimum scales below which our
galaxy 2PCF model fails and thus justifies our choice of setting
$s_{min}=15\,h^{-1}\,\mathrm{Mpc}$.

###  7.2 3PCF

We have repeated the likelihood analysis to compare the 3PCF measurements in
the G1 and G2 sample with the 3PCF model described in Sect. 5.4 and the
$25\times 25$ covariance matrix $C_{\zeta,\zeta}$ shown in Fig. 4.

This analysis improves upon the one performed by Moresco et al. (2017) in
different ways. First, we use a more accurate 3PCF model and a covariance
matrix estimated from realistic mock catalogues to perform a full likelihood
analysis. Secondly, here we use the final, PDR2 release of the VIPERS
catalogue. Finally, we consider all triangle configurations and not just a few
subsets. We stress the importance of this latter aspect since, in our
estimator, we use a triangle representation that is flexible enough to allow
us to select subsets of triangles characterized by specific configurations and
side lengths, greatly simplifying the comparison with model predictions.

The plots in Fig. 6 are analogous to those of Fig. 5. They show the 2D and 1D
marginalized posterior distributions for the parameters $b_{1}\sigma_{8}$,
$\beta$ and $\gamma$ obtained from the 3PCF measured in the G1 (left) and G2
(right) samples. The best-fitting parameters and their uncertainties are
listed in Table 3. The 3PCF analysis allows us to estimate the nonlinear bias
parameter $\gamma$, which turns out to be consistent with zero, while 2PCF is
completely insensitive to it, as shown by the green, dot-dashed curves in the
bottom right panels. This result is in agreement with Moresco et al. (2017)
and marginally consistent with those obtained from counts-in-cells analyses of
Cappi et al. (2015) and Di Porto et al. (2016) where a non-zero value of
$\gamma$ was measured but on scales much smaller than those considered here.

The constraints on the distortion parameter $\beta$ are, on the contrary,
quite weak. The best-fit value of this parameter is consistent with zero
within the errors. This is not unexpected since in our analysis we considered
only the 3PCF monopole moment, which is insensitive to the RSD effects.
Constraints on $\beta$ could be obtained by considering the 3PCF multipoles.
This is, however, beyond the scope of this work since a multipole analysis
would require the use of a significantly larger covariance matrix, which
cannot be accurately estimated with the 153 mock catalogues at our disposal.

Another interesting feature is the mild degeneracy between $\beta$ and the
clustering amplitude $b_{1}\sigma_{8}$ also seen for the 2PCF. The difference
is that here the error on $b_{1}\sigma_{8}$ from the 3PCF analysis is about
$3$ times larger than from the 2PCF one.

To test the sensitivity to the choice of $s_{min}$ we repeated the likelihood
analysis using $s_{min}=5,\,10,\,15\,\mbox{and}\,20\,h^{-1}\,\mathrm{Mpc}$. In
Figure 6 we only show the cases of $s_{min}=10$ and 15 $h^{-1}\,\mathrm{Mpc}$
to avoid overcrowding. When $s_{min}$ decreases, the maximum of the 1D
posterior distributions for $b\sigma_{8}$ shifts to smaller values. However,
due to the large errors, the best-fit values of $b\sigma_{8}$ (but also the
ones of $\gamma$ and $\beta$) are in fact consistent with each other for all
choices of $s_{min}$. In Appendix C we present a more detailed study of the
impact of $s_{min}$ on the parameters’ estimates using the VIPERS mock
catalogues.

Finally, despite their different shapes, the 1D posteriors of all parameters
obtained from the 3PCF analysis for all values of $s_{min}$ largely overlap
those obtained from the 2PCF.

To further assess the impact of nonlinear effects, we compared the measured
3PCF with a model including the damping term of Eq. (26) and varying its
strength. The result is shown in Fig. 7, with the reference case with
$\alpha=0$ corresponding to no damping. We set the pairwise velocity
dispersion, $\sigma_{12}$, equal to the best fit value obtained from the 2PCF
analysis. The plot compares model predictions with the 3PCF measured for
various triangle shapes, whose side lengths are indicated on the x-axis. The
results show that incoherent motions do not significantly affect model
predictions on the scales considered in our analysis. They may, however,
become relevant for next-generation surveys, in which statistical errors will
be significantly smaller than those considered here. We plan to further
investigate this aspect and introduce a novel nonlinear 3PCF model to tackle
the problem.

Figure 7: The impact of small-scale damping on our model for the 3PCF
monopole. The estimated 3PCF from the VIPERS G1 and G2 samples (top and bottom
respectively, data and error bars as in Fig. 3) are compared to the
theoretical predictions, varying the amplitude of nonlinear motions, through
the parameter $\alpha$. The case $\alpha=0$ corresponds to the reference case,
in which small-scales nonlinear motions are ignored.

###  7.3 Joint analysis of 2PCF and 3PCF

|
---|---
Figure 8: Marginalized 1D and 2D posterior distributions for the parameters
$\sigma_{8}(z)$, $f(z)$ and $b_{1}$, from the joint 2PCF+3PCF analysis of the
G1 and G2 VIPERS samples (left and right panels, respectively). All estimates
use $s_{min}=15\,h^{-1}\,\mathrm{Mpc}$. The meaning of contours and shaded
areas are as in Fig. 6. The vertical and horizontal black dashed lines show
for reference the values measured for the same parameters by Planck
Collaboration et al. (2020), self-consistently scaled to the corresponding
VIPERS redshift.

Finally, we perform for the first time a joint 2PCF and 3PCF correlation
analysis on the VIPERS data. The main goal is to break the parameter
degeneracy and obtain individual estimates for $\sigma_{8}$, for the linear
growth rate $f$, and for the two bias parameters $b_{1}$ and $\gamma$. Fig. 8
shows the 1D and 2D posterior probability distributions (blue regions and
solid curves) of these parameters from the likelihood analysis of the G1 and
G2 samples and using the full $45\times 45$ covariance matrix. The same range
of scales $[15,40]\,\mathrm{Mpc}\,h^{-1}$ was considered in both the 2- and
3-point statistics. The 1D posterior distributions obtained from the 2PCF-only
(green, dot-dashed) and 3PCF-only (red, dashed) analyses are also shown for
reference. The best-fit values of the joint analysis and their uncertainties
are listed in Table 3.

The joint analysis successfully breaks the $\sigma_{8}-b_{1}-f$ degeneracy,
and the parameters can be measured individually, although with different
uncertainties. The value of $\sigma_{8}$ is measured with a $\sim 20$ %
uncertainty in both G1 and G2 samples. The error on $b_{1}$ is $\sim 20$ % in
the G1 sample, increasing to $\sim 30$% in G2. The error on the growth rate
$f$ is significantly larger ($65-75$%), reflecting the fact that the
$f-\sigma_{8}$ degeneracy is only partially broken by our analysis. Though
significantly less precise, it is reassuring that our best-fit $\sigma_{8}$
and $f$ values are in agreement with those of the Planck Collaboration et al.
(2020) cosmology, scaled to the redshift of the survey (vertical and
horizontal dashed lines). The sensitivity of these results to the choice of
$s_{min}$ is thoroughly discussed in Appendix C.

The joint analysis does not improve the estimate of $\gamma$, whose precision
is driven by the 3PCF signal, as shown in Fig. 8. It does instead improve the
estimate of $\sigma_{12}$ whose value, in agreement with the one obtained from
the 2PCF analysis, is estimated with a significantly higher precision.

Overall, the results of the joint analysis clearly show that the degeneracy of
some parameters are successfully broken, thanks to the ability to extract
information from intermediate to small scales where nonlinear effects are
relevant, both in the evolution of the density fluctuation and in the biasing
relation.

Table 3: Summary of the best-fit parameter values and their 1$\sigma$ uncertainties, as obtained from the G1/G2 VIPERS samples using the 2PCF-only, 3PCF-only and joint 2PCF+3PCF analyses. The parameters common to all analyses are $b_{1}$, $\gamma$, $\beta$, $\sigma_{8}$. For the 2PCF and joint analyses, the pairwise velocity dispersion $\sigma_{12}$ is also considered. To account for parameter degeneracy, we also list some of the relevant parameter combinations constrained by the individual 2PCF and 3PCF analyses.  Sample | Probe | $\sigma_{8}(z)$ | $b_{1}$ | $\gamma$ | $\beta$ | $\sigma_{12}$ | $f$ | $b_{1}\sigma_{8}(z)$ | $f\sigma_{8}(z)$
---|---|---|---|---|---|---|---|---|---
| 2PCF | – | – | - | $1.2\pm 0.6$ | $2.4\pm 2.1$ | – | $0.56^{+0.12}_{-0.11}$ | $0.7\pm 0.2$
G1 | 3PCF | – | – | $0.1^{+0.8}_{-1.3}$ | $0.02^{+0.70}_{-0.0}$ | – | – | $1.0\pm 0.3$ | $0.0\pm 0.5$
| Joint | $0.50\pm 0.12$ | $1.60\pm 0.43$ | $-0.1^{+0.8}_{-1.3}$ | $0.4^{+0.3}_{-0.2}$ | $2.4\pm 2.0$ | $0.64^{+0.55}_{-0.37}$ | $0.84^{+0.09}_{-0.14}$ | $0.36^{+0.17}_{-0.12}$
| 2PCF | – | – | – | $0.6^{+0.9}_{-0.3}$ | $1.0^{+1.6}_{-1.0}$ | – | $0.67^{+0.15}_{-0.17}$ | $0.49^{+0.42}_{-0.16}$
G2 | 3PCF | – | – | $0.2^{+3.5}_{-1.3}$ | $0.0^{+1.4}_{-0.0}$ | – | – | $0.4^{+0.5}_{-0.2}$ | $0.0^{+0.7}_{-0.0}$
| Joint | $0.39^{+0.11}_{-0.13}$ | $1.9^{+0.8}_{-0.5}$ | $0.5^{+1.3}_{-1.2}$ | $0.49^{+0.31}_{-0.23}$ | $1.0^{+1.6}_{-1.0}$ | $1.0\pm 1.0$ | $0.74^{+0.07}_{-0.08}$ | $0.43^{+0.16}_{-0.15}$

##  8 Discussion and Conclusions

Figure 9: Comparison of our new estimates of the growth rate of structure,
parameterised by the product $f\sigma_{8}$ (light blue circles), with previous
measurements from the literature: 6dFGS (Beutler et al., 2012; Huterer et al.,
2017; Adams & Blake, 2017); 2MASS (Davis et al., 2011); GAMA (Blake et al.,
2013); WiggleZ (Blake et al., 2012); VVDS (Guzzo et al., 2008); VIPERS (de la
Torre et al., 2013, 2017; Pezzotta et al., 2017; Hawken et al., 2017; Mohammad
et al., 2018); FastSound (Okumura et al., 2016); SDSS+BOSS+eBOSS (eBOSS
Collaboration et al., 2020), "$C^{3}$ \- Cluster clustering cosmology"
(Marulli et al., 2020). The black solid line shows the $\Lambda$CDM + GR
Planck Collaboration et al. (2020) model prediction. Figure 10:  Values for
$f(z)$ (left panel) and $\sigma_{8}(z)$, as estimated from our joint 2PCF and
3PCF analysis (light blue circles), again compared to other recent
measurements using different techniques: VIPERS joint lensing-RSD analysis (de
la Torre et al., 2017), BOSS-DR12 (Gil-Marín et al., 2017), Wiggle-z (Marín et
al., 2013). For completeness, we also include results at smaller redshifts by
Shi et al. (2018), Singh et al. (2019) and Jullo et al. (2019).

In this work, we have measured the 2PCF and 3PCF of the galaxies in the final
data release of the VIPERS survey and used their joint information to break
parameter degeneracies and estimate the galaxy bias parameters, $b_{1}$ and
$\gamma$, the clustering amplitude, $\sigma_{8}$, and the linear growth rate
of density fluctuations, $f$, at two redshifts $z\simeq 0.61$ and $z\simeq
0.8$.

Because of the survey footprint, separated in two patches highly elongated in
one direction, we performed the analysis in configuration space focusing on
scales smaller than $40\,h^{-1}\,\mathrm{Mpc}$. To do so, we adopted a 1-loop
anisotropic 2PCF model to account for non-linear effects, galaxy bias and
redshift distortions. For the 3PCF we considered the tree-level model of its
monopole moment in redshift space proposed by Slepian et al. (2017b). To
estimate the 3PCF we have used the SHD method of Slepian & Eisenstein (2015)
whose efficiency allowed us to measure the 3PCF of all 153 mock VIPERS
catalogs and thus to estimate the covariance matrix used in the likelihood
analysis. Because of all these aspects, we believe that our analysis improves
over the previous joint 2PCF and 3PCF analysis of the WiggleZ galaxy sample
that was performed by Marín et al. (2013) using the brute-force triplet
counting estimator simplified analytic models for 2PCF, 3PCF and their
covariance matrix.

The main results of this work can be summarised as follow:

  * •

The results of our 2PCF analysis agree, as expected, with those obtained by
Pezzotta et al. (2017) in the range $[5,50]\,h^{-1}\,\mathrm{Mpc}$ despite
having used a different covariance matrix and a diffeent 2PCF model. Indeed,
our estimates of

$\displaystyle f\sigma_{8}(z=0.61)$ $\displaystyle=0.58\pm 0.12$
$\displaystyle f\sigma_{8}(z=0.87)$ $\displaystyle=0.43\pm 0.12$

and

$\displaystyle b_{1}\sigma_{8}(z=0.61)$ $\displaystyle=0.69\pm 0.04$
$\displaystyle b_{1}\sigma_{8}(z=0.87)$ $\displaystyle=0.70\pm 0.04$

match those of Pezzotta et al. (2017) ($f\sigma_{8}(z=0.61)=0.55\pm
0.12,\,f\sigma_{8}(z=0.87)=0.40\pm 0.11$; $b_{1}\sigma_{8}(z=0.61)=0.73\pm
0.03,\,b_{1}\sigma_{8}(z=0.87)=0.74\pm 0.04$).

When we consider pairs with $s_{min}\geq 10$ $h^{-1}\,\mathrm{Mpc}$ the 1D
posterior probability distribution of $\gamma\equiv b_{2}/b_{1}$ is flat,
whereas it peaks at negative values when pairs at separations as small as 5
$h^{-1}\,\mathrm{Mpc}$ are included. Since a negative value is neither
confirmed by previous VIPERS analyses (Cappi et al., 2015; Di Porto et al.,
2016) nor by our joint 2PCF and 3PCF analysis, we consider it as a piece of
evidence that the 2PCF model is inadequate to describe galaxy clustering on
such small scales. Moreover, we found that including pairs with separations
larger than 40 $h^{-1}\,\mathrm{Mpc}$, where shot noise starts to be
significant, increases the computational cost of a less accurate covariance
matrix and does not improve the quality of the fit.

For all these reasons, we focused on a smaller scale range
$[15,40]\,h^{-1}\,\mathrm{Mpc}$, and then assessed the robustness of our
results to this choice. We then repeated the 2PCF analysis of G1 and G2
samples using this range and found $f\sigma_{8}(z=0.61)=0.47^{+0.17}_{-0.15}$.
This is in good agreement with the estimates of de la Torre et al. (2017),
$f\sigma_{8}(z=0.6)=0.48\pm 0.11$, and of Gil-Marín et al. (2017),
$f\sigma_{8}(z=0.57)=0.432\pm 0.022$. Our measurements are shown in Fig. 9
together with existing estimates of $f\sigma_{8}$ obtained at different
redshifts. As evident, there is an excellent agreement with previous results
and with the theoretical expectation of a Planck cosmology (Planck
Collaboration et al., 2020).

  * •

We measured the 3PCF of the G1 and G2 samples in the same range as the 2PCF,
$[15,40]$ $h^{-1}\,\mathrm{Mpc}$, for all triangle configurations. We compared
these measurements with the model predictions to estimate the bias parameters
$b_{1}$ and $\gamma$ (for the remaining parameters we have the assumed a local
Lagrangian bias model). We found $\gamma(z=0.61)=0.1^{+0.8}_{-1.3}$ and
$\gamma(z=0.8)=0.2^{+3.5}_{-1.3}$, in tension with the result of Moresco et
al. (2017) ($\gamma(z\simeq 0.6)=-0.47\pm 0.144$). The mismatch probably
reflects the different scales probed by the two cases, since Moresco et al.
(2017) pushed their analysis to scales as small as $1h^{-1}\,\mathrm{Mpc}$
(and fixed all parameters of the model except $b_{1}$ and $b_{2}$).

Di Porto et al. (2016) estimated the biasing function from a count-in-cells
analysis of the first data release of VIPERS. Among the shape parameters are
used to characterize the biasing function in that analysis the one dubbed
$B\equiv[1-\tilde{b}/\hat{b}]$, can be directly compared to $\gamma$. The
measured $B$ values at the redshifts of the G1 and G2 samples and on the scale
of $8$ $h^{-1}\,\mathrm{Mpc}$ ($B(z=0.6)=0.007\pm 0.006$ and
$B(z=0.8)=0.005\pm 0.005$) agree with our estimates of $\gamma$,

Cappi et al. (2015) also performed counts in cells to derive volume-averaged
higher-order correlation functions from which they inferred the bias
parameters $b_{1}$ and $b_{2}$. The scales considered in their analyses and
the magnitude cuts are, however, different from ours. Nevertheless, it is
reassuring that, similarly to our case, they did not detect significant
deviations from a linear biasing.

  * •

We have shown that a joint 2- and 3-point correlation analysis of the VIPERS
samples successfully breaks the degeneracy among $\sigma_{8}$, $f$ and $b_{1}$
(and, consequently, $b_{2}$) and reduces the errors in the estimate of their
combinations $f\sigma_{8}$, $b_{1}\sigma_{8}$ and $f/b_{1}$. It is instructive
to compare our results to those of similar analyses, aimed at breaking
parameter degeneracies by combining 2- and 3-point statistics in configuration
(Marín et al., 2013) and Fourier space (Gil-Marín et al., 2017) and by
combining clustering and gravitational lensing measurements using the VIPERS
PDR2 data (de la Torre et al., 2017).

Our estimates of the clustering amplitude, $\sigma_{8}(z=0.61)=0.50\pm 0.12$
and $\sigma_{8}(z=0.8)=0.39^{+0.11}_{-0.13}$, are in good agreement with those
obtained by de la Torre et al. (2017)($\sigma_{8}(z=0.6)=0.52\pm 0.06$ and
$\sigma_{8}(z=0.86)=0.48\pm 0.04$), Marín et al. (2013)
($\sigma_{8}(z=0.55)=0.61^{+0.08}_{-0.09}$), Gil-Marín et al. (2017)
($\sigma_{8}(z=0.57)=0.66\pm 0.067$) and with the Planck $\Lambda$CDM
predictions, as shown in the right panel of Fig. 10.

Similarly, our estimates of the growth rate $f(z=0.61)=0.64^{+0.55}_{-0.37}$
and $f(z=0.8)=1.0\pm 1.0$ agree with those of de la Torre et al. (2017)
($f(z=0.6)=0.93\pm 0.22$ and $f(z=0.86)=0.99\pm 0.19$) and Gil-Marín et al.
(2017) ($f(z=0.57)=0.649\pm 0.076$) within the error that is much larger in
our case. All these measurements are shown in the left panel of Fig. 10
together with the Planck $\Lambda$CDM predictions (Planck Collaboration et
al., 2020).

We notice that the errors on both $f$ and $\sigma_{8}$ in our case are
significantly larger than in de la Torre et al. (2017), despite having used
similar datasets. The reason for this is twofold. First of all, we have
adopted a more conservative approach and considered $s_{min}=15$
$h^{-1}\,\mathrm{Mpc}$, whereas de la Torre et al. (2017) used $s_{min}=5$
$h^{-1}\,\mathrm{Mpc}$. This choice, motivated by the discrepant $\gamma$
value found in the 2PCF analysis, prevents us from accessing information
encoded on smaller scales. If these are included, then the errors in the
measured $f$ and $\sigma_{8}$ significantly decrease towards values comparable
to those of de la Torre et al. (2017), as illustrated in Appendix C. The
second reason is that the clustering-lensing analysis is more efficient in
breaking parameter degeneracy than the joint 2PCF-3PCF. In fact, the latter
combines the galaxy 2PCF, which has a $b_{1}\sigma_{8}$ degeneracy, to the
galaxy 3PCF, which has a similar degeneracy, $b_{1}^{3}\sigma_{8}^{4}$. On the
contrary, the lensing analysis has a $b_{1}\sigma_{8}^{2}$ degeneracy, which
is effectively broken when combined with the galaxy 2PCF.

  * •

We have performed many tests to check the robustness of our results against
non-linear effects, which are expected to be relevant on the scales probed
here and are the reason for the conservative cut at $s_{min}=15$
$h^{-1}\,\mathrm{Mpc}$. Decreasing $s_{min}$ reduces parameters errors from
the 2PCF only analysis but also picks up a negative value for $\gamma$, in
mild tension (given the considerable uncertainty) with the joint 2PCF-3PCF
analysis and previous VIPERS analyses. On the contrary, reducing $s_{min}$
does not have a significant impact on the 3PCF and the joint analyses. In
Table 4 we list the best fit values of $\sigma_{8}$, $f$, $b_{1}$ and $\gamma$
obtained for different choices of $s_{min}$. Errors increase when $s_{min}$
increases, as a result of the reduction of the number of pairs/triplets and
scales included in the analysis. In particular, with
$s_{min}=20\,h^{-1}\,\mathrm{Mpc}$ no significant constraints can be set on
$\gamma$. A detailed analysis of the sensitivity of our results to the choice
of $s_{min}$ is presented in Appendix C.

Table 4: Best-fit values of the parameters obtained in the joint analysis of 2PCF and 3PCF, using different values of $s_{min}$. The values of $s_{min}$ are in $h^{-1}\,\mathrm{Mpc}$. All the constraints are compatible, within 1 $\sigma$. Sample | $s_{min}$ | $\sigma_{8}(z)$ | $f$ | $b_{1}$ | $\gamma$
---|---|---|---|---|---
| $10$ | $0.4\pm 0.1$ | $1.0^{+0.7}_{-0.5}$ | $1.8^{+0.4}_{-0.3}$ | $-0.1\pm 0.6$
G1 | 15 | $0.5\pm 0.1$ | $0.6^{+0.6}_{-0.4}$ | $1.6\pm 0.4$ | $-0.1^{+0.9}_{-1.3}$
| 20 | $0.4^{+0.3}_{-0.2}$ | $1.1^{+1.6}_{-1.0}$ | $0.5^{+0.7}_{-0.4}$ | unconstrained
| 10 | $0.3\pm 0.1$ | $1.4^{+1.2}_{-0.8}$ | $2.0^{+1.0}_{-0.6}$ | $-1.1^{+1.4}_{-2.2}$
G2 | 15 | $0.4\pm 0.1$ | $1.0\pm 1.0$ | $1.9^{+0.8}_{-0.5}$ | $0.5^{+1.3}_{-1.2}$
| 20 | $0.2^{+0.2}_{-0.1}$ | $0.9^{+1.4}_{-0.8}$ | $1.8^{+1.9}_{-0.9}$ | $-0.2^{+1.5}_{-3.9}$

Our work confirms the importance of clustering analyses beyond 2-point
statistics. Here, we were able to break parameter degeneracies using a
relatively modest number of objects (23352 and 13046 for the G1 and G2
samples, respectively). In this respect, this work represents a successful
pilot study in the preparation for the next generation spectroscopic surveys,
such as the DESI project, the Euclid (Laureijs et al., 2011) and the Roman
(Akeson et al., 2019) space telescope missions. These surveys will be able to
perform clustering analyses on scales much larger than those considered here.
However, small scales will still encode an even larger amount of information.
Its extraction will require measuring higher-order statistics and comparing
results to non-linear models, as we have done in this work. In this respect,
one of the main lessons learned from our analysis is the need to develop a
full 1-loop model for the 3PCF, matching the 2PCF one, which we plan to
present in future work.

The availability of such a model, along with that of a new generation of
efficient 3PCF estimators (Slepian & Eisenstein, 2015, 2016), would make
higher-order clustering analyses in configuration space more palatable and a
serious contender to more traditional Fourier-space methods. Bispectrum
analyses enjoy the availability of fast estimators and non-linear models.
However, they suffer from mode coupling induced by complex survey geometries,
which is difficult to account for. Comparing the performances of joint 2-point
and 3-point clustering analyses in configuration and Fourier space for next-
generation spectroscopic surveys is another open issue, which we plan to
investigate in the future.

Finally, in this analysis, we were able to break parameter degeneracies using
higher-order statistics. Alternatively, it can be broken by combining 2-point
clustering and gravitational lensing, as in de la Torre et al. (2017). Both
approaches have advantages and disadvantages. The 2- and 3-point clustering
analysis requires a single (spectroscopic) dataset to be performed. However,
it is less efficient in removing parameters’ degeneracy. The clustering-
lensing analysis is more effective in breaking degeneracy but needs both a
photometric and a spectroscopic survey to be performed. How to best combine
these two types of analyses is another interesting issue that deserves a
dedicated future study.

## Data Availability

The VIPERS PDR2 data, as well as the mock samples used here are publicly
available from the VIPERS web site (http://vipers.inaf.it). The clustering
measurements and covariance matrices from this paper are available from the
authors, upon request.

## Acknowledgements

AV and EB thank Emiliano Sefusatti, Alexander Eggemeier, Elena Sarpa, Massimo
Guidi for useful discussions. This paper uses data from the VIMOS Public
Extragalactic Redshift Survey (VIPERS). VIPERS has been performed using the
ESO Very Large Telescope, under the “Large Programme" 182.A-0886. The
participating institutions and funding agencies are listed at
http://vipers.inaf.it. AV, EB, LG and LM are supported by ASI/INAF agreement
n. 2018-23-HH.0 “Scientific activity for Euclid mission, Phase D" and INFN
project “InDark". EB and LG are further supported by MIUR/PRIN 2017 “From
Darklight to Dark Matter: understanding the galaxy-matter connection to
measure the Universe". EB is also supported by ASI/INAF agreement n.
2017-14-H.O “Unveiling Dark Matter and Missing Baryons in the high-energy
sky". MM acknowledges the grants ASI n.I/023/12/0, ASI n.2018-23-HH.0, and
support from MIUR, PRIN 2017 (grant 20179ZF5KS). This research was supported
by the Munich Institute for Astro- and Particle Physics (MIAPP) which is
funded by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) under Germany’s Excellence Strategy- EXC - 2094 - 390783311. LM
also acknowledge the support from the grant PRIN-MIUR 2017 WSCC32.

## References

  * Adams & Blake (2017) Adams C., Blake C., 2017, MNRAS, 471, 839
  * Akeson et al. (2019) Akeson R., et al., 2019, arXiv e-prints, p. arXiv:1902.05569
  * Alam et al. (2021) Alam S., et al., 2021, Phys. Rev. D, 103, 083533
  * Alcock & Paczynski (1979) Alcock C., Paczynski B., 1979, Nature, 281, 358
  * Baldauf et al. (2012) Baldauf T., Seljak U., Desjacques V., McDonald P., 2012, Phys. Rev. D, 86, 083540
  * Bel et al. (2019) Bel J., Pezzotta A., Carbone C., Sefusatti E., Guzzo L., 2019, A&A, 622, A109
  * Bernardeau et al. (2002) Bernardeau F., Colombi S., Gaztañaga E., Scoccimarro R., 2002, Phys. Rep., 367, 1
  * Beutler et al. (2012) Beutler F., et al., 2012, MNRAS, 423, 3430
  * Beutler et al. (2014) Beutler F., et al., 2014, MNRAS, 443, 1065
  * Blake et al. (2012) Blake C., et al., 2012, MNRAS, 425, 405
  * Blake et al. (2013) Blake C., et al., 2013, MNRAS, 436, 3089
  * Cappi et al. (2015) Cappi A., et al., 2015, A&A, 579, A70
  * Chan et al. (2012) Chan K. C., Scoccimarro R., Sheth R. K., 2012, Phys. Rev. D, 85, 083509
  * Cole et al. (2005) Cole S., et al., 2005, MNRAS, 362, 505
  * DESI Collaboration et al. (2016a) DESI Collaboration et al., 2016a, arXiv e-prints, p. arXiv:1611.00036
  * DESI Collaboration et al. (2016b) DESI Collaboration et al., 2016b, arXiv e-prints, p. arXiv:1611.00037
  * Davidzon et al. (2016) Davidzon I., et al., 2016, A&A, 586, A23
  * Davis et al. (2011) Davis M., Nusser A., Masters K. L., Springob C., Huchra J. P., Lemson G., 2011, MNRAS, 413, 2906
  * Di Porto et al. (2016) Di Porto C., et al., 2016, A&A, 594, A62
  * Eisenstein et al. (2005) Eisenstein D. J., et al., 2005, ApJ, 633, 560
  * Fry (1994) Fry J. N., 1994, Phys. Rev. Lett., 73, 215
  * Garilli et al. (2014) Garilli B., et al., 2014, A&A, 562, A23
  * Gaztañaga et al. (2009) Gaztañaga E., Cabré A., Castand er F., Crocce M., Fosalba P., 2009, MNRAS, 399, 801
  * Gil-Marín et al. (2014) Gil-Marín H., Wagner C., Noreña J., Verde L., Percival W., 2014, J. Cosmology Astropart. Phys., 2014, 029
  * Gil-Marín et al. (2015) Gil-Marín H., Noreña J., Verde L., Percival W. J., Wagner C., Manera M., Schneider D. P., 2015, MNRAS, 451, 539
  * Gil-Marín et al. (2017) Gil-Marín H., Percival W. J., Verde L., Brownstein J. R., Chuang C.-H., Kitaura F.-S., Rodríguez-Torres S. A., Olmstead M. D., 2017, MNRAS, 465, 1757
  * Granett et al. (2019) Granett B. R., Favole G., Montero-Dorta A. D., Branchini E., Guzzo L., de la Torre S., 2019, MNRAS, 489, 653
  * Groth & Peebles (1977) Groth E. J., Peebles P. J. E., 1977, ApJ, 217, 385
  * Guo et al. (2014) Guo H., Li C., Jing Y. P., Börner G., 2014, ApJ, 780, 139
  * Guo et al. (2015) Guo H., et al., 2015, MNRAS, 449, L95
  * Guzzo et al. (2008) Guzzo L., et al., 2008, Nature, 451, 541
  * Guzzo et al. (2014) Guzzo L., et al., 2014, A&A, 566, A108
  * Hamilton (2000) Hamilton A. J. S., 2000, MNRAS, 312, 257
  * Hartlap et al. (2007) Hartlap J., Simon P., Schneider P., 2007, A&A, 464, 399
  * Hawken et al. (2017) Hawken A. J., et al., 2017, A&A, 607, A54
  * Hoffmann et al. (2018) Hoffmann K., Gaztañaga E., Scoccimarro R., Crocce M., 2018, MNRAS, 476, 814
  * Hui & Gaztañaga (1999) Hui L., Gaztañaga E., 1999, ApJ, 519, 622
  * Huterer & Linder (2007) Huterer D., Linder E. V., 2007, Phys. Rev. D, 75, 023519
  * Huterer et al. (2017) Huterer D., Shafer D. L., Scolnic D. M., Schmidt F., 2017, J. Cosmology Astropart. Phys., 2017, 015
  * Jing & Börner (2004) Jing Y. P., Börner G., 2004, ApJ, 607, 140
  * Jullo et al. (2019) Jullo E., et al., 2019, A&A, 627, A137
  * Kaiser (1987) Kaiser N., 1987, MNRAS, 227, 1
  * Kayo et al. (2004) Kayo I., et al., 2004, PASJ, 56, 415
  * Kazin et al. (2013) Kazin E. A., et al., 2013, MNRAS, 435, 64
  * Klypin et al. (2016) Klypin A., Yepes G., Gottlöber S., Prada F., Heß S., 2016, MNRAS, 457, 4340
  * Kulkarni et al. (2007) Kulkarni G. V., Nichol R. C., Sheth R. K., Seo H.-J., Eisenstein D. J., Gray A., 2007, MNRAS, 378, 1196
  * Landy & Szalay (1993) Landy S. D., Szalay A. S., 1993, ApJ, 412, 64
  * Laureijs et al. (2011) Laureijs R., et al., 2011, arXiv e-prints, p. arXiv:1110.3193
  * Lesgourgues (2011) Lesgourgues J., 2011, arXiv e-prints, p. arXiv:1104.2932
  * Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473
  * Maartens et al. (2015) Maartens R., Abdalla F. B., Jarvis M., Santos M. G., 2015, arXiv e-prints, p. arXiv:1501.04076
  * Marín (2011) Marín F., 2011, ApJ, 737, 97
  * Marín et al. (2013) Marín F. A., et al., 2013, MNRAS, 432, 2654
  * Marulli et al. (2012) Marulli F., Bianchi D., Branchini E., Guzzo L., Moscardini L., Angulo R. E., 2012, MNRAS, 426, 2566
  * Marulli et al. (2013) Marulli F., et al., 2013, A&A, 557, A17
  * Marulli et al. (2016) Marulli F., Veropalumbo A., Moresco M., 2016, Astronomy and Computing, 14, 35
  * Marulli et al. (2020) Marulli F., Veropalumbo A., García-Farieta J. E., Moresco M., Moscardini L., Cimatti A., 2020, arXiv e-prints, p. arXiv:2010.11206
  * Matarrese et al. (1997) Matarrese S., Verde L., Heavens A. F., 1997, MNRAS, 290, 651
  * McBride et al. (2011a) McBride C. K., Connolly A. J., Gardner J. P., Scranton R., Newman J. A., Scoccimarro R., Zehavi I., Schneider D. P., 2011a, ApJ, 726, 13
  * McBride et al. (2011b) McBride C. K., Connolly A. J., Gardner J. P., Scranton R., Scoccimarro R., Berlind A. A., Marín F., Schneider D. P., 2011b, ApJ, 739, 85
  * McDonald & Roy (2009) McDonald P., Roy A., 2009, J. Cosmology Astropart. Phys., 2009, 020
  * Mohammad et al. (2018) Mohammad F. G., et al., 2018, A&A, 610, A59
  * Moresco et al. (2014) Moresco M., Marulli F., Baldi M., Moscardini L., Cimatti A., 2014, MNRAS, 443, 2874
  * Moresco et al. (2017) Moresco M., et al., 2017, A&A, 604, A133
  * Moresco et al. (2020) Moresco M., Veropalumbo A., Marulli F., Moscardini L., Cimatti A., 2020, arXiv e-prints, p. arXiv:2011.04665
  * Nichol et al. (2006) Nichol R. C., et al., 2006, MNRAS, 368, 1507
  * Okumura et al. (2016) Okumura T., et al., 2016, PASJ, 68, 38
  * Peacock et al. (2001) Peacock J. A., et al., 2001, Nature, 410, 169
  * Pearson & Samushia (2019) Pearson D. W., Samushia L., 2019, MNRAS, 486, L105
  * Percival et al. (2014) Percival W. J., et al., 2014, MNRAS, 439, 2531
  * Pezzotta et al. (2017) Pezzotta A., et al., 2017, A&A, 604, A33
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Rota et al. (2017) Rota S., et al., 2017, A&A, 601, A144
  * Saito et al. (2014) Saito S., Baldauf T., Vlah Z., Seljak U., Okumura T., McDonald P., 2014, Phys. Rev. D, 90, 123522
  * Scoccimarro (2004) Scoccimarro R., 2004, Phys. Rev. D, 70, 083007
  * Scoccimarro et al. (1999) Scoccimarro R., Zaldarriaga M., Hui L., 1999, ApJ, 527, 1
  * Scodeggio et al. (2018) Scodeggio M., et al., 2018, A&A, 609, A84
  * Sefusatti et al. (2006) Sefusatti E., Crocce M., Pueblas S., Scoccimarro R., 2006, Phys. Rev. D, 74, 023522
  * Sellentin & Heavens (2016) Sellentin E., Heavens A. F., 2016, MNRAS, 456, L132
  * Sereno et al. (2015) Sereno M., Veropalumbo A., Marulli F., Covone G., Moscardini L., Cimatti A., 2015, MNRAS, 449, 4147
  * Shi et al. (2018) Shi F., et al., 2018, ApJ, 861, 137
  * Singh et al. (2019) Singh S., Alam S., Mandelbaum R., Seljak U., Rodriguez-Torres S., Ho S., 2019, MNRAS, 482, 785
  * Slepian & Eisenstein (2015) Slepian Z., Eisenstein D. J., 2015, MNRAS, 454, 4142
  * Slepian & Eisenstein (2016) Slepian Z., Eisenstein D. J., 2016, MNRAS, 455, L31
  * Slepian & Eisenstein (2017) Slepian Z., Eisenstein D. J., 2017, MNRAS, 469, 2059
  * Slepian & Eisenstein (2018) Slepian Z., Eisenstein D. J., 2018, MNRAS, 478, 1468
  * Slepian et al. (2017a) Slepian Z., et al., 2017a, MNRAS, 468, 1070
  * Slepian et al. (2017b) Slepian Z., et al., 2017b, MNRAS, 469, 1738
  * Sugiyama et al. (2020) Sugiyama N. S., Saito S., Beutler F., Seo H.-J., 2020, MNRAS, 497, 1684
  * Szapudi (2004) Szapudi I., 2004, ApJ, 605, L89
  * Szapudi & Szalay (1998) Szapudi I., Szalay A. S., 1998, ApJ, 494, L41
  * Takahashi et al. (2012) Takahashi R., Sato M., Nishimichi T., Taruya A., Oguri M., 2012, ApJ, 761, 152
  * Taruya et al. (2010) Taruya A., Nishimichi T., Saito S., 2010, Phys. Rev. D, 82, 063522
  * Wang & Steinhardt (1998) Wang L., Steinhardt P. J., 1998, ApJ, 508, 483
  * Yankelevich & Porciani (2019) Yankelevich V., Porciani C., 2019, MNRAS, 483, 2078
  * de la Torre et al. (2013) de la Torre S., et al., 2013, A&A, 557, A54
  * de la Torre et al. (2017) de la Torre S., et al., 2017, A&A, 608, A44
  * eBOSS Collaboration et al. (2020) eBOSS Collaboration et al., 2020, arXiv e-prints, p. arXiv:2007.08991

##  Appendix A 3PCF estimator binning

The algorithm of Slepian & Eisenstein (2015) used in this work estimates the
multipole coefficients of the 3PCF Legendre expansion in radial bins
$\Delta_{r_{12}}$ and $\Delta_{r_{13}}$,
$\zeta_{l}(\Delta_{r_{12}},\Delta_{r_{13}})$. To obtain an unbiased estimate
of $\zeta(\Delta_{r_{12}},\Delta_{r_{13}},\Delta_{\mu})$ in bins
$\Delta_{\mu}$, one needs to use bin-averaged Legendre polynomials
$\overline{\mathcal{P}}_{l}(\Delta_{\mu})$, i.e.

$\zeta(\Delta_{r_{12}},\Delta_{r_{13}},\Delta_{\mu})\,=\sum_{l=0}^{l_{max}}\zeta_{l}(\Delta_{r_{12}},\Delta_{r_{13}})\overline{\mathcal{P}}_{l}(\Delta_{\mu})\,,$
(32)

where

$\overline{\mathcal{P}}_{l}(\Delta_{\mu})\equiv\overline{\mathcal{P}}_{l}(\mu_{min}\leq\mu\leq\mu_{max})=\frac{1}{\Delta\mu}\left[\frac{\mathcal{P}_{l+1}(\mu)-\mathcal{P}_{l-1}(\mu)}{2l+1}\right]_{\mu_{min}}^{\mu_{max}}\,,$
(33)

in the bin $\Delta_{\mu}=[\mu_{min},\mu_{max}]$; when $l=0$,
$\overline{\mathcal{P}}_{l=0}(\Delta_{\mu})=1$.

In this work, we express the 3PCF as a function of $r_{23}$, rather than
$\mu$, with the two quantities related through the relation

$\mu=\frac{r^{2}_{12}+r^{2}_{13}-r^{2}_{23}}{2r_{12}r_{13}}\,,$ (34)

so that, given $r_{12}$ and $r_{13}$, the cosine angle $\mu$ varies in the
range $[0,1]$, whereas $r_{23}$ varies between $|r_{12}-r_{13}|$ and
$|r_{12}+r_{13}|$. One can show that in this case the binned Legendre
polynomials are of the form

$\overline{\mathcal{P}}_{l}(\Delta_{r_{12}},\Delta_{r_{13}},\Delta_{r_{23}})\,=\frac{32\pi}{V_{12}V_{13}}\int
k^{2}I_{l}(k;\Delta_{r_{12}})I_{l}(k;\Delta_{r_{13}})I_{0}(k;\Delta_{r_{23}})\mathrm{d}k\,,$
(35)

where $V_{12}$ and $V_{13}$ represent the volume of the spherical shells of
width $\Delta_{r_{12}}$, $\Delta_{r_{13}}$ respectively, and

$I_{l}(k;\Delta_{r})=\int_{r_{min}}^{r_{max}}r^{2}j_{l}(kr)\mathrm{d}r\,,$
(36)

where $j_{l}(kr)$ are spherical Bessel functions.

##  Appendix B Explicit expressions for the 3PCF model

The explicit expression of the $A_{l}$ terms that appear in Eq. (25) is

$\displaystyle A_{0}=$ $\displaystyle
b_{1}^{3}\left\\{\frac{34}{21}\left[1+\frac{4}{3}\beta+\frac{1154}{1275}\beta^{2}+\frac{936}{2975}\beta^{3}+\frac{21}{425}\beta^{4}\right]\right.$
(37)
$\displaystyle\left.+\gamma\left[1+\frac{2}{3}\beta+\frac{1}{9}\beta^{2}\right]+\frac{16}{675}\beta^{2}\gamma^{\prime}\right\\},$
$\displaystyle A_{1}=$ $\displaystyle-
b_{1}^{3}\left[1+\frac{4}{3}\beta+\frac{82}{75}\beta^{2}+\frac{12}{25}\beta^{3}+\frac{3}{35}\beta^{4}\right],$
$\displaystyle A_{2}=$ $\displaystyle
b_{1}^{3}\left\\{\frac{8}{21}\left[1+\frac{4}{3}\beta+\frac{52}{21}\beta^{2}+\frac{81}{49}\beta^{3}+\frac{12}{35}\beta^{4}\right]\right.$
$\displaystyle\left.+\frac{32\gamma}{945}\beta^{2}+\frac{5}{2}\left(\frac{8}{15}+\frac{16\beta}{45}+\frac{344\beta^{2}}{4725}\right)\gamma^{\prime}\right\\},$
$\displaystyle A_{3}=$ $\displaystyle-
b_{1}^{3}\left[\frac{8}{75}\beta^{2}+\frac{16}{175}\beta^{3}+\frac{8}{315}\beta^{4}\right],$
$\displaystyle A_{4}=$ $\displaystyle
b_{1}^{3}\left[-\frac{32}{3675}\beta^{2}+\frac{32}{8575}\beta^{3}+\frac{128}{11025}\right].$

In the same equation, the $f_{l}$ terms are:

$\displaystyle
f_{l}(r_{12},r_{13},r_{23})=\left\\{\begin{array}[]{lr}\xi^{[l]}(r_{12})\xi^{[l]}(r_{13})\mathcal{L}_{l}\left(\mu_{23}\right)+\text{cyc.}&\text{
if $l$ is even}\\\ \\\ \left[\xi^{[l+]}(r_{12})\xi^{[l-]}(r_{13})+\right.&\\\
\left.\xi^{[l+]}(r_{13})\xi^{[l-]}(r_{12})\right]\mathcal{L}_{l}\left(\mu_{23}\right)+\text{cyc.}&\text{
if $l$ is odd}\,,\end{array}\right.$ (38)

where

$\begin{split}\xi_{i}^{[l]}&=\int\frac{k^{2}\mathrm{d}k}{2\pi^{2}}{2\pi^{2}}P_{\rm
lin}(k)j_{l}\left(kr_{i}\right)\\\
\xi_{i}^{\left[l\pm\right]}&=\int\frac{k^{2}\mathrm{d}k}{2\pi^{2}}k^{\pm
1}P_{\rm lin}(k)j_{l}\left(kr_{i}\right).\end{split}$ (39)

and $P_{\rm lin}(k)$ is the linear matter power spectrum.

##  Appendix C Sensitivity to nonlinear effects

Figure 11: Average estimates over the 153 mock samples of $f\sigma_{8}$ using
the 2PCF (top, filled points) and the 3PCF (bottom, filled squares), as a
function of the minimum scale included in the analysis, $s_{min}$. The black
dashed line in both cases gives the reference value, corresponding to the
cosmology of the simulated mocks, while error bars are given by the scatter
among the mocks. Note how the error bars for the 3PCF estimates are $\sim 10$
times larger than those from the 2PCF, due to the fact that for the latter
only the monopole (i.e. the isotropic information), has been considered.
Figure 12: Average estimates of $f(z)$ and $\sigma_{8}(z)$ from the joint
2PCF-3PCF analysis of the 153 mocks (filled blue squares), compared to the
reference cosmology (black dashed line). As usual, error bars are the standard
deviation over the mocks. While the estimates of $f(z)$ are robust, recovering
the input value of the mocks, there is a small bias for $\sigma_{8}$ (still
within 1$\sigma$ at the reference $s_{min}=15h^{-1}\,\mathrm{Mpc}$.

Since our clustering analysis includes scales smaller than
$40h^{-1}\,\mathrm{Mpc}$, nonlinear effects in both dynamics and galaxy bias
cannot be neglected. To minimise the impact of potential systematic errors
derived from an incorrect model of these effects, we have adopted a
conservative approach and excluded scales below $15h^{-1}\,\mathrm{Mpc}$ from
the analysis. These are scales, however, where a significant amount of
cosmological information is stored. In this Appendix, we use the realistic
mock VIPERS catalogues to perform several tests to assess the sensitivity of
our results to the choice of this minimum scale, $s_{min}$. We are confident
that these tests provide useful indications for the real data analysis since,
as we have discussed in Sect. 3, the 2- and 3-point correlation functions
measured in the mocks agree well with those measured on the real data on all
scales considered here.

Let us first consider the parameter combination $f\sigma_{8}$, as obtained
from the 2PCF and the 3PCF. We show here the results from the G2 sample only,
since they are representative of the generality of the results. In Fig. 11 we
show the average and rms scatter over the 153 mocks of the value of
$f\sigma_{8}$ estimated using the 2PCF and the 3PCF (top and bottom panels,
respectively), as a function of $s_{min}$. The results show that including
scales as small as $5\,h^{-1}\,\mathrm{Mpc}$ significantly reduces the
statistical error without compromising the accuracy. This conclusion is true
for both the 2PCF and the 3PCF cases, indicating that the models we have used
are adequate for measuring this parameter combination.

Things are different when one tries to break the parameter degeneracy through
the joint 2PCF and 3PCF analysis, as shown in Fig. 12. In the two panels, the
average estimates of $\sigma_{8}$ (top) and $f$ (bottom) are plotted as a
function of $s_{min}$. Note how the size of the error bars decreases, when
smaller scales are progressively included in the analysis. For
$s_{min}=5\,h^{-1}\,\mathrm{Mpc}$ they are significantly smaller than in our
baseline case $(s_{min}=15h^{-1}\,\mathrm{Mpc})$. However, while random errors
decrease, systematic errors increase: the value of $\sigma_{8}$ is
systematically underestimated in general (although at the $1\sigma$ level
only) and, to compensate, the measured growth rate is larger than the true
one. These results justify our conservative choice to set
$s_{min}=15\,h^{-1}\,\mathrm{Mpc}$ and indicate that a better modelling is
required to push a joint 2PCF and 3PCF analysis to smaller scales.