# New constraints on cosmological modified gravity theories from anisotropic
three-point correlation functions of BOSS DR12 galaxies

Naonori S. Sugiyama1
1 National Astronomical Observatory of Japan, Mitaka, Tokyo 181-8588, Japan
Daisuke Yamauchi2, Tsutomu Kobayashi3, Tomohiro Fujita4,5, Shun Arai6, and
Shin’ichi Hirano7
2 Faculty of Engineering, Kanagawa University, Kanagawa, 221-8686, Japan
3 Department of Physics, Rikkyo University, Toshima, Tokyo 171-8501, Japan
4 Waseda Institute for Advanced Study, Shinjuku, Tokyo 169-8050, Japan
5 Research Center for the Early Universe, The University of Tokyo, Bunkyo,
Tokyo 113-0033, Japan
6 Kobayashi-Maskawa Institute, Nagoya University, Nagoya 464-8602, Japan
7Department of Physics, Tokyo Institute of Technology, 2-12-1 Ookayama,
Meguro-ku, Tokyo 152-8551, Japan
Shun Saito8,9, Florian Beutler10, and Hee-Jong Seo11,12,13
8 Institute for Multi-messenger Astrophysics and Cosmology, Department of
Physics,
Missouri University of Science and Technology, 1315 N. Pine St., Rolla MO
65409, USA
9 Kavli Institute for the Physics and Mathematics of the Universe (WPI),
Todai Institutes for Advanced Study, The University of Tokyo, Chiba 277-8582,
Japan
10 Institute for Astronomy, University of Edinburgh, Royal Observatory,
Blackford Hill, Edinburgh EH9 3HJ, UK
11 Department of Physics and Astronomy, Ohio University, Clippinger Labs,
Athens, OH 45701, USA
12 Physics Division, Lawrence Berkeley National Laboratory, 1 Cyclotron Road,
Berkeley, CA 94720, USA
13 Berkeley Center for Cosmological Physics, Department of Physics, University
of California, Berkeley, CA 94720, USA E-mail<EMAIL_ADDRESS>
###### Abstract

We report a new test of modified gravity theories using the large-scale
structure of the Universe. This paper is the first attempt to (1) apply a
joint analysis of the anisotropic components of galaxy two- and three-point
correlation functions (2 and 3PCFs) to actual galaxy data and (2) constrain
the nonlinear effects of degenerate higher-order scalar-tensor (DHOST)
theories on cosmological scales. Applying this analysis to the Baryon
Oscillation Spectroscopic Survey (BOSS) data release 12, we obtain the lower
bounds of $-1.655<\xi_{\rm t}$ and $-0.504<\xi_{\rm s}$ at the $95\%$
confidence level on the parameters characterising the time evolution of the
tidal and shift terms of the second-order velocity field. These constraints
are consistent with GR predictions of $\xi_{\rm t}=15/1144$ and $\xi_{\rm
s}=0$. Moreover, they represent a $35$-fold and $20$-fold improvement,
respectively, over the joint analysis with only the isotropic 3PCF. We ensure
the validity of our results by investigating various quantities, including
theoretical models of the 3PCF, window function corrections, cumulative ${\rm
S/N}$, Fisher matrices, and statistical scattering effects of mock simulation
data. We also find statistically significant discrepancies between the BOSS
data and the Patchy mocks for the 3PCF measurement. Finally, we package all of
our 3PCF analysis codes under the name HITOMI and make them publicly available
so that readers can reproduce all the results of this paper and easily apply
them to ongoing future galaxy surveys.

###### keywords:

cosmology: large-scale structure of Universe – cosmology: dark matter –
cosmology: observations – cosmology: theory

††pubyear: ††pagerange: New constraints on cosmological modified gravity
theories from anisotropic three-point correlation functions of BOSS DR12
galaxies–A

##  1 Introduction

###  1.1 Outline and summary

This paper presents a comprehensive study of the joint analysis of galaxy two-
and three-point correlation functions (2 and 3PCFs) with isotropic and
anisotropic components to constrain the non-linear effects of modified gravity
theories on a cosmological scale. Section 1 outlines the theoretical
development and the present constraints for scalar-tensor theories. We also
outline the development of the measurement and analysis of the three-point
correlation function of galaxies. We organize this paper such that readers
unfamiliar with both or one of the two areas follow the recent developments
and understand how they fit together.

Readers interested in the theoretical aspects may read Sections 2, 3, and 7.
Section 2 reviews the non-linear evolution of the large-scale structure of the
Universe in scalar-tensor theories. Section 3 presents detailed calculations
of the theoretical model of the 3PCF and, in particular, investigates the
dependence of the parameters that characterise the effect of scalar-tensor
theories on the 3PCF model. Finally, Section 7 discusses the extent to which
the 3PCF contains information on the non-linear effects of scalar-tensor
theories through Fisher analysis.

Readers interested in the analysis method of the 3PCF may read Sections 4, 5,
6, and 8. Section 4 reviews how to measure the 3PCF from galaxy data and
examines the effect of the window function on the measured 3PCFs. Section 5
presents the results of the 2PCF and 3PCF covariance matrices computed from
mock simulations. Section 6 describes the setup for the data analysis in this
paper. Finally, Section 8 discusses in detail whether the 3PCFs measured from
the galaxy data in this paper can be fitted using the corresponding
theoretical model in terms of $p$-values.

For readers familiar with the two areas in the literature and interested in
the final results, we suggest they jump directly to Section 9. The novel
aspect of this paper is to focus on observationally constraining the second-
order velocity field, which is a key to seeking a deviation from GR in scalar-
tensor theories. We also show that the second-order velocity field imprints a
unique signature in the anisotropic 3PCF on large scales. Following Yamauchi &
Sugiyama (2022) and Section 3.4, we then parameterise the effects of scalar-
tensor theories in the time evolution of the shift and tidal terms of the
second-order velocity field using parameters $\xi_{\rm s}$ and $\xi_{\rm t}$
defined in Eq. (3.27). Constraining these parameters using Baryon Oscillation
Spectroscopic Survey Data Release 12 galaxies (Eisenstein et al., 2011; Bolton
et al., 2012; Dawson et al., 2013; Alam et al., 2015), we obtain the following
lower bounds given in Eqs. (9.8) and (9.10):

$\displaystyle-1.655<\xi_{\rm t}\ \ \mbox{and}\ \ -0.504<\xi_{\rm
s}\quad\mbox{($95\%$CL)}.$

Since $\xi_{\rm t}=15/1144$ and $\xi_{\rm s}=0$ in GR, these results are
consistent with GR. Finally, we summarise the final results and the various
findings leading up to them in Section 10, which concludes this paper.

We package all the code used to complete this paper under the name HITOMI
111https://github.com/naonori/hitomi.git and make it publicly available.
Appendix A summarises the structure and usage of HITOMI.

###  1.2 General motivation

The greatest mystery in current cosmology is the cause of the accelerated
expansions that have presumably occurred twice in the cosmic expansion
history: i.e., inflation and late-time acceleration.

Scalar-tensor theories, modified gravity theories that add a single scalar
field degree of freedom to General Relativity (GR), have been actively studied
as a promising candidate to explain these accelerated expansions (for reviews,
see Langlois, 2019; Kase & Tsujikawa, 2019; Kobayashi, 2019; Amendola et al.,
2020; Frusciante & Perenon, 2020)222For reviews of modified gravity theories,
including other theories than scalar-tensor theories, see Nojiri & Odintsov
2011; Sebastiani et al. 2017; Nojiri et al. 2017; Cataneo & Rapetti 2018;
Ishak 2019; Ferreira 2019; Baker et al. 2021; Arai et al. 2022..

The accelerated expansion in the very early Universe, called inflation
(Starobinsky, 1980; Guth, 1981; Sato, 1981; Linde, 1982; Albrecht &
Steinhardt, 1982), is thought to be caused by a single scalar field in the
simplest model, generating the seeds of the cosmic fluctuations currently
observed. Furthermore, the statistical properties of these fluctuations are in
excellent agreement with current observations of the cosmic microwave
background (CMB; Aghanim et al., 2020) and the large-scale structure (LSS;
Alam et al., 2021a). On the other hand, the cosmological constant may explain
the late-time accelerated expansion (Riess et al., 1998; Perlmutter et al.,
1999). However, its smallness implies a serious fine-tuning problem in
fundamental physics (Weinberg, 1989; Martin, 2012), and in order to avoid this
problem, it is preferable to adopt a scalar field that varies with time.

In order to test scalar-tensor theories in the late-time Universe, it is
crucial to follow the time evolution of the large-scale structure in detail.
Examples of already completed galaxy surveys are the Baryon Oscillation
Spectroscopic Survey (BOSS Eisenstein et al., 2011; Bolton et al., 2012;
Dawson et al., 2013; Alam et al.,
2015)333https://www.sdss3.org/science/boss_publications.php and the Extended
BOSS (eBOSS; Dawson et al., 2016; Alam et al.,
2021a)444https://www.sdss.org/surveys/eboss/. Furthermore, The next-generation
galaxy surveys, such as the Dark Energy Spectroscopic Instrument (DESI; DESI
Collaboration et al., 2016)555http://desi.lbl.gov/, Euclid (Laureijs et al.,
2011)666www.euclid-ec.org, and the Subaru Prime Focus Spectrograph (PFS;
Takada et al., 2014)777https://pfs.ipmu.jp/index.html, will provide
unprecedented accuracy in testing scalar-tensor theories.

###  1.3 DHOST theories and their constraints

In this paper, we pay particular attention to the behaviour in the late-time
Universe of Degenerate Higher-Order Scalar-Tensor (DHOST) theory (for reviews,
see Langlois, 2019; Kobayashi, 2019), which are a quite general theoretical
framework of scalar-tensor theories that can evade the Ostrogradsky
instability (Ostrogradsky, 1850; Woodard, 2015; Ganz & Noui, 2020). Scalar-
tensor theories have been developing rapidly over the last decade. In 2011,
Deffayet et al. (2011); Kobayashi et al. (2011) rediscovered the most general
theory with second-order equations of motion for metric tensor and scalar
fields, Horndeski theories (Horndeski, 1974). To go beyond Horndeski theories,
Gleyzes et al. (2015a, b) found a class of healthy theories having higher-
order field equations that reduce to a second-order system by combining
different components (see also Zumalacárregui & García-Bellido, 2014, for
examples beyond Horndeski). This discovery results from a degeneracy between
the kinetic terms of the scalar field and the metric. This class of theories
has been extended to reach DHOST theories (Langlois & Noui, 2016; Crisostomi
et al., 2016; Ben Achour et al., 2016a; Ben Achour et al., 2016b; Langlois,
2017; Langlois et al., 2020), encompassing Horndeski and Beyond Horndeski
theories888Hereafter, we do not distinguish between Beyond Horndeski theories
and DHOST theories.. So far, DHOST theories have been constrained primarily by
three observations999As other probes of DHOST theories, for example, Babichev
& Lehébel (2018) shows that the scalar field in DHOST theories can
significantly modify the speed of sound in the atmosphere of the Earth;
Beltran Jimenez et al. (2016); Dima & Vernizzi (2018) strongly constrain DHOST
models using Hulse-Taylor pulsar observations; Saltas & Lopes (2019) proposes
helioseismology as a precise way to test DHOST theories on astrophysical
scales.: gravitational waves (GW), celestial objects, and cosmological data
that is the subject of this paper.

Since GW170817 was observed by LIGO and Virgo (Abbott et al., 2017a), the
situation surrounding the observational constraints of modified gravity has
changed dramatically. The simultaneous observation of GRB170817 (Abbott et
al., 2017b), a Gamma Ray burst, confirmed that the speed of GWs matches the
speed of electromagnetic waves with high accuracy, ruling out various scalar-
tensor theories that change the speed of GWs at low redshifts (Lombriser &
Taylor, 2016; Lombriser & Lima, 2017; Creminelli & Vernizzi, 2017; Sakstein &
Jain, 2017; Ezquiaga & Zumalacárregui, 2017; Baker et al., 2017; Langlois et
al., 2018). Creminelli et al. (2018, 2019) pointed out that a subset of DHOST
theories leads to the decay of gravitational waves, resulting in further tight
constraints on DHOST theories. However, the theory of gravity considered in
that paper, the class I DHOST theory (Langlois & Noui, 2016; Crisostomi et
al., 2016; Ben Achour et al., 2016b), still survives and can modify gravity in
cosmology without pathological instability (de Rham & Matas, 2016; Langlois et
al., 2017; Amendola et al., 2018). Furthermore, de Rham & Melville (2018)
showed that such cosmological scalar-tensor theories, which predict the speed
of GWs to be different from the speed of light, break down on high energy
scales ($\sim 10^{2}\,{\rm Hz}$) seen in neutron star mergers, indicating that
the constraints from GW observations may not necessarily apply to cosmological
scales. Therefore, it is essential to test modified gravity theories
independently at various energy scales, such as the GW and cosmological
scales.

DHOST theories generally have characteristic non-linear effects that violate
the Vainshtein screening mechanism inside any gravitational source (Kobayashi
et al., 2015; Koyama & Sakstein, 2015; Crisostomi & Koyama, 2018; Langlois et
al., 2018; Dima & Vernizzi, 2018; Hirano et al., 2019b; Crisostomi et al.,
2019). As an alternative to the cosmological constant, scalar-tensor theories
must give an ${\cal O}(1)$ modification from GR at cosmological scales, but at
small scales, they must satisfy tests in weakly gravitational regions such as
the solar system. The Vainshtein screening mechanism (for a review, Babichev &
Deffayet, 2013), universally found in scalar-tensor theories, is a typical
mechanism that satisfies these requirements, suppressing scalar interactions
and restoring standard gravity through non-linear effects. While Horndeski
theories allow for a natural implementation of the Vainshtein mechanism
(Kimura et al., 2012; Narikawa et al., 2013; Koyama, 2016), DHOST theories
partially violate it, allowing one to test DHOST theories by examining the
internal structure of objects such as Newtonian stars (Saito et al., 2015;
Sakstein, 2015a, b; Jain et al., 2016; Sakstein et al., 2017a; Saltas et al.,
2018; Saltas & Christensen-Dalsgaard, 2022), Neutron stars (Babichev et al.,
2016; Sakstein et al., 2017b), and galaxy clusters (Sakstein et al., 2016;
Salzano et al., 2017). The Vainshtein radius, the maximum scale at which the
Vainshtein mechanism works, is estimated to be ${\cal O}(100)\,\textrm{pc}$
for the Sun and ${\cal O}(1)\,\textrm{Mpc}$ for a galaxy cluster.

DHOST theories predict a characteristic gravitational non-linear effect on
even cosmological scales exceeding tens of Mpc. That is, DHOST theories
violate the consistency relation for LSS (Crisostomi et al., 2020;
Lewandowski, 2020) (see also Hirano et al., 2018). The LSS consistency
relation (Peloso & Pietroni, 2013; Kehagias & Riotto, 2013; Creminelli et al.,
2013) is an analogue of the consistency relation originally proposed for
single-field inflation models (Maldacena, 2003; Creminelli & Zaldarriaga,
2004), which relates $n$-point statistics of cosmological fluctuations to
$(n-1)$-point statistics in a non-perturbative matter. It is valid in the
limit where the wavenumber of one of the $n$-points is hugely smaller than the
others. This consistency relation is because the equations that the
fluctuations obey are invariant under a Galilean transformation (Scoccimarro &
Frieman, 1996; Creminelli et al., 2013). In particular, in the so-called
equal-time consistency relation, the Galilean transformation eliminates the
large-scale flow of matter and thus cancels all non-linear contributions when
calculating the $n$-point statistics. This behavior is also known as infrared
(IR) cancellation (Jain & Bertschinger, 1996; Scoccimarro & Frieman, 1996;
Kehagias & Riotto, 2013; Peloso & Pietroni, 2013; Sugiyama & Futamase, 2013;
Sugiyama & Spergel, 2014; Blas et al., 2013, 2016b; Lewandowski & Senatore,
2017). On the other hand, the LSS consistency relation breaks down when
considering multiple fluids (Tseliakhovich & Hirata, 2010; Yoo et al., 2011;
Bernardeau et al., 2012, 2013; Peloso & Pietroni, 2014; Creminelli et al.,
2014a; Lewandowski et al., 2015; Slepian & Eisenstein, 2017) or primordial
non-Gaussianities (Berezhiani & Khoury, 2014; Valageas et al., 2017; Esposito
et al., 2019; Goldstein et al., 2022), or when the equivalence principle
breaks (Creminelli et al., 2014b). DHOST theories have a structure similar to
that of multiple fluids, and on large scales, the Galilean transformation
cannot make the relative velocity between the scalar field and matter vanish
(for details, see Crisostomi et al., 2020; Lewandowski, 2020). As a result,
DHOST theories violate the LSS consistency relation.

Our interest in this paper is to constrain DHOST theories on cosmological
scales, i.e., ${\cal O}(10-100)\,\textrm{Mpc}$ scales. However, studies using
cosmological data to constrain DHOST theories are still limited (Hirano et
al., 2019a; Traykova et al., 2019; Peirone et al., 2019; Hiramatsu, 2022). On
the other hand, many papers on Horndeski theories have used cosmological data
to constrain the model (Okada et al., 2013; Barreira et al., 2014; Bellini et
al., 2016; Mueller et al., 2018; Kreisch & Komatsu, 2018; Arai & Nishizawa,
2018; Noller & Nicola, 2020, 2019; Raveri, 2020; Melville & Noller, 2020;
Perenon et al., 2019; Noller, 2020). Therefore, exploring new cosmological
methods for constraining DHOST theories is of great significance.

###  1.4 Constraints on modified gravity theories using galaxy two-point
statistics

The logarithmic growth rate function $f$ of dark matter fluctuations, measured
using redshift-space distortions (RSD; Kaiser, 1987), plays an important role
in constraining modified gravity theories in the late-time Universe. In the
power spectrum analysis, we cannot measure the growth rate function by itself,
but usually, by the combination $f\sigma_{8}=d\ln\sigma_{8}/d\ln a$ (Song &
Percival, 2009; Percival & White, 2009) using $\sigma_{8}$ representing the
rms of matter fluctuations on the $8\,h^{-1}\,{\rm Mpc}$ scale. For example,
the most recent observations, BOSS and eBOSS, measured $f\sigma_{8}$ with a
precision of $\sim 5\%$ in the redshift range $0.2<z<1.0$ (Alam et al.,
2021a).

One concern is to test modified gravity theories directly using existing
$f\sigma_{8}$ measurements. The standard practice is constructing a model of
the non-linear galaxy power spectrum assuming GR, then using that model to
measure $f\sigma_{8}$ from data up to the mildly non-linear region $(k\sim
0.2\,h\,{\rm Mpc^{-1}})$ (for recent studies, e.g., d’Amico et al., 2020;
Ivanov et al., 2020; Lange et al., 2022; Kobayashi et al., 2022; Yuan et al.,
2022). Therefore, it is worth noting that many existing analysis results using
galaxy data up to the non-linear region only verify the consistency of GR.
Thus, to test the gravity theory by consistently considering both linear and
non-linear effects, a power spectrum model that considers non-linear effects
specific to the modified gravity theory of interest is necessary. Several
studies have been done on this for various modified gravity theories (Koyama
et al., 2009; Taruya et al., 2014a; Taruya et al., 2014b; Takushima et al.,
2015; Bellini & Zumalacarregui, 2015; Taruya, 2016; Barreira et al., 2016;
Bose & Koyama, 2016; Cusin et al., 2018a, b; Bose et al., 2017; Bose et al.,
2018; Aviles et al., 2018; Hernández-Aguayo et al., 2019; Cataneo et al.,
2019; Valogiannis et al., 2020; Valogiannis & Bean, 2019; Bose et al., 2020b).
However, only one study constrained the theory from actual galaxy data using a
galaxy power spectrum model that consistently includes the non-linear effects
arising from modified gravity (Song et al., 2015a), where the authors focused
on $f(R)$ gravity (Hu & Sawicki, 2007).

In particular, Hirano et al. (2020) pointed out that in DHOST theories, even
the next-order solutions of the power spectrum in perturbation theory, the so-
called one-loop solutions, are challenging to perform physically meaningful
theoretical calculations due to the divergence of the wavenumber integral in
the ultraviolet (UV) region. Therefore, the modelling of non-linear power
spectra in DHOST theories is still highly uncertain.

###  1.5 Developments in the study of galaxy three-point statistics

A more straightforward way to investigate the non-linearity of scalar-tensor
theories is to use three-point statistics of cosmological fluctuations, i.e.,
the 3PCF or the bispectrum. The reason is that, on large scales, the three-
point statistics consist of a combination of second-order and linear-order
dark matter fluctuations. The second-order fluctuations depend on two wave
vectors in Fourier space and can be decomposed into three components using the
angle between the two wave vectors: monopole (growth), dipole (shift), and
quadrupole (tidal force) (Schmittfull et al., 2015). For example, Horndeski
theories deviate only the coefficient of the tidal term from GR while keeping
the shift term among these three components (Bernardeau & Brax, 2011;
Takushima et al., 2014; Bartolo et al., 2013; Bellini et al., 2015; Burrage et
al., 2019). On the other hand, DHOST theories change both the shift and tidal
terms (Hirano et al., 2018; Crisostomi et al., 2020; Lewandowski, 2020), and
this change in the shift term leads to a violation of the LSS consistency
relation (Crisostomi et al., 2020; Lewandowski, 2020). In addition to scalar-
tensor theories, there has been much researches on higher-order statistics in,
for example, $f(R)$ gravity theory (Gil-Marín et al., 2011; Borisov & Jain,
2009; Hellwing et al., 2013; Bose & Taruya, 2018; Bose et al., 2020a). Several
observational proposals have been made to test modified gravity theories using
cosmological three-point statistics, such as galaxy clustering (Yamauchi et
al., 2017b; Yamauchi & Sugiyama, 2022), weak lensing (Dinda, 2018; Munshi et
al., 2020b, a; Munshi & McEwen, 2020), and CMB lensing (Namikawa et al., 2018;
Namikawa et al., 2019), but none have been applied to actual observational
data yet.

In the context of the galaxy three-point statistics, 3PCF resolves the
degeneracy between the linear bias $b_{1}$ and $\sigma_{8}$ and allows us to
directly study the evolution of dark matter density fluctuations apart from
the RSD effect (Fry, 1994; Frieman & Gaztanaga, 1994; Matarrese et al., 1997;
Verde et al., 1998; Gaztañaga & Scoccimarro, 2005; Sefusatti et al., 2006;
Greig et al., 2013; Hoffmann et al., 2015; Samushia et al., 2021).
Furthermore, many previous studies have proposed to constrain primordial non-
Gaussianities from the galaxy three-point statistics (Fry & Scherrer, 1994;
Verde et al., 2000; Scoccimarro et al., 2004; Sefusatti & Komatsu, 2007;
Sefusatti, 2009; Sefusatti et al., 2010; Liguori et al., 2010; Desjacques &
Seljak, 2010; Sefusatti et al., 2012; Scoccimarro et al., 2012; Alvarez et
al., 2014; Tellarini et al., 2015, 2016; Welling et al., 2016; Yamauchi et
al., 2017a; Karagiannis et al., 2018; Bharadwaj et al., 2020; Moradinezhad
Dizgah et al., 2021; Shirasaki et al., 2021; Coulton et al., 2023; Karagiannis
et al., 2022). Recently, as in the case of the galaxy two-point statistics
(e.g., Matsubara, 2004; Taruya et al., 2011), the anisotropic component of the
galaxy three-point statistics induced by the RSD effect and the Alcock-
Paczyński (AP) effect (Alcock & Paczyński, 1979) has attracted much attention,
and its cosmological utility has been actively studied (Song et al., 2015b;
Gagrani & Samushia, 2017; Yankelevich & Porciani, 2019; Gualdi & Verde, 2020;
Mazumdar et al., 2020; Sugiyama et al., 2021; Agarwal et al., 2021; Rizzo et
al., 2023; Tsedrik et al., 2023).

Based on standard perturbation theory (SPT), many theoretical studies of the
galaxy three-point statistics have been conducted to calculate higher-order
non-linearities, redshift-space distortions, and bias effects, and the results
of these calculations have been tested for validity by comparing them with
measurements from N-body simulations (Peebles, 1980; Fry, 1984; Goroff et al.,
1986; Hivon et al., 1995; Scoccimarro, 1997; Scoccimarro et al., 1998; Jing &
Boerner, 1997; Scoccimarro et al., 1999; Scoccimarro, 2000; Barriga &
Gaztanaga, 2002; Barriga & Gaztañaga, 2002; Gaztañaga & Scoccimarro, 2005; Pan
et al., 2007; Marín et al., 2008; Guo & Jing, 2009; Pollack et al., 2012;
Lazanu et al., 2016; McCullagh et al., 2016; Lazanu et al., 2016; Lazanu &
Liguori, 2018; Hoffmann et al., 2018; Desjacques et al., 2018c; Child et al.,
2018; Eggemeier et al., 2019; Oddo et al., 2020; Eggemeier et al., 2021; Oddo
et al., 2021; Philcox et al., 2022). Other approaches have been widely used in
research, such as the halo models (Ma & Fry, 2000; Scoccimarro et al., 2001b;
Takada & Jain, 2003; Fosalba et al., 2005; Smith et al., 2008; Yamamoto et
al., 2017; Nan et al., 2018) and fitting formulas (Scoccimarro & Frieman,
1999; Scoccimarro & Couchman, 2001; Gil-Marín et al., 2012, 2014; Takahashi et
al., 2020). Beyond SPT, several improved perturbation theories have been
proposed. Rampf & Wong (2012) used a resummation method based on Lagrangian
perturbation theory. Baldauf et al. (2015b); Munshi & Regan (2017); Ivanov et
al. (2022) discussed some correction terms for SPT based on the effective
field theory of large-scale structure. Hashimoto et al. (2017) applied a
resummation method similar to the TNS model of the power spectrum (Taruya et
al., 2010). Kuruvilla & Porciani (2020) generalised the so-called streaming
model to higher-order statistics. Blas et al. (2016a); Ivanov & Sibiryakov
(2018) developed the time-sliced perturbation theory (TSPT) to resum the IR
modes of the bulk flow and describe the non-linear damping of Baryon Acoustic
Oscillations (BAOs; Peebles & Yu, 1970; Sunyaev & Zeldovich, 1970). Sugiyama
et al. (2021) constructed a new IR-resummed bispectrum model by adding a term
to the model proposed by TSPT.

The measurement of three-point statistics for galaxies, galaxy clusters, and
quasars has a long history. As a simple method, two-dimensional three-point
angular statistics have been observed from the dawn of the study of
cosmological three-point statistics to the present (Peebles & Groth, 1975;
Peebles, 1975; Groth & Peebles, 1977; Fry & Peebles, 1980; Fry & Seldner,
1982; Sharp et al., 1984; Jing & Zhang, 1989; Jing et al., 1991; Toth et al.,
1989; Frieman & Gaztañaga, 1999; Szapudi et al., 2001; de Carvalho et al.,
2020). Eventually, with the development of spectroscopic observations of
galaxies, three-dimensional three-point statistics have become the primary
targets observed in configuration space (Bean et al., 1983; Efstathiou &
Jedrzejewski, 1994; Hale-Sutton et al., 1989; Gott et al., 1991; Jing &
Borner, 1998; Jing & Börner, 2004; Kayo et al., 2004; Wang et al., 2004;
Gaztañaga et al., 2005; Pan & Szapudi, 2005; Nichol et al., 2006; Kulkarni et
al., 2007; Gaztañaga et al., 2009; Marin, 2011; McBride et al., 2011a, b;
Marin et al., 2013; Guo et al., 2014; Moresco et al., 2017a; Slepian et al.,
2017a, b; Moresco et al., 2017b, 2021) or in Fourier space (Baumgart & Fry,
1991; Scoccimarro et al., 2001a; Feldman et al., 2001; Verde et al., 2002;
Nishimichi et al., 2007; Gil-Marín et al., 2015a, b; Gil-Marín et al., 2017a;
Pearson & Samushia, 2018; Sugiyama et al., 2019; Philcox & Ivanov, 2022;
Cabass et al., 2022b; D’Amico et al., 2022a; Cabass et al., 2022a; D’Amico et
al., 2022b). As another approach, Chiang et al. (2015) measured the squeezed
limit bispectrum by splitting the observing region and measuring the position-
dependent power spectrum. Since the first measurement of the galaxy three-
point statistics by Peebles & Groth (1975), the three-point statistic
measurement has long been limited to measuring only certain scale-dependence
of the three-point statistics. However, it is now possible to perform
cosmological analysis using the information on the full shape of galaxy three-
point statistics at cosmological scales $(\sim 100\,h^{-1}\,{\rm Mpc})$.

In recent years, cosmological analysis of the three-point statistics of
galaxies has made remarkable progress, mostly focusing on the isotropic
component, i.e., _monopole_ , of the three-point statistics. Slepian et al.
(2017b) and Pearson & Samushia (2018) reported the detection of the BAO signal
through the monopole 3PCF and the monopole bispectrum, respectively. Gil-Marín
et al. (2017a); d’Amico et al. (2020); Philcox & Ivanov (2022) performed a
joint analysis of the monopole and quadrupole power spectra and the monopole
bispectrum to constrain the cosmological parameters of interest. Cabass et al.
(2022b); D’Amico et al. (2022a); Cabass et al. (2022a) constrained primordial
non-Gaussianities using the monopole bispectrum.

The anisotropic components, i.e., _quadrupole_ and _hexadecapole_ , of the
galaxy three-point statistics have been the subject of pretty limited studies
of measurements and cosmological analyses from actual galaxy data. Sugiyama et
al. (2019) reported the first detection of the quadrupole bispectrum signal at
the $14\sigma$ level from the BOSS DR12 galaxies. Sugiyama et al. (2021)
performed an anisotropic BAO analysis using the monopole and quadrupole
components of the 2PCF and 3PCF for the MultiDark-Patchy mock catalogues
(Patchy mocks; Kitaura et al., 2016) reproducing the BOSS galaxy distribution,
showing the improvement of the Hubble parameter constraint by $\sim 30\%$
compared to the 2PCF-only analysis result. D’Amico et al. (2022b) performed
the first joint analysis of the monopole and quadrupole components of the
power and bispectra measured from the BOSS DR12 galaxy data. More recently,
Ivanov et al. (2023) presented the results of an anisotropic bispectrum
analysis including quadrupole and hexadecapole components measured from the
BOSS DR12 data.

###  1.6 Goal of this paper

The primary goal of this paper is to use the 3PCF of galaxies to perform a
consistent cosmological analysis that constrains DHOST theories and their
subclass, Horndeski theories, while accounting for linear and non-linear
effects. To this end, Yamauchi & Sugiyama (2022) pointed out that the
parameters characterising non-linear density fluctuations in DHOST theories
degenerate with the non-linear bias parameter, so measuring the non-linear
velocity field due to the RSD effect is essential. In addition, the authors
proposed a simple parameterisation scheme that characterises the time
evolution of the scale dependence of the non-linear velocity field to
facilitate the combined analysis of galaxy samples at different redshifts.
Specifically, the time evolution of the shift and tidal terms of the second-
order velocity field is represented by $\xi_{\rm s}$ and $\xi_{\rm t}$,
respectively, where $\xi_{\rm s}=0$ and $\xi_{\rm t}=15/1144$ in GR. Following
the suggestion of Yamauchi & Sugiyama (2022), we apply the joint analysis
method of the anisotropic 2PCF and 3PCF of galaxies established by Sugiyama et
al. (2021) to BOSS Data Release 12 galaxies (Eisenstein et al., 2011; Bolton
et al., 2012; Dawson et al., 2013; Alam et al., 2015) to constrain these
$\xi_{\rm s}$ and $\xi_{\rm t}$ parameters.

When we need to use values of fiducial cosmology parameters in our analysis,
we adopt a flat $\Lambda$CDM model with the following parameters: matter
density $\Omega_{\rm m0}=0.31$, Hubble constant $h\equiv H_{0}/(100\,{\rm
km\,s^{-1}\,Mpc^{-1}})=0.676$, baryon density $\Omega_{\rm b0}h^{2}=0.022$,
and spectral tilt $n_{\rm s}=0.97$, which are the same as those used in the
final cosmological analysis in the BOSS project (Alam et al., 2017) and
consistent with the best-fit values in Planck 2018 (Aghanim et al., 2020). We
adopt a value for the total neutrino mass of $\sum m_{\nu}=0.06\,{\rm eV}$
close to the minimum allowed by neutrino oscillation experiments. We use these
fiducial parameters to estimate the distance to galaxies from the observed
redshift of each galaxy and to calculate the shape of the linear matter power
spectrum at the redshifts of interest with CLASS (Blas et al., 2011).

##  2 DHOST theories

In this section, we briefly review the analytic expressions of DHOST theories.
Section 2.1 introduces the class I DHOST theory and the perturbative solutions
of the density and velocity fields of dark matter and galaxies solved up to
the second-order in that theory. In Eqs. (2.1)-(2.14) of this subsection, we
adopt the expressions and notations given by Hirano et al. (2018). Section 2.2
discusses the limitation of the assumptions adopted to derive the perturbative
solutions used in this paper.

###  2.1 Density and velocity fluctuations in DHOST theories

We begin by summarising the theoretical models we will investigate in this
paper and the assumptions used to derive those models.

  1. 1.

Gravity theory is a subclass of quadratic DHOST theories, the class I DHOST
theory (Crisostomi et al., 2016), which encompasses Horndeski and Beyond
Horndeski theories and is free from the instabilities of a cosmological
background (de Rham & Matas, 2016; Langlois et al., 2017).

  2. 2.

Matter is cold dark matter (CDM) that can be described as a pressureless
perfect fluid without vorticity (Bernardeau et al., 2002).

  3. 3.

Matter is minimally coupled to gravity, and the effects of the DHOST gravity
appear only through the gravitational potential.

  4. 4.

When solving the equations of motion of metric tensor and scalar fields in
DHOST theories, the quasi-static approximation (e.g., Pace et al., 2021) is
used. Then, the gravitational potential is determined by a modified Poisson
equation (Hirano et al., 2018; Crisostomi et al., 2020; Lewandowski, 2020;
Hirano et al., 2020).

  5. 5.

Statistical properties of the CDM fluctuations are those derived in the
standard theory of inflation, which satisfy the following properties:
adiabaticity, negligibly weak non-Gaussianity, nearly scale-free, statistical
homogeneity, statistical isotropy, and statistical parity symmetry.

  6. 6.

Galaxy biases are assumed to be present only in the density field, and three
biases are considered: linear bias $b_{1}$, second-order local bias $b_{2}$,
and second-order non-local bias (tidal bias) $b_{s^{2}}$ (for a review, see
e.g., Saito et al., 2014; Desjacques et al., 2018b). Any bias effects related
to higher-order derivatives and the velocity field of the galaxy are ignored.

The action of quadratic DHOST theories is given by (Langlois & Noui, 2016;
Crisostomi et al., 2016)

$\displaystyle S_{\rm DHOST}$ $\displaystyle=$ $\displaystyle\int
d^{4}x\sqrt{-g}\Big{[}{\cal G}_{2}(\phi,X)-{\cal G}_{3}(\phi,X)\Box\phi+{\cal
F}(\phi,X)R$ (2.1) $\displaystyle+$ $\displaystyle
a_{1}\phi_{\mu\nu}\phi^{\mu\nu}+a_{2}(\Box\phi)^{2}+a_{3}(\Box\phi)\phi^{\mu}\phi_{\mu\nu}\phi^{\nu}$
$\displaystyle+$ $\displaystyle
a_{4}\phi^{\mu}\phi_{\mu\rho}\phi^{\rho\nu}\phi_{\nu}+a_{5}(\phi^{\mu}\phi_{\mu\nu}\phi^{\nu})^{2}\Big{]},$

where $\phi_{\mu}=\nabla_{\mu}\phi$,
$\phi_{\mu\nu}=\nabla_{\mu}\nabla_{\nu}\phi$, $X=-\phi_{\mu}\phi^{\mu}/2$, and
$a_{i}=a_{i}(\phi,X)$ for $i=1,\dots,5$. The functions $a_{i}$ ($i=1,\dots,5$)
satisfy the degeneracy condition given by (Crisostomi et al., 2016) to avoid
the Ostrogradsky ghost (Ostrogradsky, 1850; Woodard, 2015).

The density perturbation $\delta$ and velocity field v of dark matter follow
the equations of a pressureless perfect fluid without vorticity:

$\displaystyle\dot{\delta}(\textbf{{x}})+a^{-1}\partial_{i}\left((1+\delta(\textbf{{x}}))v^{i}(\textbf{{x}})\right)$
$\displaystyle=$ $\displaystyle 0,$
$\displaystyle\dot{\theta}(\textbf{{x}})+H\theta(\textbf{{x}})+a^{-1}\partial_{i}\left(v^{j}(\textbf{{x}})\partial_{j}v^{i}(\textbf{{x}})\right)$
$\displaystyle=$ $\displaystyle-a^{-1}\partial^{2}\Phi(\textbf{{x}}),$ (2.2)

where $a$ and $H=\dot{a}/a$ respectively denote the scale factor and the
Hubble parameter, and $\theta=\partial_{i}v^{i}$ is the divergence of the
velocity field. Because of no vorticity, the velocity field is represented as
$v^{i}=(\partial_{i}/\partial^{2})\theta$. The gravitational potential $\Phi$
is determined by the following modified Poisson equation (Hirano et al.,
2018):

$\displaystyle\frac{\partial^{2}\Phi(\textbf{{x}})}{a^{2}H^{2}}=\kappa\delta(\textbf{{x}})+\nu\frac{\dot{\delta}(\textbf{{x}})}{H}+\mu\frac{\ddot{\delta}(\textbf{{x}})}{H^{2}}+\frac{\partial^{2}S_{\Phi}^{\rm
NL}(\textbf{{x}})}{a^{2}H^{2}},$ (2.3)

where $\kappa$, $\nu$, and $\mu$ are functions that depend only on time, and
$S_{\Phi}^{\rm NL}$ is a non-linear source term obtained from the equation of
motion of the scalar field.

To solve the above equations, we expand all the fluctuations as follows:
$X=\sum_{n}X_{n}$, where $X=\\{\delta,\theta,\Phi,S_{\Phi}^{\rm NL}\\}$, and
$X_{n}={\cal O}(\delta_{1}^{n})$. Then, the non-linear source $S_{\Phi}^{\rm
NL}$ up to the second-order is given by

$\displaystyle\frac{\partial^{2}S_{\Phi,1}^{\rm
NL}(\textbf{{x}})}{a^{2}H^{2}}$ $\displaystyle=$ $\displaystyle 0,$
$\displaystyle\frac{\partial^{2}S_{\Phi,2}^{\rm
NL}(\textbf{{x}})}{a^{2}H^{2}}$ $\displaystyle=$
$\displaystyle\tau_{\alpha}W_{\alpha}(\textbf{{x}})-\tau_{\gamma}W_{\gamma}(\textbf{{x}}),$
(2.4)

where

$\displaystyle W_{\alpha}(\textbf{{x}})$ $\displaystyle=$
$\displaystyle\left[\delta_{1}(\textbf{{x}})\right]^{2}+\left[\frac{\partial_{i}}{\partial^{2}}\delta_{1}(\textbf{{x}})\right]\left[\partial_{i}\delta_{1}(\textbf{{x}})\right],$
$\displaystyle W_{\gamma}(\textbf{{x}})$ $\displaystyle=$
$\displaystyle\left[\delta_{1}(\textbf{{x}})\right]^{2}-\left[\frac{\partial_{i}\partial_{j}}{\partial^{2}}\delta_{1}(\textbf{{x}})\right]^{2}.$
(2.5)

The evolution of the density perturbation follows

$\displaystyle\ddot{\delta}(\textbf{{x}})+(2+\varsigma)H\dot{\delta}(\textbf{{x}})-\frac{3}{2}\Omega_{\rm
m}\Xi H^{2}\delta(\textbf{{x}})=H^{2}S_{\delta}^{\rm NL}(\textbf{{x}}),$ (2.6)

where $\varsigma=(2\mu-\nu)/(1-\mu)$, $(3/2)\Omega_{\rm m}\Xi=\kappa/(1-\mu)$,
and $S_{\delta}^{\rm NL}$ is a non-linear source of the density perturbation,
vanishing at linear order and given at second-order by

$\displaystyle S_{\delta,2}^{\rm
NL}=S_{\alpha}W_{\alpha}(\textbf{{x}})-S_{\gamma}W_{\gamma}(\textbf{{x}})$
(2.7)

with

$\displaystyle(1-\mu)S_{\alpha}$ $\displaystyle=$ $\displaystyle
2f^{2}+\frac{3}{2}\Omega_{\rm m}\Xi-\varsigma f+\tau_{\alpha},$
$\displaystyle(1-\mu)S_{\gamma}$ $\displaystyle=$ $\displaystyle
f^{2}+\tau_{\gamma}.$ (2.8)

Once the solution of $\delta$ is obtained, the solution of $\theta$ is also
derived from the continuity equation in Eq. (2.2). In Fourier space101010Our
convention for the Fourier transform is
$\displaystyle\widetilde{f}(\textbf{{k}})=\int
d^{3}xe^{-i\textbf{{k}}\cdot\textbf{{x}}}f(\textbf{{x}}).$ , Eqs. (2.2) and
(2.3) determine $\delta_{n}(\textbf{{k}})$ and $\theta_{n}(\textbf{{k}})$ in
terms of the linear density fluctuations to be:

$\displaystyle\widetilde{\delta}_{n}(\textbf{{k}})$ $\displaystyle=$
$\displaystyle\int\frac{d^{3}p_{1}}{(2\pi)^{3}}\cdots\int\frac{d^{3}p_{n}}{(2\pi)^{3}}(2\pi)^{3}\delta_{\rm
D}(\textbf{{k}}-\textbf{{p}}_{[1n]})$ $\displaystyle\times$ $\displaystyle
F^{(\rm
m)}_{n}(\textbf{{p}}_{1},\dots,\textbf{{p}}_{2})\delta_{1}(\textbf{{p}}_{1})\cdots\delta_{1}(\textbf{{p}}_{n}),$
$\displaystyle\widetilde{\theta}_{n}(\textbf{{k}})$ $\displaystyle=$
$\displaystyle-
aHf\int\frac{d^{3}p_{1}}{(2\pi)^{3}}\cdots\int\frac{d^{3}p_{n}}{(2\pi)^{3}}(2\pi)^{3}\delta_{\rm
D}(\textbf{{k}}-\textbf{{p}}_{[1n]})$ (2.9) $\displaystyle\times$
$\displaystyle G^{(\rm
m)}_{n}(\textbf{{p}}_{1},\dots,\textbf{{p}}_{2})\delta_{1}(\textbf{{p}}_{1})\cdots\delta_{1}(\textbf{{p}}_{n}),$

where $\textbf{{p}}_{[1n]}=\textbf{{p}}_{1}+\dots+\textbf{{p}}_{n}$, and
$\delta_{\rm D}$ is the delta function. The functions $F^{(\rm m)}_{2}$ and
$G^{(\rm m)}_{2}$ are kernel functions that characterise the gravitational
non-linear effects, and the superscript $(\rm m)$ stands for “matter”. In the
second-order, $F^{(\rm m)}_{2}$ and $G^{(\rm m)}_{2}$ are given by

$\displaystyle F^{(\rm
m)}_{2}(\textbf{{p}}_{1},\textbf{{p}}_{2})=\kappa_{\delta}\alpha_{\rm
s}(\textbf{{p}}_{1},\textbf{{p}}_{2})-\frac{2}{7}\lambda_{\delta}\gamma(\textbf{{k}}_{1},\textbf{{k}}_{2})$
$\displaystyle G^{(\rm
m)}_{2}(\textbf{{p}}_{1},\textbf{{p}}_{2})=\kappa_{\theta}\alpha_{\rm
s}(\textbf{{p}}_{1},\textbf{{p}}_{2})-\frac{4}{7}\lambda_{\theta}\gamma(\textbf{{k}}_{1},\textbf{{k}}_{2}),$
(2.10)

where

$\displaystyle\alpha_{\rm s}(\textbf{{k}}_{1},\textbf{{k}}_{2})$
$\displaystyle=$ $\displaystyle
1+(\hat{k}_{1}\cdot\hat{k}_{2})\frac{(k_{1}^{2}+k_{2}^{2})}{2k_{1}k_{2}},$
$\displaystyle\gamma(\textbf{{k}}_{1},\textbf{{k}}_{2})$ $\displaystyle=$
$\displaystyle 1-(\hat{k}_{1}\cdot\hat{k}_{2})^{2},$ (2.11)

and

$\displaystyle\kappa_{\theta}$ $\displaystyle=$ $\displaystyle
2\kappa_{\delta}\left[1+\frac{1}{2f}\frac{d\ln\kappa_{\delta}}{d\ln
a}\right]-1,$ $\displaystyle\lambda_{\theta}$ $\displaystyle=$
$\displaystyle\lambda_{\delta}\left[1+\frac{1}{2f}\frac{d\ln\lambda_{\delta}}{d\ln
a}\right].$ (2.12)

The evolutions of $\kappa_{\delta}$ and $\lambda_{\delta}$ follow

$\displaystyle\ddot{\kappa}_{\delta}+[4f+(2+\varsigma)]H\dot{\kappa}_{\delta}+H^{2}\left(2f^{2}+\frac{3}{2}\Omega_{\rm
m}\Xi\right)\kappa_{\delta}$ $\displaystyle=H^{2}S_{\alpha},$ (2.13)
$\displaystyle\ddot{\lambda}_{\delta}+[4f+(2+\varsigma)]H\dot{\lambda}_{\delta}+H^{2}\left(2f^{2}+\frac{3}{2}\Omega_{\rm
m}\Xi\right)\lambda_{\delta}$ $\displaystyle=\frac{7}{2}H^{2}S_{\gamma}.$
(2.14)

Since the galaxy density field is a biased quantity, we assume the linear bias
$b_{1}$, the second-order local bias $b_{2}$, and the second-order tidal bias
$b_{s^{2}}$ as the bias parameters that describe the galaxy density
fluctuation up to second order (e.g., Desjacques et al., 2018b):

$\displaystyle\delta^{(\rm
g)}(\textbf{{x}})=b_{1}\delta(\textbf{{x}})+\frac{b_{2}}{2}\left[\delta(\textbf{{x}})\right]^{2}+b_{s^{2}}[s_{ij}]^{2},$
(2.15)

where the superscript $(\rm g)$ stands for “galaxy”, and $[s_{ij}]^{2}$ is
given by

$\displaystyle[s_{ij}]^{2}=\left[\frac{\partial_{i}\partial_{j}}{\partial^{2}}\delta(\textbf{{x}})\right]^{2}-\frac{1}{3}[\delta(\textbf{{x}})]^{2}.$
(2.16)

Then, the second-order kernel functions for galaxies are given by

$\displaystyle F_{2}^{\rm(g)}$ $\displaystyle=$ $\displaystyle
b_{1}F_{2}^{\rm(m)}(\textbf{{p}}_{1},\textbf{{p}}_{2})+\frac{1}{2}b_{2}+b_{s^{2}}\left[(\hat{p}_{1}\cdot\hat{p}_{2})^{2}-\frac{1}{3}\right],$
$\displaystyle G_{2}^{\rm(g)}$ $\displaystyle=$ $\displaystyle
G_{2}^{\rm(m)}(\textbf{{p}}_{1},\textbf{{p}}_{2}).$ (2.17)

The RSD effect shifts the observed position of galaxies $\textbf{{x}}_{\rm
red}$ from their real-space position $\textbf{{x}}^{\prime}$ due to the
peculiar velocity of galaxies along the line-of-sight (LOS) direction:

$\displaystyle\textbf{{x}}_{\rm
red}(\textbf{{x}}^{\prime})=\textbf{{x}}^{\prime}+\frac{\textbf{{v}}(\textbf{{x}}^{\prime})\cdot\hat{n}}{aH}\hat{n},$
(2.18)

where $\hat{n}$ is a unit vector pointing to the galaxy from the origin. The
observed galaxy density fluctuation is then distorted along the LOS direction
as follows:

$\displaystyle\delta_{\rm s}^{(\rm g)}(\textbf{{x}})=\int
d^{3}x^{\prime}\left(1+\delta^{(\rm
g)}(\textbf{{x}}^{\prime})\right)\delta_{\rm
D}\left(\textbf{{x}}-\textbf{{x}}_{\rm red}(\textbf{{x}}^{\prime})\right)-1.$
(2.19)

In Fourier space, the $n$-th order solution of $\delta_{\rm s}^{(\rm g)}$ is
represented as

$\displaystyle\widetilde{\delta}_{{\rm s},\,n}^{(\rm g)}(\textbf{{k}})$
$\displaystyle=$
$\displaystyle\int\frac{d^{3}p_{1}}{(2\pi)^{3}}\cdots\int\frac{d^{3}p_{n}}{(2\pi)^{3}}(2\pi)^{3}\delta_{\rm
D}(\textbf{{k}}-\textbf{{p}}_{[1n]})$ (2.20) $\displaystyle\times$
$\displaystyle
Z_{n}(\textbf{{p}}_{1},\dots,\textbf{{p}}_{2})\delta_{1}(\textbf{{p}}_{1})\cdots\delta_{1}(\textbf{{p}}_{n}).$

The first and second-order kernel functions are given by (Scoccimarro et al.,
1999)

$\displaystyle Z_{1}$ $\displaystyle=$ $\displaystyle
b_{1}+f(\hat{p}\cdot\hat{n})^{2},$ $\displaystyle Z_{2}$ $\displaystyle=$
$\displaystyle
F_{2}^{(g)}(\textbf{{p}}_{1},\textbf{{p}}_{2})+f(\hat{k}\cdot\hat{n})^{2}G^{(g)}_{2}(\textbf{{p}}_{1},\textbf{{p}}_{2})$
(2.21) $\displaystyle+$
$\displaystyle\frac{f(\textbf{{k}}\cdot\hat{n})}{2}\left[\frac{(\hat{p}_{1}\cdot\hat{n})}{p_{1}}Z_{1}(\textbf{{p}}_{2})+\frac{(\hat{p}_{2}\cdot\hat{n})}{p_{2}}Z_{1}(\textbf{{p}}_{1})\right],$

where $\textbf{{k}}=\textbf{{p}}_{1}+\textbf{{p}}_{2}$. In the rest of this
paper, we focus only on the galaxy density fluctuation with RSDs, so for
simplicity of notation, we refer to it simply as $\delta$ instead of
$\delta_{\rm s}^{(\rm g)}$. We also omit the angle-dependence $\hat{n}$ of any
function that includes RSDs.

At the leading-order in perturbation theory, the redshift-space power spectrum
and bispectrum are represented as

$\displaystyle P(\textbf{{k}})$ $\displaystyle=$
$\displaystyle[Z_{1}(\textbf{{k}})]^{2}P_{\rm lin}(k),$ $\displaystyle
B(\textbf{{k}}_{1},\textbf{{k}}_{2},\textbf{{k}}_{3})$ $\displaystyle=$
$\displaystyle
2Z_{2}(\textbf{{k}}_{1},\textbf{{k}}_{2})Z_{1}(\textbf{{k}}_{1})Z_{1}(\textbf{{k}}_{2})P_{\rm
lin}(k_{1})P_{\rm lin}(k_{2})$ (2.22) $\displaystyle+$ $\displaystyle\mbox{2
perms.},$

where $\textbf{{k}}_{1}+\textbf{{k}}_{2}+\textbf{{k}}_{3}=0$, and $P_{\rm
lin}$ is the linear matter power spectrum. In what follows, we omit the
$\textbf{{k}}_{3}$-dependence of the bispectrum for notational simplicity:
$B(\textbf{{k}}_{1},\textbf{{k}}_{2})=B(\textbf{{k}}_{1},\textbf{{k}}_{2},\textbf{{k}}_{3}=-\textbf{{k}}_{1}-\textbf{{k}}_{2})$.

Finally, we conclude this subsection by summarising the key points about
galaxy fluctuations from a theoretical point of view. First, in the case of
$\nu=\mu=\tau_{\alpha}=0$ in Eqs. (2.3) and (2.4), Horndeski theories are
recovered; a $\Lambda$CDM model additionally has $\kappa=(3/2)\Omega_{\rm
m}(z)$ and $\tau_{\gamma}=0$; in both Horndeski theories and $\Lambda$CDM,
$\kappa_{\delta}=\kappa_{\theta}=1$ from Eq. (2.13), and $\lambda_{\delta}$
and $\lambda_{\theta}$ are still time-dependent from Eq. (2.14); for the
approximation $f^{2}=\Omega_{\rm m}$ in $\Lambda$CDM,
$\lambda_{\delta}=\lambda_{\theta}=1$. Second, since the linear equation of
the density fluctuation (2.6) omits space-dependence as in the $\Lambda$CDM
case, under the assumption that the scalar field becomes to prevail during the
accelerated Universe, the shape of the linear matter power spectrum can be the
usual $\Lambda$CDM one determined in the matter-dominant era. In other words,
the characteristic scale-dependences in $\delta$ and $v$ due to scalar-tensor
theories appear only through the non-linear kernel functions $F^{(m)}_{n\geq
2}$ and $G^{(m)}_{n\geq 2}$. Third, the non-linear terms that appear in the
fluid equation in Eq. (2.2) and the Poisson equation in Eq. (2.3), such as
$\partial_{i}(\delta v^{i})$, $\partial_{i}(v^{j}\partial_{j}v^{i})$,
$W_{\alpha}$ and $W_{\gamma}$, become zero when the volume average or ensemble
average is calculated. Therefore, the resulting non-linear solutions satisfy
$\int d^{3}x\delta_{n}=\langle\delta_{n}\rangle=0$ and $\int
d^{3}x\theta_{n}=\langle\theta_{n}\rangle=0$ for $n\geq 2$, and the
corresponding kernel functions satisfy $F^{(\rm m)}_{n\geq 2}=G^{(\rm
m)}_{n\geq 2}=0$ when $\textbf{{p}}_{1}+\dots+\textbf{{p}}_{n}=0$ as known in
the case of $\Lambda$CDM. This condition partially breaks when the non-linear
bias effect is taken into account, resulting in $F_{2}^{(\rm
g)}(\textbf{{p}},-\textbf{{p}})\neq 0$ and $G_{2}^{(\rm
g)}(\textbf{{p}},-\textbf{{p}})=0$ (see Eq. (2.17)).

###  2.2 Limitation of our assumptions

This subsection discusses the possible cases where the assumptions adopted in
building the theoretical model in the previous subsection are violated,
introducing some previous studies. The following bullet labels correspond to
those in Section 2.1.

  1. 1.

Besides scalar-tensor theories, two other examples of modified gravity
theories have been widely studied in cosmology: the Hu-Sawicki model (Hu &
Sawicki, 2007) of $f(R)$ gravity (see Capozziello & Francaviglia, 2008;
Sotiriou & Faraoni, 2010, for reviews) and the normal branch of the 5D brane-
world Dvali-Gabadadze-Porrati model (nDGP; Dvali et al., 2000). These two
models have been investigated in detail by Alam et al. (2021b) as
representative targets in DESI. Focusing on the non-linear effects, the nDGP
model generates a scale dependence of the same form as Horndeski theories,
characterised by the function $\gamma(\textbf{{p}}_{1},\textbf{{p}}_{2})$
(2.11). On the other hand, the Hu-Sawicki $f(R)$ model produces a kernel
function different from the one predicted by scalar-tensor theories.
Specifically, in the modified Poisson equation of Eq. (2.3), $\kappa$ is
scale-dependent, resulting in the linear growth function that depends on the
wavenumber. In addition, the non-linear source $S_{\Phi}^{\rm NL}$ for the Hu-
Sawicki model also appears as a form that cannot be described by $W_{\alpha}$
and $W_{\gamma}$, unlike Eq. (2.4). Such non-linearities in the density field
specific to the Hu-Sawicki model have been studied by (Koyama et al., 2009;
Taruya, 2016) in the context of perturbation theory, and the model has been
tested by applying the theory to BOSS galaxy data (Song et al., 2015a).

  2. 2.

The effect of the relative velocity of baryons and CDM enters the galaxy
density fluctuation quadratically together with the corresponding bias
parameter (Dalal et al., 2010), thus modifying the shape of the measured
bispectrum. In particular, as in the case of the $\kappa_{\delta}$ parameter
in DHOST theories, it corrects the term in $F_{2}^{(\rm
g)}(\textbf{{p}}_{1},\textbf{{p}}_{2})$ that depends on
$(\hat{p}_{1}\cdot\hat{p}_{2})$ called the shift term (Yoo et al., 2011). The
relative velocity effect on galaxy clustering has been measured using the
galaxy power spectrum (Yoo & Seljak, 2013; Beutler et al., 2016) and 3PCF
(Slepian et al., 2018), but any signature has not yet been detected.

Although massive neutrinos can also change the shape of the bispectrum, the
results of simulations performed by Ruggeri et al. (2018) confirm that the CDM
component in the bispectrum is dominant; Interestingly, Kamalinejad & Slepian
(2020) has shown that the effect of neutrino corrections appears in the shift
term as well as the growth and tidal terms in the second-order velocity field
(3.16). Hence, the anisotropic 3PCF (or bispectrum) may help to constrain the
neutrino masses (see e.g., Saito et al., 2009; Levi & Vlah, 2016; Yoshikawa et
al., 2020).

  3. 3.

The case of non-minimally coupled scalar fields with CDM has already been the
subject of several studies in the context of cosmology (Kimura et al., 2018;
Chibana et al., 2019; Kase & Tsujikawa, 2020b; Chiba et al., 2020; Kase &
Tsujikawa, 2020a). For example, Kimura et al. (2018); Chibana et al. (2019)
have shown that in this case, the continuity equation (2.2) is modified, and
thus the relation between the density fluctuations in real and redshift
spaces, i.e. the Kaiser formula in linear theory (Kaiser, 1987), is also
modified.

  4. 4.

The quasi-static approximation breaks when the scale of interest is close to
the sound horizon scale. Even in GR, it is known that there are relativistic
corrections to $F_{2}^{(\rm g)}$ when approaching the horizon size (Tram et
al., 2016; Jolicoeur et al., 2017, 2018; Koyama et al., 2018; Castiblanco et
al., 2019; Umeh et al., 2019; Calles et al., 2020; de Weerd et al., 2020).

  5. 5.

Various possibilities have been proposed for how the initial conditions of
cosmic fluctuations predicted by inflation theory could affect observables.
One of the most critical examples relevant to this paper is the existence of
primordial non-Gaussianity, which breaks the LSS consistency relation
(Berezhiani & Khoury, 2014; Valageas et al., 2017; Esposito et al., 2019).

  6. 6.

Fujita & Vlah (2020) proposed a bias expansion formalism dubbed “Monkey bias”
based on the LSS consistency relation and showed that it is equivalent to the
existing bias expansion framework. In other words, in DHOST theories, which
violate the LSS consistency relation, the existing bias expansion we adopted
(2.15) may not be valid, and a new bias in the shift term of non-linear galaxy
density fluctuations, i.e., the shift bias parameter, may appear. Moreover,
the shift bias may also induce velocity bias effects.

In Section 9.12, we will discuss and clarify which parts of theories can be
tested with the anisotropic 3PCF, even in the presence of the shift and
velocity biases.

##  3 Theoretical models

This section describes how to calculate the theoretical models of multipole
2PCFs and 3PCFs. Section 3.1 summarises the decomposition formalism for the
anisotropic three-point statistics (bispectra and 3PCFs). Section 3.2
introduces the power and bispectrum models used to compute the 2PCF and 3PCF.
Section 3.3 discusses what parameters should be varied to perform the
cosmological analysis and shows the specific parameter dependence of the
bispectrum model we use. Section 3.4 reviews new parameters helpful in testing
DHOST theories proposed by Yamauchi & Sugiyama (2022) and their time
evolution. Section 3.5 discusses the limits of applying our theoretical models
of the 2PCF and 3PCF to the data analysis.

###  3.1 Decomposition formalisms of the 2PCF and 3PCF

We follow the decomposition formalism of redshift-space bispectra proposed by
Sugiyama et al. (2019) using the tri-polar spherical harmonics (TripoSH) as a
basis function. In that formalism, under statistical homogeneity, isotropy,
and parity-symmetry assumptions, we define the base function to expand the
bispectrum using three spherical harmonics $Y_{\ell m}$ as

$\displaystyle{\cal
S}_{\ell_{1}\ell_{2}\ell}(\hat{k}_{1},\hat{k}_{2},\hat{n})$ $\displaystyle=$
$\displaystyle\frac{4\pi}{h_{\ell_{1}\ell_{2}\ell}}\sum_{m_{1}m_{2}m}\left(\begin{smallmatrix}\ell_{1}&\ell_{2}&\ell\\\
m_{1}&m_{2}&m\end{smallmatrix}\right)$ (3.1) $\displaystyle\times$
$\displaystyle
Y_{\ell_{1}m_{1}}(\hat{k}_{1})Y_{\ell_{2}m_{2}}(\hat{k}_{2})Y_{\ell
m}(\hat{n}),$

where

$\displaystyle
h_{\ell_{1}\ell_{2}\ell}=\sqrt{\frac{(2\ell_{1}+1)(2\ell_{2}+1)(2\ell+1)}{4\pi}}\left(\begin{smallmatrix}\ell_{1}&\ell_{2}&\ell\\\
0&0&0\end{smallmatrix}\right),$ (3.2)

and the circle bracket with $6$ multipole indices, $(\dots)$, denotes the
Wigner-3j symbol. The bispectrum is then expanded as

$\displaystyle
B(\textbf{{k}}_{1},\textbf{{k}}_{2},\hat{n})=\sum_{\ell_{1}+\ell_{2}+\ell={\rm
even}}B_{\ell_{1}\ell_{2}\ell}(k_{1},k_{2}){\cal
S}_{\ell_{1}\ell_{2}\ell}(\hat{k}_{1},\hat{k}_{2},\hat{n}),$ (3.3)

and the corresponding multipole components are given by

$\displaystyle B_{\ell_{1}\ell_{2}\ell}(k_{1},k_{2})$ $\displaystyle=$
$\displaystyle 4\pi
h_{\ell_{1}\ell_{2}\ell}^{2}\int\frac{d^{2}\hat{k}_{1}}{4\pi}\int\frac{d^{2}\hat{k}_{2}}{4\pi}\int\frac{d^{2}\hat{n}}{4\pi}$
(3.4) $\displaystyle\times$ $\displaystyle{\cal
S}^{*}_{\ell_{1}\ell_{2}\ell}(\hat{k}_{1},\hat{k}_{2})B(\textbf{{k}}_{1},\textbf{{k}}_{2}).$

Since the bispectrum multipoles defined here are independent of the coordinate
system in which they are calculated, it is possible to compare theoretical
calculations with observations in different coordinate systems. Specifically,
we use the following coordinate system with $\hat{k}_{1}$ as the $z$-axis for
theoretical calculations:

$\displaystyle\hat{k}_{1}$ $\displaystyle=$ $\displaystyle\\{0,0,1\\}$
$\displaystyle\hat{k}_{2}$ $\displaystyle=$
$\displaystyle\\{\sin\theta_{k_{2}},0,\cos\theta_{k_{2}}\\}$
$\displaystyle\hat{n}$ $\displaystyle=$
$\displaystyle\\{\sin\theta\cos\varphi,\sin\theta\sin\varphi,\cos\theta\\}.$
(3.5)

On the other hand, when measuring the bispectrum from galaxy data, we use the
Cartesian coordinate and take the north pole as our $z$-axis (see Section
4.2).

We perform the expansion of the 3PCF in the same way as for the bispectrum.
The resulting 3PCF multipoles are related to $B_{\ell_{1}\ell_{2}\ell}$
through a two-dimensional Hankel transform:

$\displaystyle\zeta_{\ell_{1}\ell_{2}\ell}(r_{1},r_{2})$ $\displaystyle=$
$\displaystyle
i^{\ell_{1}+\ell_{2}}\int\frac{dk_{1}k_{1}^{2}}{2\pi^{2}}\int\frac{dk_{2}k_{2}^{2}}{2\pi^{2}}$
(3.6) $\displaystyle\times$ $\displaystyle
j_{\ell_{1}}(r_{1}k_{1})j_{\ell_{2}}(r_{2}k_{2})B_{\ell_{1}\ell_{2}\ell}(k_{1},k_{2}),$

where $j_{\ell}$ is the spherical Bessel function at the $\ell$-th order. This
relation means that $\zeta_{\ell_{1}\ell_{2}\ell}$ have in principle the same
information as $B_{\ell_{1}\ell_{2}\ell}$, facilitating the comparison of the
configuration-space and Fourier-space analyses.

Note that
$B_{\ell_{1}\ell_{2}\ell}(k_{1},k_{2})=B_{\ell_{2}\ell_{1}\ell}(k_{2},k_{1})$
and
$\zeta_{\ell_{1}\ell_{2}\ell}(r_{1},r_{2})=\zeta_{\ell_{2}\ell_{1}\ell}(r_{2},r_{1})$.
From this relation, when $\ell_{1}=\ell_{2}$, only $k_{1}\geq k_{2}$ and
$r_{1}\geq r_{2}$ need to be computed for the bispectrum and 3PCF,
respectively. Also, when $\ell>0$, only $\ell_{1}\geq\ell_{2}$ should be
considered.

In the case of the power spectrum, it is common to expand the power spectrum
using Legendre polynomial functions ${\cal L}_{\ell}$ (e.g., Hamilton, 1997):

$\displaystyle P(\textbf{{k}})=\sum_{\ell}P_{\ell}(k){\cal
L}_{\ell}(\hat{k}\cdot\hat{n}),$ (3.7)

and the corresponding multipole components of the 2PCF are given by

$\displaystyle\xi_{\ell}(r)$ $\displaystyle=$ $\displaystyle
i^{\ell}\int\frac{dkk^{2}}{2\pi^{2}}j_{\ell}(rk)P_{\ell}(k).$ (3.8)

This paper tests DHOST theories by measuring $\xi_{\ell}$ and
$\zeta_{\ell_{1}\ell_{2}\ell}$ from the BOSS galaxy data and comparing them
with the corresponding theoretical models. The index $\ell$ that is common for
both $\xi_{\ell}$ and $\zeta_{\ell_{1}\ell_{2}\ell}$ represents the
decomposition related to the RSD or AP effect, where $\ell=0$ means monopole,
$\ell=2$ quadrupole, and $\ell=4$ hexadecapole. Relativistic effects can
generate $\ell={\rm odd}$ components (e.g., McDonald, 2009; Desjacques et al.,
2018a; Clarkson et al., 2019), but we ignore them here. Furthermore, we also
ignore the $\ell=4$ modes; although the signal of the $\ell=4$ modes is too
small to be detected in the BOSS data, it should be taken into account in the
future as it helps to improve the constraints on the cosmological parameters
(Beutler et al., 2017; Sugiyama et al., 2019). Therefore, in this paper, we
focus on only two modes, $\ell=0$ and $\ell=2$. In particular, for the 3PCF,
we consider the first two terms of the monopole ($\zeta_{000}$ and
$\zeta_{110}$) and the first two terms of the quadrupole ($\zeta_{202}$ and
$\zeta_{112}$).

Finally, we discuss the relation with the widely used decomposition formalism
of the bispectrum proposed by Scoccimarro et al. (1999). As in Eq. (3.5), this
formalism decomposes the bispectrum by choosing the coordinate system with
$k_{1}$ as the $z$-axis and using the spherical harmonic function for the LOS
direction:
$B(\textbf{{k}}_{1},\textbf{{k}}_{2},\hat{n})=\sum_{LM}B_{LM}(\textbf{{k}}_{1},\textbf{{k}}_{2})Y_{LM}(\hat{n})$.
The relation between Scoccimarro et al. (1999)’s decomposition formalism and
our TripoSH decomposition has already been shown in Eq. (25) of Sugiyama et
al. (2019). According to the relation, $\zeta_{202}$ contains only $M=0$ mode
in Scoccimarro et al. (1999)’s formalism, while $\zeta_{112}$ further contains
the $M\neq 0$ modes in addition to the $M=0$ mode. The ability to handle the
$M\neq 0$ modes, including window function corrections (see Section 4.3), is
one advantage of our TripoSH decomposition formalism. For example, studies of
the quadrupole bispectrum using Scoccimarro et al. (1999)’s method have mainly
dealt only with the $M=0$ mode (D’Amico et al., 2022b). One reason is that the
correction formula for the window function effect is only given for the $M=0$
case (Pardede et al., 2022). Moreover, we show in Section 7 that $\zeta_{112}$
gives additional cosmological information to $\zeta_{202}$, pointing out the
importance of the $M\neq 0$ modes.

###  3.2 IR-resummed power spectrum and bispectrum models

In this paper, we focus on the 2PCF and 3PCF at scales above $80\,h^{-1}\,{\rm
Mpc}$ (Section 9), where we can ignore loop corrections arising from higher-
order non-linear effects. The power spectrum and bispectrum shapes can be
described at those scales by their leading solutions, the so-called tree-level
solutions (2.22). However, we need to consider the non-linear damping effect
of BAOs due to the linear gravity that shifts the position of galaxies.

The non-linear damping of BAO can be described by a large-scale bulk flow that
is position-independent in a given observed region (Eisenstein et al., 2007a;
Crocce & Scoccimarro, 2008; Matsubara, 2008; Sugiyama & Spergel, 2014; Baldauf
et al., 2015a), called the infra-red (IR) flow. In the limit where the IR flow
does not correlate with small-scale density fluctuations, based on the
Galilean invariance of the system of equations in the IR limit, all the
effects of the IR flow are cancelled out in equal-time $n$-point statistics
(Jain & Bertschinger, 1996; Scoccimarro & Frieman, 1996; Kehagias & Riotto,
2013; Peloso & Pietroni, 2013; Sugiyama & Futamase, 2013; Sugiyama & Spergel,
2014; Blas et al., 2013, 2016b; Lewandowski & Senatore, 2017). However, when
we deviate from such an extreme situation, we find a correlation between the
IR flow and the small-scale density field. By extracting this correlation in
the full perturbative order only for the BAO signal, it becomes possible to
describe the non-linear effects of BAOs. This kind of construction of
$n$-point statistics models is called the IR resummation method (Crocce &
Scoccimarro, 2008; Matsubara, 2008; Sugiyama & Spergel, 2014; Senatore &
Zaldarriaga, 2015; Baldauf et al., 2015a; Blas et al., 2016a; Senatore &
Trevisan, 2018; Ivanov & Sibiryakov, 2018; Lewandowski & Senatore, 2020;
Sugiyama et al., 2021). In this paper, we will use the IR resummed power and
bispectrum models given in Eqs. (3.9) and (3.12), even in DHOST theories that
break the IR cancellation, but we will mention the issues that may arise in
this case in Section 3.5.

For the power spectrum, we adopt the following IR-resummed model:

$\displaystyle P(\textbf{{k}})=\left[Z_{1}(\textbf{{k}})\right]^{2}\left[{\cal
D}^{2}(\textbf{{k}})P_{\rm w}(k)+P_{\rm nw}(k)\right],$ (3.9)

where $P_{\rm lin}$ is decomposed into two parts: the "no-wiggle (nw)" part
$P_{\rm nw}$ that is a smooth version of $P_{\rm lin}$ with the baryon
oscillations removed (Eisenstein & Hu, 1998), and the "wiggle (w)" part
defined as $P_{\rm w}=P_{\rm lin}-P_{\rm nw}$. The non-linear BAO degradation
is represented by the two-dimensional Gaussian damping factor derived from a
differential motions of Lagrangian displacements (Eisenstein et al., 2007a;
Crocce & Scoccimarro, 2008; Matsubara, 2008):

$\displaystyle{\cal
D}(\textbf{{k}})=\exp\left(-\frac{k^{2}(1-\mu^{2})\sigma_{\perp}^{2}+k^{2}\mu^{2}\sigma^{2}_{\parallel}}{2}\right),$
(3.10)

where $\mu=\hat{k}\cdot\hat{n}$. We compute the radial and transverse
components of smoothing parameters, $\sigma_{\perp}$ and $\sigma_{\parallel}$,
using the Zel’dovich approximation (Zel’Dovich, 1970; Crocce & Scoccimarro,
2008; Matsubara, 2008):

$\displaystyle\sigma^{2}_{\perp}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\int\frac{dp}{2\pi^{2}}P_{\rm lin}(p),$
$\displaystyle\sigma^{2}_{\parallel}$ $\displaystyle=$
$\displaystyle(1+f)^{2}\,\sigma^{2}_{\perp}.$ (3.11)

The power spectrum model in Eq. (3.9) was first proposed empirically by
Eisenstein et al. (2007a). Subsequently, the damping factor ${\cal D}^{2}$ in
front of $P_{\rm lin}$ was derived in the context of perturbation theory by
Crocce & Scoccimarro (2008); Matsubara (2008); an additional term to recover a
smooth linear power spectrum without BAOs, $(1-{\cal D}^{2})P_{\rm nw}$, was
derived using the IR resummation method (Sugiyama & Spergel, 2014; Baldauf et
al., 2015a; Blas et al., 2016a; Ivanov & Sibiryakov, 2018; Sugiyama et al.,
2021).

For the bispectrum, we adopt the following IR-resummed model (Sugiyama et al.,
2021):

$\displaystyle B(\textbf{{k}}_{1},\textbf{{k}}_{2})$ $\displaystyle=$
$\displaystyle
2\,Z_{2}(\textbf{{k}}_{1},\textbf{{k}}_{2})Z_{1}(\textbf{{k}}_{1})Z_{1}(\textbf{{k}}_{2})$
(3.12) $\displaystyle\times$ $\displaystyle\Big{\\{}{\cal
D}(\textbf{{k}}_{1}){\cal D}(\textbf{{k}}_{2}){\cal D}(\textbf{{k}}_{3})P_{\rm
w}(k_{1})P_{\rm w}(k_{2})$ $\displaystyle+$ $\displaystyle{\cal
D}^{2}(\textbf{{k}}_{1})P_{\rm w}(k_{1})P_{\rm nw}(k_{2})+{\cal
D}^{2}(\textbf{{k}}_{2})P_{\rm nw}(k_{1})P_{\rm w}(k_{2})$ $\displaystyle+$
$\displaystyle P_{\rm nw}(k_{1})P_{\rm nw}(k_{2})\Big{\\}}+\mbox{2 perms.},$

where $\textbf{{k}}_{1}+\textbf{{k}}_{2}+\textbf{{k}}_{3}=0$. As in the case
of the power spectrum, this bispectrum model restores the tree-level solution
(2.22) consisting of a smooth version (without BAOs) of the linear power
spectrum after degrading the BAO signature 111111Blas et al. (2016a); Ivanov &
Sibiryakov (2018) proposed a bispectrum model similar to Eq. (3.12). However,
the authors ignore the ${\cal O}(P^{2}_{\rm w}/P^{2}_{\rm nw})$ term, so their
model does not include the second line term,
$D(\textbf{{k}}_{1})D(\textbf{{k}}_{2})D(\textbf{{k}}_{3})P_{\rm
w}(k_{1})P_{\rm w}(k_{2})$, in Eq. (3.12). This term added by Sugiyama et al.
(2021) contains the full tree-level solution..

###  3.3 Parameterization method for the bispectrum

The non-linear kernel functions $F^{(\rm m)}_{2}$ and $G^{(\rm m)}_{2}$ can be
decomposed into three terms using Legendre polynomial functions ${\cal
L}_{\ell}(\hat{p}_{1}\cdot\hat{p}_{2})$: i.e., monopole, dipole, and
quadrupole components (Schmittfull et al., 2015). They are called the growth,
shift, and tidal terms, and are understood in $\Lambda$CDM as follows: the
growth term represents the spherical collapse of density fluctuations (Fosalba
& Gaztanaga, 1998); the shift term appears in the form
$\Psi_{1}^{i}\partial_{i}\delta_{1}$ or $\Psi_{1}^{i}\partial_{i}\theta_{1}$
as a coordinate transformation of $\delta$ or $\theta$ by the displacement
vector $\mathbf{\Psi}$; the last term represents the tidal force (2.16). Then,
$F^{(\rm m)}_{2}$ and $G^{(\rm m)}_{2}$ (2.10) are rewritten as (e.g., Bouchet
et al., 1992; Sherwin & Zaldarriaga, 2012; Baldauf et al., 2012; Schmittfull
et al., 2015)

$\displaystyle F^{(\rm m)}_{2}$ $\displaystyle=$
$\displaystyle\left(\kappa_{\delta}-\frac{4}{21}\lambda_{\delta}\right)+\kappa_{\delta}S(\textbf{{k}}_{1},\textbf{{k}}_{2})+\frac{2}{7}\lambda_{\delta}T(\textbf{{k}}_{1},\textbf{{k}}_{2}),$
$\displaystyle G^{(\rm m)}_{2}$ $\displaystyle=$
$\displaystyle\left(\kappa_{\theta}-\frac{8}{21}\lambda_{\theta}\right)+\kappa_{\theta}S(\textbf{{k}}_{1},\textbf{{k}}_{2})+\frac{4}{7}\lambda_{\theta}T(\textbf{{k}}_{1},\textbf{{k}}_{2}),$
(3.13)

where $S$ and $T$ are the scale-dependent functions characterising the shift
and tidal terms:

$\displaystyle S(\textbf{{k}}_{1},\textbf{{k}}_{2})$ $\displaystyle=$
$\displaystyle\frac{1}{2}(\hat{k}_{1}\cdot\hat{k}_{2})\left(\frac{k_{1}}{k_{2}}+\frac{k_{2}}{k_{1}}\right),$
$\displaystyle T(\textbf{{k}}_{1},\textbf{{k}}_{2})$ $\displaystyle=$
$\displaystyle(\hat{k}_{1}\cdot\hat{k}_{2})^{2}-\frac{1}{3}.$ (3.14)

As mentioned in Section 2.1, the coefficients of the growth, shift, and tidal
terms are not independent of each other but are related to under the condition
that $F_{2}^{(\rm m)}(\textbf{{p}},-\textbf{{p}})=G_{2}^{(\rm
m)}(\textbf{{p}},-\textbf{{p}})=0$. Therefore, the coefficient of the growth
term is determined from the coefficients of the shift and tidal terms.

Considering the linear and non-linear bias effects, that the second-order
fluctuations are proportional to $\sigma_{8}^{2}$, and that $G_{2}^{(\rm m)}$
always appears with $f$, we introduce the following parameterisation,

$\displaystyle F^{\rm(g)}_{2}\sigma_{8}^{2}$ $\displaystyle=$
$\displaystyle(b_{1}\sigma_{8})\left[(F_{\rm g}\sigma_{8})+(F_{\rm
s}\sigma_{8})S+(F_{\rm t}\sigma_{8})T\right],$ $\displaystyle
fG^{\rm(g)}_{2}\sigma_{8}^{2}$ $\displaystyle=$
$\displaystyle(f\sigma_{8})\left[(G_{\rm g}\sigma_{8})+(G_{\rm
s}\sigma_{8})S+(G_{\rm t}\sigma_{8})T\right].$ (3.15)

DHOST theories have $G_{\rm g}=G_{\rm s}-(2/3)G_{\rm t}$ from the condition
$G_{2}^{(\rm g)}(\textbf{{p}},-\textbf{{p}})=0$; Horndeski theories further
have $F_{\rm s}=G_{\rm s}=1$. The specific form of each coefficient in DHOST
theories is given by

$\displaystyle F_{\rm g}$ $\displaystyle=$
$\displaystyle\kappa_{\delta}-\frac{4}{21}\lambda_{\delta}+\frac{1}{2}\frac{b_{2}}{b_{1}},$
$\displaystyle F_{\rm s}$ $\displaystyle=$ $\displaystyle\kappa_{\delta},$
$\displaystyle F_{\rm t}$ $\displaystyle=$
$\displaystyle\frac{2}{7}\lambda_{\delta}+\frac{b_{s^{2}}}{b_{1}},$
$\displaystyle G_{\rm g}$ $\displaystyle=$
$\displaystyle\kappa_{\theta}-\frac{8}{21}\lambda_{\theta},$ $\displaystyle
G_{\rm s}$ $\displaystyle=$ $\displaystyle\kappa_{\theta},$ $\displaystyle
G_{\rm t}$ $\displaystyle=$ $\displaystyle\frac{4}{7}\lambda_{\theta}.$ (3.16)

In Eq. (3.16), $F_{\rm g}$ and $F_{\rm t}$ do not contain any cosmological
information because they are degenerate with the non-linear bias parameters,
and $G_{\rm g}$ is determined from $G_{\rm t}$ and $G_{\rm s}$. Thus,
cosmologically meaningful parameters are $F_{\rm s}$, $G_{\rm s}$, and $G_{\rm
t}$.

Following the method proposed by Sugiyama et al. (2021), we decompose the IR-
resummed bispectrum model into

$\displaystyle
B(\textbf{{k}}_{1},\textbf{{k}}_{2})=\sum_{p=1}^{22}X^{(p)}B^{(p)}(\textbf{{k}}_{1},\textbf{{k}}_{2}),$
(3.17)

with

$\displaystyle B^{(p)}(\textbf{{k}}_{1},\textbf{{k}}_{2})$ $\displaystyle=$
$\displaystyle 2\,H^{(p)}(\textbf{{k}}_{1},\textbf{{k}}_{2})$ (3.18)
$\displaystyle\times$ $\displaystyle\Big{\\{}{\cal D}(\textbf{{k}}_{1}){\cal
D}(\textbf{{k}}_{2}){\cal D}(\textbf{{k}}_{3})P^{\rm(n)}_{\rm
w}(k_{1})P^{\rm(n)}_{\rm w}(k_{2})$ $\displaystyle+$ $\displaystyle{\cal
D}^{2}(\textbf{{k}}_{1})P^{\rm(n)}_{\rm w}(k_{1})P^{\rm(n)}_{\rm nw}(k_{2})$
$\displaystyle+$ $\displaystyle{\cal D}^{2}(\textbf{{k}}_{2})P^{\rm(n)}_{\rm
nw}(k_{1})P^{\rm(n)}_{\rm w}(k_{2})$ $\displaystyle+$ $\displaystyle
P^{\rm(n)}_{\rm nw}(k_{1})P^{\rm(n)}_{\rm nw}(k_{2})\Big{\\}}+\mbox{2
perms.},$

where $P_{\rm w}^{(\rm n)}$ and $P_{\rm nw}^{(\rm n)}$ are respectively the
wiggle and no-wiggle linear matter power spectra normalized by
$\sigma^{2}_{8}$: $P_{\rm w}^{(\rm n)}=P_{\rm w}/\sigma_{8}^{2}$ and $P_{\rm
nw}^{(\rm n)}=P_{\rm nw}/\sigma_{8}^{2}$. The functions $X^{(p)}$ ($p=1-22$)
represent the combinations of the parameters of interest and are given by

$\displaystyle X^{(1)}$ $\displaystyle=$ $\displaystyle(F_{\rm
g}\sigma_{8})(b_{1}\sigma_{8})^{3},$ $\displaystyle X^{(2)}$ $\displaystyle=$
$\displaystyle(F_{\rm s}\sigma_{8})(b_{1}\sigma_{8})^{3},$ $\displaystyle
X^{(3)}$ $\displaystyle=$ $\displaystyle(F_{\rm
t}\sigma_{8})(b_{1}\sigma_{8})^{3},$ $\displaystyle X^{(4)}$ $\displaystyle=$
$\displaystyle(F_{\rm g}\sigma_{8})(b_{1}\sigma_{8})^{2}(f\sigma_{8}),$
$\displaystyle X^{(5)}$ $\displaystyle=$ $\displaystyle(F_{\rm
s}\sigma_{8})(b_{1}\sigma_{8})^{2}(f\sigma_{8}),$ $\displaystyle X^{(6)}$
$\displaystyle=$ $\displaystyle(F_{\rm
t}\sigma_{8})(b_{1}\sigma_{8})^{2}(f\sigma_{8}),$ $\displaystyle X^{(7)}$
$\displaystyle=$ $\displaystyle(F_{\rm
g}\sigma_{8})(b_{1}\sigma_{8})(f\sigma_{8})^{2},$ $\displaystyle X^{(8)}$
$\displaystyle=$ $\displaystyle(F_{\rm
s}\sigma_{8})(b_{1}\sigma_{8})(f\sigma_{8})^{2},$ $\displaystyle X^{(9)}$
$\displaystyle=$ $\displaystyle(F_{\rm
t}\sigma_{8})(b_{1}\sigma_{8})(f\sigma_{8})^{2},$ $\displaystyle X^{(10)}$
$\displaystyle=$ $\displaystyle(G_{\rm
g}\sigma_{8})(b_{1}\sigma_{8})^{2}(f\sigma_{8}),$ $\displaystyle X^{(11)}$
$\displaystyle=$ $\displaystyle(G_{\rm
s}\sigma_{8})(b_{1}\sigma_{8})^{2}(f\sigma_{8}),$ $\displaystyle X^{(12)}$
$\displaystyle=$ $\displaystyle(G_{\rm
t}\sigma_{8})(b_{1}\sigma_{8})^{2}(f\sigma_{8}),$ $\displaystyle X^{(13)}$
$\displaystyle=$ $\displaystyle(G_{\rm
g}\sigma_{8})(b_{1}\sigma_{8})(f\sigma_{8})^{2},$ $\displaystyle X^{(14)}$
$\displaystyle=$ $\displaystyle(G_{\rm
s}\sigma_{8})(b_{1}\sigma_{8})(f\sigma_{8})^{2},$ $\displaystyle X^{(15)}$
$\displaystyle=$ $\displaystyle(G_{\rm
t}\sigma_{8})(b_{1}\sigma_{8})(f\sigma_{8})^{2},$ $\displaystyle X^{(16)}$
$\displaystyle=$ $\displaystyle(G_{\rm g}\sigma_{8})(f\sigma_{8})^{3},$
$\displaystyle X^{(17)}$ $\displaystyle=$ $\displaystyle(G_{\rm
s}\sigma_{8})(f\sigma_{8})^{3},$ $\displaystyle X^{(18)}$ $\displaystyle=$
$\displaystyle(G_{\rm t}\sigma_{8})(f\sigma_{8})^{3},$ $\displaystyle
X^{(19)}$ $\displaystyle=$ $\displaystyle(b_{1}\sigma_{8})^{3}(f\sigma_{8}),$
$\displaystyle X^{(20)}$ $\displaystyle=$
$\displaystyle(b_{1}\sigma_{8})^{2}(f\sigma_{8})^{2},$ $\displaystyle
X^{(21)}$ $\displaystyle=$ $\displaystyle(b_{1}\sigma_{8})(f\sigma_{8})^{3},$
$\displaystyle X^{(22)}$ $\displaystyle=$ $\displaystyle(f\sigma_{8})^{4}.$
(3.19)

The scale-dependent functions $H^{(p)}$ ($p=1-22$) are derived by decomposing
the non-linear kernel functions $Z_{1}Z_{1}Z_{2}$ in terms of the parameters,
given by

$\displaystyle H^{(1)}$ $\displaystyle=$ $\displaystyle 1,$ $\displaystyle
H^{(2)}$ $\displaystyle=$ $\displaystyle
S(\textbf{{k}}_{1},\textbf{{k}}_{2}),$ $\displaystyle H^{(3)}$
$\displaystyle=$ $\displaystyle T(\textbf{{k}}_{1},\textbf{{k}}_{2}),$
$\displaystyle H^{(4)}$ $\displaystyle=$
$\displaystyle(\mu_{1}^{2}+\mu_{2}^{2}),$ $\displaystyle H^{(5)}$
$\displaystyle=$ $\displaystyle
S(\textbf{{k}}_{1},\textbf{{k}}_{2})(\mu_{1}^{2}+\mu_{2}^{2}),$ $\displaystyle
H^{(6)}$ $\displaystyle=$ $\displaystyle
T(\textbf{{k}}_{1},\textbf{{k}}_{2})(\mu_{1}^{2}+\mu_{2}^{2}),$ $\displaystyle
H^{(7)}$ $\displaystyle=$ $\displaystyle(\mu_{1}^{2}\mu_{2}^{2}),$
$\displaystyle H^{(8)}$ $\displaystyle=$ $\displaystyle
S(\textbf{{k}}_{1},\textbf{{k}}_{2})(\mu_{1}^{2}\mu_{2}^{2}),$ $\displaystyle
H^{(9)}$ $\displaystyle=$ $\displaystyle
T(\textbf{{k}}_{1},\textbf{{k}}_{2})(\mu_{1}^{2}\mu_{2}^{2}),$ $\displaystyle
H^{(10)}$ $\displaystyle=$ $\displaystyle(\mu^{2}),$ $\displaystyle H^{(11)}$
$\displaystyle=$ $\displaystyle
S(\textbf{{k}}_{1},\textbf{{k}}_{2})(\mu^{2}),$ $\displaystyle H^{(12)}$
$\displaystyle=$ $\displaystyle
T(\textbf{{k}}_{1},\textbf{{k}}_{2})(\mu^{2}),$ $\displaystyle H^{(13)}$
$\displaystyle=$ $\displaystyle(\mu^{2})(\mu_{1}^{2}+\mu_{2}^{2}),$
$\displaystyle H^{(14)}$ $\displaystyle=$ $\displaystyle
S(\textbf{{k}}_{1},\textbf{{k}}_{2})(\mu^{2})(\mu_{1}^{2}+\mu_{2}^{2}),$
$\displaystyle H^{(15)}$ $\displaystyle=$ $\displaystyle
T(\textbf{{k}}_{1},\textbf{{k}}_{2})(\mu^{2})(\mu_{1}^{2}+\mu_{2}^{2}),$
$\displaystyle H^{(16)}$ $\displaystyle=$
$\displaystyle(\mu^{2})(\mu_{1}^{2}\mu_{2}^{2}),$ $\displaystyle H^{(17)}$
$\displaystyle=$ $\displaystyle
S(\textbf{{k}}_{1},\textbf{{k}}_{2})(\mu^{2})(\mu_{1}^{2}\mu_{2}^{2}),$
$\displaystyle H^{(18)}$ $\displaystyle=$ $\displaystyle
T(\textbf{{k}}_{1},\textbf{{k}}_{2})(\mu^{2})(\mu_{1}^{2}\mu_{2}^{2}),$
$\displaystyle H^{(19)}$ $\displaystyle=$ $\displaystyle
DV(\textbf{{k}}_{1},\textbf{{k}}_{2}),$ $\displaystyle H^{(20)}$
$\displaystyle=$ $\displaystyle
DV(\textbf{{k}}_{1},\textbf{{k}}_{2})(\mu_{1}^{2}+\mu_{2}^{2})+V(\textbf{{k}}_{1},\textbf{{k}}_{2}),$
$\displaystyle H^{(21)}$ $\displaystyle=$ $\displaystyle
DV(\textbf{{k}}_{1},\textbf{{k}}_{2})(\mu_{1}^{2}\mu_{2}^{2})+V(\textbf{{k}}_{1},\textbf{{k}}_{2})(\mu_{1}^{2}+\mu_{2}^{2}),$
$\displaystyle H^{(22)}$ $\displaystyle=$ $\displaystyle
V(\textbf{{k}}_{1},\textbf{{k}}_{2})(\mu_{1}^{2}\mu_{2}^{2}),$ (3.20)

where $\textbf{{k}}=\textbf{{k}}_{1}+\textbf{{k}}_{2}$,
$\mu=\hat{k}\cdot\hat{n}$, $\mu_{1}=\hat{k}_{1}\cdot\hat{n}$,
$\mu_{2}=\hat{k}_{2}\cdot\hat{n}$, and

$\displaystyle V(\textbf{{k}}_{1},\textbf{{k}}_{2})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\frac{k^{2}}{k_{1}k_{2}}\mu^{2}\mu_{1}\mu_{2},$
$\displaystyle DV(\textbf{{k}}_{1},\textbf{{k}}_{2})$ $\displaystyle=$
$\displaystyle\frac{1}{2}k\mu\left[\frac{\mu_{1}}{k_{1}}+\frac{\mu_{2}}{k_{2}}\right].$
(3.21)

We pre-compute $B^{(p)}(\textbf{{k}}_{1},\textbf{{k}}_{2})$ using the fiducial
cosmology introduced in Section 1 and save the resulting data in a file. In
this way, when constraining $X^{(p)}$ from the BOSS data, we can quickly
calculate the bispectrum by loading the data file containing $B^{(p)}$ and
substituting them into Eq. (3.17) along with $X^{(p)}$.

Here we demonstrate how the growth, shift, and tidal terms of the second-order
density and velocity fields affect the multipole components of the 3PCF. To do
so, we consider the following seven bispectra:

$\displaystyle B_{\rm FG}(\textbf{{k}}_{1},\textbf{{k}}_{2})$ $\displaystyle=$
$\displaystyle\sum_{p=1,4,7}X^{(p)}B^{(p)}(\textbf{{k}}_{1},\textbf{{k}}_{2}),$
$\displaystyle B_{\rm FS}(\textbf{{k}}_{1},\textbf{{k}}_{2})$ $\displaystyle=$
$\displaystyle\sum_{p=2,5,8}X^{(p)}B^{(p)}(\textbf{{k}}_{1},\textbf{{k}}_{2}),$
$\displaystyle B_{\rm FT}(\textbf{{k}}_{1},\textbf{{k}}_{2})$ $\displaystyle=$
$\displaystyle\sum_{p=3,6,9}X^{(p)}B^{(p)}(\textbf{{k}}_{1},\textbf{{k}}_{2}),$
$\displaystyle B_{\rm GG}(\textbf{{k}}_{1},\textbf{{k}}_{2})$ $\displaystyle=$
$\displaystyle\sum_{p=10,13,16}X^{(p)}B^{(p)}(\textbf{{k}}_{1},\textbf{{k}}_{2}),$
$\displaystyle B_{\rm GS}(\textbf{{k}}_{1},\textbf{{k}}_{2})$ $\displaystyle=$
$\displaystyle\sum_{p=11,14,17}X^{(p)}B^{(p)}(\textbf{{k}}_{1},\textbf{{k}}_{2}),$
$\displaystyle B_{\rm GT}(\textbf{{k}}_{1},\textbf{{k}}_{2})$ $\displaystyle=$
$\displaystyle\sum_{p=12,15,18}X^{(p)}B^{(p)}(\textbf{{k}}_{1},\textbf{{k}}_{2}),$
$\displaystyle B_{\rm BF}(\textbf{{k}}_{1},\textbf{{k}}_{2})$ $\displaystyle=$
$\displaystyle\sum_{p=19,20,21,22}X^{(p)}B^{(p)}(\textbf{{k}}_{1},\textbf{{k}}_{2}).$
(3.22)

where $B_{\rm FG}$, $B_{\rm FS}$, $B_{\rm FT}$, $B_{\rm GG}$, $B_{\rm GS}$,
and $B_{\rm GT}$ are proportional to $(F_{\rm g}\sigma_{8})$, $(F_{\rm
s}\sigma_{8})$, $(F_{\rm t}\sigma_{8})$, $(G_{\rm g}\sigma_{8})$, $(G_{\rm
s}\sigma_{8})$, and $(G_{\rm t}\sigma_{8})$, respectively, and $B_{\rm BF}$
depends only on $(b_{1}\sigma_{8})$ and $(f\sigma_{8})$. When computing the
above seven bispectra, we assume the cosmological parameters in $\Lambda$CDM
given in Section 1, the linear bias parameter $b_{1}=2$, no non-linear bias,
i.e. $b_{2}=b_{s^{2}}=0$, and the redshift $z=0.61$. Next, we decompose the
seven bispectra using TripoSHs according to Section 3.1 and compute the 3PCF
multipoles via the 2D Hankel transform (3.6). We plot the resulting 3PCF
multipoles in Figures 1 and 2 as a function of $r_{2}$ after fixing $r_{1}$ to
$50$, $80$, $90$, $100$, and $130\,h^{-1}\,{\rm Mpc}$.

As shown in Sugiyama et al. (2021), in the monopole component (Figure 1), the
growth term (“FG”) is positive for scales smaller than $\sim 130\,h^{-1}\,{\rm
Mpc}$ and has a peak at $r_{1}=r_{2}$, while it goes from positive to negative
and behaves like a trough for scales above $\sim 130\,h^{-1}\,{\rm Mpc}$. On
the other hand, the shift (“FS”) and tidal (“FT”) terms have troughs for any
scale. Depending on the scale of interest, the shift term dominates for scales
above $\sim 30\,h^{-1}\,{\rm Mpc}$, and the total 3PCF (“total”), which is the
sum of all components, is found to have a trough. To illustrate the trough-
like behavior of the 3PCF at $r_{1}=r_{2}$, we have drawn vertical black lines
representing $r_{1}=r_{2}$ in Figures 1 and 2. It can be seen that the bottom
of the trough of the black curve representing the total 3PCF is always on the
line $r_{1}=r_{2}$. Around $r_{1}\sim 100\,h^{-1}\,{\rm Mpc}$, the BAO peak
appears and has a wavy shape as it cancels out the trough due to non-linear
gravity effects (e.g., see the middle panels). At $r_{1}=130\,h^{-1}\,{\rm
Mpc}$ (the bottom panels), almost all the components have troughs, so the 3PCF
has a more significant trough at $r_{1}=r_{2}$.

The quadrupole component (Figure 2) of the 3PCF only shows an overall trough
behaviour because the BAO signal is sufficiently non-linearly damped. The most
dominant term in the quadrupole 3PCF is the “BF” term, which does not depend
on any non-linear coefficients such as $F_{\rm g}$ or $G_{\rm g}$. This “BF”
term consists of two effects: first, a term expressed as the product of a
linear density field and a linear velocity field, and second, a term expressed
as the square of the linear velocity field. In particular, the former can be
interpreted as a new shift term resulting from the coordinate transformation
from real to redshift space (2.18), and it dominates the “BF” term. Therefore,
it behaves similarly to the shift term in the monopole 3PCF and explains most
of the trough structure in the quadrupole 3PCF. The growth (“GG”), shift
(“GS”), and tidal (“GT”) terms in the non-linear velocity field contribute to
the quadrupole 3PCF comparably to those in the non-linear density field, and
thus we can use the quadrupole 3PCF to determine the “GG”, “GS”, and “GT”
terms. In contrast to the monopole case, the growth terms (“FG” and “GG”) are
negative and behave as troughs, while the shift terms (“FS” and “GS”) are
positive.

Figure 1: The monopole 3PCFs, $\zeta_{000}$ (left) and $\zeta_{110}$ (right),
calculated from the decomposed bispectra (3.22) according to the parameter
dependence, are shown as a function of $r_{2}$ after fixing $r_{1}$ to $50$,
$80$, $90$, $100$, and $130\,h^{-1}\,{\rm Mpc}$. The “FG”, “FS”, and “FT”
terms arise from the growth, shift and tidal effects of the non-linear density
fluctuation; the “GG”, “GS”, and “GT” terms arise from those of the non-linear
velocity field; the “BF” term consists only of linear density and linear
velocity fields; the “total” term is the sum of all the decomposed components.
For these calculations, the $\Lambda$CDM model at $z=0.61$, the linear bias
$b_{1}=2.0$, and no non-linear bias are assumed.

Figure 2: Same as Figure 1, except that the quadrupole 3PCFs, $\zeta_{202}$
and $\zeta_{112}$, are shown.

###  3.4 Time-dependences of parameters

We review the discussion by Yamauchi & Sugiyama (2022) on introducing new
parameters to test DHOST theories and their time-dependences.

Note that some previous works predict that constraining $\sigma_{8}$ alone
from the 3PCF can break the degeneracy between $f\sigma_{8}$ and $\sigma_{8}$,
but this no longer happens in the framework of DHOST theories. To illustrate
this fact in the context of our parameterisation, we can see from Eq. (3.15)
that the coefficient of the shift term in the second-order density fluctuation
in $\Lambda$CDM ($F_{\rm s}=1$) determines $\sigma_{8}$ because both the
growth and tidal terms are degenerate with the non-linear bias parameters
(Schmittfull et al., 2015). However, in the case of DHOST theories, there
appears the parameter $\kappa_{\delta}$ in the coefficient of the shift term,
which makes it impossible to measure $\sigma_{8}$ alone. Therefore, we
introduce three new parameters that are not degenerate with $\sigma_{8}$
following Yamauchi & Sugiyama (2022):

$\displaystyle E_{f}$ $\displaystyle=$
$\displaystyle\frac{f}{\kappa_{\delta}}=\frac{f\sigma_{8}}{F_{\rm
s}\sigma_{8}},$ $\displaystyle E_{\rm s}$ $\displaystyle=$
$\displaystyle\frac{\kappa_{\theta}}{\kappa_{\delta}}=\frac{G_{\rm
s}\sigma_{8}}{F_{\rm s}\sigma_{8}},$ $\displaystyle E_{\rm t}$
$\displaystyle=$
$\displaystyle\frac{\lambda_{\theta}}{\kappa_{\delta}}=\frac{7}{4}\frac{G_{\rm
t}\sigma_{8}}{F_{\rm s}\sigma_{8}}.$ (3.23)

In GR or Horndeski theories, $E_{f}=f$, $E_{\rm s}=1$ and $E_{\rm
t}=\lambda_{\theta}$, because $\kappa_{\delta}=\kappa_{\theta}=1$. Horndeski
theories differ from $\Lambda$CDM only in $f$ and $E_{\rm t}$ while keeping
$E_{\rm s}=1$. If $E_{\rm s}\neq 1$, then the signal is specific to DHOST
theories; $E_{\rm s}\neq 1$ is a sufficient condition for detecting DHOST
theories because there can be DHOST theories satisfying $E_{\rm s}=1$.

It has been known for a long time that the coefficient of the tidal term in
the non-linear density field, $\lambda_{\delta}$, is time-dependent in GR
(e.g., Bouchet et al., 1992), and in the case of $\Lambda$CDM, the following
approximation holds well with an precision better than $0.6\%$ (Bouchet et
al., 1995; Yamauchi et al., 2017b)121212The original derivation of the
equation was calculated in the Lagrangian picture and is given in the form
(Bouchet et al., 1995)
$\displaystyle\frac{2}{7}\lambda_{\delta}=\frac{1}{2}\left[1-\frac{3}{7}\Omega_{\rm
m}^{-1/143}\right].$ (3.24) This equation can be rewritten to Eq. (3.25) under
the condition $(1-\Omega_{\rm m})\ll 1$. :

$\displaystyle\lambda_{\delta}\sim\Omega_{\rm m}^{3/572}.$ (3.25)

Through Eq. (2.12), the coefficient of the tidal term in the non-linear
velocity field, $\lambda_{\theta}$, is also given by (Yamauchi & Sugiyama,
2022)

$\displaystyle\lambda_{\theta}\sim\Omega_{\rm m}^{15/1144}.$ (3.26)

Yamauchi et al. (2017b) extended the above discussion to Horndeski theories
and showed that $\lambda_{\delta}$ is parameterised as a power of $\Omega_{\rm
m}$ in Horndeski theories. In addition, Yamauchi & Sugiyama (2022) performed
similar calculations for DHOST theories and showed that the coefficient of the
shift term, $\kappa_{\delta}$, is also described by a power of $\Omega_{\rm
m}$. The coefficients of the shift and tidal terms in the non-linear velocity
field can be calculated through Eq. (2.12), and they also follow the powers of
$\Omega_{\rm m}$. Therefore, we can parameterise the time-dependences of
$E_{\rm f}$, $E_{\rm s}$, and $E_{\rm t}$ as follows:

$\displaystyle E_{f}\sim\Omega_{\rm m}^{\xi_{f}},\quad E_{\rm
s}\sim\Omega_{\rm m}^{\xi_{\rm s}},\quad E_{\rm t}\sim\Omega_{\rm m}^{\xi_{\rm
t}}.$ (3.27)

In GR, we have

$\displaystyle\xi_{f}=\frac{6}{11},\quad\xi_{\rm s}=0,\quad\xi_{\rm
t}=\frac{15}{1144}.$ (3.28)

In summary, we parameterise the second-order kernel function of the velocity
field (3.15) as

$\displaystyle fG_{2}^{\rm(g)}\sigma_{8}^{2}=\Omega_{\rm m}^{\xi_{f}}\,(F_{\rm
s}\sigma_{8})^{2}\Bigg{[}\left(G_{\rm g}\right)+\Omega_{\rm m}^{\xi_{\rm
s}}\,S+\frac{4}{7}\Omega_{\rm m}^{\xi_{\rm t}}\,T\Bigg{]},$ (3.29)

where $(G_{\rm g})=\Omega_{\rm m}^{\xi_{\rm s}}-(8/21)\Omega_{\rm m}^{\xi_{\rm
t}}$, and the functions $S$ and $T$ are given in Eq. (3.14). We will test the
theory of gravity by measuring the above three parameters, $\xi_{f}$,
$\xi_{\rm s}$, and $\xi_{\rm t}$, from the BOSS data in Section 9.

In DHOST theories, the Planck mass is time-varying, and the time variation of
the Hubble parameter is different from GR. Therefore, one may be concerned
that the time dependence of $\Omega_{\rm m}$ is different from $\Omega_{\rm
m}^{\rm GR}$ that is calculated assuming GR. However, Appendix $C$ in Yamauchi
& Sugiyama (2022) showed that the difference between DHOST theories and GR is
suppressed by $(1-\Omega_{\rm m}^{\rm GR})$. Hence, we can replace
$\Omega_{\rm m}$ in Eq. (3.27) with $\Omega_{\rm m}^{\rm GR}$ as an
approximation and perform the analysis to constrain $\xi_{f}$, $\xi_{\rm s}$,
and $\xi_{\rm t}$.

###  3.5 Limitations of our theoretical approach to the 2PCF and 3PCF

In this subsection, we discuss the validity of the calculation methods of the
2PCF and 3PCF models described so far and the limitations of their
application.

First, we can use the TripoSH decomposed 3PCF (3.6) to constrain all the scale
dependencies in the 3PCF, such as the shift and tidal terms, as shown Figures
1 and 2, because it does not focus only on specific scale dependencies such as
the squeezed limit. However, our analysis that uses only some multipoles of
the TripoSH decomposition does not fully utilize the information on the scale
dependence of the 3PCF. The reason for restricting the multipole components
used in this work is to keep the number of data bins much smaller than the
number of mock simulations used to compute the covariance matrix (Section 5).
Therefore, increasing the number of multipoles in the 3PCF to be considered
will improve the results of this work when more mock catalogues are created in
the future.

Second, note that the power spectrum and bispectrum models in Eqs. (3.9) and
(3.12) are valid for any theory in which the IR cancellation occurs based on
the Galilean invariance of the system of equations in the IR limit: i.e.,
these models hold not only for $\Lambda$CDM but also for Horndeski theories
(Crisostomi et al., 2020). On the other hand, as Lewandowski (2020) pointed
out in the power spectrum case, additional terms arise when performing the IR
resummation in DHOST theories because of the violation of the IR cancellation.
Specifically, when one applies the IR limit to the one-loop solution of the
power spectrum in DHOST theories, a term proportional to $k^{2}P_{\rm lin}(k)$
appears, changing the shape of the power spectrum (Crisostomi et al., 2020;
Lewandowski, 2020; Hirano et al., 2020). Note also that this additional term
is degenerate with the higher-order derivative bias, which is ignored in this
paper. However, since this additional term in the IR limit is proportional to
$k^{2}P_{\rm lin}(k)$, it is considered to be negligible at the large scales
of interest in this paper ($\geq 80\,h^{-1}\,{\rm Mpc}$). Assuming that the
same should happen in the bispectrum, we directly use the power and bispectrum
models in Eqs. (3.9) and (3.12) in the present analysis. In addition, it
should be noted that Hirano et al. (2020) showed that in DHOST theories, a
term consisting of the product of first- and third-order fluctuations in the
one-loop power spectrum causes UV divergence. Further model development is
thus needed to take advantage of smaller-scale information by solving these
problems.

Third, since the linear equation for density fluctuations is scale-independent
(2.6), we assume that we can use the shape of the linear matter power spectrum
determined in the high-$z$ region, where the scalar field is expected to be
sub-dominant. Thus, we can pre-compute the $\sigma_{8}^{2}$-normalized wiggle
and no-wiggle power spectra, $P_{\rm w}^{(\rm n)}$ and $P_{\rm nw}^{(\rm n)}$,
appearing in the $B^{(p)}$ terms (3.18), using a $\Lambda$CDM model.

Fourth, there is a concern about the pre-computation of ${\cal
D}(\textbf{{k}})$ (3.10) appearing in the $B^{(p)}$ terms (3.18). It is known
that $\sigma_{\perp}$ and $\sigma_{\parallel}$, which characterise ${\cal
D}(\textbf{{k}})$, can be calculated successfully using linear displacement
vectors (e.g., Matsubara, 2008), and we adopt the same calculation in this
paper (3.11). Since $\sigma_{\perp}$ and $\sigma_{\parallel}$ in the linear
theory depend on $f$ and $\sigma_{8}$, their values should differ for
different gravity theories. For this reason, it is desirable to vary
$\sigma_{\perp}$ and $\sigma_{\parallel}$ as free parameters in the data
analysis. However, to do so, the bispectrum decomposition method in Eq. (3.17)
cannot be applied, and the computation time of the bispectrum model increases
significantly, making it challenging to perform cosmological analysis.
Fortunately, the BAO signal does not significantly impact the shape of the
3PCF. The reason is that the BAO signal is maximized when $r_{1}\sim r_{2}\sim
100\,h^{-1}\,{\rm Mpc}$, while $r_{1}$ and $r_{2}$ can take various
combinations in the 3PCF (Sugiyama et al., 2021). Therefore, in this paper, we
ignore the concern about ${\cal D}(\textbf{{k}})$ and pre-compute
$\sigma_{\perp}$ and $\sigma_{\parallel}$ using the linear theory in
$\Lambda$CDM. Furthermore, to keep consistency with the 3PCF calculation, we
fix $\sigma_{\perp}$ and $\sigma_{\parallel}$ to those calculated using the
$\Lambda$CDM model in the 2PCF calculation as well.

Finally, to simplify the analysis, we ignore the AP effect (Alcock &
Paczyński, 1979), which can directly measure the Hubble parameter and angular
radial distance at the redshift of the galaxy distribution of interest.
Ignoring the AP effect means that the values of the angular diameter distance
and the Hubble parameter, which should be constrained by the AP effect, are
given by the fiducial $\Lambda$CDM. In this sense, this paper assumes the
fiducial $\Lambda$CDM for the expansion of the Universe in the background
spacetime. However, the AP effect can be determined by the 2PCF by a few
percent and is not expected to significantly affect the constraint results for
the parameters that characterize the nonlinear fluctuations of interest in
this paper, such as $\xi_{\rm t}$ and $\xi_{\rm s}$. Since DHOST theories vary
these parameter values, the AP effect is expected to provide further
information into the constraint on DHOST theories. Sugiyama et al. (2021) have
performed a joint analysis of the anisotropic 2PCF and 3PCF to constrain the
AP parameters under the GR assumption. Combining that method with the analysis
method developed in this paper allows for consistent DHOST theory constraints
that simultaneously account for the AP and nonlinear gravity effects, which is
left as future work.

##  4 Measurements

This section summarises how to measure multipole 2PCFs and 3PCFs from BOSS
galaxy data according to the method proposed by Sugiyama et al. (2021). First,
Section 4.1 introduces the BOSS galaxy data used in this paper and the mock
simulation data designed to reproduce it. Then, Section 4.2 describes the
measurements of the multipole 2PCFs and 3PCFs. Finally, Section 4.3 explains
how to correct for the window function effects on the measured 2PCF and 3PCF.

###  4.1 Data

We use the final galaxy clustering data set, Data Release 12 (DR12; Alam et
al., 2015), from the Baryon Oscillation Spectroscopic Survey (BOSS; Dawson et
al., 2013). The BOSS survey is part of the Sloan Digital Sky Survey III (SDSS
III Eisenstein et al., 2011), selected galaxies from multicolour SDSS imaging
(Fukugita et al., 1996; Gunn et al., 1998; Smith et al., 2002; Gunn et al.,
2006; Doi et al., 2010) and used the SDSS multi-fibre spectrograph (Bolton et
al., 2012; Smee et al., 2013) to measure spectroscopic redshifts of the
galaxies. As detailed in Reid et al. (2016), the BOSS survey has four samples,
CMASS, LOWZ, LOWZ2 and LOWZ3, and those four samples are combined into one
sample. In brief, the survey footprint, veto masks and survey-related
systematics (such as fibre collisions and redshift failures) are considered to
construct data and random catalogues for the DR12 BOSS galaxies. This DR12
combined sample comprises $1.2$ million massive galaxies over an effective
area of $9329\,{\rm deg}^{2}$ and covers a redshift range of $0.2-0.75$. In
our analysis, we split this redshift range into two redshift bins defined by
$0.2<z<0.5$ and $0.5<z<0.75$ with the effective redshifts $z_{\rm eff}=0.38$
and $0.61$, respectively, where the effective redshifts are calculated as the
weighted average over all galaxies (see e.g. Eq. (67) in Beutler et al.,
2014). The DR12 combined sample is observed across the two Galactic
hemispheres, referred to as the Northern and Southern galactic caps (NGC and
SGC, respectively), and the NGC and SGC samples probe slightly different
galaxy populations in the low-redshift part of the combined sample (see
Appendix A in Alam et al., 2015).

To derive the covariance matrices of the 2PCF and 3PCF and test the validity
of the 2PCF and 3PCF models given in Eqs. (3.9) and (3.12), we use the
MultiDark-Patchy mock catalogues (Patchy mocks; Kitaura et al., 2016). The
Patchy mocks have been calibrated to an $N$-body simulation-based reference
sample using approximate galaxy solvers and analytical-statistical biasing
models and incorporate observational effects including the survey geometry,
veto mask and fiber collisions. The reference catalogue is extracted from one
of the BigMultiDark simulations (Klypin et al., 2016), which was performed
using GADGET-2 (Springel, 2005) with $3840^{3}$ particles on a volume of
$(2.5\,h^{-1}\,{\rm Mpc})^{3}$. Halo abundance matching is used to reproduce
the observed BOSS two and three-point clustering measurements (Rodríguez-
Torres et al., 2016). There are $2048$ catalogues available for each the NGC
and SGC over the redshift range $z=0.2-0.75$. The fiducial cosmology for these
mocks assumes a $\Lambda$CDM cosmology with $(\Omega_{\rm\Lambda},\Omega_{\rm
m},\Omega_{\rm b},\sigma_{8},h)$
$=(0.692885,0.307115,0.048206,0.8288,0.6777)$. These fiducial parameters are
slightly different from those used in our analysis of the BOSS galaxy data
introduced in the introduction (Section 1), but we expect that such
differences do not significantly affect the covariance matrix estimations of
the 2PCF and 3PCF.

We include three different incompleteness weights to account for shortcomings
of the BOSS dataset: a fiber collision weight, $w_{\rm cp}$, a redshift
failure weight, $w_{\rm noz}$, and a systematics weight, $w_{\rm sys}$, which
is a combination of a stellar density weight and a seeing condition weight.
Each galaxy observed at position x is counted with the following weight (Ross
et al., 2012; Anderson et al., 2014; Reid et al., 2016):

$\displaystyle w_{\rm c}(\textbf{{x}})=w_{\rm sys}(\textbf{{x}})\left(w_{\rm
cp}(\textbf{{x}})+w_{\rm noz}(\textbf{{x}})-1\right).$ (4.1)

In addition, we use a signal-to-noise weight, the so-called FKP weight,
proposed by Feldman et al. (1994), $w_{\rm
FKP}(\textbf{{x}})=1/[1+\bar{n}_{0}(\textbf{{x}})\bar{P}]$, where
$\bar{P}=10^{4}\,(h^{-1}\,{\rm Mpc})^{3}$. The FKP weight function is
effective not only for the power spectrum but also for the bispectrum when
assuming Gaussian errors (Scoccimarro, 2000), and bispectrum measurements from
the Patchy mock catalogue confirm that the FKP weight improves the bispectrum
signal-to-noise ratio even when including non-Gaussian errors (see Appendix D
in Sugiyama et al., 2019). We expect the validity of the FKP weight to hold
for the 2PCF and 3PCF in configuration space because we measure the 2PCF and
3PCF as Fourier transforms of the power spectrum and bispectrum, respectively
(Section 4.2). For the galaxy data, multiplying the completeness weights by
the FKP weights yields the local weight function that is used in our analysis,
while the random catalogues have only the FKP weights:

$\displaystyle w^{(\rm gal)}(\textbf{{x}})$ $\displaystyle=$ $\displaystyle
w_{\rm c}(\textbf{{x}})\,w_{\rm FKP}(\textbf{{x}}),$ $\displaystyle w^{(\rm
ran)}(\textbf{{x}})$ $\displaystyle=$ $\displaystyle w_{\rm
FKP}(\textbf{{x}}),$ (4.2)

where the superscripts, “(gal)” and “(ran)”, stand for “galaxy” and “random”.

###  4.2 Estimators of the 2PCF and 3PCF

We measure the number densities of both real and random galaxies weighted by
the spherical harmonic function $Y_{\ell m}$:

$\displaystyle D_{\ell m}(\textbf{{x}})$ $\displaystyle=$
$\displaystyle\sum_{i}^{N_{\rm gal}}w^{(\rm gal)}(\textbf{{x}}_{i})Y_{\ell
m}^{*}\left(\hat{x}^{(\rm gal)}_{i}\right)\delta_{\rm
D}\left(\textbf{{x}}-\textbf{{x}}^{(\rm gal)}_{i}\right),$ $\displaystyle
R_{\ell m}(\textbf{{x}})$ $\displaystyle=$ $\displaystyle\sum_{j}^{N_{\rm
ran}}w^{(\rm ran)}(\textbf{{x}}_{j})Y_{\ell m}^{*}\left(\hat{x}^{(\rm
ran)}_{j}\right)\delta_{\rm D}\left(\textbf{{x}}-\textbf{{x}}^{(\rm
ran)}_{j}\right),$ (4.3)

where $N_{\rm gal}$ and $N_{\rm ran}$ are the total number of real and random
galaxies, respectively, and the normal number densities are given by
$D(\textbf{{x}})=\sqrt{4\pi}D_{00}(\textbf{{x}})$ and
$R(\textbf{{x}})=\sqrt{4\pi}R_{00}(\textbf{{x}})$. Defining $N^{\prime}_{\rm
gal}\equiv\int d^{3}xD(\textbf{{x}})$ and $N^{\prime}_{\rm ran}\equiv\int
d^{3}xR(\textbf{{x}})$, we can estimate the survey volume as

$\displaystyle V=\frac{N_{\rm ran}^{{}^{\prime}2}}{\int
d^{3}x\,[R(\textbf{{x}})]^{2}}.$ (4.4)

Then, the observed density fluctuation weighted by $Y_{\ell m}$ is

$\displaystyle\delta_{\rm obs,\ell m}(\textbf{{x}})=V\left[D_{\ell
m}(\textbf{{x}})/N^{\prime}_{\rm gal}-R_{\ell m}(\textbf{{x}})/N^{\prime}_{\rm
ran}\right],$ (4.5)

and

$\displaystyle\delta_{\rm obs}(\textbf{{x}})=\sqrt{4\pi}\delta_{\rm
obs,\,00}(\textbf{{x}}).$ (4.6)

We use the fast Fourier transform (FFT) algorithm 131313http://fftw.org/ to
calculate

$\displaystyle\widetilde{\delta}_{\rm obs,\ell
m}(\textbf{{k}})=\frac{1}{W_{\rm mass}(\textbf{{k}})}\int
d^{3}xe^{-i\textbf{{k}}\cdot\textbf{{x}}}\delta_{\rm obs,\ell
m}(\textbf{{x}}),$ (4.7)

where the Fourier transform of the normal density fluctuation is given by
$\widetilde{\delta}_{\rm
obs}(\textbf{{k}})=\sqrt{4\pi}\,\widetilde{\delta}_{\rm
obs,00}(\textbf{{k}})$, and $W_{\rm mass}(\textbf{{k}})$ is the mass
assignment function that corrects for the effect when arising assign particles
on a regular grid in position space (Jing, 2005). The most popular mass
assignment function is given by (Hockney & Eastwood, 1981)

$\displaystyle W_{\rm mass}(\textbf{{k}})=\prod_{i=x,y,z}\left[{\rm
sinc}\left(\frac{\pi k_{i}}{2k_{\rm N,i}}\right)\right]^{p},$ (4.8)

where $k_{\rm N,i}$ is the Nyquist frequency of $i$-axis with the grid spacing
$H_{i}$ on the axis. The indexes $p=1$, $p=2$, and $p=3$ correspond to the
nearest grid point (NGP), cloud-in-cell (CIC), and triangular-shaped cloud
(TSC) assignment functions, respectively.

The FFT-based estimator of the multipole 2PCFs is given by (Hand et al., 2017;
Sugiyama et al., 2018) (see also Bianchi et al., 2015; Scoccimarro, 2015)

$\displaystyle\widehat{\xi}_{\ell}(r)$ $\displaystyle=$
$\displaystyle\frac{(4\pi)}{V}\sum_{m}\int\frac{d^{2}\hat{r}}{4\pi}Y_{\ell
m}(\hat{r})\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\textbf{{k}}\cdot\textbf{{r}}}$
(4.9) $\displaystyle\times$ $\displaystyle\bigg{[}\widetilde{\delta}_{\rm
obs,\ell m}(\textbf{{k}})\widetilde{\delta}^{*}_{\rm
obs}(\textbf{{k}})-S_{\ell m}(\textbf{{k}})\bigg{]}.$

The shot-noise term $S_{\ell m}(\textbf{{k}})$ is given by

$\displaystyle S_{\ell m}(\textbf{{k}})$ $\displaystyle=$
$\displaystyle\frac{C_{\rm shot}(\textbf{{k}})}{W_{\rm
mass}^{2}(\textbf{{k}})}\left(\frac{V}{N^{\prime}_{\rm gal}}\right)^{2}$
(4.10) $\displaystyle\times$ $\displaystyle\Bigg{[}\sum_{i}^{N_{\rm
gal}}\left[w^{(\rm gal)}(\textbf{{x}}_{i})\right]^{2}Y_{\ell
m}^{*}\left(\hat{x}^{(\rm gal)}_{i}\right)$ $\displaystyle+$
$\displaystyle\left(\frac{N^{\prime}_{\rm gal}}{N^{\prime}_{\rm
ran}}\right)^{2}\sum_{j}^{N_{\rm ran}}\left[w^{(\rm
ran)}(\textbf{{x}}_{j})\right]^{2}Y_{\ell m}^{*}\left(\hat{x}^{(\rm
ran)}_{j}\right)\Bigg{]}.$

where $C_{\rm shot}(\textbf{{k}})$ represents the correction for the
assignment effect to the shot-noise term, given by (Eq. (20) in Jing, 2005)

$\displaystyle C_{\rm shot}(\textbf{{k}})$ (4.11) $\displaystyle=$
$\displaystyle\begin{cases}1,\hskip 152.22241pt\mbox{NGP};\\\
\prod_{i}\left[1-\frac{2}{3}\sin^{2}\left(\frac{\pi k_{i}}{2k_{\rm
N,i}}\right)\right],\hskip 62.59596pt\mbox{CIC};\\\
\prod_{i}\left[1-\sin^{2}\left(\frac{\pi k_{i}}{2k_{\rm
N,i}}\right)+\frac{2}{15}\sin^{4}\left(\frac{\pi k_{i}}{2k_{\rm
N,i}}\right)\right],\mbox{TSC}.\end{cases}$

The angle integral $\int d^{2}\hat{r}/(4\pi)$ in Eq. (4.9) can be rewritten as

$\displaystyle\int\frac{d^{2}\hat{r}}{4\pi}=\frac{1}{N_{r}(r)}\sum_{r-\Delta
r/2<r<r+\Delta r/2},$ (4.12)

where $\Delta r$ is the width of the $r$-bins, and $N_{r}(r)$ is the number of
three-dimensional data contianed in each $r$-bin width. From the expression of
the shot noise term in the 2PCF given in Eq. (4.10), we compute the weighted
mean number density as

$\displaystyle\bar{n}=\left\\{\frac{V}{N^{{}^{\prime}2}_{\rm
gal}}\sum_{i}^{N_{\rm gal}}\left[w^{(\rm
gal)}(\textbf{{x}}_{i})\right]^{2}\right\\}^{-1}.$ (4.13)

The FFT-based estimator of the multipole 3PCFs is given by (Sugiyama et al.,
2019) (see also Scoccimarro, 2015; Slepian & Eisenstein, 2016)

$\displaystyle\widehat{\zeta}_{\ell_{1}\ell_{2}\ell}(r_{1},r_{2})$
$\displaystyle=$
$\displaystyle\frac{(4\pi)^{2}h_{\ell_{1}\ell_{2}\ell}}{V}\sum_{m_{1}m_{2}m}\left(\begin{smallmatrix}\ell_{1}&\ell_{2}&\ell\\\
m_{1}&m_{2}&m\end{smallmatrix}\right)$ (4.14) $\displaystyle\times$
$\displaystyle\Bigg{[}\int
d^{3}xF_{\ell_{1}m_{1}}(\textbf{{x}};r_{1})F_{\ell_{2}m_{2}}(\textbf{{x}};r_{2})G_{\ell
m}(\textbf{{x}})$
$\displaystyle-\delta_{r_{1}r_{2}}^{\rm(K)}S_{\ell_{1}m_{1};\ell_{2}m_{2};\ell
m}(r_{1})\Bigg{]},$

where

$\displaystyle F_{\ell m}(\textbf{{x}};r)$ $\displaystyle=$ $\displaystyle
i^{\ell}\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\textbf{{k}}\cdot\textbf{{x}}}j_{\ell}(rk)Y_{\ell
m}^{*}(\hat{k})\widetilde{\delta}_{{\rm obs}}(\textbf{{k}}),$ $\displaystyle
G_{\ell m}(\textbf{{x}})$ $\displaystyle=$
$\displaystyle\int\frac{d^{3}k}{(2\pi)^{3}}\,e^{i\textbf{{k}}\cdot\textbf{{x}}}\,\widetilde{\delta}_{{\rm
obs},\ell m}(\textbf{{k}}).$ (4.15)

Note that the shot-noise term only contributes to the 3PCF measurement for the
$r_{1}=r_{2}$ bins, represented by the Kronecker delta $\delta^{(\rm
K)}_{r_{1}r_{2}}$ in Eq. (4.14). To specifically calculate the shot-noise term
in the 3PCF, we first measure the following density field

$\displaystyle N(\textbf{{x}})$ $\displaystyle=$
$\displaystyle\sum_{i}^{N_{\rm
gal}}\left[w_{i}^{\rm(gal)}(\textbf{{x}}_{i})\right]^{2}\delta_{\rm
D}\left(\textbf{{x}}-\textbf{{x}}_{i}^{\rm(gal)}\right)$ (4.16)
$\displaystyle+$ $\displaystyle\left(\frac{N^{\prime}_{\rm
gal}}{N^{\prime}_{\rm ran}}\right)^{2}\sum_{i}^{N_{\rm
ran}}\left[w_{i}^{\rm(ran)}(\textbf{{x}}_{i})\right]^{2}\delta_{\rm
D}\left(\textbf{{x}}-\textbf{{x}}_{i}^{\rm(ran)}\right),$

and divide it by $(N^{\prime}_{\rm gal}/V)$ to have

$\displaystyle\delta_{\rm N}(\textbf{{x}})=(V/N^{\prime}_{\rm
gal})\,N(\textbf{{x}}).$ (4.17)

Then, we calculate the Fourier transform of $\delta_{\rm N}(\textbf{{x}})$ in
the same manner as in Eq. (4.7) and denote it as $\widetilde{\delta}_{\rm
N}(\textbf{{k}})$. Finally, we derive $S_{\ell_{1}m_{1};\ell_{2}m_{2};\ell
m}(r)$ by substituting $\widetilde{\delta}_{\rm N}(\textbf{{k}})$ into the
following equation

$\displaystyle S_{\ell_{1}m_{1};\ell_{2}m_{2};\ell m}(r)$ $\displaystyle=$
$\displaystyle\left(\frac{1}{4\pi r^{2}\Delta
r}\right)\left(\frac{V}{N^{\prime}_{\rm gal}}\right)(-1)^{\ell_{1}+\ell_{2}}$
(4.18) $\displaystyle\times$
$\displaystyle\int\frac{d^{2}\hat{r}}{4\pi}Y^{*}_{\ell_{1}m_{1}}(\hat{r})Y^{*}_{\ell_{2}m_{2}}(\hat{r})\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\textbf{{k}}\cdot\textbf{{r}}}$
$\displaystyle\times$ $\displaystyle\left[\widetilde{\delta}_{\ell
m}(\textbf{{k}})\widetilde{\delta}_{\rm N}^{*}(\textbf{{k}})-S^{(\rm
3PCF)}_{\ell m}(\textbf{{k}})\right],$

where

$\displaystyle S^{\rm(3PCF)}_{\ell m}(\textbf{{k}})$ $\displaystyle=$
$\displaystyle\frac{C_{\rm shot}(\textbf{{k}})}{W_{\rm
mass}^{2}(\textbf{{k}})}\left(\frac{V}{N^{\prime}_{\rm gal}}\right)^{2}$
(4.19) $\displaystyle\times$ $\displaystyle\Bigg{[}\sum_{i}^{N_{\rm
gal}}\left[w^{(\rm gal)}(\textbf{{x}}_{i})\right]^{3}Y_{\ell
m}^{*}\left(\hat{x}^{(\rm gal)}_{i}\right)$ $\displaystyle-$
$\displaystyle\left(\frac{N^{\prime}_{\rm gal}}{N^{\prime}_{\rm
ran}}\right)^{3}\sum_{j}^{N_{\rm ran}}\left[w^{(\rm
ran)}(\textbf{{x}}_{j})\right]^{3}Y_{\ell m}^{*}\left(\hat{x}^{(\rm
ran)}_{j}\right)\Bigg{]}.$

The factor $(1/(4\pi r^{2}\Delta r))$ can be rewritten as

$\displaystyle\frac{1}{4\pi r^{2}\Delta r}=\frac{1}{N_{r}(r)}\frac{N_{\rm
grid}}{V_{\rm FFT}},$ (4.20)

where $V_{\rm FFT}$ is the volume of the Cartesian box in which the galaxies
are placed before the FFT is performed, and $N_{\rm grid}$ is the number of
FFT grid cells.

In the scale range of $80\leq r\leq 150\,h^{-1}\,{\rm Mpc}$, we choose $\Delta
r=5\,h^{-1}\,{\rm Mpc}$ for the 2PCF and $\Delta r=10\,h^{-1}\,{\rm Mpc}$ for
the 3PCF. Considering
$\zeta_{\ell_{1}\ell_{2}\ell}(r_{1},r_{2})=\zeta_{\ell_{2}\ell_{1}\ell}(r_{2},r_{1})$,
the numbers of data bins for the 2PCF and 3PCF multipoles are $15$, $15$,
$36$, $36$, $64$, and $36$ for $\xi_{0}$, $\xi_{2}$, $\zeta_{000}$,
$\zeta_{110}$, $\zeta_{202}$, and $\zeta_{112}$, respectively.

We use the Cartesian coordinates $\textbf{{x}}=\\{x,y,z\\}$ with the $z$-axis
pointing to the north pole to define a cuboid of dimension
$\textbf{{L}}[\,h^{-1}\,{\rm Mpc}]=(L_{x},L_{y},L_{z})$ containing the galaxy
sample; to perform the FFT, each axis of this cuboid is delimited into
$\textbf{{N}}=(N_{x},N_{y},N_{z})$ grids. We then distribute the galaxies on
the FFT grid using the TSC assignment function. We adopt the same values for L
and N that were used by the Fourier space analysis of the two-point statistics
performed by Beutler et al. (2017). They are chosen so that the width of each
grid is $\sim 5\,h^{-1}\,{\rm Mpc}$, which is well below the scales $r\geq
80\,h^{-1}\,{\rm Mpc}$ that we are interested in. We summarise the specific
values of L and N, as well as the survey volume (4.4) and the weighted mean
number density (4.13) computed using these values of L and N in Table 1.

| $(L_{x},L_{y},L_{z})\,[\,h^{-1}\,{\rm Mpc}]$ | $(N_{x},N_{y},N_{z})$ | $V\,[(\,h^{-1}\,{\rm Gpc})^{3}]$ | $\bar{n}/10^{-4}\,[(\,h^{-1}\,{\rm Mpc})^{-3}]$
---|---|---|---|---
NGC at $z_{\rm eff}=0.38$ ($0.2<z<0.5$) | (1350, 2450, 1400) | (250, 460, 260) | 1.51 | 2.65
SGC at $z_{\rm eff}=0.38$ ($0.2<z<0.5$) | (1000, 1900, 1100) | (190, 360, 210) | 0.56 | 2.88
NGC at $z_{\rm eff}=0.61$ ($0.5<z<0.75$) | (1800, 3400, 1900) | (340, 650, 360) | 2.35 | 1.37
SGC at $z_{\rm eff}=0.61$ ($0.5<z<0.75$) | (1000, 2600, 1500) | (190, 500, 280) | 0.87 | 1.29
Table 1:  The length of each side of the cube containing the observed
galaxies, defined for performing FFTs, and the number of grids on which the
cube is delimited are shown for the four BOSS samples (Section 4.1). Also
shown are the survey volume (4.4) and the mean galaxy number density (4.13),
calculated using the values of these parameters.

###  4.3 Window function corrections

When measuring 2PCFs and 3PCFs in configuration space from galaxy data, if we
directly measure their angle-averaged multipole components, we can not
eliminate the effect of the window function (Appendix A in Sugiyama et al.,
2021). The FFT-based estimators introduced in Section (4.2) are a typical
example of this, but even when measuring multipole 3PCFs without using the
FFT, we need to be aware of the window function effect (Slepian & Eisenstein,
2015, 2018). Since the window function $W(\textbf{{x}})$ characterising the
geometry of the observed region can be estimated as
$W(\textbf{{x}})=(V/N^{\prime}_{\rm ran})\,R(\textbf{{x}})$, we can
quantitatively estimate the corrections due to the window function by
measuring the multipole 2PCFs and 3PCFs from the random catalogue.

For the 2PCF, we compute

$\displaystyle Q_{\ell}(r)$ $\displaystyle=$
$\displaystyle\frac{(4\pi)}{V}\sum_{m}\int\frac{d^{2}\hat{r}}{4\pi}Y_{\ell
m}(\hat{r})\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\textbf{{k}}\cdot\textbf{{r}}}$
(4.21) $\displaystyle\times$ $\displaystyle\bigg{[}\widetilde{W}_{\ell
m}(\textbf{{k}})\widetilde{W}^{*}(\textbf{{k}})-S^{\rm(w)}_{\ell
m}(\textbf{{k}})\bigg{]}.$

where $\widetilde{W}_{\ell m}(\textbf{{k}})$ is the Fourier transform of
$(V/N^{\prime}_{\rm ran})R_{\ell m}(\textbf{{x}})$ computed in the same manner
as in Eq. (4.7),
$\widetilde{W}(\textbf{{k}})=\sqrt{4\pi}\,\widetilde{W}_{00}(\textbf{{k}})$,
and the shot-noise term is given by

$\displaystyle S^{\rm(w)}_{\ell m}(\textbf{{k}})$ $\displaystyle=$
$\displaystyle\frac{C_{\rm shot}(\textbf{{k}})}{W_{\rm
mass}^{2}(\textbf{{k}})}\left(\frac{V}{N^{\prime}_{\rm ran}}\right)^{2}$
(4.22) $\displaystyle\times$ $\displaystyle\Bigg{[}\sum_{i}^{N_{\rm
ran}}\left[w^{(\rm ran)}(\textbf{{x}}_{i})\right]^{2}Y_{\ell
m}^{*}\left(\hat{x}^{(\rm ran)}_{i}\right)\Bigg{]}.$

Then, we have the theoretical model of $\xi_{\ell}(r)$ taking the survey
window effect into account as follows (Wilson et al., 2017; Beutler et al.,
2017):

$\displaystyle\xi^{(\rm
w)}_{\ell}(r)=(2\ell+1)\sum_{\ell_{1}\ell_{2}}\left(\begin{smallmatrix}\ell_{1}&\ell_{2}&\ell\\\
0&0&0\end{smallmatrix}\right)^{2}Q_{\ell_{1}}(r)\,\xi_{\ell_{2}}(r).$ (4.23)

For the 3PCF, we compute

$\displaystyle Q_{\ell_{1}\ell_{2}\ell}(r_{1},r_{2})$ $\displaystyle=$
$\displaystyle\frac{(4\pi)^{2}h_{\ell_{1}\ell_{2}\ell}}{V}\sum_{m_{1}m_{2}m}\left(\begin{smallmatrix}\ell_{1}&\ell_{2}&\ell\\\
m_{1}&m_{2}&m\end{smallmatrix}\right)$ (4.24) $\displaystyle\times$
$\displaystyle\Bigg{[}\int d^{3}xF_{\ell_{1}m_{1}}^{({\rm
w})}(\textbf{{x}};r_{1})F_{\ell_{2}m_{2}}^{({\rm
w})}(\textbf{{x}};r_{2})G_{\ell m}^{({\rm w})}(\textbf{{x}})$
$\displaystyle-\delta_{r_{1}r_{2}}^{\rm(K)}S^{({\rm
w})}_{\ell_{1}m_{1};\ell_{2}m_{2};\ell m}(r_{1})\Bigg{]},$

where

$\displaystyle F_{\ell m}^{({\rm w})}(\textbf{{x}};r)$ $\displaystyle=$
$\displaystyle
i^{\ell}\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\textbf{{k}}\cdot\textbf{{x}}}j_{\ell}(rk)Y_{\ell
m}^{*}(\hat{k})\widetilde{W}(\textbf{{k}}),$ $\displaystyle G_{\ell m}^{({\rm
w})}(\textbf{{x}})$ $\displaystyle=$
$\displaystyle\int\frac{d^{3}k}{(2\pi)^{3}}\,e^{i\textbf{{k}}\cdot\textbf{{x}}}\,\widetilde{W}_{\ell
m}(\textbf{{k}}).$ (4.25)

The shot-noise term is given by

$\displaystyle S^{{\rm(w)}}_{\ell_{1}m_{1};\ell_{2}m_{2};\ell m}(r)$
$\displaystyle=$ $\displaystyle\left(\frac{1}{4\pi r^{2}\Delta
r}\right)\left(\frac{V}{N^{\prime}_{\rm ran}}\right)(-1)^{\ell_{1}+\ell_{2}}$
(4.26) $\displaystyle\times$
$\displaystyle\int\frac{d^{2}\hat{r}}{4\pi}Y^{*}_{\ell_{1}m_{1}}(\hat{r})Y^{*}_{\ell_{2}m_{2}}(\hat{r})\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\textbf{{k}}\cdot\textbf{{r}}}$
$\displaystyle\times$ $\displaystyle\left[\widetilde{W}_{\ell
m}(\textbf{{k}})\widetilde{\delta}_{\rm N}^{{(\rm w)*}}(\textbf{{k}})-S^{(\rm
3PCF,w)}_{\ell m}(\textbf{{k}})\right],$

where

$\displaystyle S^{\rm(3PCF,w)}_{\ell m}(\textbf{{k}})$ $\displaystyle=$
$\displaystyle\frac{C_{\rm shot}(\textbf{{k}})}{W_{\rm
mass}^{2}(\textbf{{k}})}\left(\frac{V}{N^{\prime}_{\rm ran}}\right)^{2}$
(4.27) $\displaystyle\times$ $\displaystyle\Bigg{[}\sum_{i}^{N_{\rm
ran}}\left[w^{(\rm ran)}(\textbf{{x}}_{i})\right]^{3}Y_{\ell
m}^{*}\left(\hat{x}^{(\rm ran)}_{i}\right)\Bigg{]},$

and $\widetilde{\delta}_{\rm N}^{\rm(w)}(\textbf{{k}})$ is the Fourier
transform of

$\displaystyle\delta_{\rm N}^{(\rm
w)}(\textbf{{x}})=\left(\frac{V}{N^{\prime}_{\rm ran}}\right)\sum_{i}^{N_{\rm
ran}}\left[w_{i}^{\rm(ran)}\right]^{2}\delta_{\rm
D}\left(\textbf{{x}}-\textbf{{x}}_{i}^{(\rm ran)}\right).$ (4.28)

Then, we have the theoretical model of
$\zeta_{\ell_{1}\ell_{2}\ell}(r_{1},r_{2})$ taking the survey window effect
into account as follows (Sugiyama et al., 2019; Sugiyama et al., 2021):

$\displaystyle\zeta^{(\rm w)}_{\ell_{1}\ell_{2}\ell}(r_{1},r_{2})$
$\displaystyle=$
$\displaystyle(4\pi)\sum_{\ell^{\prime}_{1}+\ell^{\prime}_{2}+\ell^{\prime}={\rm
even}}\ \
\sum_{\ell^{\prime\prime}_{1}+\ell^{\prime\prime}_{2}+\ell^{\prime\prime}={\rm
even}}$ (4.29) $\displaystyle\times$
$\displaystyle\left\\{\begin{smallmatrix}\ell^{\prime\prime}_{1}&\ell^{\prime\prime}_{2}&\ell^{\prime\prime}\\\
\ell^{\prime}_{1}&\ell^{\prime}_{2}&\ell^{\prime}\\\
\ell_{1}&\ell_{2}&\ell\end{smallmatrix}\right\\}\left[\frac{h_{\ell_{1}\ell_{2}\ell}h_{\ell_{1}\ell^{\prime}_{1}\ell^{\prime\prime}_{1}}h_{\ell_{2}\ell^{\prime}_{2}\ell^{\prime\prime}_{2}}h_{\ell\ell^{\prime}\ell^{\prime\prime}}}{h_{\ell^{\prime}_{1}\ell^{\prime}_{2}\ell^{\prime}}h_{\ell^{\prime\prime}_{1}\ell^{\prime\prime}_{2}\ell^{\prime\prime}}}\right]$
$\displaystyle\times$ $\displaystyle
Q_{\ell^{\prime\prime}_{1}\ell^{\prime\prime}_{2}\ell^{\prime\prime}}(r_{1},r_{2})\,\zeta_{\ell^{\prime}_{1}\ell^{\prime}_{2}\ell^{\prime}}(r_{1},r_{2}),$

where the bracket with $9$ multipole indices, $\\{\dots\\}$, denotes the
Wigner-$9$j symbol. In the likelihood fitting performed in Section 9, we use
$\xi_{\ell}^{(\rm w)}$ and $\zeta_{\ell_{1}\ell_{2}\ell}^{(\rm w)}$ to compare
the measured multipole 2PCF and 3PCF estimators with the theoretical models
given in Eqs. (3.9) and (3.12). In this paper, we ignore the contribution from
the integral constraint (Peacock & Nicholson, 1991) for both the 2PCF and the
3PCF.

In the 2PCF case, the correction equation for the window function effect shown
in Eq. (4.23) calculates only the three multipole components for both
$Q_{\ell_{1}}$ and $\xi_{\ell_{2}}$, i.e., $\ell_{1},\ell_{2}=0,2,4$. The
reason is that our analysis focuses only on large scales above
$80\,h^{-1}\,{\rm Mpc}$, where the linear theory is dominant, and the linear
Kaiser effect gives only up to the hexadecapole $\ell=4$. For the window
correction formula of the 3PCF (4.29), Sugiyama et al. (2021) examined in
detail which multipole components contribute to the observed estimator (4.14)
and to what extent, for the NGC sample at $0.4<z<0.6$, and showed that a
finite number of multipole components can correct for the window effect on the
3PCF with sufficiently good accuracy. Assuming that this result is not
significantly different for the other BOSS samples, we calculate a total of
$14$ multipole components for both
$Q_{\ell_{1}^{\prime\prime}\ell_{2}^{\prime\prime}\ell^{\prime\prime}}$ and
$\zeta_{\ell_{1}^{\prime}\ell_{2}^{\prime}\ell^{\prime}}$ as follows:
$(\ell_{1},\ell_{2},\ell)=(0,0,0)$, $(1,1,0)$, $(2,2,0)$, $(3,3,0)$ and
$(4,4,0)$ for the monopole 3PCF ($\ell=0$), and
$(\ell_{1},\ell_{2},\ell)=(0,2,2)$, $(1,1,2)$, $(2,0,2)$, $(1,3,2)$,
$(2,2,2)$, $(3,1,2)$, $(2,4,2)$, $(3,3,2)$ and $(4,2,2)$ for the quadrupole
3PCF ($\ell=2$).

Figures 3 and 4 plot the $13$ window 3PCF multipoles normalized by $Q_{000}$
as a function of $r_{2}$ after fixing $r_{1}$ to $60$ and $120\,h^{-1}\,{\rm
Mpc}$. For the monopole components ($Q_{110}$, $Q_{220}$, $Q_{330}$, and
$Q_{440}$), we find that the window 3PCF multipoles measured at different
redshift bins in each sky region (NGC or SGC) behave similarly (see, for
example, the solid blue and dashed orange lines). On the other hand, for the
quadrupole component, we see that the four BOSS samples may behave
differently. The first few terms of the monopole and quadrupole components,
such as $Q_{110}$, $Q_{220}$, $Q_{202}$, $Q_{112}$, and $Q_{022}$, have values
of ${\cal O}(0.01)-{\cal O}(0.1)$, while the higher-order terms have values of
${\cal O}(0.01)$ or less. Therefore, we can conclude that the higher-order
window 3PCF multipoles have no significant effect on the final $\zeta^{(\rm
w)}_{\ell_{1}\ell_{2}\ell}(r_{1},r_{2})$, as long as we measure the first few
terms of the monopole and quadrupole components, i.e., $\zeta^{(\rm
w)}_{000}(r_{1},r_{2})$, $\zeta^{(\rm w)}_{110}(r_{1},r_{2})$, $\zeta^{(\rm
w)}_{202}(r_{1},r_{2})$, and $\zeta^{(\rm w)}_{112}(r_{1},r_{2})$.

Figures 5 and 6 plot the theoretical predictions for the 3PCF multipoles,
including window function effects, corresponding to the four BOSS samples.
These calculations assume the $\Lambda$CDM and linear bias as in Figures 1 and
2, with redshifts of $0.38$ and $0.61$. As the value of $r_{1}$ increases, the
difference between NGC and SGC due to the window function effect becomes more
considerable.

To quantitatively estimate the extent to which the multipole component of
interest, $\zeta_{\ell_{1}\ell_{2}\ell}^{(\rm w)}$, is affected by the other
multipole components,
$\zeta_{\ell_{1}^{\prime}\ell_{2}^{\prime}\ell^{\prime}}$, through window
function effects, we compute the following quantities (Sugiyama et al., 2021):

$\displaystyle\Delta\bar{\zeta}_{\ell^{\prime}_{1}\ell^{\prime}_{2}\ell^{\prime}}$
$\displaystyle=$
$\displaystyle\frac{\mbox{Sum}\Big{[}\Delta\zeta_{\ell^{\prime}_{1}\ell^{\prime}_{2}\ell^{\prime}}^{\ell_{1}\ell_{2}\ell}/Q_{000}\Big{]}}{\mbox{Sum}\left[\zeta_{\ell_{1}\ell_{2}\ell}^{(\rm
w)}/Q_{000}\right]}$

with

$\displaystyle\Delta\zeta_{\ell^{\prime}_{1}\ell^{\prime}_{2}\ell^{\prime}}^{\ell_{1}\ell_{2}\ell}(r_{1},r_{2})$
$\displaystyle=$ $\displaystyle(4\pi)\ \
\sum_{\ell^{\prime\prime}_{1}+\ell^{\prime\prime}_{2}+\ell^{\prime\prime}={\rm
even}}\left\\{\begin{smallmatrix}\ell^{\prime\prime}_{1}&\ell^{\prime\prime}_{2}&\ell^{\prime\prime}\\\
\ell^{\prime}_{1}&\ell^{\prime}_{2}&\ell^{\prime}\\\
\ell_{1}&\ell_{2}&\ell\end{smallmatrix}\right\\}$ (4.31) $\displaystyle\times$
$\displaystyle\left[\frac{h_{\ell_{1}\ell_{2}\ell}h_{\ell_{1}\ell^{\prime}_{1}\ell^{\prime\prime}_{1}}h_{\ell_{2}\ell^{\prime}_{2}\ell^{\prime\prime}_{2}}h_{\ell\ell^{\prime}\ell^{\prime\prime}}}{h_{\ell^{\prime}_{1}\ell^{\prime}_{2}\ell^{\prime}}h_{\ell^{\prime\prime}_{1}\ell^{\prime\prime}_{2}\ell^{\prime\prime}}}\right]$
$\displaystyle\times$ $\displaystyle
Q_{\ell^{\prime\prime}_{1}\ell^{\prime\prime}_{2}\ell^{\prime\prime}}(r_{1},r_{2})\,\zeta_{\ell^{\prime}_{1}\ell^{\prime}_{2}\ell^{\prime}}(r_{1},r_{2})$

and

$\displaystyle\mbox{Sum}\left[\zeta_{\ell_{1}\ell_{2}\ell}\right]=\begin{cases}\sum_{r_{1}\geq
r_{2}}\zeta_{\ell_{1}\ell_{2}\ell}(r_{1},r_{2})\quad\mbox{for
$\ell_{1}=\ell_{2}$};\\\ \sum_{r_{1},r_{2}}\
\zeta_{\ell_{1}\ell_{2}\ell}(r_{1},r_{2})\quad\mbox{for
$\ell_{1}\neq\ell_{2}$},\end{cases}$ (4.32)

where $\Delta\bar{\zeta}_{\ell^{\prime}_{1}\ell^{\prime}_{2}\ell^{\prime}}$
satisfies
$\sum_{\ell_{1}^{\prime}\ell_{2}^{\prime}\ell^{\prime}}\Delta\bar{\zeta}_{\ell^{\prime}_{1}\ell^{\prime}_{2}\ell^{\prime}}=1$,
and the summation is performed in the range of $80\leq r\leq 150\,h^{-1}\,{\rm
Mpc}$ which we use for our data analysis.

Table 2 summarises the
$\Delta\bar{\zeta}_{\ell^{\prime}_{1}\ell^{\prime}_{2}\ell^{\prime}}$ results
calculated from Eq. (LABEL:Eq:Delta_zeta) for the four BOSS samples.
Naturally, the multipole component that is the same as the target one has the
largest contribution. For example, for $\zeta_{000}^{(\rm w)}$ at $z_{\rm
eff}=0.38$ in NGC, $95.34\%$ of the contribution comes from $\zeta_{000}$. For
all four samples, multipole components other than the measured one have
positive or negative values, and their overall contribution is about $5-10\%$.
As expected, the contributions of higher-order components such as
$\zeta_{330}$, $\zeta_{440}$, and $\zeta_{332}$ are mostly below $0.5\%$.
Therefore, we conclude that the window function correction equation in Eq.
(4.29) can account for the window function effect on the 3PCF in BOSS with
sufficient accuracy, even if it is truncated at a finite number of 14
multipole components used in this work.

We note here the importance of $\Delta\bar{\zeta}_{112}$, which includes the
$M\neq 0$ modes of Scoccimarro et al. (1999)’ decomposition method in the
correction for window function effects: it gives a contribution comparable to
$\Delta\bar{\zeta}_{202}$ and $\Delta\bar{\zeta}_{022}$, which include only
the $M=0$ mode, and tends to have the opposite sign to that of
$\Delta\bar{\zeta}_{202}$ and $\Delta\bar{\zeta}_{022}$. Therefore, failure to
properly account for effects such as $\Delta\bar{\zeta}_{112}$ that include
the $M\neq 0$ modes may result in an error of $\sim 5\%$ in the correction for
the window function effect.

$z_{\rm eff}=0.38$ ($0.2<z<0.5$)
---
| NGC |  | SGC
| $\Delta\bar{\zeta}_{\ell^{\prime}_{1}\ell^{\prime}_{2}\ell^{\prime}}\,[\%]$ | $\zeta^{(\rm w)}_{000}$ | $\zeta^{(\rm w)}_{110}$ | $\zeta^{(\rm w)}_{202}$ | $\zeta^{(\rm w)}_{112}$ |  | $\zeta^{(\rm w)}_{000}$ | $\zeta^{(\rm w)}_{110}$ | $\zeta^{(\rm w)}_{202}$ | $\zeta^{(\rm w)}_{112}$
monopole ($\ell=0$) | $\Delta\bar{\zeta}_{000}$ | $\mathbf{95.34}$ | $\mathbf{1.59}$ | $\mathbf{-1.18}$ | $\mathbf{0.97}$ |  | $\mathbf{86.41}$ | $\mathbf{2.43}$ | $-0.15$ | $0.19$
| $\Delta\bar{\zeta}_{110}$ | $\mathbf{9.39}$ | $\mathbf{102.60}$ | $\mathbf{1.41}$ | $\mathbf{-3.77}$ |  | $\mathbf{13.66}$ | $\mathbf{98.60}$ | $0.25$ | $\mathbf{-0.50}$
| $\Delta\bar{\zeta}_{220}$ | $-0.09$ | $\mathbf{-0.53}$ | $\mathbf{1.14}$ | $0.09$ |  | $-0.12$ | $\mathbf{-0.79}$ | $0.13$ | $0.01$
| $\Delta\bar{\zeta}_{330}$ | $0.26$ | $0.21$ | $0.08$ | $-0.01$ |  | $0.41$ | $0.35$ | $0.02$ | $-0.02$
| $\Delta\bar{\zeta}_{440}$ | $0.06$ | $0.00$ | $0.06$ | $-0.00$ |  | $0.09$ | $0.01$ | $-0.00$ | $-0.00$
quadrupole ($\ell=2$) | $\Delta\bar{\zeta}_{202}$ | $\mathbf{-3.75}$ | $\mathbf{0.59}$ | $\mathbf{90.14}$ | $\mathbf{2.07}$ |  | $-0.43$ | $0.10$ | $\mathbf{89.28}$ | $\mathbf{2.38}$
| $\Delta\bar{\zeta}_{112}$ | $\mathbf{4.36}$ | $\mathbf{-2.89}$ | $\mathbf{4.21}$ | $\mathbf{96.78}$ |  | $\mathbf{0.86}$ | $-0.38$ | $\mathbf{5.06}$ | $\mathbf{92.64}$
| $\Delta\bar{\zeta}_{022}$ | $\mathbf{-4.57}$ | $\mathbf{0.81}$ | $\mathbf{0.50}$ | $\mathbf{2.75}$ |  | $-0.48$ | $0.14$ | $\mathbf{0.51}$ | $\mathbf{3.16}$
| $\Delta\bar{\zeta}_{312}$ | $-0.10$ | $-0.38$ | $\mathbf{1.95}$ | $0.13$ |  | $-0.05$ | $-0.05$ | $\mathbf{2.87}$ | $0.26$
| $\Delta\bar{\zeta}_{222}$ | $0.04$ | $0.04$ | $0.38$ | $0.07$ |  | $-0.02$ | $-0.01$ | $-0.04$ | $0.12$
| $\Delta\bar{\zeta}_{132}$ | $\mathbf{-0.56}$ | $\mathbf{-1.54}$ | $0.14$ | $0.23$ |  | $-0.19$ | $-0.19$ | $0.18$ | $\mathbf{0.66}$
| $\Delta\bar{\zeta}_{422}$ | $-0.21$ | $-0.23$ | $\mathbf{1.02}$ | $0.21$ |  | $-0.02$ | $-0.08$ | $\mathbf{1.53}$ | $0.31$
| $\Delta\bar{\zeta}_{332}$ | $0.10$ | $0.08$ | $0.02$ | $0.10$ |  | $-0.03$ | $-0.02$ | $0.10$ | $0.24$
| $\Delta\bar{\zeta}_{242}$ | $-0.29$ | $-0.34$ | $0.12$ | $0.38$ |  | $-0.08$ | $-0.14$ | $0.26$ | $\mathbf{0.55}$
|  |  |  |  |  |  |  |  |  |
$z_{\rm eff}=0.61$ ($0.5<z<0.75$)
| NGC |  | SGC
| $\Delta\bar{\zeta}_{\ell^{\prime}_{1}\ell^{\prime}_{2}\ell^{\prime}}\,[\%]$ | $\zeta^{(\rm w)}_{000}$ | $\zeta^{(\rm w)}_{110}$ | $\zeta^{(\rm w)}_{202}$ | $\zeta^{(\rm w)}_{112}$ |  | $\zeta^{(\rm w)}_{000}$ | $\zeta^{(\rm w)}_{110}$ | $\zeta^{(\rm w)}_{202}$ | $\zeta^{(\rm w)}_{112}$
monopole ($\ell=0$) | $\Delta\bar{\zeta}_{000}$ | $\mathbf{96.00}$ | $\mathbf{1.77}$ | $\mathbf{-2.00}$ | $\mathbf{1.78}$ |  | $\mathbf{89.49}$ | $\mathbf{2.40}$ | $\mathbf{-1.29}$ | $\mathbf{1.19}$
| $\Delta\bar{\zeta}_{110}$ | $\mathbf{10.43}$ | $\mathbf{104.18}$ | $\mathbf{2.92}$ | $\mathbf{-6.81}$ |  | $\mathbf{13.59}$ | $\mathbf{101.41}$ | $\mathbf{2.10}$ | $\mathbf{-4.36}$
| $\Delta\bar{\zeta}_{220}$ | $-0.08$ | $\mathbf{-0.60}$ | $\mathbf{1.99}$ | $0.01$ |  | $-0.06$ | $\mathbf{-0.80}$ | $\mathbf{1.33}$ | $-0.04$
| $\Delta\bar{\zeta}_{330}$ | $0.25$ | $0.21$ | $0.04$ | $0.01$ |  | $0.35$ | $0.31$ | $0.00$ | $-0.01$
| $\Delta\bar{\zeta}_{440}$ | $0.06$ | $0.01$ | $0.05$ | $-0.00$ |  | $0.08$ | $0.01$ | $0.02$ | $-0.00$
quadrupole ($\ell=2$) | $\Delta\bar{\zeta}_{202}$ | $\mathbf{-6.26}$ | $\mathbf{1.21}$ | $\mathbf{86.94}$ | $\mathbf{2.72}$ |  | $\mathbf{-3.73}$ | $\mathbf{0.84}$ | $\mathbf{86.45}$ | $\mathbf{2.94}$
| $\Delta\bar{\zeta}_{112}$ | $\mathbf{7.84}$ | $\mathbf{-5.12}$ | $\mathbf{5.36}$ | $\mathbf{97.84}$ |  | $\mathbf{5.09}$ | $\mathbf{-3.31}$ | $\mathbf{5.98}$ | $\mathbf{94.88}$
| $\Delta\bar{\zeta}_{022}$ | $\mathbf{-7.65}$ | $\mathbf{1.65}$ | $0.48$ | $\mathbf{3.64}$ |  | $\mathbf{-4.54}$ | $\mathbf{1.14}$ | $0.44$ | $\mathbf{3.93}$
| $\Delta\bar{\zeta}_{312}$ | $-0.01$ | $\mathbf{-0.54}$ | $\mathbf{2.25}$ | $0.09$ |  | $0.01$ | $-0.30$ | $\mathbf{2.89}$ | $0.17$
| $\Delta\bar{\zeta}_{222}$ | $0.04$ | $0.04$ | $\mathbf{0.69}$ | $0.07$ |  | $-0.00$ | $0.01$ | $0.42$ | $0.11$
| $\Delta\bar{\zeta}_{132}$ | $-0.28$ | $\mathbf{-2.63}$ | $0.13$ | $-0.02$ |  | $-0.06$ | $\mathbf{-1.66}$ | $0.15$ | $0.25$
| $\Delta\bar{\zeta}_{422}$ | $-0.18$ | $-0.13$ | $\mathbf{0.98}$ | $0.20$ |  | $-0.07$ | $-0.03$ | $\mathbf{1.21}$ | $0.27$
| $\Delta\bar{\zeta}_{332}$ | $0.08$ | $0.09$ | $0.03$ | $0.10$ |  | $-0.02$ | $0.01$ | $0.09$ | $0.19$
| $\Delta\bar{\zeta}_{242}$ | $-0.24$ | $-0.15$ | $0.13$ | $0.37$ |  | $-0.12$ | $-0.03$ | $0.21$ | $0.48$
Table 2:  Contributions of other 3PCF multipole components to the observed
3PCF multipole components, as manifested through the effect of the window
function, are shown for the four BOSS samples. When the contribution to the
final result exceeds $0.5\%$, it is written in bold. The value of the same
multipole component $\Delta\bar{\zeta}_{\ell_{1}\ell_{2}\ell}$
(LABEL:Eq:Delta_zeta) as the measured $\zeta_{\ell_{1}\ell_{2}\ell}^{(\rm w)}$