(4.29) is larger (smaller) than $100\%$, when the total contribution from all
the other multipole components is negative (positive).

Figure 3: The monopole and quadrupole components of the window 3PCF (4.24),
$Q_{000}$, $Q_{110}$, $Q_{220}$, $Q_{330}$, $Q_{440}$, $Q_{202}$, $Q_{112}$,
and $Q_{022}$, measured from the four BOSS samples are shown as a function of
$r_{2}$ after fixing $r_{1}$ to $60\,h^{-1}\,{\rm Mpc}$ (left) and
$120\,h^{-1}\,{\rm Mpc}$ (right). Figure 4: Same as Figure 3, except that the
higher-order quadrupole components of the window 3PCF, $Q_{312}$, $Q_{222}$,
$Q_{132}$, $Q_{422}$, $Q_{332}$, and $Q_{242}$, are shown.

Figure 5: The monopole 3PCFs, $\zeta_{000}$ (left) and $\zeta_{110}$ (right),
that include the window function effect (4.29) are shown for the four BOSS
samples. The results are plotted as a function of $r_{2}$ after fixing $r_{1}$
to $50$, $80$, $90$, $100$, $130\,h^{-1}\,{\rm Mpc}$ from top to bottom
panels. 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 6: Same as Figure 5, except that the quadrupole 3PCFs, $\zeta_{202}$
(left) and $\zeta_{112}$ (right), are shown.

##  5 Covariance matrix

We estimate the covariance matrix from the $2048$ Patchy mock catalogues
described in Section 4.1. Let $\textbf{{d}}^{(r)}$ be the data vector measured
from the $r$-th catalogue, and $\overline{\textbf{{d}}}=(1/N_{\rm
s})\sum_{r=1}^{N_{\rm s}}\textbf{{d}}^{(r)}$ be its mean value, then the
covariance matrix of the data vector is given by

$\displaystyle\textbf{{C}}=\frac{1}{N_{\rm s}-1}\sum_{r=1}^{N_{\rm
s}}\left(\textbf{{d}}^{(r)}-\overline{\textbf{{d}}}\right)\left(\textbf{{d}}^{(r)}-\overline{\textbf{{d}}}\right)^{T}.$
(5.1)

where $N_{\rm s}=2048$ is the number of the Patchy mock catalogues.

###  5.1 Effects of a finite number of mocks

The covariance matrix C inferred from the mock catalogues suffers from noise
due to the finite number of mocks, which directly leads to an increase in the
uncertainty of the cosmological parameters (Hartlap et al., 2007; Taylor et
al., 2013; Dodelson & Schneider, 2013; Percival et al., 2014; Taylor &
Joachimi, 2014). This effect is decomposed into two factors. First, the
inverse covariance matrix, $\textbf{{C}}^{-1}$, provides a biased estimate of
the true inverse covariance matrix. To correct this bias, we rescale the
inverse covariance matrix as (Hartlap et al., 2007)

$\displaystyle\textbf{{C}}_{\rm Hartlap}^{-1}=\left(\frac{N_{\rm s}-N_{\rm
b}-2}{N_{\rm s}-1}\right)\textbf{{C}}^{-1},$ (5.2)

where the pre-factor on the right-hand side, $(N_{\rm s}-N_{\rm b}-2)/(N_{\rm
s}-1)$, is the so-called “Hartlap” factor, and $N_{\rm b}$ is the number of
data bins. Second, we need to consider the propagation of the error in the
covariance matrix to the error on the estimated parameters. This effect is
corrected by multiplying the final result of the parameter errors by the
following factor (Percival et al., 2014)

$\displaystyle M_{1}=\sqrt{\frac{1+B(N_{\rm b}-N_{\rm p})}{1+A+B(N_{\rm
p}+1)}}$ (5.3)

with

$\displaystyle A$ $\displaystyle=$ $\displaystyle\frac{2}{(N_{\rm s}-N_{\rm
b}-1)(N_{\rm s}-N_{\rm b}-4)}$ $\displaystyle B$ $\displaystyle=$
$\displaystyle\frac{N_{\rm s}-N_{\rm b}-2}{(N_{\rm s}-N_{\rm b}-1)(N_{\rm
s}-N_{\rm b}-4)},$ (5.4)

where $N_{\rm p}$ is the number of parameters.

The derivation of the Hartlap factor (5.2) assumes that the data vector
follows a Gaussian distribution. On the other hand, Sellentin & Heavens (2016)
shows that in covariance matrix estimates from simulations, the data vector
follows a multivariate $t$-distribution. When the number of simulations is
sufficiently larger than the number of data bins, this $t$-distribution
approaches a Gaussian distribution (Heavens et al., 2017), and the present
analysis satisfies this condition. The reason is that the number of the Patchy
mocks we use to estimate the covariance matrix is $2048$, while the maximum
number of data in our analysis is $202$ (Section 6.2). In addition, the
derivation of the $M_{1}$ factor (5.3) also assumes the Gaussian distribution
of the data vector, but there is no known value for the correction factor that
corresponds to $M_{1}$ in the Sellentin & Heavens (2016)’s method. Therefore,
in this paper, we have decided to use Eqs. (5.2) and (5.3) to correct the
uncertainty in parameter estimation due to a finite number of simulations (see
also e.g., Percival et al., 2022).

We can therefore evaluate the effect of a finite number of mocks on the final
error estimation using the square root of the Hartlap factor multiplied by the
$M_{1}$ factor (Percival et al., 2014),

$\displaystyle M_{2}=\sqrt{\frac{N_{\rm s}-1}{N_{\rm s}-N_{\rm b}-2}}M_{1}.$
(5.5)

Note that this $M_{2}$ factor is not used in the actual analysis. It is
essential to increase the number of simulations and reduce the number of data
bins to keep the value of $M_{2}$ as close to $1$ as possible for a
conservative analysis. The reason is that the Hartlap and $M_{1}$ factors
cannot be always accurately correct for parameter errors for any number of
simulations. For example, using both monopole and quadrupole components of the
2PCF and 3PCF, as in this paper, Sugiyama et al. (2021) performed an
anisotropic BAO analysis with the AP effect. The result showed that the error
in the angular diameter distance for $M_{2}=1.32$ is underestimated by about
$10\%$ compared to the case for $M_{2}=1.06$ by changing the number of
simulations. We will calculate the $M_{2}$ factor in Section 6.4 and summarise
the results in Table 3, where $M_{2}\sim 1.1$, indicating that our analysis
achieves $M_{2}$ values sufficiently close to $1$.

###  5.2 Correlation matrix

Figure 7: The correlation matrices of the monopole and quadrupole 2PCFs
($\xi_{0}$ and $\xi_{2}$), the first two monopole 3PCFs ($\zeta_{000}$ and
$\zeta_{110}$), and the first two quadrupole 3PCFs ($\zeta_{202}$ and
$\zeta_{112}$) are shown for the four BOSS samples. For simplicity of the
figure, only the results for the $r_{1}=r_{2}$ case multipole 3PCFs, i.e.
$\zeta(r_{1},r_{2}=r_{1})$, are plotted, but in the actual analysis (Section
9), the data bins for the $r_{1}\neq r_{2}$ case are also used. The plotted
scale range is $80\leq r\leq 150\,h^{-1}\,{\rm Mpc}$, and the $r$-bin width is
$\Delta r=10\,h^{-1}\,{\rm Mpc}$.

The $(i,j)$ elements of the correlation matrix is computed from the covariance
matrix as

$\displaystyle r_{ij}=\frac{C_{ij}}{\sqrt{C_{ii}C_{jj}}}.$ (5.6)

Considering the data vector
$\textbf{{d}}=\\{\xi_{0},\xi_{2},\zeta_{000},\zeta_{110},\zeta_{202},\zeta_{112}\\}$,
we show the results of the correlation matrix for the four BOSS samples in
Figure 7. To simplify the figure, we only plot the results for the diagonal
component of the 3PCF multipoles, i.e.,
$\zeta_{\ell_{1}\ell_{2}\ell}(r_{1},r_{2}=r_{1})$. The range of scales shown
in the figure is $80\leq r\leq 150\,h^{-1}\,{\rm Mpc}$, and the width of the
$r$-bin is $\Delta r=10\,h^{-1}\,{\rm Mpc}$. The four samples show similar
results, and we summarise the overall features below. First, the monopole 2PCF
and the monopole 3PCFs have a moderate correlation ($0.25<r_{ij}<0.5$); the
same is true for the quadrupole 2PCF and the quadrupole 3PCFs. Next, the first
two terms of the monopole 3PCFs ($\zeta_{000}$ and $\zeta_{110}$) are strongly
correlated with each other ($0.5<r_{ij}<0.75$); on the other hand, the first
two quadrupole 3PCFs ($\zeta_{202}$ and $\zeta_{112}$) are weakly correlated
($0.0<r_{ij}<0.25$). This result indicates that $\zeta_{202}$ and
$\zeta_{112}$ have independent information from each other. These results are
consistent with the results in the bispectrum presented by Sugiyama et al.
(2019).

###  5.3 Standard deviation

Figure 8: The mean values and standard deviations of $\xi_{0}$, $\xi_{2}$,
$\zeta_{000}$, $\zeta_{110}$, $\zeta_{202}$, and $\zeta_{112}$ calculated from
the $2048$ Patchy mock catalogues. The results are plotted at the two
redshifts, $z=0.38$ (magenta) and $0.61$ (blue), for the NGC sample. The
magenta dashed lines are the standard deviation at $z=0.61$ multiplied by
$\sqrt{V_{z=0.61}/V_{z=0.38}}$ and normalized to have the same survey volume
as the sample at $z=0.38$, where the survey volumes at $z=0.38$ and $z=0.61$,
$V_{z=0.38}$ and $V_{z=0.61}$, respectively, are given in Table 1. For
simplicity of the figure, only the results in the $r_{1}=r_{2}$ case for the
3PCF are plotted.

The standard deviation is given by the square root of the diagonal components
of the covariance matrix: i.e., $\sqrt{C_{ii}}$. Figure 8 shows the mean and
standard deviation of $\xi_{\ell}(r)$ and
$\zeta_{\ell_{1}\ell_{2}\ell}(r_{1},r_{2})$ calculated from the $2048$ Patchy
mock simulations. The mock data used are the NGC samples at $z=0.38$ and
$0.61$. For $\zeta_{\ell_{1}\ell_{2}\ell}(r_{1},r_{2})$, the measured values
and the standard deviations are plotted as a function of the scale variable
$r_{1}=r_{2}=r$ to simplify the figure.

From this figure, it can be seen that the mean values of $\xi_{\ell}$ and
$\zeta_{\ell_{1}\ell_{2}\ell}$ do not differ much for the different redshifts,
i.e., $z=0.38$ and $0.61$ (compare the magenta and blue lines). One may expect
the amplitudes of $\xi_{\ell}$ and $\zeta_{\ell_{1}\ell_{2}\ell}$ to be larger
at lower redshifts because the tree-level solutions (2.22) of $\xi_{\ell}$ and
$\zeta_{\ell_{1}\ell_{2}\ell}$ are proportional to $D^{2}$ and $D^{4}$,
respectively, with $D$ being the linear growth function. However, this is not
the case in Figure 8. There are two possible reasons for this. The first is
the bias effect. For halos with similar mass, the lower the redshift, the
smaller the value of the linear bias $b_{1}$ tends to be. Therefore,
$b_{1}D(z)$ is less time-dependent and does not show significant differences
at the different redshifts, especially for the monopole components of
$\xi_{\ell}$ and $\zeta_{\ell_{1}\ell_{2}\ell}$. A similar effect to the
linear bias is likely to occur for the non-linear bias included in the 3PCF.
Second, the product of the linear growth rate $f$ and the linear growth
function $D$ is also a less time-dependent function. Therefore, the redshift
dependence of $\xi_{\ell}$ and $\zeta_{\ell_{1}\ell_{2}\ell}$ is not
pronounced, even for the quadrupole component.

On the other hand, the standard deviations of $\xi_{\ell}$ and
$\zeta_{\ell_{1}\ell_{2}\ell}$ are significantly different for the different
redshifts. In general, the so-called Gaussian terms in the covariance matrix
depend only on the two-point statistic, while higher-order statistics such as
the three-point statistic appear in the non-Gaussian terms. It is also known
that the covariance matrix is inversely proportional to the survey volume, and
that the higher the number density of observed galaxies, the smaller the
covariance matrix. Therefore, the fact that the $\xi_{\ell}$ and
$\zeta_{\ell_{1}\ell_{2}\ell}$ signals measured from the Patchy mock do not
differ significantly at the different redshifts suggests that the redshift
dependence in the standard deviation may be due to the survey volume and
galaxy number density.

In Figure 8, the standard deviation at $z=0.61$ (blue) multiplied by
$\sqrt{V_{z=0.61}/V_{z=0.38}}$ is plotted as magenta dashed lines, with the
survey volumes at $z=0.38$ and $z=0.61$ denoted as $V_{z=0.38}$ and
$V_{z=0.61}$, respectively. In the case of the 2PCF, the magenta dashed line
is similar to the result at $z=0.38$ (magenta), indicating that the difference
in the standard deviation of the 2PCF due to differences in redshift can be
explained mainly by differences in the survey volume. However, this is not the
case for the 3PCF, where the standard deviation of the 3PCF at $z=0.38$ is
smaller than the magenta dashed line. This fact suggests that the effect of
the galaxy number density on the covariance matrix is more significant for the
3PCF than for the 2PCF. In other words, as can be seen from Table 1, the
sample at $z=0.38$ has a higher galaxy number density than the sample at
$z=0.61$, even though the survey volume is smaller. Therefore, the standard
deviation at $z=0.38$ is smaller than the standard deviation at $z=0.61$
normalised to the survey volume at $z=0.38$. This result is consistent with
the finding of Sugiyama et al. (2020) that the galaxy number density plays an
essential role in the covariance matrix of the bispectrum, even on large
scales.

###  5.4 Cumulative signal-to-noise ratio

Figure 9: Cumulative signal-to-noise ratios $({\rm S/N}s)$ for the multipole
components of the 2PCF and 3PCF, where both signal and noise (covariance
matrix) are computed from the $2048$ Patchy mock catalogues. The maximum scale
used for the ${\rm S/N}$ calculation is fixed at $r_{\rm
max}=150\,h^{-1}\,{\rm Mpc}$ and the ${\rm S/N}s$ are plotted as a function of
the minimum scale $r_{\rm min}$. The blue and magenta solid lines show the
results for the samples at $z=0.61$ and $z=0.38$, respectively. The magenta
dashed lines are the ${\rm S/N}$ values in the sample at $z=0.61$ multiplied
by $\sqrt{V_{z=0.38}/V_{z=0.61}}$.

The covariance matrix is a two-dimensional quantity in the 2PCF case and a
four-dimensional quantity in the 3PCF case. Therefore, a useful way
compressing and quantifying this multi-dimensional information in the
covariance matrix is to estimate the cumulative signal-to-noise $({\rm S/N})$
ratios, given by

$\displaystyle\left(\rm{\frac{S}{N}}\right)=\left(\overline{\textbf{{d}}}^{T}\cdot\textbf{{C}}_{\rm
Hartlap}^{-1}\cdot\overline{\textbf{{d}}}\right)^{1/2}.$ (5.7)

We calculate the cumulative ${\rm S/N}$ for each multipole component of the
2PCF and 3PCF: i.e., $\overline{d}=\xi_{0}$, $\xi_{2}$, $\zeta_{000}$,
$\zeta_{110}$, $\zeta_{202}$, or $\zeta_{112}$. We also fix the maximum scale
$r_{\rm max}=150\,h^{-1}\,{\rm Mpc}$, vary the minimum scale $r_{\rm min}$
from $150\,h^{-1}\,{\rm Mpc}$ to $30\,h^{-1}\,{\rm Mpc}$, and calculate the
${\rm S/N}$ as a function of $r_{\rm min}$. In Figure 9, we plot the ${\rm
S/N}$ for the four BOSS samples, NGC and SGC at $z=0.38$ and $0.61$. Note that
we do not consider cross-covariance matrices between different multipole
components, e.g., between $\xi_{0}$ and $\zeta_{000}$. How the information in
the covariance matrix, including all cross-covariance matrices, ultimately
propagates to the errors in the cosmological parameters of interest will be
discussed through the Fisher analysis in Section 7.

The top two panels of Figure 9 show the ${\rm S/N}$ of $\xi_{0}$ and
$\xi_{2}$. In all cases shown in the panels, the ${\rm S/N}$ at $z=0.61$ (blue
line) is larger than the ${\rm S/N}$ at $z=0.38$ (magenta line). The
difference is because the ${\rm S/N}$ of the 2PCF is proportional to the
square root of the survey volume $V$, and the survey volume at $z=0.61$,
denoted $V_{\rm z=0.61}$, is larger than the survey volume at $z=0.38$,
$V_{\rm z=0.38}$. Therefore, multiplying the ${\rm S/N}$ at $z=0.61$ by
$\sqrt{V_{\rm z=0.38}/V_{\rm z=0.61}}$ approximately reproduces the ${\rm
S/N}$ at $z=0.38$ (see magenta dashed lines). This result is consistent with
the findings in the signal and standard deviation of the 2PCF in Figure 8.

The middle and bottom results are for $\zeta_{000}$, $\zeta_{110}$,
$\zeta_{202}$, and $\zeta_{112}$. These results show that, in contrast to the
2PCF case, the ${\rm S/N}$ at $z=0.38$ is comparable to the ${\rm S/N}$ at
$z=0.61$. The difference in the ${\rm S/N}$ at $z=0.38$ and $z=0.61$ in the
3PCF case cannot be explained by the difference in the survey volumes (see
magenta dashed lines). This behaviour of the ${\rm S/N}$ of the 3PCF can be
explained by the finding shown in Figure 8 that the galaxy number density
strongly influences the standard deviation of the 3PCF. In particular, in the
present case, the effect of the galaxy number density is more pronounced when
considering correlations between different scales, resulting in the ${\rm
S/N}$ at $z=0.38$ that is comparable to the ${\rm S/N}$ at $z=0.61$. This
result shows that a higher galaxy number density is as important for obtaining
cosmological information from the 3PCF as increasing the survey volume.

##  6 Analysis settings

###  6.1 Likelihoods

We assume that the likelihood of the data compared to the model predictions
follows a multivariate Gaussian distribution:

$\displaystyle\ln{\cal
L}(\textbf{{d}}|\boldsymbol{\theta})=-\frac{1}{2}\left[\textbf{{d}}-\textbf{{t}}(\boldsymbol{\theta})\right]\textbf{{C}}^{-1}_{\rm
Hartlap}\left[\textbf{{d}}-\textbf{{t}}(\boldsymbol{\theta})\right]^{T},$
(6.1)

where d is the data vector, t is the model prediction of the data vector given
the model parameters $\boldsymbol{\theta}$, and $\textbf{{C}}^{-1}_{\rm
Hartlap}$ is the inverse of the covariance matrix after correction by the
Hartlap factor (5.2). We can then obtain the posterior distribution of the
model parameters given the data by performing Bayesian inference:

$\displaystyle{\cal P}(\boldsymbol{\theta}|\textbf{{d}})\propto{\cal
L}(\textbf{{d}}|\boldsymbol{\theta})\Pi(\boldsymbol{\theta})$ (6.2)

where ${\cal P}(\boldsymbol{\theta}|\textbf{{d}})$ is the posterior
distribution of $\boldsymbol{\theta}$ given the data vector, d, and
$\Pi(\boldsymbol{\theta})$ is the prior distribution.

We assume that the four BOSS galaxy samples (Table 1) are far enough apart
that they each have independent cosmological information. Then, when
constraining the model parameters common to each galaxy sample, we add up the
likelihood functions of each galaxy sample. For example, when using all four
galaxy samples, the total likelihood function is given by

$\displaystyle\ln{\cal L}_{\rm total}$ $\displaystyle=$ $\displaystyle\ln{\cal
L}_{{\rm NGC\,at}\,z=0.38}+\ln{\cal L}_{{\rm SGC\,at}\,z=0.38}$ (6.3)
$\displaystyle+$ $\displaystyle\ln{\cal L}_{{\rm NGC\,at}\,z=0.61}+\ln{\cal
L}_{{\rm SGC\,at}\,z=0.61}.$

###  6.2 Multipoles used, scale range, and number of bins

To repeat what was explained in Section 4.2, the scale range used for
parameter estimation in Section 9 is $80\,h^{-1}\,{\rm Mpc}\leq r\leq
150\,h^{-1}\,{\rm Mpc}$, and we choose $\Delta r=5\,h^{-1}\,{\rm Mpc}$ and
$10\,h^{-1}\,{\rm Mpc}$ for the 2PCF and 3PCF bin widths, respectively. The
2PCF and 3PCF multipoles used are the monopole and quadrupole 2PCFs ($\xi_{0}$
and $\xi_{2}$), the two monopole 3PCFs ($\zeta_{000}$ and $\zeta_{202}$), and
the two quadrupole 3PCFs ($\zeta_{202}$ and $\zeta_{112}$). 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. The reason why
the bin width for the 3PCF is wider than for the 2PCF is to reduce the number
of data bins and to conservatively estimate the inverse covariance matrix for
the 2PCF and 3PCF. The total number of data bins is then $202$, which is small
enough compared to the $2048$ Patchy mock simulations (Section 5).

###  6.3 Parameter setting

The parameters we constrain are as follows:

$\displaystyle\boldsymbol{\theta}=\boldsymbol{\theta}_{\rm
phys}+\boldsymbol{\theta}_{\rm bias},$ (6.4)

where

$\displaystyle\boldsymbol{\theta}_{\rm bias}=\\{(b_{1}\sigma_{8}),(F_{\rm
g}\sigma_{8}),(F_{\rm t}\sigma_{8})\\},$ (6.5)

and

$\displaystyle\boldsymbol{\theta}_{\rm
phys}=\begin{cases}\\{f\sigma_{8},\sigma_{8}\\},\quad\mbox{GR};\\\
\\{\sigma_{8},\xi_{f},\xi_{\rm t}\\},\quad\mbox{Horndeski};\\\ \\{F_{\rm
s}\sigma_{8},\xi_{f},\xi_{\rm s},\xi_{\rm t}\\},\quad\mbox{DHOST}.\\\
\end{cases}$ (6.6)

The parameters $F_{\rm g,s,t}$ and $\xi_{f,\rm s,t}$ are given in Eq. (3.15)
and (3.27), respectively; the reason $F_{\rm s}$ does not appear in GR and
Horndeski theorfies is that $F_{\rm s}=1$ in those theories. We assume that
the bias parameters $\boldsymbol{\theta}_{\rm bias}$ take different values in
all four BOSS samples. $f\sigma_{8}$, $\sigma_{8}$, and $F_{\rm s}\sigma_{8}$
have common values in NGC and SGC. $\xi_{f,{\rm s},{\rm t}}$ are common to all
four BOSS samples. For the 2PCF analysis, we only consider $b_{1}\sigma_{8}$
and $f\sigma_{8}$. For example, if all four BOSS samples are used to constrain
DHOST theories, the total number of parameters is $17$. Once again, note that
the AP parameters are not varied in this analysis.

###  6.4 $M_{1}$ and $M_{2}$ factors

| $M_{1}$ | $M_{2}$ |  $\\#$ of params.
---|---|---|---
2PCF only ($z_{\rm eff}=0.38$) | $1.006$ | $1.013$ | $3$
2PCF only ($z_{\rm eff}=0.61$) | $1.006$ | $1.013$ | $3$
GR ($z_{\rm eff}=0.38$) | $1.049$ | $1.105$ | $8$
GR ($z_{\rm eff}=0.61$) | $1.049$ | $1.105$ | $8$
Horndeski ($z_{\rm eff}=0.38$) | $1.048$ | $1.104$ | $9$
Horndeski ($z_{\rm eff}=0.61$) | $1.048$ | $1.104$ | $9$
Horndeski ($z_{\rm eff}=0.38,\,0.61$) | $1.044$ | $1.100$ | $16$
DHOST ($z_{\rm eff}=0.38$) | $1.048$ | $1.104$ | $10$
DHOST ($z_{\rm eff}=0.61$) | $1.048$ | $1.104$ | $10$
DHOST ($z_{\rm eff}=0.38,\,0.61$) | $1.044$ | $1.100$ | $17$
Table 3:  A summary of the $M_{1}$ (5.3) and $M_{2}$ (5.5) factor values used
in our analysis. These values are calculated from the number of the Patchy
mock simulations, $2048$ (Section 4.1), the number of data bins, $30$ for the
2PCF only and $202$ for the 2+3PCF (Section 6.2), and the number of parameters
summarised in the rightmost column (Section 6.3).

As mentioned in Section 5, the number of the Patchy mock simulations used to
calculate the covariance matrices for the 2PCF and 3PCF is finite, so the
inverse of the covariance matrix must be multiplied by the Hartlap factor and
the final parameter error by $M_{1}$.

The $M_{1}$ factor (5.3) is derived assuming that all parameters are
constrained from a single data set. However, when constraining the common
parameters $\xi_{f}$, $\xi_{\rm s}$, and $\xi_{\rm t}$ (3.27) from the four
independent BOSS samples, as in the present analysis, the $M_{1}$ factor is
expected to take a different form, but we do not know the correct correction
factor corresponding to the $M_{1}$ factor in such as case. Therefore, when we
use different galaxy samples simultaneously, we first count up all common and
non-common parameters in the galaxy samples. Then, we calculate the $M_{1}$
factor using the number of data bins computed from a single galaxy sample and
the number of the Patchy mocks, $2048$, corresponding to that galaxy sample,
and multiply it by the final parameter error. Specifically, the multipole
components of the 2PCF and 3PCF measured from a single galaxy sample are
$\xi_{0}$, $\xi_{2}$, $\zeta_{000}$, $\zeta_{110}$, $\zeta_{202}$, and
$\zeta_{112}$, for a total data bin number of $202$ (Section 6.2). The number
of parameters depends on the type of analysis; for example, we need $17$
parameters to test DHOST theories using all four galaxy samples (Section 6.3).

Table 3 summarises the values of the $M_{1}$ and $M_{2}$ factors calculated in
our analysis, leading to $M_{2}\sim 1.1$ for all the 2+3PCF joint analyses.
Thus, even without considering the Hartlap and $M_{1}$ factors in our
analysis, the effect of finite mocks is at most $10\%$. In other words, since
our analysis correctly considers these factors, the error due to the finite
mock effect in the estimated parameter error is guaranteed to be $\lesssim
10\%$.

###  6.5 MCMC

We apply the Metropolis-Hastings (MH) algorithm, an Markov Chain Monte Carlo
(MCMC) method, implemented in the publicly available software package Monte
Python (Audren et al., 2013; Brinckmann & Lesgourgues, 2019) to estimate the
posterior distribution of parameters in a multi-dimensional parameter space.
In doing so, we set the super-update parameter to $20$, as recommended by the
developers. In order to improve the convergence of the posterior distributions
obtained by MCMC, we first perform an MCMC analysis with the number of steps
set to $N_{\rm step}=200\,000$ and calculate the best-fit values and
covariance matrix of the parameters. Then, we add the information of the best-
fit values and covariance matrix and perform an MCMC analysis again with the
same number of steps .

We ensure convergence of each MCMC chain, imposing $R-1\lesssim{\cal
O}(10^{-4})$ where $R$ is the standard Gelman-Rubin criteria (Gelman & Rubin,
1992). Furthermore, the convergence of the results is also checked through the
following method. First, we create eight independent MCMC chains and compute
the mean and standard deviation of the parameters from each chain. Then, from
the eight means and standard deviations, we compute the standard deviation of
the mean and the mean of the standard deviation and check that the ratio of
them is less than about $20\%$ for all the results.

###  6.6 Mock tests

We perform MCMC analyses on $100$ Patchy mock catalogues (Kitaura et al.,
2016) using the same set of cosmological and nuisance parameters as in the
actual BOSS galaxy data analysis. We then verify that our analysis can
correctly return the values of the non-linear parameters predicted by GR for
the Patchy mock catalogues designed under the assumption of a $\Lambda$CDM
model. We also discuss the statistical scatter of the $100$ values of the
parameters to be estimated.

##  7 Fisher analysis

| $(b_{1}\sigma_{8})_{\rm fid}$ | $(f\sigma_{8})_{\rm fid}$ | $(F_{\rm g}\sigma_{8})_{\rm fid}$ | $(F_{\rm s}\sigma_{8})_{\rm fid}$ | $(F_{\rm t}\sigma_{8})_{\rm fid}$ | $(G_{\rm s}\sigma_{8})_{\rm fid}$ | $(G_{\rm t}\sigma_{8})_{\rm fid}$
---|---|---|---|---|---|---|---
| $1.362$ | $0.485$ | $0.552$ | $0.681$ | $0.194$ | $0.681$ | $0.386$
| $\sigma_{\rm fisher}(b_{1}\sigma_{8})$ | $\sigma_{\rm fisher}(f\sigma_{8})$ | $\sigma_{\rm fisher}(F_{\rm g}\sigma_{8})$ | $\sigma_{\rm fisher}(F_{\rm s}\sigma_{8})$ | $\sigma_{\rm fisher}(F_{\rm t}\sigma_{8})$ | $\sigma_{\rm fisher}(G_{\rm s}\sigma_{8})$ | $\sigma_{\rm fisher}(G_{\rm t}\sigma_{8})$
GR
Case 1 | $0.159$ | $0.093$ | $-$ | $-$ | $-$ | $-$ | $-$
Case 2 | $0.154$ | $0.093$ | $0.418$ | $0.472$ | $0.312$ | $-$ | $-$
Case 3 | $0.159$ | $0.092$ | $0.802$ | $1.170$ | $2.387$ | $-$ | $-$
Case 4 | $0.155$ | $0.091$ | $0.420$ | $0.643$ | $0.426$ | $-$ | $-$
Case 5 | $0.151$ | $0.090$ | $0.330$ | $0.450$ | $0.291$ | $-$ | $-$
Case 6 | $0.155$ | $0.089$ | $0.396$ | $0.619$ | $0.411$ | $-$ | $-$
Case 7 | $0.151$ | $0.089$ | $0.315$ | $0.442$ | $0.283$ | $-$ | $-$
Case 8 | $0.153$ | $0.091$ | $0.363$ | $0.492$ | $0.324$ | $-$ | $-$
Horndeski
Case 2 | $0.154$ | $0.093$ | $2.361$ | $0.479$ | $3.156$ | $-$ | $29.25$
Case 3 | $0.159$ | $0.092$ | $0.816$ | $2.433$ | $2.633$ | $-$ | $1.312$
Case 4 | $0.155$ | $0.092$ | $0.430$ | $0.726$ | $0.431$ | $-$ | $1.111$
Case 5 | $0.151$ | $0.091$ | $0.331$ | $0.463$ | $0.295$ | $-$ | $1.044$
Case 6 | $0.155$ | $0.091$ | $0.409$ | $0.715$ | $0.419$ | $-$ | $1.015$
Case 7 | $0.151$ | $0.091$ | $0.316$ | $0.459$ | $0.285$ | $-$ | $0.946$
Case 8 | $0.153$ | $0.091$ | $0.366$ | $0.498$ | $0.344$ | $-$ | $1.383$
DHOST
Case 2 | $0.154$ | $0.093$ | $3.100$ | $3.840$ | $4.774$ | $35.22$ | $38.88$
Case 3 | $0.159$ | $0.093$ | $2.211$ | $3.652$ | $2.633$ | $2.976$ | $1.981$
Case 4 | $0.157$ | $0.092$ | $0.726$ | $0.748$ | $0.446$ | $1.896$ | $1.112$
Case 5 | $0.153$ | $0.092$ | $0.499$ | $0.476$ | $0.303$ | $1.636$ | $1.051$
Case 6 | $0.156$ | $0.091$ | $0.710$ | $0.739$ | $0.431$ | $1.757$ | $1.017$
Case 7 | $0.153$ | $0.091$ | $0.487$ | $0.468$ | $0.293$ | $1.519$ | $0.953$
Case 8 | $0.154$ | $0.092$ | $0.595$ | $0.528$ | $0.344$ | $2.292$ | $1.580$

Table 4:  The standard deviations of the parameters as predicted by the Fisher
analysis, denoted as $\sigma_{\rm fisher}(\boldsymbol{\theta})$, are shown.
These results are for the NGC at $z_{\rm eff}=0.38$. The parameter vectors of
interest, $\boldsymbol{\theta}$, are different for each of the three gravity
theories, GR, Horndeski, and DHOST theories (Eq. 7.3). The classification of
the data vectors used is as shown in Eq. (7.4). The fiducial values of the
parameters are calculated under the assumption of GR and are denoted as
$(\boldsymbol{\theta})_{\rm fid}$. The scale range used for this Fisher
analysis is $80\,h^{-1}\,{\rm Mpc}\leq r\leq 150\,h^{-1}\,{\rm Mpc}$.

| $(E_{f})_{\rm fid}$ | $(E_{\rm s})_{\rm fid}$ | $(E_{\rm t})_{\rm fid}$ |  | $(\xi_{f})_{\rm fid}$ | $(\xi_{\rm s})_{\rm fid}$ | $(\xi_{\rm t})_{\rm fid}$
---|---|---|---|---|---|---|---
| $0.713$ | $1.000$ | $0.992$ |  | $0.545$ | $0.000$ | $0.013$
| $\sigma_{\rm fisher}(E_{f})$ | $\sigma_{\rm fisher}(E_{\rm s})$ | $\sigma_{\rm fisher}(E_{\rm t})$ |  | $\sigma_{\rm fisher}(\xi_{f})$ | $\sigma_{\rm fisher}(\xi_{\rm s})$ | $\sigma_{\rm fisher}(\xi_{\rm t})$
Horndeski
Case 2 | $0.502$ | $-$ | $75.28$ |  | $1.146$ | $-$ | $124.3$
Case 3 | $2.558$ | $-$ | $6.018$ |  | $5.844$ | $-$ | $9.934$
Case 4 | $0.760$ | $-$ | $3.189$ |  | $1.737$ | $-$ | $5.264$
Case 5 | $0.488$ | $-$ | $2.767$ |  | $1.114$ | $-$ | $4.567$
Case 6 | $0.750$ | $-$ | $2.969$ |  | $1.714$ | $-$ | $4.901$
Case 7 | $0.483$ | $-$ | $2.522$ |  | $1.103$ | $-$ | $4.163$
Case 8 | $0.523$ | $-$ | $3.597$ |  | $1.194$ | $-$ | $5.938$
DHOST
Case 2 | $4.018$ | $57.31$ | $103.8$ |  | $9.180$ | $93.85$ | $171.3$
Case 3 | $3.824$ | $6.519$ | $9.837$ |  | $8.737$ | $10.68$ | $16.24$
Case 4 | $0.785$ | $2.835$ | $3.213$ |  | $1.794$ | $4.642$ | $5.304$
Case 5 | $0.503$ | $2.460$ | $2.771$ |  | $1.150$ | $4.029$ | $4.574$
Case 6 | $0.777$ | $2.633$ | $2.998$ |  | $1.776$ | $4.311$ | $4.950$
Case 7 | $0.495$ | $2.260$ | $2.527$ |  | $1.130$ | $3.701$ | $4.172$
Case 8 | $0.558$ | $3.542$ | $3.981$ |  | $1.274$ | $5.801$ | $6.572$

Table 5:  Same as Table 5, but the standard deviations of the parameters
defined in Eqs. (3.23) and (3.27), calculated through variable
transformations, are shown.

Before proceeding to the MCMC analysis using actual galaxies in Sections 8 and
9, in this section, we will understand how the 3PCF contains cosmological
information through the Fisher analysis.

There are several reasons for performing the Fisher analysis before the MCMC
analysis. First, calculating the Fisher matrix in Section 7.1 is less
computationally intensive than performing the MCMC analysis, making it easier
to compare the analysis results in various settings that take too much time in
the MCMC analysis. Taking advantage of this, Section 7.2 examines how the
constraints on the parameters of interest change for various combinations of
the multipole components of the 3PCF; in doing so, we focus only on the NGC
sample at $z=0.38$ as a representative example. Section 7.3 also discusses the
relation among the values of the predicted parameter errors for the four BOSS
samples, NGC and SGC at $z=0.38,\,0.61$, and how the combination of the four
BOSS samples affects the final results. Finally, in Section 9.11, we compare
the results obtained from the above Fisher analysis with those obtained from
the MCMC parameter estimation and check their consistency to confirm the
validity of the final results in this paper.

In Section 7.4, the Fisher analysis also allows us to estimate cosmological
information at scales smaller than the scale range used in the MCMC analysis.
The results are expected to motivate the construction of theoretical models
applicable to smaller scales.

Finally, in Section 7.5, we use the results of the Fisher analysis to
determine the range of a flat prior used when performing the MCMC analysis.

###  7.1 Fisher matrix

From the likelihood function given in Eq. (6.1), we calculate the Fisher
matrix as

$\displaystyle F_{ij}$ $\displaystyle=$
$\displaystyle-\left\langle\frac{\partial}{\partial\theta_{i}}\frac{\partial}{\partial\theta_{j}}\ln{\cal
L}\right\rangle$ (7.1) $\displaystyle=$
$\displaystyle\frac{\partial\,\textbf{{t}}(\boldsymbol{\theta})}{\partial\theta_{i}}\textbf{{C}}_{\rm
Hartlap}^{-1}\frac{\partial\,\textbf{{t}}^{\rm
T}(\boldsymbol{\theta})}{\partial\theta_{j}},$

where we assumed that the covariance matrix C is independent of the
parameters. The indices $i$ and $j$ run over parameters of interest. In the
limit of the Gaussian likelihood surface, the Cramer-Rao inequality shows that
the Fisher matrix provides the minimum standard deviation on parameters,
marginalized over all the other parameters: $\sigma(\theta_{i})\geq\sigma_{\rm
fisher}(\theta_{i})=\left(F^{-1}\right)_{ii}^{1/2}$. We note that we adopt the
inverse covariance matrix, $\textbf{{C}}_{\rm Hartlap}^{-1}$, that is non-
Gaussian estimated from the Patchy mock simulations.

We consider three parameter vectors, depending on the gravity theory of
interest:

$\displaystyle\boldsymbol{\theta}=\\{(b_{1}\sigma_{8}),(f\sigma_{8})\\}+\boldsymbol{\theta}_{\rm
3PCF},$ (7.2)

with

$\displaystyle\boldsymbol{\theta}_{\rm 3PCF}$ $\displaystyle=$
$\displaystyle\begin{cases}\\{(F_{\rm g}\sigma_{8}),(F_{\rm
s}\sigma_{8}),(F_{\rm t}\sigma_{8})\\},\hskip 85.07388pt\mbox{GR};\\\
\\{(F_{\rm g}\sigma_{8}),(F_{\rm s}\sigma_{8}),(F_{\rm t}\sigma_{8}),(G_{\rm
t}\sigma_{8})\\},\hskip 19.91684pt\quad\mbox{Horndeski};\\\ \\{(F_{\rm
g}\sigma_{8}),(F_{\rm s}\sigma_{8}),(F_{\rm t}\sigma_{8}),(G_{\rm
s}\sigma_{8}),(G_{\rm t}\sigma_{8})\\},\hskip 5.69046pt\mbox{DHOST},\\\
\end{cases}$ (7.3)

where $F_{\rm s}\sigma_{8}=\sigma_{8}$ for GR and Horndeski theories. We
obtain the results for $E_{f,\rm s,t}$ and $\xi_{f,\rm s,t}$ using the
variable transformations in Eqs. (3.23) and (3.27). In particular, the results
including $\xi_{f,\rm s,t}$ correspond to the parameter set (6.6) used in the
MCMC analysis performed in Section 9.

The fiducial values of the cosmological parameters needed to compute the
Fisher matrix are the values in the $\Lambda$CDM model presented in Section 1.
In doing so, we assume that the linear bias is $b_{1}=2$, and the values of
the non-linear biases are zero: i.e., $b_{2}=b_{\rm s^{2}}=0$.

###  7.2 Information contained in 3PCF multipoles

For the NGC sample at $z=0.38$, we perform Fisher analyses on the following
eight data vectors consisting of combinations of the 2PCF and 3PCF multipole
components, using the same settings as the MCMC analysis performed in Section
9 to investigate which components and how they affect parameter estimates.

$\displaystyle\mbox{Case 1}\quad\textbf{{d}}=\\{\xi_{0},\xi_{2}\\};$
$\displaystyle\mbox{Case
2}\quad\textbf{{d}}=\\{\xi_{0},\xi_{2},\zeta_{000},\zeta_{110}\\};$
$\displaystyle\mbox{Case
3}\quad\textbf{{d}}=\\{\xi_{0},\xi_{2},\zeta_{202},\zeta_{112}\\};$
$\displaystyle\mbox{Case
4}\quad\textbf{{d}}=\\{\xi_{0},\xi_{2},\zeta_{000},\zeta_{202}\\};$
$\displaystyle\mbox{Case
5}\quad\textbf{{d}}=\\{\xi_{0},\xi_{2},\zeta_{000},\zeta_{202},\zeta_{110}\\};$
$\displaystyle\mbox{Case
6}\quad\textbf{{d}}=\\{\xi_{0},\xi_{2},\zeta_{000},\zeta_{202},\zeta_{112}\\};$
$\displaystyle\mbox{Case
7}\quad\textbf{{d}}=\\{\xi_{0},\xi_{2},\zeta_{000},\zeta_{202},\zeta_{110},\zeta_{112}\\};$
$\displaystyle\mbox{Case
8}\quad\textbf{{d}}=\\{\xi_{0},\xi_{2},\zeta_{110},\zeta_{112}\\}.$ (7.4)

Case $1$ constrains $f\sigma_{8}$ using only the monopole and quadrupole
2PCFs. Cases $2$, $3$, and $4$ add to Case $1$ the two monopole 3PCFs
($\zeta_{000}$ and $\zeta_{110}$), the two quadrupole 3PCFs ($\zeta_{202}$ and
$\zeta_{112}$), and the first terms of the monopole and quadrupole 3PCFs
($\zeta_{000}$ and $\zeta_{202}$), respectively. These three cases will
highlight the importance of simultaneously considering both monopoles and
quadrupoles in the 3PCF. Moreover, Cases $5$, $6$, and $7$ reveal the extent
to which the final results can be improved by adding higher-order multipole
components to Case $4$. Finally, Case $8$ only uses the higher-order
multipoles, $\zeta_{110}$ and $\zeta_{112}$, for the monopole and quadrupole
components.

We summarise the results of the Fisher analysis in Table 5. In Horndeski and
DHOST theories, the case $2$ results show that using only the monopole 3PCFs
very weakly constrains the non-linear velocity parameters $G_{\rm
s}\sigma_{8}$ and $G_{\rm t}\sigma_{8}$. On the other hand, in Case $3$, using
only the quadrupole 3PCFs, we can mildly constrain the non-linear coefficients
of both the density field and the velocity field. The reason is that the
density and velocity fluctuations contribute to the quadrupole 3PCFs to the
same extent (Figure 2). Moreover, Cases $4$, $5$, $6$, $7$, and $8$, using
both the monopole and quadrupole components, can constrain the non-linear
coefficients more strongly than Cases $2$ and $3$. In particular, for the
$G_{\rm s}\sigma_{8}$ and $G_{\rm t}\sigma_{8}$ constraints in DHOST theories,
Case $7$ is $\sim 20$ and $\sim 40$ times better than Case $2$, respectively:

$\displaystyle\sigma_{\rm fisher}(G_{\rm s}\sigma_{8})$ $\displaystyle=$
$\displaystyle 35.22\quad\mbox{for Case 2},$ $\displaystyle\sigma_{\rm
fisher}(G_{\rm t}\sigma_{8})$ $\displaystyle=$ $\displaystyle
38.88\quad\mbox{for Case 2},$ $\displaystyle\sigma_{\rm fisher}(G_{\rm
s}\sigma_{8})$ $\displaystyle=$ $\displaystyle 1.519\quad\mbox{for Case 7},$
$\displaystyle\sigma_{\rm fisher}(G_{\rm t}\sigma_{8})$ $\displaystyle=$
$\displaystyle 0.953\quad\mbox{for Case 7}.$ (7.5)

These results support the argument of this paper that we should use both
monopole and quadrupole 3PCFs to study the non-linearity of the velocity
field.

Case $7$, which uses all components of $\zeta_{000}$, $\zeta_{110}$,
$\zeta_{202}$, and $\zeta_{112}$, provides the best constraints on $F_{\rm
s}\sigma_{8}$, $G_{\rm s}\sigma_{8}$, and $G_{\rm t}\sigma_{8}$, as expected.
Therefore, we can conclude that all these multipole components should be used
in the MCMC analysis in Section 9.

Case $7$ yields results that are about $10\%$ better than Case $5$, which uses
$\zeta_{000}$, $\zeta_{110}$, and $\zeta_{202}$. This result indicates that
while $\zeta_{202}$ contains the main cosmological information, $\zeta_{112}$
contains other information in addition to $\zeta_{202}$. Existing studies
using Scoccimarro et al. (1999)’ decomposition method of the bispectrum tend
to ignore the $M\neq 0$ mode of the quadrupole component as not containing
much cosmological information (e.g., Gagrani & Samushia, 2017; Rizzo et al.,
2023; D’Amico et al., 2022b). However, our results show the importance of the
$M\neq 0$ modes because $\zeta_{202}$ contains only the $M=0$ mode, while
$\zeta_{112}$ further contains the $M\neq 0$ modes in addition to the $M=0$
mode (see also Section 3.1).

By comparing the results of Case $4$, consisting of $\zeta_{000}$ and
$\zeta_{202}$, with those of Case $8$, consisting of $\zeta_{110}$ and
$\zeta_{112}$, we can find another viewpoint on the importance of higher-order
multipole components. For example, the $(F_{\rm s}\sigma_{8})$ constraint is
better in Case $8$, and the $(G_{\rm s}\sigma_{8})$ and $(G_{\rm
t}\sigma_{8})$ constraints are better in Case $4$. Also, the $(G_{\rm
s,t}\sigma_{8})$ result in Case $4$ is only about $30\%$ better than Case $8$.
Thus, although $\zeta_{202}$ is more informative than $\zeta_{112}$, we
interpret the information on both sides as overlapping to some extent.

We further calculate $\sigma_{\rm fisher}(\boldsymbol{\theta})$ for
$\theta=E_{f}$, $E_{\rm s}$, $E_{\rm t}$, $\xi_{f}$, $\xi_{\rm s}$, and
$\xi_{\rm t}$ through the variable transformations in Eqs. (3.23) and (3.27),
summarising the results in Table 5. We find that both the monopole and
quadrupole components of the 3PCF are needed to constrain $E_{\rm s}$, $E_{\rm
t}$, $\xi_{\rm s}$, and $\xi_{\rm t}$ better. In Case $7$, the standard
deviations of $E_{\rm s}$ and $E_{\rm t}$ are more than twice larger than the
fiducial values of $E_{\rm s}$ and $E_{\rm t}$, i.e., $\sigma_{\rm
fisher}(E_{\rm s,t})/(E_{\rm s,t})_{\rm fid}>2$, indicating that it is
impossible to detect the $E_{\rm s}$ and $E_{\rm t}$ signals in the BOSS data.
We can also confirm that for each of the $\xi_{\rm s}$ and $\xi_{\rm t}$
constraints in DHOST theories, the results of Case $7$ are $\sim 25$ and $\sim
40$ times stronger than those of Case $2$, respectively:

$\displaystyle\sigma_{\rm fisher}(\xi_{\rm s})$ $\displaystyle=$
$\displaystyle 93.85\quad\mbox{for Case 2},$ $\displaystyle\sigma_{\rm
fisher}(\xi_{\rm t})$ $\displaystyle=$ $\displaystyle 171.3\quad\mbox{for Case
2},$ $\displaystyle\sigma_{\rm fisher}(\xi_{\rm s})$ $\displaystyle=$
$\displaystyle 3.701\quad\mbox{for Case 7},$ $\displaystyle\sigma_{\rm
fisher}(\xi_{\rm t})$ $\displaystyle=$ $\displaystyle 4.172\quad\mbox{for Case
7}.$ (7.6)

In GR, adding any multipole component of the 3PCF can only improve the
$b_{1}\sigma_{8}$ and $f\sigma_{8}$ constraints by a few per cent. This result
is consistent with the MCMC analysis of Sugiyama et al. (2021) on the Patchy
mock catalogues. Furthermore, the 3PCF-specific information, $\sigma_{8}$, is
also uninformative compared to $f\sigma_{8}$. Specifically, $f\sigma_{8}$ can
be determined with a precision of $\sim 20\%$, while $\sigma_{8}$ can only
reach a precision of $\sim 60\%$. These results are for large scales ($r\geq
80\,h^{-1}\,{\rm Mpc}$); what happens when even smaller scales are used will
be discussed in Section 7.4.

###  7.3 Fisher forecasts with all four BOSS samples

DHOST
---
$\sigma_{\rm fisher}(\theta)/(\theta)_{\rm fid}$ | $(b_{1}\sigma_{8})$ | $(f\sigma_{8})$ | $(F_{\rm g}\sigma_{8})$ | $(F_{\rm s}\sigma_{8})$ | $(F_{\rm t}\sigma_{8})$ | $(G_{\rm s}\sigma_{8})$ | $(G_{\rm t}\sigma_{8})$
NGC at $z=0.38$ | $0.112$ | $0.188$ | $0.883$ | $0.687$ | $1.511$ | $2.231$ | $2.469$
NGC at $z=0.61$ | $0.110$ | $0.181$ | $1.201$ | $0.880$ | $2.077$ | $2.789$ | $3.265$
SGC at $z=0.38$ | $0.188$ | $0.311$ | $1.446$ | $1.142$ | $2.609$ | $3.642$ | $4.123$
SGC at $z=0.61$ | $0.188$ | $0.304$ | $2.055$ | $1.444$ | $3.609$ | $4.721$ | $5.460$

Table 6:  The standard deviations computed by the Fisher analysis of Case $7$
in DHOST theories divided by the fiducial values of the parameters,
$\sigma_{\rm fisher}(\boldsymbol{\theta})/(\boldsymbol{\theta})_{\rm fid}$,
are shown, where
$\boldsymbol{\theta}=\\{(b_{1}\sigma_{8}),(f\sigma_{8}),(F_{\rm
g}\sigma_{8}),(F_{\rm s}\sigma_{8}),(F_{\rm t}\sigma_{8}),(G_{\rm
s}\sigma_{8}),(G_{\rm t}\sigma_{8})\\}$. These results are for the NGC and SGC
at $z=0.38,\,0.61$.

DHOST
---
| $\sigma_{\rm fisher}(E_{f})/(E_{f})_{\rm fid}$ | $\sigma_{\rm fisher}(E_{\rm s})/(E_{\rm s})_{\rm fid}$ | $\sigma_{\rm fisher}(E_{\rm t})/(E_{\rm t})_{\rm fid}$ |  | $\sigma_{\rm fisher}(\xi_{f})$ | $\sigma_{\rm fisher}(\xi_{\rm s})$ | $\sigma_{\rm fisher}(\xi_{\rm t})$
NGC at $z=0.38$ | $0.626$ | $2.260$ | $2.541$ |  | $1.130$ | $3.701$ | $4.172$
NGC at $z=0.61$ | $0.890$ | $2.931$ | $3.384$ |  | $2.084$ | $6.890$ | $7.955$
SGC at $z=0.38$ | $1.035$ | $3.653$ | $4.234$ |  | $1.867$ | $5.983$ | $6.950$
SGC at $z=0.61$ | $1.469$ | $4.878$ | $5.677$ |  | $3.439$ | $11.468$ | $13.345$

Table 7:  Same as Figure 7, but $\sigma_{\rm
fisher}(\boldsymbol{\theta})/(\boldsymbol{\theta})_{\rm fid}$ for
$\boldsymbol{\theta}=\\{E_{f},E_{\rm s},E_{\rm t}\\}$ and $\sigma_{\rm
fisher}(\boldsymbol{\theta})$ for $\boldsymbol{\theta}=\\{\xi_{f},\xi_{\rm
s},\xi_{\rm t}\\}$.

In this subsection, we repeat the analysis of Case $7$ in DHOST theories,
performed in Section 7.2, for the other three BOSS samples, NGC at $z=0.61$
and SGC at $z=0.38,\,0.61$, and summarise the results in Tables 7 and 7.

Table 7 shows that the results for $(b_{1}\sigma_{8})$ and $(f\sigma_{8})$,
which are mainly determined by the 2PCF, are slightly better for the sample at
$z=0.61$ than for the sample at $z=0.38$ for both NGC and SGC. On the other
hand, for the 3PCF-specific parameters, $(F_{\rm g,s,t}\sigma_{8})$ and
$(G_{\rm s,t}\sigma_{8})$, the error is smaller for the $z=0.38$ sample than
for the $z=0.61$ sample. This result reflects the different characteristics of
the cumulative ${\rm S/N}$ between the 2PCF and the 3PCF, as discussed in
Section 5.4. In other words, it suggests that higher number densities are more
favourable than larger survey volumes for constraining the non-linear
parameters, $(F_{\rm g,s,t}\sigma_{8})$ and $(G_{\rm s,t}\sigma_{8})$, using
3PCF measurements.

Table 7 summarises the results of $E_{f,{\rm s},{\rm t}}$ and $\xi_{f,{\rm
s},{\rm t}}$. As expected, the $z=0.38$ sample gives a smaller error than the
$z=0.61$ sample for both $E_{f,{\rm s},{\rm t}}$ and $\xi_{f,{\rm s},{\rm
t}}$. However, in the case of $\xi_{f,{\rm s},{\rm t}}$, the error at $z=0.38$
is almost twice as small as that at $z=0.61$, which is extremely favourable
for the $z=0.38$ sample. For example, the $\xi_{s}$ results for the NGC
samples are

$\displaystyle\sigma_{\rm fisher}(\xi_{\rm s})$ $\displaystyle=$
$\displaystyle 3.701\quad\mbox{for NGC at $z=0.38$},$
$\displaystyle\sigma_{\rm fisher}(\xi_{\rm s})$ $\displaystyle=$
$\displaystyle 6.890\quad\mbox{for NGC at $z=0.61$}.$ (7.7)

This is because we parameterise the time evolution of $E_{f,{\rm s},{\rm t}}$
as $E_{f,{\rm s},{\rm t}}=\Omega_{\rm m}^{\xi_{f,{\rm s},{\rm t}}}$. That is,
because $d\xi_{f,{\rm s},{\rm t}}=d\ln E_{f,{\rm s},{\rm t}}/(\ln\Omega_{\rm
m})$, the errors in $\xi_{f,{\rm s},{\rm t}}$ are smaller for lower redshifts
with smaller values of $\Omega_{\rm m}$. Specifically, in the LCDM model
introduced in Section 1, $\Omega_{\rm m}(z=0.38)=0.54$ and $\Omega_{\rm
m}(z=0.61)=0.65$, so $1/\ln\Omega_{\rm m}(z=0.38)=1.62$ and $1/\ln\Omega_{\rm
m}(z=0.61)=2.32$. Even if $\sigma_{\rm fisher}(E_{f,\rm s,t})/(E_{f,\rm
s,t})_{\rm fid}$ has the same value at the two redshifts of $z=0.38$ and
$z=0.61$, the value of $\sigma_{\rm fisher}(\xi_{f,\rm s,t})$ at $z=0.38$ is
$2.32/1.62=1.42$ times smaller than at $z=0.61$.

###  7.4 Fisher forecasts using smaller scales

So far, we have performed the Fisher analysis in the same setting as the MCMC
analysis that will be performed in Section 9. There, we have dealt with the
behaviour of only large scales, $80\,h^{-1}\,{\rm Mpc}\leq r\leq
150\,h^{-1}\,{\rm Mpc}$. However, seeing how the parameter constraints improve
when the minimum scale used, $r_{\rm min}$, is varied should be an excellent
motivation for the future development of theoretical models.

Figure 10 plots $\sigma_{\rm fisher}(\theta)/(\theta)_{\rm fid}$ as a function
of $r_{\rm min}$ for the three gravity theories, GR, Horndeski, and DHOST, at
two redshifts of $z=0.38$ (magenta lines) and $z=0.61$ (blue lines). The
multipole components of the 2PCF and 3PCF used here are Case $7$ (Eq. (7.4)).
First, even on the smaller scale, adding the 3PCF hardly improves the
$f\sigma_{8}$ constraint compared to the case where only the 2PCF is used
(compare solid and dashed lines in the top left panel of Figure 10.). On the
other hand, at $r_{\rm min}=30\,h^{-1}\,{\rm Mpc}$, the $\sigma_{8}$
constraint reaches a precision of $\sim 10\%$, from which useful cosmological
information may be extracted: e.g., $f=(f\sigma_{8})/\sigma_{8}$ can be
determined with a precision of $\sim 10\%$.

In addition, the non-linear velocity parameters, $G_{\rm s}\sigma_{8}$ and
$G_{\rm t}\sigma_{8}$, can be determined with $30\mathchar 45\relax 50\%$
precision at $r_{\rm min}=30\,h^{-1}\,{\rm Mpc}$. Thus, future galaxy surveys
with even larger volumes than the BOSS survey, such as DESI, Euclid, and PFS,
may detect such non-linear coefficients of the velocity field.

Note that we obtained the Fisher analysis results using the IR-resumed tree-
level solutions of the 2PCF and 3PCF given in Eqs. (3.9) and (3.12). Although
these models accurately describe the non-linear damping behaviour of BAO on
large scales, they cannot predict the 2PCF and 3PCF on small scales with high
accuracy. Therefore, to apply these models to smaller scales, it is necessary
to account for non-linear effects, called loop correction terms. We leave to
investigate how the results change when such a loop correction is added for
future research.

Figure 10: The standard deviations computed by the Fisher analysis divided by
the fiducial values of the parameter, $\sigma_{\rm
fisher}(\theta)/(\theta)_{\rm fid}$, are shown as a function of the minimum
scale used, $r_{\rm min}$. These results are for the NGC at $z_{\rm
eff}=0.38,\,0.61$. The solid lines show the results using the multipole
components of the 2PCF and 3PCF given in Case $7$ (7.4), and the dashed lines
are for the 2PCF-only analysis, Case $1$. The points at $r_{\rm
min}=80\,h^{-1}\,{\rm Mpc}$ in the right panels correspond to the results in
Table 7.

###  7.5 Flat priors

Prior range
---
$(b_{1}\sigma_{8})_{\rm NGC,\,z=0.38}$ | $[0.60,2.13]$
$(b_{1}\sigma_{8})_{\rm SGC,\,z=0.38}$ | $[0.08,2.64]$
$(b_{1}\sigma_{8})_{\rm NGC,\,z=0.61}$ | $[0.54,1.88]$
$(b_{1}\sigma_{8})_{\rm SGC,\,z=0.61}$ | $[0.07,2.36]$
$(f\sigma_{8})_{\rm NGC,\,z=0.38}$ | $[0.03,0.94]$
$(f\sigma_{8})_{\rm SGC,\,z=0.38}$ | $[0.00,1.24]$
$(f\sigma_{8})_{\rm NGC,\,z=0.61}$ | $[0.05,0.91]$
$(f\sigma_{8})_{\rm SGC,\,z=0.61}$ | $[0.00,1.21]$
$(F_{\rm g}\sigma_{8})_{\rm NGC,\,z=0.38}$ | $[0.00,2.99]$
$(F_{\rm g}\sigma_{8})_{\rm SGC,\,z=0.38}$ | $[0.00,4.54]$
$(F_{\rm g}\sigma_{8})_{\rm NGC,\,z=0.61}$ | $[0.00,3.44]$
$(F_{\rm g}\sigma_{8})_{\rm SGC,\,z=0.61}$ | $[0.00,5.54]$
$(F_{\rm s}\sigma_{8})_{\rm NGC,\,z=0.38}$ | $[0.00,3.02]$
$(F_{\rm s}\sigma_{8})_{\rm SGC,\,z=0.38}$ | $[0.00,4.57]$
$(F_{\rm s}\sigma_{8})_{\rm NGC,\,z=0.61}$ | $[0.00,3.27]$
$(F_{\rm s}\sigma_{8})_{\rm SGC,\,z=0.61}$ | $[0.00,4.99]$
$(F_{\rm t}\sigma_{8})_{\rm NGC,\,z=0.38}$ | $[0.00,1.66]$
$(F_{\rm t}\sigma_{8})_{\rm SGC,\,z=0.38}$ | $[0.00,2.72]$
$(F_{\rm t}\sigma_{8})_{\rm NGC,\,z=0.61}$ | $[0.00,1.97]$
$(F_{\rm t}\sigma_{8})_{\rm SGC,\,z=0.61}$ | $[0.00,3.29]$
$(\xi_{f})_{\rm NGC,\,z=0.38}$ | $[-5.11,6.20]$
$(\xi_{f})_{\rm SGC,\,z=0.38}$ | $[-8.79,9.88]$
$(\xi_{f})_{\rm NGC,\,z=0.61}$ | $[-9.88,10.97]$
$(\xi_{f})_{\rm SGC,\,z=0.61}$ | $[-16.65,17.74]$
$(\xi_{\rm s})_{\rm NGC,\,z=0.38}$ | $[-18.51,11.10]$
$(\xi_{\rm s})_{\rm SGC,\,z=0.38}$ | $[-29.91,17.95]$
$(\xi_{\rm s})_{\rm NGC,\,z=0.61}$ | $[-34.45,20.67]$
$(\xi_{\rm s})_{\rm SGC,\,z=0.61}$ | $[-57.34,34.40]$
$(\xi_{\rm t})_{\rm NGC,\,z=0.38}$ | $[-20.84,12.53]$
$(\xi_{\rm t})_{\rm SGC,\,z=0.38}$ | $[-34.74,20.86]$
$(\xi_{\rm t})_{\rm NGC,\,z=0.61}$ | $[-39.76,23.88]$
$(\xi_{\rm t})_{\rm SGC,\,z=0.61}$ | $[-66.71,40.05]$
Table 8: The flat priors for the parameters that we employ in our MCMC
analysis are shown. The results are calculated from Case $7$ in Eq. (7.4)
assuming DHOST theories.

As shown by the results of the Fisher analysis in Section 7, the constraints
on the non-linear parameters $\xi_{f,\rm s,t}$ constrained by the 3PCF
measured from BOSS are weak. Therefore, we need to set appropriate priors to
efficiently perform the MCMC analysis.

We use the Fisher analysis results of Case $7$ in DHOST theories for the four
BOSS galaxy samples, performed in Section 7.3. Then, we adopt a flat prior of
$\theta_{\rm fid}\pm 5\sigma_{\rm fisher}(\theta)$ as the base setting for all
parameters. If using several samples to constrain common parameters, we adopt
a narrower range of priors for those samples. For example, at $z_{\rm
eff}=0.38$, when constraining $f\sigma_{8}$ using both NGC and SGC samples, we
adopt the prior computed in NGC. After this basic setting, we set a stronger
prior based on further physical considerations below.

The linear bias $b_{1}$, the linear growth rate $f$, and $\sigma_{8}$ are
always positive by definition: i.e., $b_{1}\sigma_{8}\geq 0$ and
$f\sigma_{8}\geq 0$.

In the case of GR, the non-linear parameters to be constrained are $F_{\rm
g}\sigma_{8}$, $\sigma_{8}$, and $F_{\rm t}\sigma_{8}$. The non-linear local
bias parameter $(1/2)(b_{2}/b_{1})$ appearing in $F_{\rm g}$ is calculated to
be $-0.02$ for $b_{1}=2.0$ using the fitting formula given by Lazeyras et al.
(2016), which is sufficiently small compared to $17/21$. The tidal bias
parameter $(b_{s^{2}}/b_{1})$ appearing in $F_{\rm t}$ is also calculated to
be $b_{s^{2}}/b_{1}=(-2/7)(1-1/b_{1})=-0.14$ for the linear Lagrangian bias
model (e.g., Desjacques et al., 2018b), and its value is also smaller than
$2/7$. Therefore, even if the non-linear bias parameter is present, $F_{\rm
g}\sigma_{8}$ and $F_{\rm t}\sigma_{8}$ are expected to be larger than zero:
i.e., $F_{\rm g}\sigma_{8}\geq 0$ and $F_{\rm t}\sigma_{8}\geq 0$. We will
discuss the validity of the analysis results when these conditions are imposed
in Section 9.5 by comparing them with the results when $F_{\rm g}\sigma_{8}$
and $F_{\rm t}\sigma_{8}$ can take negative values.

In the cases of Horndeski and DHOST theories, the parameterisation we adopt
describes the time evolution of the coefficients of the tidal and shift terms
as powers of $\Omega_{\rm m}$ (Section 3.4), implicitly assuming that these
coefficients are always positive: i.e., $F_{\rm s}\sigma_{8}\geq 0$, $G_{\rm
s}\sigma_{8}\geq 0$, and $G_{\rm t}\sigma_{8}\geq 0$. For $F_{\rm
g}\sigma_{8}$ and $F_{\rm t}\sigma_{8}$, assuming that Horndeski and DHOST
theories are not far from GR, we adopt $F_{\rm g}\sigma_{8}\geq 0$ and $F_{\rm
t}\sigma_{8}\geq 0$, just like GR.

The Fisher analysis shows that the BOSS data cannot detect the $E_{{\rm
s},{\rm t}}$ signals and only give them an upper limit (Section 7). This fact
means that as $E_{{\rm s},{\rm t}}$ approach zero, the parameters $\xi_{{\rm
s},{\rm t}}=\log_{\Omega_{\rm m}}E_{{\rm s},{\rm t}}$ can be as large as
desired because of $\Omega_{\rm m}<1$. Therefore, in this analysis, we set the
upper limit of $\xi_{{\rm s},{\rm t}}$ to $(\xi_{\rm s,t})_{\rm
fid}+3\sigma_{\rm fisher}(\xi_{\rm s,t})$, which is narrower than the basic
setting. If $\xi_{{\rm s},{\rm t}}$ reach their upper bounds set here, we
report only the lower bounds for those parameters as the final results.

We summarise the results of the above discussion in Table 8.

##  8 Goodness of fit

In this section, we examine the extent to which our analysis can give good
fits to the 2PCF and 3PCF measurements from the BOSS data or Patchy mocks for
a variety of cases, before presenting specific parameter constraint values in
Section 9.

For this purpose, we calculate the minimum of $\chi^{2}=-2\ln{\cal L}$ (6.1),
denoted $\chi^{2}_{\rm min}$, from the best-fit parameter values obtained from
the joint analysis of the 2PCF and 3PCF. We use two multipole 2PCFs ($\xi_{0}$
and $\xi_{2}$), two monopole 3PCFs ($\zeta_{000}$ and $\zeta_{110}$), and two
quadrupole 3PCFs ($\zeta_{202}$ and $\zeta_{112}$) in this analysis; the
assumed gravity theories are GR, Horndeski, and DHOST theories. Tables 9-14
show the $\chi^{2}_{\rm min}$ divided by the degrees of freedom (DoF), i.e.,
the reduced $\chi^{2}_{\rm min}$, and the corresponding one-tailed $p$-values.
At two redshift bins, $z=0.38$ and $z=0.61$, results are presented for NGC
only, SGC only, and both NGC and SGC. In Horndeski and DHOST theories, we
constrain the common parameters $\xi_{f,{\rm s},{\rm t}}$ among different
redshift bins using the samples at both the redshift bins. Finally, we also
include the results of the analysis using only 2PCF.

If the theoretical model fits the measurements well, the $p$-value should be
close to $0.5$. A $p$-value close to $1$ does not mean that the theoretical
model is correct, but that the theoretical model can explain the measurements
within the error range, thanks to too large statistical errors in the
measurements. On the other hand, a $p$-value close to $0$ indicates that the
theoretical model cannot explain the measurements. In this paper, we decide
that if $p<0.05$, attention should be paid to the consistency between the
theoretical model and the measurements, and if $p<0.01$, there is an apparent
discrepancy between them. We write in bold the $\chi^{2}$ and $p$ values shown
in Tables 9-14 if $p<0.01$. Finally, we comment on the behaviour of $p$-values
when combining different galaxy samples. For example, suppose that the reduced
$\chi^{2}_{\rm min}$ is larger than $1$: i.e., $\chi^{2}_{\rm min}/{\rm
DoF}>1$. In this case, if we increase the values of $\chi^{2}_{\rm min}$ and
DoF by an equal factor while keeping the value of the reduced $\chi^{2}_{\rm
min}$, the resulting $p$-value will be smaller than the original value, and
conversely, if $\chi^{2}_{\rm min}/{\rm DoF}<1$, it will be larger than the
original value. Since we treat the different galaxy samples as statistically
independent, a similar situation occurs in analyses with multiple galaxy
samples. Thus, if the $p$-value obtained from each galaxy sample is small,
combining galaxy samples will yield a smaller $p$-value.

Section 8.1 reports an unexplained discrepancy between the 3PCF measured from
the BOSS galaxy data at $z=0.38$ and our theoretical model on large scales,
even considering DHOST theories, which is beyond GR. Section 8.2 shows that
this discrepancy between the data and the theoretical model appears from the
monopole 3PCF. Section 8.3 shows that the discrepancy still appears even when
the parameter prior set introduced in Section 7.5 is removed. Section 8.4
confirms that the discrepancy does not appear in the analysis using the Patchy
mock. Finally, as a temporary measure, we rescale the covariance matrix of the
3PCF at $z=0.38$ to generate acceptable $p$-values in Section 8.5. Section 9
will report the parameter estimation results with and without rescaling the
covariance matrix.

###  8.1 BOSS galaxies

BOSS DR12
---
$\chi^{2}_{\rm min}/{\rm DoF}$ ($p$-value)
| NGC + SGC | NGC | SGC
2PCF only ($z_{\rm eff}=0.38$) |  $56.04/57$ $(0.511)$ |  $32.18/28$ $(0.267)$ |  $23.36/28$ $(0.715)$
2PCF only ($z_{\rm eff}=0.61$) |  $80.24/57$ $(0.023)$ |  $42.08/28$ $(0.043)$ |  $36.94/28$ $(0.120)$
GR ($z_{\rm eff}=0.38$) |  $\mathbf{488.38/396}$ $\mathbf{(0.001)}$ |  $238.22/197$ $(0.024)$ |  $\mathbf{248.72/197}$ $\mathbf{(0.007)}$
GR ($z_{\rm eff}=0.61$) |  $428.60/396$ $(0.125)$ |  $218.48/197$ $(0.140)$ |  $209.24/197$ $(0.262)$
Horndeski ($z_{\rm eff}=0.38$) |  $\mathbf{488.24/395}$ $\mathbf{(0.001)}$ |  $236.38/196$ $(0.026)$ |  $\mathbf{251.62/196}$ $\mathbf{(0.004)}$
Horndeski ($z_{\rm eff}=0.61$) |  $427.56/395$ $(0.125)$ |  $218.06/196$ $(0.134)$ |  $209.16/196$ $(0.247)$
Horndeski ($z_{\rm eff}=0.38,\,0.61$) |  $\mathbf{918.76/792}$ $\mathbf{(0.001)}$ |  $456.50/394$ $(0.016)$ |  $458.94/394$ $(0.013)$
DHOST ($z_{\rm eff}=0.38$) |  $\mathbf{487.48/394}$ $\mathbf{(0.001)}$ |  $235.80/195$ $(0.024)$ |  $\mathbf{248.36/195}$ $\mathbf{(0.006)}$
DHOST ($z_{\rm eff}=0.61$) |  $427.62/394$ $(0.117)$ |  $217.94/195$ $(0.125)$ |  $209.06/195$ $(0.233)$
DHOST ($z_{\rm eff}=0.38,\,0.61$) |  $\mathbf{918.04/791}$ $\mathbf{(0.001)}$ |  $455.96/393$ $(0.015)$ |  $458.88/393$ $(0.012)$
Table 9: The reduced $\chi^{2}$ and $p$-values (in round brackets) obtained
from the joint analysis of the 2PCF and 3PCF are shown. These values are
written in bold if $p<0.01$. The minimum $\chi^{2}$, denoted $\chi_{\rm
min}^{2}$, is calculated from the best-fit parameters. The data used is the
BOSS DR12 galaxy, split into two sky regions, NGC and SGC, and two redshift
bins, $z=0.38$ and $z=0.61$. The joint analysis shows the results for three
gravity theories, i.e., GR, Horndeski, and DHOST theories; for Horndeski and
DHOST theories, also shown are the results using the two redshift bins to
constrain the parameters $\xi_{f,{\rm s},{\rm t}}$, which characterises the
time evolution of the linear and non-linear effects of the velocity field.
Furthermore, the results for the 2PCF-only analysis are shown. The
combinations of the 2PCF and 3PCF multipole components used in this analysis
correspond to Case $1$ and Case $7$ in Eq. 7.4.

Table 9 shows the results from the analysis method described in this paper. We
have performed the MCMC analysis (Section 6) using $\xi_{0}$, $\xi_{2}$,
$\zeta_{000}$, $\zeta_{110}$, $\zeta_{202}$, and $\zeta_{112}$ measured from
the BOSS galaxy data (Section 4), the covariance matrix computed from the
$2048$ Patchy mocks (Section 5), and the flat prior of the parameter range
(Section 7.5).

First, we focus on the analysis case using only the 2PCF. For the NGC+SGC
sample, the obtained $p$-values are $p=0.511$ at $z=0.38$ and $p=0.023$ at
$z=0.61$. This $p=0.023$ indicates a small amount of a poor fit between the
model and the measurements, but we consider it not problematic.

Next, turning to the joint analysis results of the 2PCF and 3PCF assuming GR,
we find that the $p$-value at $z=0.38$ obtained for the NGC+SGC sample is
extremely small, $0.001$. At $z=0.38$, the results for only NGC and only SGC
are $p=0.024$ and $p=0.007$, indicating that the SGC sample is more
problematic than the NGC. On the other hand, the $p$-value at $z=0.61$ for the
NGC+SGC sample is $p=0.125$, indicating that our model explains the measured
values without problems.

Finally, for Horndeski and DHOST theories, we find results similar to the GR
case: the $p$-value is $0.001$ at $z=0.38$ and $p\sim 0.1$ at $z=0.61$ for the
NGC+SGC sample.

Thus, we conclude that there is an unexplained discrepancy between the 3PCF
measurement from the BOSS sample at $z=0.38$ and the theoretical model we are
using. Even Horndeski and DHOST theories, which are modified gravity theories
beyond GR, cannot explain this discrepancy.

###  8.2 Monopole- or Quadrupole-only 3PCF

Joint analysis with monopole 3PCFs ($\zeta_{000}$ and $\zeta_{110}$) only
---
$\chi^{2}_{\rm min}/{\rm DoF}$ ($p$-value)
| NGC + SGC | NGC | SGC
GR ($z_{\rm eff}=0.38$) |  $\mathbf{252.08/196}$ $\mathbf{(0.004)}$ |  $122.56/97$ $(0.041)$ |  $127.32/97$ $(0.021)$
GR ($z_{\rm eff}=0.61$) |  $216.54/196$ $(0.150)$ |  $109.00/97$ $(0.191)$ |  $106.94/97$ $(0.230)$
Horndeski ($z_{\rm eff}=0.38$) |  $\mathbf{252.02/195}$ $\mathbf{(0.004)}$ |  $122.54/96$ $(0.035)$ |  $127.30/96$ $(0.018)$
Horndeski ($z_{\rm eff}=0.61$) |  $216.40/195$ $(0.140)$ |  $109.02/96$ $(0.172)$ |  $106.90/96$ $(0.210)$
Horndeski ($z_{\rm eff}=0.38,\,0.61$) |  $\mathbf{469.58/392}$ $\mathbf{(0.004)}$ |  $231.70/194$ $(0.033)$ |  $234.68/194$ $(0.024)$
DHOST ($z_{\rm eff}=0.38$) |  $\mathbf{251.78/194}$ $\mathbf{(0.003)}$ |  $122.48/95$ $(0.030)$ |  $127.28/95$ $(0.015)$
DHOST ($z_{\rm eff}=0.61$) |  $216.40/194$ $(0.129)$ |  $108.96/95$ $(0.155)$ |  $106.90/95$ $(0.190)$
DHOST ($z_{\rm eff}=0.38,\,0.61$) |  $\mathbf{469.86/391}$ $\mathbf{(0.004)}$ |  $231.78/193$ $(0.029)$ |  $234.94/193$ $(0.021)$
Table 10:  Same as Table 9, except that only the monopole component of the
3PCF is used in the joint analysis of the 2PCF and 3PCF. The combination of
the 2PCF and 3PCF multipole components used in this analysis corresponds to
Case $2$ in Eq. (7.4).  Joint analysis with quadrupole 3PCFs ($\zeta_{202}$
and $\zeta_{112}$) only
---
$\chi^{2}_{\rm min}/{\rm DoF}$ ($p$-value)
| NGC + SGC | NGC | SGC
GR ($z_{\rm eff}=0.38$) |  $280.18/252$ $(0.107)$ |  $146.56/125$ $(0.091)$ |  $132.68/125$ $(0.302)$
GR ($z_{\rm eff}=0.61$) |  $281.58/252$ $(0.097)$ |  $148.04/125$ $(0.078)$ |  $132.48/125$ $(0.306)$
Horndeski ($z_{\rm eff}=0.38$) |  $279.66/251$ $(0.103)$ |  $144.68/124$ $(0.099)$ |  $132.38/124$ $(0.287)$
Horndeski ($z_{\rm eff}=0.61$) |  $281.68/251$ $(0.089)$ |  $148.14/124$ $(0.069)$ |  $132.48/124$ $(0.285)$
Horndeski ($z_{\rm eff}=0.38,\,0.61$) |  $563.66/504$ $(0.034)$ |  $294.48/250$ $(0.028)$ |  $265.08/250$ $(0.245)$
DHOST ($z_{\rm eff}=0.38$) |  $279.18/250$ $(0.099)$ |  $144.52/123$ $(0.090)$ |  $132.28/123$ $(0.268)$
DHOST ($z_{\rm eff}=0.61$) |  $281.68/250$ $(0.082)$ |  $148.04/123$ $(0.062)$ |  $132.50/123$ $(0.263)$
DHOST ($z_{\rm eff}=0.38,\,0.61$) |  $563.26/503$ $(0.032)$ |  $294.66/249$ $(0.025)$ |  $265.24/249$ $(0.229)$
Table 11:  Same as Table 9, except that only the quadrupole component of the
3PCF is used in the joint analysis of the 2PCF and 3PCF. The combination of
the 2PCF and 3PCF multipole components used in this analysis corresponds to
Case $3$ in Eq. (7.4).

We investigate whether the discrepancy between the 3PCF measurement from the
galaxy sample at $z=0.38$ and the theoretical model, shown in Table 9,
originates from the monopole or quadrupole component.

For this purpose, Tables 10 and 11 show the joint analysis results using only
monopole 3PCFs or only quadrupole 3PCFs in addition to the monopole and
quadrupole 2PCFs. For a fair comparison with Table 9, the prior distributions
of the parameters used here are those given in Table 8. For the NGC+SGC at
$z=0.38$, the $p$-value obtained using the monopole 3PCFs is less than $0.01$
in all three gravity theories, whereas the $p$-value obtained using the
quadrupole 3PCFs is $p\sim 0.1$. Therefore, we can conclude that the monopole
component of the 3PCF measurement at $z=0.38$ is inconsistent with the
theoretical model.

###  8.3 No prior in DHOST theories

No prior in DHOST
---
$\chi^{2}_{\rm min}/{\rm DoF}$ ($p$-value)
| NGC + SGC | NGC | SGC
DHOST ($z_{\rm eff}=0.38$) |  $\mathbf{486.14/394}$ $\mathbf{(0.001)}$ |  $235.80/195$ $(0.024)$ |  $\mathbf{246.28/195}$ $\mathbf{(0.008)}$
DHOST ($z_{\rm eff}=0.61$) |  $428.16/394$ $(0.114)$ |  $218.58/195$ $(0.119)$ |  $209.58/195$ $(0.225)$
DHOST ($z_{\rm eff}=0.38,\,0.61$) |  $\mathbf{916.04/791}$ $\mathbf{(0.001)}$ |  $458.02/393$ $(0.013)$ |  $456.78/393$ $(0.014)$
Table 12:  The analysis is repeated as in Table 9, except that the prior is
removed for all parameters related to non-linear effects in DHOST theories,
$F_{\rm g}\sigma_{8}$, $F_{\rm s}\sigma_{8}$, $F_{\rm t}\sigma_{8}$,
$\xi_{f}$, $\xi_{\rm s}$, and $\xi_{\rm t}$, allowing them to vary from
$-\infty$ to $+\infty$.

As an attempt to explain the discrepancy between the 3PCF measurement from the
galaxy sample at $z=0.38$ and the theoretical model, we remove all flat prior
for the non-linear parameters, $F_{\rm g}\sigma_{8}$, $F_{\rm s}\sigma_{8}$,
$F_{\rm t}\sigma_{8}$, $\xi_{f}$, $\xi_{\rm s}$, and $\xi_{\rm t}$, set in
Table 8 and perform parameter fitting without imposing any prior. In
particular, we investigate the possibility that imposing the conditions
$F_{\rm g}\geq 0$ and $F_{\rm t}\geq 0$ on the parameters with the non-linear
bias may have caused some problems fitting the monopole 3PCF. This subsection
focuses on DHOST theories because they have the largest number of parameters
to be varied.

Table 12 summarises the results of the calculations and confirms that the
$p$-value obtained from the NGC+SGC sample at $z=0.38$ is $0.001$, even if we
assume no prior for the non-linear parameters. Therefore, we can conclude that
the discrepancy between the galaxy data and the theoretical model at $z=0.38$
is not due to the prior imposed in Table 8.

###  8.4 Patchy mocks

MultiDark-Patchy mocks
---
$\chi^{2}_{\rm min}/{\rm DoF}$ ($p$-value)
| NGC + SGC | NGC | SGC
2PCF only ($z_{\rm eff}=0.38$) |  $(57.76\pm 12.25)/57$ $(0.447^{+0.416}_{-0.331})$ |  $(28.84\pm 7.92)/28$ $(0.421^{+0.408}_{-0.297})$ |  $(27.75\pm 8.15)/28$ $(0.478^{+0.401}_{-0.332})$
2PCF only ($z_{\rm eff}=0.61$) |  $(56.50\pm 9.54)/57$ $(0.494^{+0.332}_{-0.301})$ |  $(27.99\pm 6.79)/28$ $(0.465^{+0.352}_{-0.289})$ |  $(27.52\pm 6.76)/28$ $(0.490^{+0.345}_{-0.298})$
GR ($z_{\rm eff}=0.38$) |  $(364.08\pm 29.31)/396$ $(0.873^{+0.115}_{-0.346})$ |  $(181.04\pm 20.14)/197$ $(0.786^{+0.186}_{-0.382})$ |  $(180.91\pm 21.55)/197$ $(0.788^{+0.189}_{-0.408})$
GR ($z_{\rm eff}=0.61$) |  $(362.41\pm 24.97)/396$ $(0.886^{+0.099}_{-0.274})$ |  $(179.10\pm 18.41)/197$ $(0.815^{+0.158}_{-0.339})$ |  $(181.45\pm 18.02)/197$ $(0.780^{+0.181}_{-0.342})$
Horndeski ($z_{\rm eff}=0.38$) |  $(363.90\pm 29.45)/395$ $(0.867^{+0.121}_{-0.353})$ |  $(180.68\pm 20.20)/196$ $(0.777^{+0.193}_{-0.386})$ |  $(180.62\pm 21.48)/196$ $(0.778^{+0.197}_{-0.410})$
Horndeski ($z_{\rm eff}=0.61$) |  $(361.99\pm 25.00)/395$ $(0.882^{+0.102}_{-0.278})$ |  $(178.48\pm 18.49)/196$ $(0.810^{+0.162}_{-0.343})$ |  $(181.21\pm 17.92)/196$ $(0.768^{+0.189}_{-0.344})$
Horndeski ($z_{\rm eff}=0.38,\,0.61$) |  $(728.47\pm 36.96)/792$ $(0.948^{+0.048}_{-0.203})$ |  $(360.88\pm 25.46)/394$ $(0.883^{+0.102}_{-0.284})$ |  $(363.77\pm 26.65)/394$ $(0.860^{+0.122}_{-0.319})$
DHOST ($z_{\rm eff}=0.38$) |  $(363.26\pm 29.37)/394$ $(0.865^{+0.123}_{-0.355})$ |  $(180.14\pm 20.18)/195$ $(0.770^{+0.199}_{-0.388})$ |  $(180.39\pm 21.43)/195$ $(0.766^{+0.207}_{-0.412})$
DHOST ($z_{\rm eff}=0.61$) |  $(361.46\pm 25.03)/394$ $(0.879^{+0.105}_{-0.282})$ |  $(178.00\pm 18.40)/195$ $(0.803^{+0.167}_{-0.345})$ |  $(180.95\pm 17.98)/195$ $(0.757^{+0.197}_{-0.348})$
DHOST ($z_{\rm eff}=0.38,\,0.61$) |  $(728.23\pm 37.03)/791$ $(0.946^{+0.050}_{-0.208})$ |  $(360.45\pm 25.33)/393$ $(0.879^{+0.105}_{-0.286})$ |  $(363.47\pm 26.59)/393$ $(0.855^{+0.127}_{-0.322})$
Table 13: Same as Table 9, except that the $\chi_{\rm min}^{2}$ is calculated
from each of $100$ Patchy mocks, showing their mean values and standard
deviations, and the means and $1\sigma$ errors of the corresponding
$p$-values.

Table 13 shows the means and standard deviations of the $\chi_{\rm min}^{2}$
and the corresponding means and $1\sigma$ errors of the $p$-values obtained
from the $100$ Patchy mock catalogues. The setup for the data analysis is the
same as that performed in Table 9.

In all $30$ cases shown in Table 13, the mean $p$-values obtained are almost
always $p\gtrsim 0.5$, both in the analysis using only the 2PCF and in the
joint analysis with the 3PCF. This result means that our 2PCF and 3PCF
theoretical templates fit well with the Patchy mock simulation data,
indicating that the small $p$-values found in Table 9 are a peculiar property
of the BOSS galaxies.

As two representative examples, the rest of this subsection focuses on the
DHOST theory analyses using only the SGC sample at $z=0.38$ and all four
galaxy samples (NGC+SGC at $z=0.38,0.61$). The reasons are as follows: (1) our
primary goal is to test DHOST theories; (2) the analysis of the SGC at
$z=0.38$ in the BOSS data gives a p-value of $0.006$, which is the most
significant discrepancy from the theoretical model among the four galaxy
samples; (3) the analysis using all four BOSS galaxy samples gives our final
results in Section 9.

For the SGC sample at $z=0.38$, the $\chi_{\rm min}^{2}$ values for the BOSS
samples and the Patchy mocks are

$\displaystyle\chi_{\rm min}^{2}\left({\rm BOSS}\right)$ $\displaystyle=$
$\displaystyle 248.36,$ $\displaystyle\chi_{\rm min}^{2}\left({\rm
Patchy\,mocks}\right)$ $\displaystyle=$ $\displaystyle 180.39\pm 21.43,$ (8.1)

where ${\rm DoF}=195$. The above result means that assuming that the
$\chi_{\rm min}^{2}$ follows a Gaussian distribution, the BOSS galaxy sample
deviates from the Patchy mocks at the $3.2\sigma$ significance level.

For the NGC+SGC sample at $z=0.38,0.61$, we have

$\displaystyle\chi_{\rm min}^{2}\left({\rm BOSS}\right)$ $\displaystyle=$
$\displaystyle 918.04,$ $\displaystyle\chi_{\rm min}^{2}\left({\rm
Patchy\,mocks}\right)$ $\displaystyle=$ $\displaystyle 728.23\pm 37.03$ (8.2)

where ${\rm DoF}=791$. This result implies a discrepancy between the BOSS
galaxy sample and the Patchy mocks at the $5.1\sigma$ level. Thus, we conclude
that the discrepancy with the theoretical model in the BOSS galaxies cannot be
explained by the statistical scatter of the Patchy mocks.

Although Table 13 has shown the results obtained from $100$ Patchy mocks, for
a more detailed exploration, we perform MCMC analysis on all $2048$ publicly
available Patchy mocks for the two examples above to see if it is possible to
find realizations that return the similar $p$-values to the BOSS galaxy
sample. For the SGC sample of $z=0.38$, only one Patchy mock catalogue gives
$p=0.005$, close to the BOSS result. In this case, the Patchy mocks have a
probability of $100\times(1/2048)=0.0488\%$ to reproduce the BOSS galaxy
results. On the other hand, using all four galaxy samples, not a single
catalogue among the $2048$ Patchy mocks reproduced the BOSS results. This
result means that the BOSS result has less than a $0.0488\%$ probability of
appearing in the Patchy mocks. These results are consistent with the
$3.2\sigma$ and $5.1\sigma$ discrepancies between the BOSS and Patchy mock
data presented in Eqs. (8.1) and (8.2).

Figure 11 visualizes the results for the DHOST theory analysis in Tables 9 and
13. As expected, the histogram of $\chi_{\rm min}^{2}$ computed from the
Patchy mocks (blue bars) can be well approximated by a Gaussian function
(orange line) with input values for the mean and standard deviation of
$\chi_{\rm min}^{2}$ computed from the Patchy mocks. In the cases of SGC at
$z=0.38$ (top right panel) and NGC+SGC at $z=0.38,0.61$ (bottom left panel),
we compute the histograms from the $2048$ Patchy mocks; otherwise, we compute
them from $100$ Patchy mocks. Also, we plot the $\chi_{\rm min}^{2}$ values
obtained from the BOSS data in magenta.

Figure 11:  Visualizations of the DHOST theory analysis results in Tables 9,
13 and 14. The histograms of $\chi_{\rm min}^{2}$ computed from the Patchy
mocks are shown. In the cases of SGC at $z=0.38$ (top right panel) and NGC+SGC
at $z=0.38,0.61$ (bottom left panel), the histograms are computed from $2048$
Patchy mocks; otherwise, they are computed from $100$ Patchy mocks. Also,
Gaussian functions (orange lines) with input values for the mean and standard
deviation of $\chi_{\rm min}^{2}$ computed from the Patchy mocks are shown.
The $\chi_{\rm min}^{2}$ values obtained from the BOSS data are plotted in
magenta. Also shown are the results for the BOSS data using the rescaled 3PCF
covariance matrix at $z=0.38$ (Section 8.5) in dashed red lines.

###  8.5 Rescaling of the covariance matrix

$\chi^{2}_{\rm min}/{\rm DoF}$ ($p$-value)
---
| NGC + SGC | NGC | SGC
GR ($z_{\rm eff}=0.38\,[{\rm rescaled}]$) |  $413.86/396$ $(0.258)$ |  $209.98/197$ $(0.250)$ |  $202.56/197$ $(0.378)$
Horndeski ($z_{\rm eff}=0.38\,[{\rm rescaled}]$) |  $413.66/395$ $(0.249)$ |  $208.24/196$ $(0.261)$ |  $202.36/196$ $(0.363)$
Horndeski ($z_{\rm eff}=0.38\,[{\rm rescaled}],\,0.61$) |  $844.18/792$ $(0.097)$ |  $428.14/394$ $(0.114)$ |  $412.46/394$ $(0.251)$
DHOST ($z_{\rm eff}=0.38\,[{\rm rescaled}]$) |  $412.88/394$ $(0.246)$ |  $207.82/195$ $(0.252)$ |  $202.34/195$ $(0.344)$
DHOST ($z_{\rm eff}=0.38\,[{\rm rescaled}],\,0.61$) |  $843.66/791$ $(0.095)$ |  $428.00/393$ $(0.108)$ |  $412.36/393$ $(0.241)$
Table 14:  Same as Table 9, except that the covariance matrix of the 3PCF at
$z=0.38$ is rescaled as in Eq. (8.3).

We have discussed the discrepancy between the 3PCF measured from the BOSS data
at $z=0.38$ and the corresponding theoretical model. Unfortunately, this paper
cannot provide a definitive answer to this question.

There are three possible reasons for this discrepancy. The first concern is
about the calculation of the covariance matrix. There may be physical effects
that the Patchy mock used to calculate the covariance matrix needs to account
for fully. For example, it is necessary to verify to what extent non-linear
galaxy bias effects (Desjacques et al., 2018b) and super-sample covariance
effects (Takada & Hu, 2013) are correctly included in the Patchy mock. The
second concern is about the theoretical model. For example, the theoretical
model may have new physical effects dominating large scales at low redshifts.
If so, we also need to account for that effect in the covariance matrix
simultaneously. Finally, we are concerned with the observed galaxy data. There
may be unknown observational effects that the weight function in Eq. (4.1)
cannot explain. In any case, the findings in this section indicate the
importance of discussing the validity of cosmological analyses that consider
the 2PCF and 3PCF simultaneously.

This paper assumes that the discrepancy between the BOSS galaxy sample and the
theoretical model is due to an improper covariance matrix for the 3PCF
calculated with the Patchy mock. Therefore, as a temporary measure, we decided
to rescale the 3PCF covariance matrix at $z=0.38$ to increase the obtained
$p$-value to an acceptable value. Specifically, we rescale the 3PCF covariance
matrix at $z=0.38$ as follows:

$\displaystyle{\rm Cov}[{\rm 3PCF}]_{\rm rescaled}=A\,{\rm Cov}[{\rm 3PCF}].$
(8.3)

where the rescaling factor $A$ is $A=1.15$ and $A=1.25$ for NGC and SGC,
respectively. The values of $A$ are determined so that the resulting
$p$-values at $z=0.38$ become similar to those at $z=0.61$.

Table 14 summarises the results of repeating the same analysis as Table 9
using the rescaled covariance matrix; Figure 11 visualizes the results for the
DHOST theory analysis in Tables 14. As expected, $p\gtrsim 0.1$ for the
NGC+SGC sample at $z=0.38$. Thus, if the discrepancy between the galaxy data
and the theoretical model in the 3PCF measurement is due to the covariance
matrix computed by the Patchy mocks, we find that we can solve this problem by
increasing the resulting covariance matrix by $15-25\%$. We will give the
results using this rescaled covariance matrix as the final result of this
paper when we perform parameter estimation in Section 9 using the galaxy data
at $z=0.38$.

##  9 Results

This section calculates the mean, standard deviation, $\pm 1\sigma$ errors,
and $95\%$ upper and lower bounds for the parameters computed from the
likelihoods, where we perform parameter estimation for each BOSS DR12 galaxy
and Patchy mock data. When using the Patchy mock data, we compute the mean,
standard deviation, $\pm 1\sigma$ errors, and $95\%$ limits for the parameters
from each of the $100$ Patchy mocks; then, we calculate the means and standard
deviations of them. All results here take into account both the NGC and SGC
samples. We have already given the $\chi_{\rm min}^{2}$ and $p$-values
calculated from the best-fit values of the parameters in the NGC+SGC columns
of Tables 9, 13, and 14.

The main results of this paper are Eqs. (9.7)-(9.10), which provide
constraints in $\xi_{\rm t}$ and $\xi_{\rm s}$. In Figure 21, we plot the one-
and two-dimensional likelihood distributions corresponding to these results.
Finally, we summarise the measurement results for the 3PCF multipole
components ($\zeta_{000}$, $\zeta_{110}$, $\zeta_{202}$, and $\zeta_{112}$)
from the BOSS galaxies used in this analysis in Figures 12-19.

The combination of the 2PCF and 3PCF multipoles used in the joint analysis
performed in this section corresponds to Case $7$ in Eq. (7.4); the analysis
using only the 2PCF corresponds to Case $1$. In Section 9.10, the results of
the joint analysis with only the monopole 3PCF, which corresponds to Case $2$,
are also presented and compared with the final results obtained from Case $7$.

###  9.1 Measurements

Figures 12-19 plot the measurement results of the monopole 3PCFs
($\zeta_{000}$ and $\zeta_{110}$) and the quadrupole 3PCFs ($\zeta_{202}$ and
$\zeta_{112}$) from the BOSS galaxies as a function of $r_{2}$ with $r_{1}$
fixed at $50\,h^{-1}\,{\rm Mpc}$, $80\,h^{-1}\,{\rm Mpc}$, $90\,h^{-1}\,{\rm
Mpc}$, $100\,h^{-1}\,{\rm Mpc}$, and $130\,h^{-1}\,{\rm Mpc}$ from top to
bottom; they are shown by blue circled points with $1\sigma$ error bars. Also
plotted are the 3PCF measurements from $100$ Patchy mocks (grey) and the mean
from the 3PCF measurements from $2048$ Patchy mocks (black). Finally, the
theoretical models computed from the best-fit parameter values obtained from
the DHOST theory analysis using all four BOSS samples are plotted with magenta
lines; they are shown as solid lines on the scales $r_{1},r_{2}\geq
80\,h^{-1}\,{\rm Mpc}$ used in the MCMC analysis and as dashed lines on
smaller scales. Note that the theoretical model shown by the magenta dashed
line does not need to explain the measurements from the galaxy data.

As can be seen from the lower left of Figure 14, the $\zeta_{000}$ measured
from the SGC sample at $z=0.38$ shows a significant discrepancy with the
theoretical model on large scales, which is to be expected from the results
presented in Section 8.2.

Theoretical predictions from Figures 1 and 2 indicate that the monopole and
quadrupole 3PCFs have trough-shaped signals at $r_{1}=r_{2}$. For example,
this characteristic trough signal is seen in the blue data points for
$\zeta_{112}$ measured in the NGC sample at $z=0.38$, shown in the first and
second panels from the top in the right panel of Figure 13. However, due to
the significant statistical scattering in the galaxy data, the trough signal
is not necessarily found in the blue points of all panels in Figures 12-19.

In particular, for the monopole 3PCF, the BAO peak appears at $r_{1}\simeq
r_{2}\simeq 100\,h^{-1}\,{\rm Mpc}$. Therefore, it is expected to cancel out
the trough signal, resulting in a smooth line with no irregularities when
plotting the 3PCF as a function of $r_{2}$ after fixing
$r_{1}=100\,h^{-1}\,{\rm Mpc}$. For example, as seen from the second panel
from the bottom in the right panel of Figure 12, the $\zeta_{110}$ measured
from the NGC sample at $z=0.38$ shows that the trough-shaped signal disappears
from the data points. Conversely, this is evidence of a BAO signal in the
monopole 3PCF. Although plotting the 3PCF as a function of $r_{1}=r_{2}=r$
makes it easier to see the BAO signal from the galaxy data points (e.g., see
Figure 8 and Figure $11$ in Sugiyama et al. 2021), we do not plot such a
figure because the subject of this paper is not the BAO signal.

Figure 12:  Monopole 3PCFs ($\zeta_{000}$ and $\zeta_{110}$) measured from the
NGC sample at $z=0.38$ (blue points). These plots are shown as a function of
$r_{2}$, with $r_{1}$ fixed from the top to $50\,h^{-1}\,{\rm Mpc}$,
$80\,h^{-1}\,{\rm Mpc}$, $90\,h^{-1}\,{\rm Mpc}$, $100\,h^{-1}\,{\rm Mpc}$,
and $130\,h^{-1}\,{\rm Mpc}$. The error bars are the standard deviation of the
3PCF measurements computed from $2048$ Patchy mocks. Also plotted are the 3PCF
measurements from $100$ Patchy mocks (grey) and the mean from the 3PCF
measurements from $2048$ Patchy mocks (black). Finally, the results of the
theoretical model calculated from the best-fit parameter values obtained from
the DHOST theory analysis using all four BOSS samples (Sections 9.7-9.9) are
shown by the magenta lines; they are shown as solid lines on the scales
$r_{1},r_{2}\geq 80\,h^{-1}\,{\rm Mpc}$ used in the analysis and as dashed
lines on smaller scales.

Figure 13: Same as Figure 12, except that the quadrupole 3PCF results
($\zeta_{202}$ and $\zeta_{112}$) measured from the NGC sample at $z=0.38$ are
shown.

Figure 14: Same as Figure 12, except that the monopole 3PCF results
($\zeta_{000}$ and $\zeta_{110}$) measured from the SGC sample at $z=0.38$ are
shown.

Figure 15: Same as Figure 12, except that the quadrupole 3PCF results
($\zeta_{202}$ and $\zeta_{112}$) measured from the SGC sample at $z=0.38$ are
shown.

Figure 16: Same as Figure 12, except that the monopole 3PCF results
($\zeta_{000}$ and $\zeta_{110}$) measured from the NGC sample at $z=0.61$ are
shown.

Figure 17: Same as Figure 12, except that the quadrupole 3PCF results
($\zeta_{202}$ and $\zeta_{112}$) measured from the NGC sample at $z=0.61$ are
shown.

Figure 18: Same as Figure 12, except that the monopole 3PCF results
($\zeta_{000}$ and $\zeta_{110}$) measured from the SGC sample at $z=0.61$ are
shown.

Figure 19: Same as Figure 12, except that the quadrupole 3PCF results
($\zeta_{202}$ and $\zeta_{112}$) measured from the SGC sample at $z=0.61$ are
shown.

###  9.2 $f\sigma_{8}$ constraints from the Patchy mocks in GR

MultiDark-Patchy mocks
---
| $\langle f\sigma_{8}\rangle_{\rm mean}$ | $\langle f\sigma_{8}\rangle_{\rm std}$ | $\langle f\sigma_{8}\rangle_{-1\sigma}$ | $\langle f\sigma_{8}\rangle_{+1\sigma}$ | $\langle f\sigma_{8}\rangle_{>95\%}$ | $\langle f\sigma_{8}\rangle_{<95\%}$
2PCF only ($z_{\rm eff}=0.38$) | $0.445\,(0.491)\pm 0.078$ | $0.094\pm 0.010$ | $-0.106\pm 0.013$ | $0.091\pm 0.009$ | $0.256\pm 0.074$ | $0.635\pm 0.090$
2PCF only ($z_{\rm eff}=0.61$) | $0.457\,(0.485)\pm 0.068$ | $0.083\pm 0.009$ | $-0.091\pm 0.011$ | $0.082\pm 0.008$ | $0.291\pm 0.066$ | $0.624\pm 0.077$
GR ($z_{\rm eff}=0.38$) | $0.498(0.491)\pm 0.085$ | $0.107\pm 0.008$ | $-0.118\pm 0.009$ | $0.105\pm 0.012$ | $0.283\pm 0.082$ | $0.718\pm 0.093$
GR ($z_{\rm eff}=0.61$) | $0.504(0.485)\pm 0.070$ | $0.095\pm 0.010$ | $-0.104\pm 0.010$ | $0.092\pm 0.011$ | $0.314\pm 0.065$ | $0.698\pm 0.082$

Table 15:  Constraint results for $f\sigma_{8}$ obtained from the $100$ Patchy
mocks. One hundred means, standard deviations, $\pm 1\sigma$ errors, and
$95\%$ upper and lower bounds are computed from the $100$ Patchy mocks; then,
the means and standard deviations of them are shown. Values in parentheses are
the input values for the Patchy mocks. Results are shown for two redshift bins
at $z=0.38$ and $0.61$ in combination with the NGC and SGC samples. Also shown
are the results for the 2PCF only analysis and the joint analysis with the
3PCF assuming GR. The $\chi_{\rm min}^{2}$ and $p$ values corresponding to
this table are shown in the NGC+SGC column of Table 13.

BOSS DR12
---
| $(f\sigma_{8})_{\rm mean}$ | $(f\sigma_{8})_{\rm std}$ | $(f\sigma_{8})_{-1\sigma}$ | $(f\sigma_{8})_{+1\sigma}$ | $(f\sigma_{8})_{>95\%}$ | $(f\sigma_{8})_{<95\%}$
2PCF only ($z_{\rm eff}=0.38$) | $0.446$ | $0.086$ | $-0.096$ | $0.084$ | $0.273$ | $0.620$
2PCF only ($z_{\rm eff}=0.61$) | $0.408$ | $0.086$ | $-0.095$ | $0.084$ | $0.236$ | $0.580$
GR ($z_{\rm eff}=0.38$) | $0.561$ | $0.108$ | $-0.122$ | $0.098$ | $0.348$ | $0.785$
GR ($z_{\rm eff}=0.61$) | $0.394$ | $0.091$ | $-0.099$ | $0.088$ | $0.208$ | $0.580$
GR ($z_{\rm eff}=0.38\,[{\rm rescaled}]$) | $0.549$ | $0.108$ | $-0.122$ | $0.097$ | $0.337$ | $0.776$

Table 16:  Means, standard deviations, $\pm 1\sigma$ errors, and $95\%$ upper
and lower bounds for $f\sigma_{8}$ obtained in the 2PCF-only analysis and the
joint analysis with the 3PCF using the BOSS DR12 galaxies, assuming GR for the
joint analysis with the 3PCF. Results are shown for two redshifts at $z=0.38$
and $0.61$ using the NGC and SGC samples. Also shown are the results at
$z=0.38$ for the analysis using the rescaled covariance matrix (8.3) to give
an acceptable $p$-value. The $\chi_{\rm min}^{2}$ and $p$ values corresponding
to this table are shown in the NGC+SGC column of Tables 9 and 14.

Table 16 shows the $f\sigma_{8}$ results obtained from the analysis of $100$
Patchy mocks, assuming GR.

The standard deviations of $f\sigma_{8}$ from the 2PCF-only analysis are
almost identical to those obtained from the joint analysis with the 3PCF. This
result is consistent with the results of the Fisher analysis in Section 7.
Therefore, we can conclude that neither monopole 3PCF nor quadrupole 3PCF
contributes to reducing the $f\sigma_{8}$ error. Nevertheless, note that we
can constrain the growth rate function $f$ using the joint analysis with the
3PCF by combining the $\sigma_{8}$ constraint in Section 9.6. Furthermore, in
the context of modified gravity theories, $f$ is extended to $E_{\rm f}$, and
Section 9.7 will constrain the parameter $\xi_{f}$ characterizing its time
evolution.

Looking at the mean of $f\sigma_{8}$, the results obtained in the joint
analysis with the 3PCF ($0.498$ at $z=0.38$ and $0.504$ at $z=0.61$) are
slightly closer to the values input to the Patchy mock ($0.491$ and $0.485$)
than those obtained with the 2PCF alone ($0.445$ and $0.457$). Thus, the 3PCF
information helps reduce the bias in the $f\sigma_{8}$ mean values.

###  9.3 $f\sigma_{8}$ constraints from the BOSS DR12 galaxies in GR

Table 16 summarises the results of the $f\sigma_{8}$ constraints obtained from
the BOSS galaxy under the assumption of GR. "GR ($z=0.38$) [rescaled]" means
the results using the rescaled covariance matrix (Section 8.5).

Note that the standard deviation result of $f\sigma_{8}$ does not change with
and without rescaling the 3PCF covariance matrix at $z=0.38$: i.e.,
$(f\sigma_{8})_{\rm std}=0.108$ in both cases. Thus, the $15-25\%$ difference
in the 3PCF covariance matrix due to rescaling (Section 8.5) does not
propagate significantly to the final $f\sigma_{8}$ error. The reason for this
may be mainly due to the effect of parameter degeneracy and other factors.
However, due to the decisively different $p$-values obtained (see Tables 9 and
14), we adopt the rescaled result at $z=0.38$ as the final result.

Comparing the results of the joint analysis with the 3PCF and the 2PCF-only
analysis, the former has a larger $f\sigma_{8}$ error: i.e.,
$(f\sigma_{8})_{\rm std}=0.108,\,0.091$ at $z=0.38,\,0.61$ for the joint
analysis with the 3PCF, and $(f\sigma_{8})_{\rm std}=0.086,\,0.086$ at
$z=0.38,\,0.61$ for the 2PCF-only analysis. Therefore, one may think that
adding the 3PCF information has weakened the constraint on $f\sigma_{8}$.
However, since Table 16 shows that the statistical uncertainty of
$(f\sigma_{8})_{\rm std}$ is $\sim 0.01$, the both results are statistically
consistent at the $\lesssim 2\sigma$ level.

Our final results for the $f\sigma_{8}$ constraints are as follows. The 2PCF-
only analysis gives, at the $1\sigma$ level,

$\displaystyle f\sigma_{8}$ $\displaystyle=$ $\displaystyle
0.446^{+0.084}_{-0.096}\quad{\rm at}\,z=0.38$ $\displaystyle f\sigma_{8}$
$\displaystyle=$ $\displaystyle 0.408^{+0.084}_{-0.095}\quad{\rm at}\,z=0.61,$
(9.1)

and the joint analysis with the 3PCF presents

$\displaystyle f\sigma_{8}$ $\displaystyle=$ $\displaystyle
0.549^{+0.097}_{-0.122}\quad{\rm at}\,z=0.38$ $\displaystyle f\sigma_{8}$
$\displaystyle=$ $\displaystyle 0.394^{+0.088}_{-0.099}\quad{\rm at}\,z=0.61.$
(9.2)

These $f\sigma_{8}$ constraints are consistent with the $f\sigma_{8}$ values
($f\sigma_{8}=0.485,\,0.479$ at $z=0.38,\,0.61$) calculated from the
cosmological parameters in a flat $\Lambda$CDM model (Section 1) given by
Planck 2018 (Aghanim et al., 2020). However, the $f\sigma_{8}$ result in this
analysis, which constrains $f\sigma_{8}$ with a $\sim 20\%$ precision, is not
as competitive as existing constraints (e.g., Alam et al., 2017; Ivanov et
al., 2020; Lange et al., 2022; Kobayashi et al., 2022) because we only use
large-scale information ($r\geq 80\,h^{-1}\,{\rm Mpc}$).

###  9.4 $\sigma_{8}$ constraints from the Patchy mocks in GR

MultiDark-Patchy mocks
---
| $\langle\sigma_{8}\rangle_{\rm mean}$ | $\langle\sigma_{8}\rangle_{\rm std}$ | $\langle\sigma_{8}\rangle_{-1\sigma}$ | $\langle\sigma_{8}\rangle_{+1\sigma}$ | $\langle\sigma_{8}\rangle_{>95\%}$ | $\langle\sigma_{8}\rangle_{<95\%}$
GR ($z_{\rm eff}=0.38$) | $0.741\,(0.691)\pm 0.347$ | $0.476\pm 0.142$ | $-0.626\pm 0.184$ | $0.235\pm 0.191$ | $0.024\pm 0.097$ | $1.668\pm 0.617$
GR ($z_{\rm eff}=0.61$) | $0.612\,(0.615)\pm 0.319$ | $0.415\pm 0.165$ | $-0.550\pm 0.220$ | $0.176\pm 0.156$ | $0.003\pm 0.023$ | $1.410\pm 0.636$

Table 17:  Constraint results for $\sigma_{8}$ obtained from the $100$ Patchy
mocks. One hundred means, standard deviations, $\pm 1\sigma$ errors, and
$95\%$ upper and lower bounds are computed from the $100$ Patchy mocks; then,
the means and standard deviations of them are shown. Values in parentheses are
the input values for the Patchy mocks. Results are shown for two redshift bins
at $z=0.38$ and $0.61$ in combination with the NGC and SGC samples. Also shown
are the results for the joint analysis with the 3PCF assuming GR. The
$\chi_{\rm min}^{2}$ and $p$ values corresponding to this table are shown in
the NGC+SGC column of Table 13.

MultiDark-Patchy mocks
---
Negative $F_{\rm g}$ and $F_{\rm t}$ allowed
| $\langle\sigma_{8}\rangle_{\rm mean}$ | $\langle\sigma_{8}\rangle_{\rm std}$ | $\langle\sigma_{8}\rangle_{-1\sigma}$ | $\langle\sigma_{8}\rangle_{+1\sigma}$ | $\langle\sigma_{8}\rangle_{>95\%}$ | $\langle\sigma_{8}\rangle_{<95\%}$
GR ($z_{\rm eff}=0.38$) | $1.204\,(0.691)\pm 0.429$ | $0.628\pm 0.109$ | $-0.771\pm 0.184$ | $0.499\pm 0.263$ | $0.151\pm 0.269$ | $2.409\pm 0.568$
GR ($z_{\rm eff}=0.61$) | $1.004\,(0.615)\pm 0.441$ | $0.584\pm 0.156$ | $-0.750\pm 0.198$ | $0.374\pm 0.260$ | $0.072\pm 0.176$ | $2.140\pm 0.719$

Table 18:  Same as Table 19, except that a prior with negative $F_{\rm g}$ and
$F_{\rm t}$ allowed is adopted.

BOSS DR12
---
| $(\sigma_{8})_{\rm mean}$ | $(\sigma_{8})_{\rm std}$ | $(\sigma_{8})_{-1\sigma}$ | $(\sigma_{8})_{+1\sigma}$ | $(\sigma_{8})_{>95\%}$ | $(\sigma_{8})_{<95\%}$
GR ($z_{\rm eff}=0.38$) | $0.702$ | $0.451$ | $-0.576$ | $0.221$ | $0.000$ | $1.563$
GR ($z_{\rm eff}=0.61$) | $0.568$ | $0.404$ | $-0.547$ | $0.144$ | $0.000$ | $1.323$
GR ($z_{\rm eff}=0.38\,[{\rm rescaled}]$) | $0.692$ | $0.459$ | $-0.591$ | $0.209$ | $0.000$ | $1.568$

Table 19:  Means, standard deviations, $\pm 1\sigma$ errors, and $95\%$ upper
and lower bounds for $\sigma_{8}$ obtained in the joint analysis of the 2PCF
and 3PCF using the BOSS DR12 galaxies, assuming GR. Results are shown for two
redshifts at $z=0.38$ and $0.61$ using the NGC and SGC samples. Also shown are
the results at $z=0.38$ for the analysis using the rescaled covariance matrix
(8.3) to give an acceptable $p$-value. The $\chi_{\rm min}^{2}$ and $p$ values
corresponding to this table are shown in the NGC+SGC column of Tables 9 and
14.

Table 19 summarises the results for $\sigma_{8}$ from $100$ Patchy mocks.

The mean values for $\sigma_{8}$ are $0.741$ and $0.612$ for $z=0.38$ and
$z=0.61$, respectively, in good agreement with the mock input values ($0.691$
and $0.615$). Specifically, they agree to an accuracy of $7\%$ and $0.5\%$,
respectively. Since $\sigma_{8}$ is the only physical parameter unique to the
3PCF in GR, the fact that we can estimate the $\sigma_{8}$ value with high
accuracy guarantees the validity of our analysis.

On the other hand, the $95\%$ lower limit of $\sigma_{8}$ is consistent with
zero, so we cannot detect a statistically significant signal for $\sigma_{8}$
in our analysis.

###  9.5 $\sigma_{8}$ constraints from the Patchy mocks in GR with negative
$F_{\rm g}\sigma_{8}$ and $F_{\rm t}\sigma_{8}$ allowed

This subsection discusses the validity of the priors set in Section 7.5 for
the parameters $F_{\rm g}\sigma_{8}$ and $F_{\rm t}\sigma_{8}$, which include
non-linear bias parameters. We impose the assumption that $F_{\rm g}$ and
$F_{\rm t}$ are positive, but there is no theoretical requirement that this
assumption is correct since the values of the non-linear biases are uncertain.
Therefore, as a test, we perform parameter estimation for $\sigma_{8}$ using a
prior with negative $F_{\rm g}$ and $F_{\rm t}$ allowed to check if it returns
the input values of the Patchy mocks. Specifically, the upper bounds of
$F_{\rm g}\sigma_{8}$ and $F_{\rm t}\sigma_{8}$ given in Table 8 are
multiplied by $(-1)$ to set the lower bounds of $F_{\rm g}\sigma_{8}$ and
$F_{\rm t}\sigma_{8}$. For example, we set $-2.99\leq(F_{\rm
g}\sigma_{8})_{\rm NGC,\,z=0.38}\leq 2.99$.

We summarise the results of this analysis in Table 19. This table shows that
the mean values for $\sigma_{8}$ are $1.204$ at $z=0.38$ and $1.004$ at
$z=0.61$, which are about $1.5$ times larger than the input values, $0.691$ at
$z=0.38$ and $0.615$ at $z=0.61$, in the Patchy mocks. Thus, if we allow
negative values of $F_{\rm g}$ and $F_{\rm t}$, we cannot estimate the correct
value of $\sigma_{8}$. We have no theoretical basis for explaining this fact,
but as a result of numerical experiments, we conclude that it is reasonable to
impose the conditions $F_{\rm g}\geq 0$ and $F_{\rm t}\geq 0$ in our analysis.

###  9.6 $\sigma_{8}$ constraints from the BOSS DR12 galaxies in GR

Table 19 summarises the results for the $\sigma_{8}$ constraints obtained from
the BOSS galaxies under the GR assumption. The “GR ($z=0.38$) [rescaled]”
refers to the results obtained using the rescaled covariance matrix (Section
8.5). Figure 20 plots the marginalized one- and two-dimensional posteriors of
$f\sigma_{8}$ and $\sigma_{8}$.

Similar to the results for the $f\sigma_{8}$ constraint in Section 9.3, the
results for the $\sigma_{8}$ constraint remain almost the same whether the
covariance matrix is rescaled or not. Adopting the result using the rescaled
covariance matrix as the final result, the $\sigma_{8}$ constraints at the
$1\sigma$ level are

$\displaystyle\sigma_{8}$ $\displaystyle=$ $\displaystyle
0.692^{+0.209}_{-0.591}\quad{\rm at}\,z=0.38,$ $\displaystyle\sigma_{8}$
$\displaystyle=$ $\displaystyle 0.568^{+0.144}_{-0.547}\quad{\rm at}\,z=0.61,$
(9.3)

Also, as expected from the results of the Patchy mocks, the $95\%$ lower
bounds for $\sigma_{8}$ reach $0$, so at the $95\%$ level, we get only the
upper bounds:

$\displaystyle\sigma_{8}$ $\displaystyle<$ $\displaystyle 1.568\,(95\%\,{\rm
CL})\quad{\rm at}\,z=0.38,$ $\displaystyle\sigma_{8}$ $\displaystyle<$
$\displaystyle 1.323\,(95\%\,{\rm CL})\quad{\rm at}\,z=0.61.$ (9.4)

These results are consistent with the $\sigma_{8}$ values,
($\sigma_{8}=0.681,\,\,0.606$ at $z=0.38,\,0.61$), calculated from the
cosmological parameters in a flat $\Lambda$CDM model given by Planck 2018
(Section 1).

The ratio of the standard deviation to the mean for $\sigma_{8}$ is
$(\sigma_{8})_{\rm std}/(\sigma_{8})_{\rm mean}=0.66$ at $z=0.38$ and $0.71$
at $z=0.61$, indicating that the galaxy sample at $z=0.38$ provides a better
constraint on $\sigma_{8}$. This result is consistent with the Fisher analysis
in Section 7.3.

###  9.7 $\xi_{f}$ constraints from the BOSS DR12 galaxies in Horndeski and
DHOST theories

BOSS DR12
---
| $(\xi_{f})_{\rm mean}$ | $(\xi_{f})_{\rm std}$ | $(\xi_{f})_{-1\sigma}$ | $(\xi_{f})_{+1\sigma}$ | $(\xi_{f})_{>95\%}$ | $(\xi_{f})_{<95\%}$
Horndeski ($z_{\rm eff}=0.38$) | $0.206$ | $1.016$ | $-0.777$ | $1.176$ | $-1.906$ | $2.201$
Horndeski ($z_{\rm eff}=0.61$) | $1.142$ | $1.671$ | $-1.431$ | $1.862$ | $-2.302$ | $4.480$
Horndeski ($z_{\rm eff}=0.38,\,0.61$) | $0.562$ | $0.818$ | $-0.703$ | $0.913$ | $-1.079$ | $2.226$
Horndeski ($z_{\rm eff}=0.38\,[{\rm rescaled}]$) | $0.202$ | $1.043$ | $-0.833$ | $1.201$ | $-1.921$ | $2.287$
Horndeski ($z_{\rm eff}=0.38\,[{\rm rescaled}],\,0.61$) | $0.485$ | $0.839$ | $-0.708$ | $0.967$ | $-1.216$ | $2.175$
DHOST ($z_{\rm eff}=0.38$) | $0.458$ | $1.013$ | $-0.790$ | $1.188$ | $-1.564$ | $2.474$
DHOST ($z_{\rm eff}=0.61$) | $1.248$ | $1.722$ | $-1.372$ | $1.981$ | $-2.318$ | $4.630$
DHOST ($z_{\rm eff}=0.38,\,0.61$) | $0.834$ | $0.829$ | $-0.686$ | $0.963$ | $-0.814$ | $2.484$
DHOST ($z_{\rm eff}=0.38\,[{\rm rescaled}]$) | $0.129$ | $1.131$ | $-0.895$ | $1.078$ | $-2.096$ | $2.473$
DHOST ($z_{\rm eff}=0.38\,[{\rm rescaled}],\,0.61$) | $0.791$ | $0.830$ | $-0.691$ | $0.963$ | $-0.907$ | $2.447$

Table 20:  Means, standard deviations, $\pm 1\sigma$ errors, and $95\%$ upper
and lower bounds for $\xi_{f}$ obtained in the joint analysis of the 2PCF and
3PCF using the BOSS DR12 galaxies, assuming Horndeski or DHOST theories. The
results for the two redshifts, $z=0.38$ and $0.61$, and their combined case
are shown. Both NGC and SGC samples are used for all cases. Also shown are the
results at $z=0.38$ for the analysis using the rescaled covariance matrix
(8.3) to give acceptable $p$-values. The $\chi_{\rm min}^{2}$ and $p$ values
corresponding to this table are shown in the NGC+SGC column of Tables 9 and
14.

BOSS DR12
---
| $(\xi_{\rm t})_{\rm mean}$ | $(\xi_{\rm t})_{\rm std}$ | $(\xi_{\rm t})_{-1\sigma}$ | $(\xi_{\rm t})_{+1\sigma}$ | $(\xi_{\rm t})_{>95\%}$ | $(\xi_{\rm t})_{<95\%}$
Horndeski ($z_{\rm eff}=0.38$) | $4.221$ | $4.693$ | $-5.982$ | $5.710$ | $-3.380$ | -
Horndeski ($z_{\rm eff}=0.61$) | $11.118$ | $7.354$ | $-7.204$ | $9.639$ | $-1.256$ | -
Horndeski ($z_{\rm eff}=0.38,\,0.61$) | $5.298$ | $4.257$ | $-4.023$ | $6.092$ | $-1.865$ | -
Horndeski ($z_{\rm eff}=0.38\,[{\rm rescaled}]$) | $4.129$ | $4.704$ | $-5.968$ | $5.268$ | $-3.485$ | -
Horndeski ($z_{\rm eff}=0.38\,[{\rm rescaled}],\,0.61$) | $5.151$ | $4.300$ | $-4.016$ | $6.112$ | $-2.098$ | -
DHOST ($z_{\rm eff}=0.38$) | $4.288$ | $4.589$ | $-6.002$ | $5.348$ | $-3.103$ | -
DHOST ($z_{\rm eff}=0.61$) | $11.361$ | $7.387$ | $-6.785$ | $10.112$ | $-1.183$ | -
DHOST ($z_{\rm eff}=0.38,\,0.61$) | $5.349$ | $4.217$ | $-3.980$ | $5.915$ | $-1.688$ | -
DHOST ($z_{\rm eff}=0.38\,[{\rm rescaled}]$) | $3.745$ | $4.732$ | $-6.165$ | $5.262$ | $-3.921$ | -
DHOST ($z_{\rm eff}=0.38\,[{\rm rescaled}],\,0.61$) | $5.414$ | $4.211$ | $-3.734$ | $6.007$ | $-1.655$ | -

Table 21: Same as Table 22, except that the results for $\xi_{\rm t}$ are
shown.

BOSS DR12
---
| $(\xi_{\rm s})_{\rm mean}$ | $(\xi_{\rm s})_{\rm std}$ | $(\xi_{\rm s})_{-1\sigma}$ | $(\xi_{\rm s})_{+1\sigma}$ | $(\xi_{\rm s})_{>95\%}$ | $(\xi_{\rm s})_{<95\%}$
DHOST ($z_{\rm eff}=0.38$) | $5.232$ | $3.495$ | $-2.948$ | $4.904$ | $-0.744$ | -
DHOST ($z_{\rm eff}=0.61$) | $8.791$ | $6.841$ | $-6.984$ | $8.889$ | $-2.653$ | -
DHOST ($z_{\rm eff}=0.38,\,0.61$) | $5.407$ | $3.360$ | $-2.988$ | $4.664$ | $-0.275$ | -
DHOST ($z_{\rm eff}=0.38\,[{\rm rescaled}]$) | $5.307$ | $3.519$ | $-2.854$ | $5.063$ | $-0.740$ | -
DHOST ($z_{\rm eff}=0.38\,[{\rm rescaled}],\,0.61$) | $5.378$ | $3.438$ | $-2.777$ | $4.993$ | $-0.504$ | -

Table 22:  Same as Table 22, except the results for $\xi_{\rm s}$ that varies
in DHOST theories are shown.

Table 22 summarises the constraint results for the parameter $\xi_{f}$,
defined as $\xi_{f}=\ln_{\Omega_{\rm m}}\left(E_{f}\right)=\ln_{\Omega_{\rm
m}}\left(f/\kappa_{\delta}\right)$ (3.27), characterising the time evolution
of the amplitude of the linear velocity field. In GR and Horndeski theories,
$\xi_{f}$ corresponds to the well-known parameter $\gamma$ since
$\kappa_{\delta}=1$; in GR, $\xi_{f}=\gamma=6/11$ (3.28).

Using all the four galaxy samples, at the $1\sigma$ level, we obtain

$\displaystyle\gamma$ $\displaystyle=$ $\displaystyle
0.485^{+0.967}_{-0.708}\quad\mbox{in Horndeski},$ $\displaystyle\xi_{f}$
$\displaystyle=$ $\displaystyle 0.791^{+0.963}_{-0.691}\quad\mbox{in DHOST},$
(9.5)

and at the $95\%$ confidence level, we have

$\displaystyle-1.216<\gamma<2.175\,(95\%{\rm CL})$ $\displaystyle\quad\mbox{in
Horndeski},$ (9.6) $\displaystyle-0.907<\xi_{f}<2.447\,(95\%{\rm CL})$
$\displaystyle\quad\mbox{in DHOST}.$

All results in Table 22 are consistent with GR within the $1\sigma$ level.

Note that the $\gamma$ constraints in Horndeski theories obtained here are not
directly comparable to those obtained from existing studies by, e.g., Gil-
Marín et al. (2017b). The reason is that we simultaneously vary the $\xi_{\rm
t}$ parameter characterising the tidal term in the non-linear velocity field
in Horndeski theories, while Gil-Marín et al. (2017b) use the bispectrum model
assuming GR.

###  9.8 $\xi_{\rm t}$ constraints from the BOSS DR12 galaxies in Horndeski
and DHOST theories

Table 22 summarises the constraint results for the $\xi_{\rm t}$ parameter,
defined as $\xi_{\rm t}=\ln_{\Omega_{\rm m}}\left(E_{\rm
t}\right)=\ln_{\Omega_{\rm m}}\left(\lambda_{\theta}/\kappa_{\delta}\right)$
(3.27), characterising the time evolution of the tidal term in the second-
order velocity field. In GR, $\xi_{\rm t}=15/1144$ (3.28), and if $\xi_{\rm
t}$ deviates from the GR value, it is evidence for Horndeski or DHOST
theories.

Using all the four galaxy samples, at the $1\sigma$ level, we obtain

$\displaystyle\xi_{\rm t}$ $\displaystyle=$ $\displaystyle
5.151^{+6.112}_{-4.016}\quad\mbox{in Horndeski},$ $\displaystyle\xi_{\rm t}$
$\displaystyle=$ $\displaystyle 5.414^{+6.007}_{-3.734}\quad\mbox{in DHOST},$
(9.7)

and at the $95\%$ confidence level, we have

$\displaystyle-2.098<\xi_{\rm t}\,(95\%{\rm CL})$ $\displaystyle\quad\mbox{in
Horndeski},$ (9.8) $\displaystyle-1.655<\xi_{\rm t}\,(95\%{\rm CL})$
$\displaystyle\quad\mbox{in DHOST}.$

Eqs. (9.7) and (9.8) are one of the main results in this paper. Since the
$95\%$ upper bounds of $\xi_{\rm t}$ obtained in this analysis reach the upper
bounds set by the flat prior distribution (Section 7.5), we present only the
$95\%$ lower bounds here.

All results in Table 22 are consistent with GR within the $95\%$ level.

###  9.9 $\xi_{\rm s}$ constraints from the BOSS DR12 galaxies in DHOST
theories

Table 22 summarises the constraint results for $\xi_{\rm s}$, defined as
$\xi_{\rm s}=\ln_{\Omega_{\rm m}}\left(E_{\rm s}\right)=\ln_{\Omega_{\rm
m}}\left(\kappa_{\theta}/\kappa_{\delta}\right)$ (3.27), characterising the
time evolution of the shift term in the second-order velocity field. In GR or
Horndeski theories, $\xi_{\rm s}=0$ (3.28) because
$\kappa_{\delta}=\kappa_{\theta}=1$. If $\xi_{\rm s}\neq 0$, then it is the
specific signal appearing in DHOST theories. Note that $\xi_{\rm s}\neq 0$ is
a sufficient condition for detecting DHOST theories because there can be DHOST
theories satisfying $\kappa_{\delta}=\kappa_{\theta}$ (see Section 3.4).

Using all the four galaxy samples, at the $1\sigma$ level, we obtain

$\displaystyle\xi_{\rm s}$ $\displaystyle=$ $\displaystyle
5.378^{+4.993}_{-2.777},$ (9.9)

and at the $95\%$ confidence level, we have

$\displaystyle-0.504<\xi_{\rm s}\,(95\%{\rm CL}),$ (9.10)

where we show only the lower limit of $\xi_{\rm s}$ for the same reason as for
$\xi_{\rm t}$. Eqs. (9.9) and (9.10) are the other main results of this paper
in addition to Eqs. (9.7) and (9.8).

All results in Table 22 are consistent with GR within the $95\%$ level.

For all the results obtained from Tables 22, 22, and 22, the standard
deviations of $\xi_{f,{\rm s},{\rm t}}$ obtained by combining the samples
$z=0.38$ and $z=0.61$ are smaller than those obtained at $z=0.38$ and
$z=0.61$, respectively. Therefore, future galaxy surveys with more redshift
bins should improve our $\xi_{f,{\rm s},{\rm t}}$ constraints.

Similar to the $f\sigma_{8}$ and $\sigma_{8}$ results in GR, we confirm that
the constraints of $\xi_{f,{\rm s},{\rm t}}$ are hardly affected by rescaling
the covariance matrix by $15\%-25\%$ at $z=0.38$. This finding indicates that
the results for $\xi_{f,{\rm s},{\rm t}}$ presented here will not change
significantly even if future re-analysis from a better mock simulation gives
acceptable $p$-values.

Figure 20:  Marginalized two- and one-dimensional posteriors of $f\sigma_{8}$
and $\sigma_{8}$ for BOSS DR12. The contours indicate $68.27\%$ and $95.45\%$
confidence levels. Asterisks indicate predictions by Planck. The left and
right panels show the cases at $z=0.38$ and $z=0.61$. The NGC and SGC samples
are always combined to obtain this result. The rescaled 3PCF covariance matrix
(Section 8.5) at $z=0.38$ is used.

Figure 21: Marginalized two- and one-dimensional posteriors of $\xi_{f}$,
$\xi_{\rm t}$, and $\xi_{\rm s}$ for BOSS DR12. DHOST theories (red) vary
these all three parameters, while Horndeski theories (blue) fix $\xi_{\rm s}$
to $\xi_{\rm s}=0$. The contours indicate $68.27\%$ and $95.45\%$ confidence
levels. Asterisks indicate predictions by GR: $\xi_{f}=6/11$, $\xi_{\rm
t}=15/1144$, and $\xi_{\rm s}=0$. The NGC and SGC samples at $z=0.38$ and
$0.61$ are combined to obtain this result. The rescaled 3PCF covariance matrix
(Section 8.5) at $z=0.38$ is used.

###  9.10 Joint analysis with the monopole 3PCF only

This subsection presents the results of a joint analysis of the monopole and
quadrupole 2PCFs ($\xi_{0}$ and $\xi_{2}$) with only the monopole 3PCFs
($\zeta_{000}$ and $\zeta_{110}$) and compares them with our main results,
revealing the importance of the information in the quadrupole 3PCFs
($\zeta_{202}$ and $\zeta_{112}$). In other words, we compare the results
corresponding to Case $2$ and Case $7$ in Eq. 7.4. For simplicity, we focus
here on the case where all four BOSS galaxy samples are used assuming DHOST
theories and present a comparison of the results for $\xi_{f}$, $\xi_{\rm t}$,
and $\xi_{\rm s}$. For Case $2$, as in Case $7$, we determine the parameter
priors according to the method described in Section 7.5, based on the results
of the Fisher analysis.

The results of the joint analysis with the monopole 3PCFs are as follows:

$\displaystyle(\xi_{\rm t})_{\rm mean}\pm(\xi_{\rm t})_{\rm std}$
$\displaystyle=$ $\displaystyle 257.451\pm 145.68,$ $\displaystyle(\xi_{\rm
s})_{\rm mean}\pm(\xi_{\rm s})_{\rm std}$ $\displaystyle=$ $\displaystyle
137.396\pm 79.215.$ (9.11)

On the other hand, adding the quadrupole 3PCFs presents $(\xi_{\rm t})_{\rm
std}=4.211$ (Table 22) and $(\xi_{\rm s})_{\rm std}=3.438$ (Table 22).

The addition of the quadrupole 3PCFs reduces the values of $(\xi_{\rm t})_{\rm
std}$ and $(\xi_{\rm s})_{\rm std}$ by a factor of $\sim 35$ and $\sim 20$,
respectively. This improvement is consistent with the Fisher analysis result
in Section 7 (see Eq. (7.6) and Table 5). Therefore, we conclude that the
quadrupole component of the 3PCF should always be used to constrain $\xi_{\rm
t}$ and $\xi_{\rm s}$. Finally, the same should hold for testing other
modified gravity theories through non-linear velocity fields.

###  9.11 Consistency check with the Fisher analysis

DHOST
---
| $(\xi_{f})_{\rm std}$ | $\sigma_{\rm fisher}(\xi_{f})$
$z_{\rm eff}=0.38$ | $1.131$ | $0.967$
$z_{\rm eff}=0.61$ | $1.722$ | $1.782$
$z_{\rm eff}=0.38,\,0.61$ | $0.830$ | $0.850$
| $(\xi_{\rm t})_{\rm std}$ | $\sigma_{\rm fisher}(\xi_{\rm t})$
$z_{\rm eff}=0.38$ | $4.732$ | $3.577$
$z_{\rm eff}=0.61$ | $7.387$ | $6.834$
$z_{\rm eff}=0.38,\,0.61$ | $4.211$ | $3.169$
| $(\xi_{\rm s})_{\rm std}$ | $\sigma_{\rm fisher}(\xi_{\rm s})$
$z_{\rm eff}=0.38$ | $3.519$ | $3.147$
$z_{\rm eff}=0.61$ | $6.841$ | $5.906$
$z_{\rm eff}=0.38,\,0.61$ | $3.438$ | $2.778$
Table 23:  Comparison of the standard deviations, $\sigma_{\rm
fisher}(\theta)$ and $(\theta)_{\rm std}$, obtained from the Fisher analysis
and MCMC for $\theta=\xi_{f},\,\xi_{\rm t},\,\xi_{\rm s}$. The results are
shown for each redshift of $z=0.38$ and $z=0.61$ and the combined case of the
two redshifts. In all cases shown here, the NGC and SGC samples are used; the
MCMC results at $z=0.38$ use the rescaled covariance matrix (Section 8.5). All
the values summarised here are those already given in Tables 7, 22, 22, and
22. When combining the results for the different galaxy samples in Table 7, we
use the standard error composition formula, assuming that each galaxy sample
is independent.

This subsection discusses the consistency between the Fisher analysis results
in Section 7 and our final results from MCMC in this section. For this
purpose, We compare the standard deviation of a parameter $\theta$ computed
from the Fisher analysis, $\sigma_{\rm fisher}(\theta)$, with that estimated
from MCMC, $(\theta)_{\rm std}$, where the parameters of interest are
$\theta=\xi_{f},\,\xi_{\rm t},\,\xi_{\rm t}$, which are the main target of
this paper.

Table 23 summarises the cases for each redshift bin of $z=0.38$ and $z=0.61$
and for using both redshift bins, assuming DHOST theories. The values shown in
this table are given from Tables 7, 22, 22, and 22. When combining the results
for the different galaxy samples in Table 7, we use the standard error
combination formula, assuming that each galaxy sample is independent.

Table 23 shows that the MCMC results satisfy $(\theta)_{\rm
std}\gtrsim\sigma_{\rm fisher}(\theta)$, indicating that the MCMC results are
consistent with the Fisher analysis results, as expected. This result
reinforces the validity of our main results shown in Tables 22, 22, and 22.

###  9.12 Comments on bias effects on shift terms

DHOST theories change the shift term of the non-linear density fluctuation
from GR, which may introduce a new bias effect in the shift term, i.e., the
shift bias parameter. Since $E_{f,\rm s,t}$ are the parameters that cancel the
$\sigma_{8}$-dependence using the coefficients of the shift term of the
density fluctuation, when the shift bias appears, $E_{f,\rm s,t}$ will also be
contaminated by the bias effect. Furthermore, the shift bias may induce bias
effects in linear and non-linear velocity fields. In such cases, we cannot use
the parameterisation $E_{f,{\rm s},{\rm t}}=\Omega_{\rm m}^{\xi_{f,{\rm
s},{\rm t}}}$ adopted in this paper to characterise the time dependence of
$E_{f,{\rm s},{\rm t}}$ because the time dependence of the bias parameter is
uncertain.

If we assume the presence of the shift bias effect, we propose simultaneously
constraining all the six parameters $(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})$ (3.15) that characterise the
growth, shift, and tidal terms in the density and velocity fields in each
galaxy sample as a more general test of modified gravity theories. In such an
analysis, we should remove the relation $G_{\rm g}=G_{\rm s}-(2/3)G_{\rm t}$
imposed in DHOST theories. In particular, the $E_{\rm s}$ parameter, which
represents the ratio of the coefficients of the shift terms of the non-linear
density and velocity fields: $E_{\rm s}=(G_{\rm s}\sigma_{8})/(F_{\rm
s}\sigma_{8})$, is always $E_{\rm s}=1$ in GR and Horndeski theories.
Therefore, testing whether $E_{\rm s}=1$ in each galaxy sample verifies the
theory of varying the shift term, such as DHOST-like theories. In other words,
it should provide a means to test the LSS consistency relation, which DHOST-
like theories violate (Section 1), using the galaxy 3PCF (or bispectrum).

##  10 Conclusions

This paper presents a joint analysis of the anisotropic two-point and three-
point correlation functions measured from the publicly available BOSS DR12
galaxy data to test cosmological modified gravity theories. This paper has two
important implications. First, it is the first work to extract cosmological
information from actual galaxy data using the anisotropic component of the
galaxy three-point correlation function induced by the RSD effect. Second,
this analysis is the first attempt to constrain the non-linear effects of
modified gravity theories from the galaxy three-point statistics.

We consider DHOST theories and their subclass, Horndeski theories, which are
the candidates for modified gravity theories (see Section 2.1). They are quite
general theoretical frameworks of scalar-tensor theories. Since the time
evolution equation of the linear density fluctuations in these theories is
scale-independent (2.6), the difference with GR appears only in the linear
growth rate $f$ in the linear theory (Hirano et al., 2019a). On the other
hand, the non-linear gravitational effect causes a difference in the scale-
dependence of the density fluctuation, which allows us to examine the
deviation from GR more clearly. Specifically, Horndeski theories change the
tidal term of the second-order density fluctuation from GR, while DHOST
theories change both the shift and tidal terms (2.10 and 3.16) (Hirano et al.,
2018). However, since non-linear bias parameters contaminate the density
fluctuations, Yamauchi & Sugiyama (2022) have pointed out that one should
investigate supposedly unbiased non-linear velocity fields induced by the RSD
effect (see Section 3.4 for a review). Specifically, they have suggested that
one should constrain the parameters $\xi_{\rm t}$ and $\xi_{\rm s}$, which
characterise the time evolution of the tidal and shift terms of the second-
order velocity field: $\xi_{\rm t}=15/1144$ in GR and $\xi_{\rm s}=0$ in GR
and Horndeski theories. Therefore, if $\xi_{\rm s}\neq 0$, then it is the
signal specific to DHOST theories; they have also pointed out that in DHOST
theories, the parameter $\gamma=\ln_{\Omega_{\rm m}}(f)$, which characterises
the time dependence of the linear growth rate $f$, is extended to
$\xi_{f}=\ln_{\Omega_{\rm m}}(f/\kappa)$ with $\kappa$ being the time-
dependent function appearing in the shift term of the density fluctuation. To
this end, we test DHOST and Horndeski theories by constraining these
parameters $\xi_{f}$, $\xi_{\rm t}$, and $\xi_{\rm s}$ using the joint
analysis method of the anisotropic 2PCF and 3PCF, established by Sugiyama et
al. (2021).

The following is a summary of the details of the analysis methodology and the
findings obtained.

  1. 1.

Following Sugiyama et al. (2019), we apply the TripoSH decomposition method to
the 3PCF to extract information about the anisotropic, i.e., _quadrupole_ ,
component of the 3PCF (see Sections 3.1 and 4.2). To simplify the data
analysis, we then use only two monopole components ($\zeta_{000}$ and
$\zeta_{110}$) and two quadrupole components ($\zeta_{202}$ and $\zeta_{112}$)
from the decomposed 3PCF. For the 2PCF, we adopt the commonly used Legendre
decomposition method and use the monopole and quadrupole components: i.e.,
$\xi_{0}$ and $\xi_{2}$. It is worth noting that $\zeta_{202}$ includes only
the $M=0$ mode that appears in Scoccimarro et al. (1999)’s decomposition
formalism, while $\zeta_{112}$ includes $M\neq 0$ modes in addition to the
$M=0$ mode. Furthermore, the TripoSH-decomposed 3PCF allows a quantitative
evaluation and detailed study of the survey window effect present in the
measured 3PCFs (see Section 4.3). Thus, this work is the first to extract
information on the $M\neq 0$ modes from actual galaxy data, taking into
account the window effect.

  2. 2.

We only use data at large scales of $80\,h^{-1}\,{\rm Mpc}\leq r\leq
150\,h^{-1}\,{\rm Mpc}$, where higher-order non-linear corrections, called
loop corrections, are not expected to contribute much to the 2PCF and 3PCF. In
order to test modified gravity theories consistently using smaller scales, it
is necessary to construct a model that includes the non-linear effects of
modified gravity theories so that they are also included in the loop
corrections. To our knowledge, only one such analysis has been performed so
far for the case of the power spectrum in $f(R)$ gravity (Song et al., 2015a).
However, it is known that various uncertainties arise in the non-linear power
spectrum in DHOST theories, such as IR cancellation breaking (Crisostomi et
al., 2020; Lewandowski, 2020) and UV divergence (Hirano et al., 2020). These
theoretical uncertainties should also appear in the bispectrum. Therefore,
focusing only on large scales is necessary to remove the theoretical
uncertainties and safely constrain the non-linear effects of modified gravity
theories. Our analysis is thus the second example of a consistent analysis
incorporating the non-linear effects of modified gravity from spectroscopic
galaxy surveys, and the first to use the galaxy three-point statistic.

  3. 3.

As a theoretical model for the 3PCF, we use the IR-resummed model (3.12)
proposed by Sugiyama et al. (2021) (see Section 3.2). This model can describe
the BAO damping effect while keeping the shape of the 3PCF in the tree-level
solution. For this model, we have investigated how the three decomposed non-
linear effects, i.e., the growth, shift, and tidal terms, affect the 3PCF
multipoles (see Figures 1 and 2 in Section 3.3). For example, in the
quadrupole components ($\zeta_{202}$ and $\zeta_{112}$), the dominant term is
the product of the linear density fluctuation and the linear velocity field
that appears during the coordinate transformation from real space to redshift
space; otherwise, the non-linear effects of the density and velocity fields
contribute to the quadrupole component to the same extent. Figures 12-19 in
Section 9.1 show the $\zeta_{000}$, $\zeta_{110}$, $\zeta_{202}$, and
$\zeta_{112}$ measured from the four BOSS galaxy samples and the corresponding
theoretical models calculated using the best-fit parameters.

  4. 4.

We have used the $2048$ publicly available Patchy mocks to compute the
covariance matrices of the 2PCF and 3PCF in Section 5. In our analysis, we
ensure that the number of data bins in the 2PCF and 3PCF is sufficiently
smaller than the number of the $2048$ mocks. In particular, the parameter
$M_{2}$ (5.5), which represents the impact of a finite number of mocks on the
final parameter error, is at most $M_{2}\sim 1.1$ (see Section 6.4).

  5. 5.

To understand the nature of the covariance matrix, we have calculated the
cumulative ${\rm S/N}$ of the 2PCF and the 3PCF in Section 5.4. The results
show that the cumulative ${\rm S/N}$ of the 3PCF has different characteristics
from that of the 2PCF. In the case of the 2PCF, the galaxy sample at $z=0.61$
with a larger volume has a smaller covariance matrix than the sample at
$z=0.38$, resulting in a larger ${\rm S/N}$ at $z=0.61$. On the other hand,
for the 3PCF, the ${\rm S/N}$ at $z=0.38$ is comparable to the ${\rm S/N}$ at
$z=0.61$. Therefore, the difference in survey volume cannot explain the
relationship between the ${\rm S/N}$ of the 3PCF at $z=0.38$ and $0.61$. A
possible explanation for this 3PCF ${\rm S/N}$ behaviour is that the
covariance matrix of the 3PCF depends strongly on the number density of the
galaxies (see Sugiyama et al., 2020): the BOSS sample at $z=0.38$ has a higher
number density than the sample at $z=0.61$, even with a smaller survey volume
(Table 1). We interpret this higher number density as why the ${\rm S/N}$ at
$z=0.38$ is as high as that at $z=0.61$.

  6. 6.

We have investigated the extent to which higher-order terms in the TripoSH
decomposition of the 3PCF contain cosmological information by Fisher analysis
(see Section 7.2). The results show that $\zeta_{202}$ is the main
cosmological information in the quadrupole 3PCF, while other information is
contained in the higher-order term $\zeta_{112}$ in addition to $\zeta_{202}$.
Since $\zeta_{112}$ contains the $M\neq 0$ modes in Scoccimarro et al.
(1999)’s decomposition formalism but not in $\zeta_{202}$, this result
indicates the importance of the $M\neq 0$ modes.

  7. 7.

In Section 8, we have reported that at large scales ($\geq 80\,h^{-1}\,{\rm
Mpc}$), there can be statistically significant differences between the 3PCFs
measured from the BOSS galaxies and the corresponding theoretical models,
regardless of whether we assume GR, Horndeski or DHOST theories. For example,
the $p$-value obtained from the SGC sample at $z=0.38$ is less than $0.01$,
and the $p$-value obtained from the combined sample of the four BOSS samples
is $0.001$ (see Section 8.1). This result means that the discrepancies between
the galaxy data and the theoretical models cannot be explained within the
framework of scalar-tensor theory, even if they are due to unknown physical
effects. Other results show that the discrepancy is mainly due to the monopole
component of the 3PCF rather than the quadrupole component (see Section 8.2)
and that this discrepancy cannot be explained even if the prior distribution
of the parameters is changed (see Section 8.3). Finally, we have repeated the
same analysis for the $100$ Patchy mocks as for the BOSS sample in Section
8.4. The results show a statistically significant difference of more than
$5\sigma$ between the $p$-values of the Patchy mocks and the BOSS galaxies.
Therefore, the statistical variability of the Patchy mock galaxies cannot
explain the low $p$-values ($p\sim 0.001$) obtained from the BOSS galaxies.

  8. 8.

In this paper, we assume that the discrepancy between the BOSS galaxy sample
and the theoretical model is due to an inappropriate 3PCF covariance matrix
computed from the Patch mocks. We then take a conservative approach by
artificially rescaling the 3PCF covariance matrix at $z=0.38$ by $15\%$ for
NGC and $25\%$ for SGC, resulting in acceptable $p$-values (see Section 8.5).
To confirm the validity of this method, we have presented in Section 9 the
results of constraining the parameters of interest with and without rescaling
the covariance matrix and have confirmed that there is no significant
difference in the final results obtained in these two cases. We interpret this
result as being due to a more significant degeneracy effect between the
parameters than the $\sim 20\%$ difference in the 3PCF covariance matrix.
Therefore, we do not expect that calculating the covariance matrix from
simulation data that better reproduces the distribution of the BOSS galaxies
will significantly change the results of the present paper.

  9. 9.

We have constrained $f\sigma_{8}$ from the BOSS galaxies assuming GR in
Section 9.3. There, we have shown that adding isotropic and anisotropic 3PCF
components ($\zeta_{000}$, $\zeta_{110}$, $\zeta_{202}$, and $\zeta_{112}$)
does little to improve the results compared to the 2PCF-only analysis.
Nevertheless, the analysis using the Patchy mocks shows that the 3PCF
information does help to reduce the bias of the mean value of $f\sigma_{8}$
(see Section 9.2). Finally, we obtain $f\sigma_{8}=0.549^{+0.097}_{-0.122}$ at
$z=0.38$ and $f\sigma_{8}=0.394^{+0.088}_{-0.099}$ at $z=0.61$ in the joint
analysis of the anisotropic 2PCF and 3PCF assuming GR (9.2). These
$f\sigma_{8}$ results are not as competitive as existing constraints (e.g.,
Alam et al., 2017; Ivanov et al., 2020; Lange et al., 2022; Kobayashi et al.,
2022) because we only use large-scale information ($r\geq 80\,h^{-1}\,{\rm
Mpc}$).

One may think that adding the 3PCF information does not improve the
$f\sigma_{8}$ results due to the focus on large scales only ($r\geq
80\,h^{-1}\,{\rm Mpc}$). To test this concern, we have performed a Fisher
analysis that includes small scales $(30\,h^{-1}\,{\rm Mpc}\leq r\leq
150\,h^{-1}\,{\rm Mpc})$ and find that even if we extend the used scales to
$30\,h^{-1}\,{\rm Mpc}$, there is no improvement in the $f\sigma_{8}$ results
(see Section 7.4). However, note that we use the IR-resummed tree-level model
of the 3PCF in this Fisher analysis. Therefore, if we use a theoretical model
with various loop corrections applicable down to small scales, parameter
degeneracy may break, and it may still be possible to obtain improved
$f\sigma_{8}$ constraints through a joint analysis of the 2PCF and 3PCF.

  10. 10.

We have constrained $\sigma_{8}$ from the BOSS galaxies assuming GR in Section
9.6. Thus, while the 3PCF information does not improve the $f\sigma_{8}$
constraints, it helps to break the degeneracy between parameters by providing
information on $\sigma_{8}$: e.g., it allows us to constrain $f$. We have
obtained $\sigma_{8}=0.692^{+0.209}_{-0.591}$ at $z=0.38$ and
$\sigma_{8}=0.568^{+0.144}_{-0.547}$ at $z=0.61$ at the $1\sigma$ level. These
results are consistent with $\sigma_{8}=0.681,\,0.606$ at $z=0.38,\,0.61$
calculated from the cosmological parameters in a flat $\Lambda$CDM model given
by Planck 2018. The ratio of the standard deviation to the mean for
$\sigma_{8}$ is $(\sigma_{8})_{\rm std}/(\sigma_{8})_{\rm mean}=0.66$ at
$z=0.38$ and $0.71$ at $z=0.61$, indicating that the galaxy sample at $z=0.38$
provides a better constraint on $\sigma_{8}$. This result can be attributed to
the higher number density of the sample at $z=0.38$ compared to that at
$z=0.61$, similar to the argument of the cumulative ${\rm S/N}$ in Section
5.4.

  11. 11.

Our main results, the constraints on the $\xi_{f}$, $\xi_{\rm t}$, and
$\xi_{\rm s}$ parameters in DHOST theories, are summarised in Sections 9.7,
9.8, and 9.9. There, we obtain $\xi_{f}=0.791^{+0.963}_{-0.691}$ (9.5),
$\xi_{\rm t}=5.414^{+6.007}_{-3.734}$ (9.7), and $\xi_{\rm
s}=5.378^{+4.993}_{-2.777}$ (9.9) at the $1\sigma$ level; we also have
$-0.907<\xi_{f}<2.447$ (9.6), $-1.655<\xi_{\rm t}$ (9.8), and $-0.504<\xi_{\rm
s}$ (9.10) at the $95\%$ confidence level. Since we cannot detect the signal
of the tidal and shift terms in the second-order velocity field in the present
analysis, we can only present the $95\%$ lower bounds of the $\xi_{\rm t}$ and
$\xi_{\rm s}$ parameters. These results are consistent with the GR predictions
$\xi_{f}=\gamma=6/11$, $\xi_{\rm t}=15/1144$, and $\xi_{\rm s}=0$ (see Figure
21). Moreover, we have checked the consistency of the estimated results from
the BOSS galaxy sample with the Fisher analysis for the constraints on the