1. Introduction
The late-time accelerated expansion of the Universe (Riess et al. Reference Riess1998; Perlmutter et al. Reference Perlmutter1999) cannot be explained within the framework of general relativity (GR) coupled with the matter content possessing a positive equation of state. One plausible solution is to introduce an exotic form of fluid with negative pressure, dubbed dark energy, that drives the accelerated expansion (Copeland, Sami, & Tsujikawa Reference Copeland, Sami and Tsujikawa2006). The vacuum energy (
$\Lambda$
) is the simplest candidate of the dark energy. However, there is a huge discrepancy between the theoretically predicted value and the observed value, and its incredibly small magnitude necessitates extreme fine-tuning (Weinberg Reference Weinberg1989). Additionally, the standard
$\Lambda$
CDM (
$\Lambda+$
cold dark matter) model gives rise to tensions in Hubble constant (
$H_0$
) and the amplitude of matter fluctuations (
$\sigma_8$
) between values obtained from direct measurements and those inferred by extrapolating the best-fit parameters from the cosmic microwave background (CMB) data (Efstathiou Reference Efstathiou2014; Abbott et al. Reference Abbott2018; Brieden, Gil-Marín, & Verde 2023). A wide variety of dynamic dark energy models have been explored in an effort to address these issues (Sahni & Starobinsky Reference Sahni and Starobinsky2000; Peebles & Ratra Reference Peebles and Ratra2003; Padmanabhan Reference Padmanabhan2003; Copeland et al. Reference Copeland, Sami and Tsujikawa2006; Jassal, Bagla, & Padmanabhan Reference Jassal, Bagla and Padmanabhan2005; Alam et al. Reference Alam, Sahni, Deep Saini and Starobinsky2004).
An alternative explanation of cosmic acceleration is that the gravity is not well understood on cosmological scales and needs to be modified (Clifton et al. Reference Clifton, Ferreira, Padilla and Skordis2012). Any such modification to the theory of gravity must incorporate new degrees of freedom, as the GR is the unique Lorentz invariant low energy theory (Weinberg Reference Weinberg1965). These additional degrees of freedom are associated with new fields, with scalar fields being the simplest candidates. To ensure the Ostrogradsky stability (Ostrogradsky Reference Ostrogradsky1850; Woodard Reference Woodard2015), the attention is restricted to the class of theories that propagate a single scalar degree of freedom (Langlois Reference Langlois2019; Kobayashi Reference Kobayashi2019; Kase & Tsujikawa Reference Kase and Tsujikawa2019). The Horndeski theory is the most general scalar-tensor theory that yields equations of motion that are second-order in scalar and metric fields (Horndeski Reference Horndeski1974). It also incorporates the Vainshtein screening mechanism (Vainshtein Reference Vainshtein1972), which suppresses the force mediated by the additional scalar field on solar system scales (Will Reference Will2014). The Horndeski theories were initially considered to be the only viable scalar-tensor formulation that avoids introducing additional ghost degrees of freedom (Deffayet et al. Reference Deffayet, Gao, Steer and Zahariade2011; Kobayashi, Yamaguchi, & Yokoyama Reference Kobayashi, Yamaguchi and Yokoyama2011). Beyond the Horndeski framework, however, the theory can be extended to a broader class of scalar-tensor theories, known as De-generate Higher-Order Scalar-Tensor (DHOST) theories (Zumalacárregui & García-Bellido 2014; Langlois & Noui Reference Langlois and Noui2016; Gleyzes et al. Reference Gleyzes, Langlois, Piazza and Vernizzi2015b; Gleyzes et al. Reference Gleyzes, Langlois, Piazza and Vernizzi2015a; Crisostomi, Koyama, & Tasinato Reference Crisostomi, Koyama and Tasinato2016; Ben Achour, Langlois, & Noui Reference Ben Achour, Langlois and Noui2016; Achour et al. Reference Achour2016; Langlois Reference Langlois2019; Kobayashi Reference Kobayashi2019). The action of the DHOST theories includes the higher-order terms of second derivatives of a scalar field (Langlois & Noui Reference Langlois and Noui2016). The degeneracy in kinetic terms of the scalar field and metric field enables the formulation of viable scalar-tensor theories despite the presence of higher-order equations of motion and avoids the Ostrogradsky instability by maintaining the one scalar propagating degree of freedom (Crisostomi et al. Reference Crisostomi, Koyama and Tasinato2016).
Modifications to gravity alter the growth history of density perturbations and the clustering properties of the large-scale structure (LSS) of the Universe. Measurements of the perturbation growth rate serve as a powerful tool for probing the nature of gravity (Gong Reference Gong2008; Bernardeau & Brax Reference Bernardeau and Brax2011; Sharma & Manohar Verma Reference Sharma and Manohar Verma2022; Brax & Valageas Reference Brax and Valageas2012). These measurements are typically obtained by observing the peculiar velocities of dark matter tracers along the line of sight (LOS) through redshift-space distortions (RSD) (Kaiser Reference Kaiser1987; Hamilton Reference Hamilton1997). A widely adopted approach to test gravity involves the introduction of the gravitational growth index,
$\unicode{x03B3}$
, which in the framework of the GR is predicted to have a value of approximately
$6/11$
(Linder & Cahn Reference Linder and Cahn2007; Peebles Reference Peebles1980; Fry Reference Fry1985; Lightman & Schechter Reference Lightman and Schechter1990). Previous studies have not reported any significant deviations of
$\unicode{x03B3}$
from the predictions of the GR (Zhao et al. Reference Zhao2019; Gil-Marín et al. Reference Gil-Marín2018; Sánchez et al. Reference Sánchez2017; Grieb et al. Reference Grieb2017). However, to explore the potential signatures of modified gravity more comprehensively, it is necessary to extend these investigations into the regime of non-linear growth, which may provide access to additional information about modified gravity that is not imprinted on
$\unicode{x03B3}$
(Yamauchi, Yokoyama, & Tashiro Reference Yamauchi, Yokoyama and Tashiro2017). Yamauchi & Sugiyama (Reference Yamauchi and Sugiyama2022) have proposed that incorporating non-linear velocity fields alongside non-linear density perturbations provides a robust framework for studying the non-linear features of gravity. They introduced second-order indices that can unveil new insights into modified gravity, capturing non-linear gravitational effects that remain inaccessible to the linear growth index.
The bispectrum is the lowest-order statistic sensitive to non-Gaussianity induced by non-linear gravitational interaction, making it an ideal tool for investigating the impact of modified gravity on the LSS (Gil-Marín et al. Reference Gil-Marín2017; Slepian et al. Reference Slepian2017; Pearson & Samushia Reference Pearson and Samushia2018; Sugiyama et al. Reference Sugiyama, Saito, Beutler and Seo2019; Sugiyama et al. Reference Sugiyama, Saito, Beutler and Seo2021). It is a function of the closed triangles formed by three
${\boldsymbol{k}}$
modes and captures the coupling between modes. The observed bispectrum of tracers is anisotropic along the LOS direction due to the effects of their peculiar velocities, offering additional sensitivity to deviations from GR. This anisotropy of the redshift space bispectrum can be studied by decomposing it into spherical harmonics
$Y_l^m(\hat{p})$
(Scoccimarro, Couchman, & Frieman Reference Scoccimarro, Couchman and Frieman1999; Bharadwaj, Mazumdar, & Sarkar Reference Bharadwaj, Mazumdar and Sarkar2020; Mazumdar, Bharadwaj, & Sarkar Reference Mazumdar, Bharadwaj and Sarkar2020). The bispectrum multipole moments
$B_l^m$
, which are the expansion coefficients of the redshift space bispectrum in
$Y_l^m(\hat{p})$
basis, can be directly measured from data. The
$B_l^m$
are expected to exhibit distinct signatures due to the altered LSS clustering and mode coupling introduced by the modified gravitational evolution.
The primary objective of this study is to examine how modifications in gravitational theory, at both linear and non-linear scales, influence the bispectrum multipole signal. While most of the earlier studies are either restricted to the spherically averaged bispectrum (
$B_0^0$
) only (Hirano et al. Reference Hirano, Kobayashi, Tashiro and Yokoyama2018; Crisostomi, Lewandowski, & Vernizzi Reference Crisostomi, Lewandowski and Vernizzi2020; Borisov & Jain Reference Borisov and Jain2009; Bernardeau & Brax Reference Bernardeau and Brax2011; Bartolo et al. Reference Bartolo, Bellini, Bertacca and Matarrese2013; Takushima, Terukina, & Yamamoto Reference Takushima, Terukina and Yamamoto2014; Bellini, Jimenez, & Verde Reference Bellini, Jimenez and Verde2015; Burrage, Dombrowski, & Saadeh Reference Burrage, Dombrowski and Saadeh2019; Lewandowski Reference Lewandowski2020; Borisov & Jain Reference Borisov and Jain2009; Gil-Marín et al. 2011; Hellwing et al. Reference Hellwing, Li, Frenk and Cole2013; Bose & Taruya Reference Bose and Taruya2018; Bose et al. Reference Bose, Byun, Lacasa, Moradinezhad Dizgah and Lombriser2020; Yamauchi et al. Reference Yamauchi, Yokoyama and Tashiro2017; Yamauchi & Sugiyama Reference Yamauchi and Sugiyama2022; Dinda Reference Dinda2018; Munshi et al. Reference Munshi2020; Namikawa, Bouchet, & Taruya Reference Namikawa, Bouchet and Taruya2018; Namikawa et al. Reference Namikawa, Bose, Bouchet, Takahashi and Taruya2019) or a few triangle configurations (Bose et al. Reference Bose, Byun, Lacasa, Moradinezhad Dizgah and Lombriser2020; Sugiyama et al. Reference Sugiyama2023a), Sugiyama et al. (Reference Sugiyama2023b) has incorporated
$l=2$
multipoles also in constraining the modified gravity. Here, for the first time, we consider all possible non-zero multipoles
$B_l^m$
across all possible unique triangle configurations to systematically explore their sensitivity to deviation from GR, with a particular focus on the DHOST theories. We investigate the distinctive signatures in
$B_l^m$
that indicate departures from GR and identify the triangle configurations most sensitive to the effects of modified gravity. The influence of tracer bias on
$B_l^m$
, which may entangle with the modified gravity signals, is also analysed. Specifically, we examine the impact of linear and quadratic bias parameters on the
$B_l^m$
and determine the configurations most affected. This comprehensive analysis seeks to identify robust observational signatures of modified gravity in
$B_l^m$
, complementing existing cosmological probes and offering new insights into the fundamental nature of gravity.
A brief outline of the paper is as follows. Section 2 presents the evolution of density and velocity fluctuations in the DHOST theory, including the modelling of perturbation growth. Section 3 presents the redshift space bispectrum of biased tracers and their multipole expansion. Section 4 presents the results, and Section 5 presents the Summary and Discussion. Results for some of the multipole moments, which are not included in the main text, have been presented in the Appendix A.
2. Density fluctuations in modified gravity theories
2.1 Matter field
The dark matter behaves as a pressure-less perfect fluid with no vorticity on large scales (Bernardeau et al. Reference Bernardeau, Colombi, Gaztanaga and Scoccimarro2002). The fluctuations in its density
${\unicode{x03B4}}(\boldsymbol{x},\eta)$
and velocity field
$\boldsymbol{v}(\boldsymbol{x},\eta)$
are described by fluid equations,


where a,
$H=\dot{a}/a=\frac{1}{a}\frac{da}{d\eta}$
, and
${\unicode{x03B8}}=\nabla.\boldsymbol{v}$
are the scale factor, Hubble parameter, and velocity divergence, respectively. The potential
$\phi(\boldsymbol{x},\eta)$
accounts for the influence of gravitational interactions under the assumption that matter and gravity are minimally coupled. We consider the DHOST theories (Langlois et al. Reference Langlois, Mancarella, Noui and Vernizzi2017; Langlois Reference Langlois2017), which is a broad class of modified gravity models. The gravitational potential satisfies the modified Poisson equation in the DHOST theories under the quasi-static approximation (e.g. Pace et al. Reference Pace2021; Hirano et al. Reference Hirano, Kobayashi, Tashiro and Yokoyama2018; Hirano et al. Reference Hirano, Kobayashi, Yamauchi and Yokoyama2020),

where
$\nu_0,\nu_1,\nu_2$
are time-dependent functions, and their relation to dark energy effective field theory (EFT) parameters is given in Yamauchi & Sugiyama (Reference Yamauchi and Sugiyama2022). The source term
$S_{\phi}^{\text{NL}}(\boldsymbol{x},\eta)$
arises from the equation of motion of the scalar field and encapsulates the contributions of non-linear interactions only to the evolution of potential perturbations. The set of these equations is solved perturbatively by expanding
${\unicode{x03B4}}, {\unicode{x03B8}}$
, and
$\phi$
fields in terms of successive orders in
${\unicode{x03B4}}_1$
(i.e.
$\{{\unicode{x03B4}},{\unicode{x03B8}},\phi\}\propto \sum_n\mathcal{O}({\unicode{x03B4}}_1^n)$
), where
${\unicode{x03B4}}_1$
represents the linear order fluctuations in the density (Goroff et al. Reference Goroff, Grinstein, Rey and Wise1986). The perturbative expansion of the source term
$S_{\phi}^{\text{NL}}$
in Equation (3) is required to solve this system of equations. Its leading contribution arises at the second order
$\mathcal{O}({\unicode{x03B4}}_1^2)$
and is given by (Sugiyama et al. Reference Sugiyama2023b),

where
$\tau_1,\tau_2$
are time-dependent functions related to the EFT parameters. The spatially dependent functions
$W_1$
and
$W_2$
are expressed in terms of linear density perturbations
${\unicode{x03B4}}_1$
as,

The time evolution of the density perturbations, obtained by combining Equations (1), (2) and eliminating
$\Phi$
using Equation (3), in the Fourier space is,

where
${\Delta}({\boldsymbol{k}},\eta)$
is the Fourier transform of
${\unicode{x03B4}}(\boldsymbol{x},\eta)$
, and parameters
$\sigma=(2\nu_2-\nu_1)/(1-\nu_2)$
and
$\Xi = 2\nu_0/[3\Omega_m(1-\nu_2)]$
are functions of time. This equation provides the evolution of density fluctuations for a particular DHOST gravity model specified by the set of parameters {
$\nu_0,\nu_1,\nu_2,\tau_1,\tau_2$
} (Sugiyama et al. Reference Sugiyama2023b). In the linear regime, the source term
$S_{\Delta}^{\text{NL}}$
vanishes, and we get the first order growing mode solution
$\Delta_1({\boldsymbol{k}},\eta)= D(\eta)\Delta_L({\boldsymbol{k}})$
, where the time-dependent function
$D(\eta)$
is the linear growth factor, and
$\Delta_L({\boldsymbol{k}})$
is the initial density fluctuations. It indicates that the perturbations grow independently for all
${\boldsymbol{k}}$
modes. The evolution equation (Equation (6)) at first order can be reinterpreted in terms of linear growth rate,
$f=\frac{d\ln D}{d\ln a}$
, as,

The linear perturbations in the velocity field, calculated using the continuity equation (Equation (1)), is
$\Theta_1({\boldsymbol{k}},\eta)=-aHfD\Delta_L({\boldsymbol{k}},\eta)$
, where
$\Theta({\boldsymbol{k}},\eta)$
is the Fourier conjugate of the velocity divergence
${\unicode{x03B8}}(\boldsymbol{x},\eta)$
. We proceed further to compute the second-order solutions, where the non-linear source term
$S_{\Delta}^{\text{NL}}$
(in Equation (6)) also contributes, as expressed by

where
${\boldsymbol{q}}={\boldsymbol{q}}_1+{\boldsymbol{q}}_2$
,
${\unicode{x03B4}}_D$
is the Dirac delta function,
$\sigma_1$
and
$\sigma_2$
are time-dependent functions (Hirano et al. Reference Hirano, Kobayashi, Yamauchi and Yokoyama2020; Sugiyama et al. Reference Sugiyama2023b),

and the coupling variables
$C^{{\unicode{x03BA}},\lambda}_{ij}$
,

The second-order perturbations in density (using Equations 6 and 8) and velocity field (using Equation 1) are given, respectively, by,


The second-order kernels
$F_2^{{}}$
and
$G_2^{{}}$
are (Hirano et al. Reference Hirano, Kobayashi, Yamauchi and Yokoyama2020; Takushima et al. Reference Takushima, Terukina and Yamamoto2014; Takushima, Terukina, & Yamamoto Reference Takushima, Terukina and Yamamoto2015; Crisostomi et al. Reference Crisostomi, Lewandowski and Vernizzi2020; Lewandowski Reference Lewandowski2020),

The evolution of the density field at the second-order can be reinterpreted in the form of the evolution of non-linear growth rates
${\unicode{x03BA}}_{\unicode{x03B4}}$
and
$\lambda_{\unicode{x03B4}}$
, using Equations (11) and (8) in Equation (6),


The non-linear growth rates of the velocity divergence field
${\unicode{x03BA}}_{\unicode{x03B8}}$
and
$\lambda_{\unicode{x03B8}}$
can be determined from
${\unicode{x03BA}}_{\unicode{x03B4}}$
and
$\lambda_{\unicode{x03B4}}$
, and the relations obtained using Equations (11) and (12) in the Equation (1) are (Hirano et al. Reference Hirano, Kobayashi, Yamauchi and Yokoyama2020; Sugiyama et al. Reference Sugiyama2023b),

2.2 Modelling growth of perturbations in DHOST theories
The Equations (7), (14), and (15) provide a complete description of first- and second-order evolution of density perturbations within the framework of the DHOST gravity theories. The specific model of gravity is characterised by the parameters
$\sigma$
,
$\Xi$
,
$\sigma_1$
, and
$\sigma_2$
. These parameters are, in turn, determined by the EFT parameters {
$\nu_0$
,
$\nu_1$
,
$\nu_2$
,
$\tau_1$
,
$\tau_2$
}, which encapsulate the modifications introduced by the DHOST model. In the standard
$\Lambda$
CDM model with GR, the parameters take values {
$\nu_0=\frac{3}{2}\Omega_m$
,
$\nu_1=0$
,
$\nu_2=0$
,
$\tau_1=0$
,
$\tau_2=0$
}, hence
$\sigma=0$
,
$\Xi=1$
,
$\sigma_1=2f^2+\frac{3}{2}\Omega_m$
, and
$\sigma_2=-f^2$
and the Equations (7), (14) and (15) lead to
$f=\Omega_m^{6/11}$
,
${\unicode{x03BA}}_{\unicode{x03B4}}=1$
and
$\lambda_{\unicode{x03B4}}=\Omega_m^{3/572}$
(Bouchet et al. Reference Bouchet, Colombi, Hivon and Juszkiewicz1995; Bernardeau et al. Reference Bernardeau, Colombi, Gaztanaga and Scoccimarro2002; Yamauchi et al. Reference Yamauchi, Yokoyama and Tashiro2017). The second-order growth rates of the velocity divergence field obtained using Equation (16) are
${\unicode{x03BA}}_{\unicode{x03B8}}=1$
and
$\lambda_{\unicode{x03B8}}=\Omega_m^{15/1\,144}$
. In the Horndeski gravity, the constraints are limited to {
$\nu_1=0$
,
$\nu_2=0$
,
$\tau_1=0$
} only, and the remaining parameters {
$\nu_0, \tau_2$
} are free to take any values, allowing deviation from the standard model of gravity. This results in
$\sigma=0$
,
$\sigma_1=2f^2+\frac{3}{2}\Omega_m$
and
$\{\Xi,\sigma_2\}$
depart from their standard
$\Lambda$
CDM values (Yamauchi et al. Reference Yamauchi, Yokoyama and Tashiro2017; Takushima et al. Reference Takushima, Terukina and Yamamoto2014; Sugiyama et al. Reference Sugiyama2023a). Consequently,
${\unicode{x03BA}}_{\unicode{x03B4}}={\unicode{x03BA}}_{\unicode{x03B8}}=1$
, while
$\lambda_{\unicode{x03B4}}$
(and
$\lambda_{\unicode{x03B8}}$
) captures the modification to gravity. In the general case of DHOST theories, all the growth rates
$f,{\unicode{x03BA}}_{\unicode{x03B4}},\lambda_{\unicode{x03B4}}$
(also
${\unicode{x03BA}}_{\unicode{x03B8}},\lambda_{\unicode{x03B8}}$
) diverge from those predicted by standard
$\Lambda$
CDM+GR (Crisostomi et al. Reference Crisostomi, Lewandowski and Vernizzi2020; Lewandowski Reference Lewandowski2020; Hirano et al. Reference Hirano, Kobayashi, Tashiro and Yokoyama2018). It suffices to model these growth rates to study the perturbations in the DHOST framework.
We adopt the widely used model of linear growth rate as (Wang & Steinhardt Reference Wang and Steinhardt1998),

where
$\unicode{x03B3}$
is the linear growth index and it takes value
$\unicode{x03B3}_\textrm{GR}=6/11$
for the standard
$\Lambda$
CDM + GR theory. It has been utilised in studying the growth of linear perturbations. Recent works by Yamauchi et al. (Reference Yamauchi, Yokoyama and Tashiro2017) and Yamauchi & Sugiyama (Reference Yamauchi and Sugiyama2022) have demonstrated that the second-order growth rates in the DHOST theories can similarly be expressed as a power of
$\Omega_m$
. Following their formulation, we proceed with the model,

where
$\unicode{x03C8}$
and
$\unicode{x03BE}$
are both the second-order growth indices and take values 0 and
$3/572$
, respectively, in the
$\Lambda$
CDM + GR case. Note that we take the background evolution of
$\Omega(a)$
assuming the standard GR, given by,

In principle, the time evolution of
$\Omega_m(a)$
in the DHOST theories differs from that in GR. Nevertheless, this deviation is suppressed by the factor
$(1-\Omega_\textrm{GR})$
(Yamauchi & Sugiyama Reference Yamauchi and Sugiyama2022).
2.3 Biased tracers in redshift space
The observables, such as galaxy distribution, Lyman-
$\alpha$
, 21-cm signal, etc, trace the dark matter with some bias. The density perturbations of these biased tracers
${\unicode{x03B4}}^{(b)}(\boldsymbol{x},\eta)$
are modelled using a Taylor series expansion of matter fluctuations
${\unicode{x03B4}}(\boldsymbol{x},\eta)$
, up to second order as follows (Desjacques et al. Reference Desjacques, Jeong and Schmidt2018),

where
$b_1$
and
$b_2$
are the local linear and quadratic bias parameters, respectively, and
$b_{s^2}$
is the second-order tidal bias. The source of tidal perturbations
$s_{ij}^2$
is expressed as,

The tracers are observed in redshift space. The observed redshift includes the contribution of peculiar velocity also in addition to the Hubble expansion, and the estimated distance
$\tilde{\boldsymbol{x}}$
is distorted along the LOS direction
$\hat{n}$
, expressed as
$\tilde{\boldsymbol{x}} = \boldsymbol{x}+\frac{\boldsymbol{v}.\, \hat{n}}{aH}\hat{n}$
(Kaiser Reference Kaiser1987). The
$n\textrm{th}$
order perturbations in the tracer density field in redshift space are given by (Scoccimarro et al. Reference Scoccimarro, Couchman and Frieman1999),

The first- and second-order redshift space kernels, assuming the peculiar velocity of the tracers follow the matter velocity
$\boldsymbol{v}^{(b)}=\boldsymbol{v}$
are,


Figure 1. Sensitivity of
$B_l^m$
to the linear growth index
$\unicode{x03B3}$
. It shows the variation of the ratio
$[\alpha_l^m]_{\unicode{x03B3}}=[B_l^m]_{\unicode{x03B3}}/[B_l^m]_\textrm{GR}$
, where
$[B_l^m]_\textrm{GR}$
denotes the bispectrum multipoles for GR and
$[B_l^m]_{\unicode{x03B3}}$
denotes the multipoles at different values of
$\unicode{x03B3}$
. The panels from left to right show the results for equilateral (
${\unicode{x03BC}}=0.5,\, t=1$
), stretched (
${\unicode{x03BC}}=1,\, t=0.5$
), and squeezed (
${\unicode{x03BC}}=1,\, t=1$
) triangle shapes, respectively, with fixed size
$k_1=0.1$
$\textrm{Mpc}^{-1}$
at
$z=0.61$
. The black solid line represents monopole (
$l=0,m=0$
), while different shades of green, blue, and orange-red correspond to
$ l=2, 4$
and 6 multipoles, respectively. Distinct line styles correspond to different m values.
3. Redshift space bispectrum of biased tracers
The bispectrum is defined as,

with the condition that
${\boldsymbol{k}}$
modes form a closed triangle,
${\boldsymbol{k}}_1+{\boldsymbol{k}}_2+{\boldsymbol{k}}_3=0$
. The leading order contribution is given by (using Equation (22)) (Scoccimarro et al. Reference Scoccimarro, Couchman and Frieman1999; Mazumdar et al. Reference Mazumdar, Bharadwaj and Sarkar2020),

where
$P_L(k)=\langle |\Delta_L(k)|^2\rangle$
, is the linear matter power spectrum. We have used the Boltzmann code CLASS (Blas, Lesgourgues, & Tram Reference Blas, Lesgourgues and Tram2011) to compute the linear power spectrum
$P_L(k)$
. The redshift space bispectrum is a function of five parameters that define the triangle’s shape, size, and orientation with respect to the LOS. We opt for the parametrisation prescribed by Bharadwaj et al. (Reference Bharadwaj, Mazumdar and Sarkar2020). For a triangle with
$k_1\geq k_2\geq k_3$
, its size is parameterised by the largest side
$k_1$
, shape by two dimensionless parameters
$({\unicode{x03BC}},t)$
,

and orientation with respect to the LOS using a unit vector
$\hat{p}$
. Assuming the fixed LOS along
$\hat{z}$
, the cosines of angles between the LOS and triangle sides are related to
$\hat{p}$
by,

where
${{\unicode{x03BC}}}_i=\hat{k}_i.\hat{z}$
. We use the Equation (23) in (25) and rewrite the expression for bispectrum as (Mazumdar et al. Reference Mazumdar, Bharadwaj and Sarkar2020),

Here, the term
$(2-3C^\lambda_{12})/3$
represents the scale-dependent tidal force (Equation (21)). We have used the compact notation for the kernels,
$G_{ij}^{{}}=G_2^{{}}({\boldsymbol{q}}_i,{\boldsymbol{q}}_j,\eta)$
, etc. and defined the function
$\Delta G^{{}}_{ij}$
as,
$\Delta G^{{}}_{ij} = F_{ij}^{{}} - G_{ij}^{{}}$
. For a specific model of gravity, the kernels depend on the configuration of triangle shapes via coupling parameters,
$C^{{\unicode{x03BA}},\lambda}_{ij}$
. These are expressed in terms of shape parameters
$({\unicode{x03BC}},t)$
as follows,

The orientation (
$\hat{p}$
) dependence of the redshift space bispectrum is fully quantified by the functions
$R(\hat{p}, b_1,f)$
,
$S_{ij}(\hat{p}, b_1,f)$
and
$T_{ij}(\hat{p}, b_1,f, {\unicode{x03BC}}, t)$
given by,



These functions encapsulate the linear growth factor f and linear bias
$b_1$
only. Sensitivity of the bispectrum to non-linear growth rates
$({\unicode{x03BA}}_{\unicode{x03B4}},\lambda_{\unicode{x03B4}})$
enters through the kernels
$G^{{}}_{ij}$
and
$\Delta G^{{}}_{ij}$
(Equation 13). To study the anisotropy of the bispectrum, we decompose it into spherical harmonic
$Y_l^m(\hat{{p}})$
basis and express the bispectrum multipole moments
${B}_l^m(k_1,{\unicode{x03BC}},t )$
as (Scoccimarro et al. Reference Scoccimarro, Couchman and Frieman1999; Bharadwaj et al. Reference Bharadwaj, Mazumdar and Sarkar2020; Mazumdar et al. Reference Mazumdar, Bharadwaj and Sarkar2020),

The multipole expansion of the functions R,
$S_{ij}$
, and
$T_{ij}$
, which arises by putting Equation (28) in (33), is defined in a similar manner. The expressions for their multipoles are given in Bharadwaj et al. (Reference Bharadwaj, Mazumdar and Sarkar2020) (see Equations 26–29, 31–33, A1–A7 therein) and Mazumdar et al. (Reference Mazumdar, Bharadwaj and Sarkar2020) (see Equations A1–A22, B1–B65 therein). Note that the definitions of R,
$S_{ij}$
, and
$T_{ij}$
adopted in this work differ slightly from those in their studies. Therefore, to ensure consistency, the expressions for R,
$S_{ij}$
, and
$T_{ij}$
in Bharadwaj et al. (Reference Bharadwaj, Mazumdar and Sarkar2020) and Mazumdar et al. (Reference Mazumdar, Bharadwaj and Sarkar2020) need to be multiplied by
$b_1^3$
,
$b_1^2$
, and
$b_1^4$
, respectively.
We compute the ratio of all non-zero bispectrum multipoles (Equation 33) for varying values of
$\unicode{x03B3},\unicode{x03C8},\unicode{x03BE}$
relative to those obtained under the assumption of standard GR (
$\unicode{x03B3}_\textrm{GR}=6/11, \unicode{x03C8}_\textrm{ GR}=0, \unicode{x03BE}_\textrm{GR}=3/572$
), considering all possible triangle configurations. We restrict the values of
$\unicode{x03B3},\unicode{x03C8},\unicode{x03BE}$
that satisfy the constraints given in (Sugiyama et al. Reference Sugiyama2023b). They have placed bounds on the modified growth parameters
$\unicode{x03BE}_f=\log_{\Omega_m}(f/{\unicode{x03BA}}_{\unicode{x03B4}})$
,
$\unicode{x03BE}_s=\log_{\Omega_m}({\unicode{x03BA}}_{\unicode{x03B8}}/{\unicode{x03BA}}_{\unicode{x03B4}})$
and
$\unicode{x03BE}_t=\log_{\Omega_m}(\lambda_{\unicode{x03B8}}/{\unicode{x03BA}}_{\unicode{x03B4}})$
using data from the BOSS survey, finding that
$-0.907\lt\unicode{x03BE}_f\lt2.447$
,
$-0.504\lt\unicode{x03BE}_s$
and
$-1.655\lt\unicode{x03BE}_t$
with 95% confidence. The results are discussed in the next Sections.
4. Results
The bispectrum
$B(k_1,{\unicode{x03BC}},t,\hat{p})$
(Equation 28) exhibits an approximate power-law dependence on
$k_1$
within the scales of interest. Consequently, the key features of the bispectrum multipole ratio between different gravitational theories remain consistent across
$k_1$
. Thus, showing the results for a single value of
$k_1$
is sufficient to capture the essential behaviour. Here, we present the results for
$k_1=0.1$
Mpc
$^{-1}$
and adopt redshift
$z=0.61$
. Unless stated otherwise, we use fiducial values for the bias parameters
$b_1=1.2$
,
$b_2=0$
, and
$b_{s^2}=0$
(Nandi et al. Reference Nandi2024). The following subsections present a detailed examination of the sensitivity of bispectrum multipoles to the linear and non-linear growth indices.
4.1 Sensitivity to linear growth
In this subsection, we investigate how the bispectrum multipoles
$B_l^m$
respond to modifications in the linear growth of perturbations. Figure 1 shows the ratio
$[\alpha_l^m]_{\unicode{x03B3}}=[B_l^m]_{\unicode{x03B3}}/[B_l^m]_\textrm{GR}$
, where
$[B_l^m]_\textrm{GR}$
denotes the bispectrum multipoles for GR and
$[B_l^m]_{\unicode{x03B3}}$
corresponds to those obtained by varying the linear growth index
$\unicode{x03B3}$
. The second-order indices are kept fixed at their GR values, i.e.
$\unicode{x03C8}=\unicode{x03C8}_\textrm{GR}$
and
$\unicode{x03BE}=\unicode{x03BE}_\textrm{GR}$
. The results are shown for three distinct triangle configurations
$-$
equilateral (
${\unicode{x03BC}}=0.5, t=1$
), stretched (
${\unicode{x03BC}}=1, t=0.5$
), and squeezed (
${\unicode{x03BC}}=1,$
$t=1$
)
$-$
shown from left to right across the panels. In the figure, the black solid curve shows the ratio of monopole moments (
$l=0$
), and the shades of green, blue, and orange-red respectively show the ratios for quadrupole (
$l=2$
), hexadecapole (
$l=4$
), and tetrahexacontapole (
$l=6$
), as indicated by the colour bar. Different line styles are used to distinguish between multipoles with different m values. The grey solid line shows the ratio of
$l=8$
multipole,
$[\alpha_8^0]_{\unicode{x03B3}}$
.
Considering the monopole ratio
$[\alpha_0^0]_{\unicode{x03B3}}$
, we see that it is relatively small, around 0.6, at low values of
$\unicode{x03B3}$
across all three triangle configurations. It increases steadily with
$\unicode{x03B3}$
, but remains
$\lesssim 1.3$
as
$\unicode{x03B3}\rightarrow 1$
. Similarly, the quadrupole (
$l=2$
) ratios
$- [\alpha_2^0]_{\unicode{x03B3}}$
,
$[\alpha_2^1]_{\unicode{x03B3}}$
, and
$[\alpha_2^2]_{\unicode{x03B3}} -$
also start at lower values for small
$\unicode{x03B3}$
and rises with
$\unicode{x03B3}$
across all triangle configurations, varying within the range of 0.7–1.2. We find that, for certain triangle shapes, all the quadrupole moments carry no independent information on the linear growth index. In equilateral and squeezed configurations,
$[\alpha_2^0]_{\unicode{x03B3}}$
and
$[\alpha_2^2]_{\unicode{x03B3}}$
coincide exactly, while in stretched configurations
$[\alpha_2^0]_{\unicode{x03B3}}$
and
$[\alpha_2^1]_{\unicode{x03B3}}$
are identical for all
$\unicode{x03B3}$
. Therefore, to fully exploit the potential of all quadrupole moments for probing the gravity, the analysis must include a broader variety of triangle shapes.
Next, we examine the behaviour of the hexadecapole (
$l=4$
) moments. In general, the ratios
$[\alpha_4^m]_{\unicode{x03B3}}$
for
$m = 0, 1, 2, 3, 4$
initially start at a lower value of around 0.8 and typically increase with
$\unicode{x03B3}$
, reaching a peak close to 1, and subsequently decrease as
$\unicode{x03B3}$
continues to grow across all triangle configurations. The transition points,
$\unicode{x03B3} = \unicode{x03B3}_t$
, where each multipole
$[\alpha_4^m]_{\unicode{x03B3}}$
shifts from rising to falling, vary with the configuration. The transition points
$\unicode{x03B3}_t$
for each m multipoles and triangle configuration are summarised in Table 1.
Table 1. Transition points
$\unicode{x03B3}_t$
for the hexadecapole (
$l=4$
) ratios
$[\alpha_4^m]_{\unicode{x03B3}}$
across different triangle configurations. For stretched triangles, no transition point exists for the
$m=2$
and
$m=3$
multipoles due to their monotonic decrease with
$\unicode{x03B3}$
.

While most hexadecapole moments exhibit a peak within the range
$0.35 \lesssim \unicode{x03B3}_t \lesssim 0.73$
, a notable exception arises in the
$[\alpha_{4}^{2}]_{\unicode{x03B3}}$
and
$[\alpha_4^3]_{\unicode{x03B3}}$
components for stretched triangles. In these cases, the ratios decrease monotonically with increasing
$\unicode{x03B3}$
, starting from approximately 2 as
$\unicode{x03B3} \to 0$
and steadily declining to around 0.25 as
$\unicode{x03B3} \to 1$
. This indicates that these multipoles are more sensitive to variations in
$\unicode{x03B3}$
. We also observe degeneracies among certain multipoles for specific triangle shapes. For equilateral triangles, we observe that
$[\alpha_4^4]_{\unicode{x03B3}} = [\alpha_4^2]_{\unicode{x03B3}} = [\alpha_4^0]_{\unicode{x03B3}}$
throughout the range of
$\unicode{x03B3}$
, while for stretched triangles,
$[\alpha_4^1]_{\unicode{x03B3}} = [\alpha_4^0]_{\unicode{x03B3}}$
consistently holds. In the case of squeezed triangles,
$[\alpha_4^1]_{\unicode{x03B3}}$
equals
$[\alpha_4^3]_{\unicode{x03B3}}$
when
$\unicode{x03B3} \geq \unicode{x03B3}_\textrm{GR}$
. As a result, the
$m=1, 2, 4$
hexadecapoles carry no additional information on the linear growth index beyond what is already present in the
$m=0$
moment for these cases.
We now turn to the behaviour of the tetrahexacontapole (
$l=6$
) moments. In general,
$[\alpha_6^m]_{\unicode{x03B3}}$
values are initially high at small
$\unicode{x03B3}$
and gradually decrease as
$\unicode{x03B3}$
increases. However, the
$[\alpha_6^1]_{\unicode{x03B3}}$
corresponding to equilateral triangles is an exception, which begins at a lower value and surges as
$\unicode{x03B3}$
becomes large. Additionally, the ratios
$[\alpha_6^0]_{\unicode{x03B3}}$
,
$[\alpha_6^1]_{\unicode{x03B3}}$
,
$[\alpha_6^2]_{\unicode{x03B3}}$
, and
$[\alpha_6^3]_{\unicode{x03B3}}$
also exhibit a distinct behaviour for the case of squeezed triangles. Their values are in the range 1–
$1.1$
at small
$\unicode{x03B3}$
, and they increase slightly to attain a peak value at around
$\unicode{x03B3}_t=0.17$
and then begin to decline gradually. The multipoles
$[\alpha_6^0]_{\unicode{x03B3}}$
,
$[\alpha_6^2]_{\unicode{x03B3}}$
and
$[\alpha_6^4]_{\unicode{x03B3}}$
for equilateral configuration are especially sensitive to the
$\unicode{x03B3}$
. They attain values in the range
$\approx 2-2.6$
at small
$\unicode{x03B3}\rightarrow 0$
and decline to
$0-0.4$
at
$\unicode{x03B3}=1$
. For the rest of the cases,
$[\alpha_6^m]_{\unicode{x03B3}}$
remains within the narrower range from
$0.5$
to
$1.4$
. In many cases, various
$[\alpha_6^m]_{\unicode{x03B3}}$
overlap with each other, and not all tetrahexacontapoles contain independent information on linear growth. For the case of stretched triangles,
$[\alpha_6^0]_{\unicode{x03B3}}=[\alpha_6^1]_{\unicode{x03B3}}$
and
$[\alpha_6^2]_{\unicode{x03B3}}\approx[\alpha_6^3]_{\unicode{x03B3}}\approx[\alpha_6^5]_{\unicode{x03B3}}$
. For squeezed triangles,
$[\alpha_6^0]_{\unicode{x03B3}}\approx[\alpha_6^1]_{\unicode{x03B3}}\approx[\alpha_6^2]_{\unicode{x03B3}}$
and
$[\alpha_6^4]_{\unicode{x03B3}}=[\alpha_6^5]_{\unicode{x03B3}}=[\alpha_6^6]_{\unicode{x03B3}}$
.
Finally, we consider the highest non-zero multipoles,
$l=8$
, where only terms with
$m\leq6$
are non-vanishing, while all higher-order multipoles are identically zero. We see that these
$l=8$
multipoles do not provide additional information about
$\unicode{x03B3}$
beyond what is already captured by the
$l=6$
multipoles. For all triangle configurations, each
$[\alpha_8^m]_{\unicode{x03B3}}$
multipole is precisely equal to
$[\alpha_6^6]_{\unicode{x03B3}}$
.
4.2 Sensitivity to second-order growth
Further, we analyse the sensitivity of bispectrum multipoles
$B_l^m(k_1,{\unicode{x03BC}},t)$
to variations in the second-order growth indices,
$\unicode{x03BE}$
and
$\unicode{x03C8}$
. The observed
$B_l^m$
is also influenced by the bias of tracers (Equation 28); therefore, quantifying the effects due to the bias is crucial to isolating the true signal indicative of deviations from GR. We study the variations of
$B_l^m$
with respect to both the second-order growth indices and bias parameters simultaneously.
To thoroughly assess the capability of
$B_l^m(k_1,{\unicode{x03BC}},t)$
in distinguishing second-order growth effects induced by modified gravity, we systematically explore all possible triangle configurations in our analysis. Here, the parameters
$({\unicode{x03BC}}, t)$
uniquely specify the shape of each triangle. Thus, to investigate the behaviour of
$B_l^m(k_1, {\unicode{x03BC}}, t)$
across the full range of triangle shapes, we present our results in the
${\unicode{x03BC}}$
-t parameter space while keeping the triangle size fixed at
$k_1 = 0.1\, \mathrm{Mpc}^{-1}$
. Here, we focus our discussion on the multipoles that exhibit significant indications of modified gravity, while the remaining multipoles are discussed in the 5. Notably, the monopole
$B_0^0$
, as shown in Figure 9, exhibits minimal variation with changes in gravity at the second order, and its corresponding results are provided in the Appendix A. The results for the higher-order multipoles are discussed in detail below.
4.2.1 Quadrupole (
$l=2$
) moments
Figure 2 shows the results for
$l=2,m=1$
multipole. In the top row, we plot the ratio
$[\alpha_l^m]_{\unicode{x03BE}}=[B_l^m]_{\unicode{x03BE}}/[B_l^m]_\textrm{GR}$
as a function of
$({\unicode{x03BC}},t)$
, with the
$k_1=0.1$
Mpc
$^{-1}$
fixed. Here,
$[B_l^m]_{\unicode{x03BE}}$
denotes the bispectrum multipoles evaluated by varying
$\unicode{x03BE}$
, while the other growth indices are held fixed at their GR values, i.e.
$\unicode{x03B3}=\unicode{x03B3}_\textrm{GR}$
and
$\unicode{x03C8}=\unicode{x03C8}_\textrm{GR}$
. Subsequent rows (second to fifth) present the ratios
$[\alpha_l^m]_{\unicode{x03C8}}$
,
$[\alpha_l^m]_{b_1}$
,
$[\alpha_l^m]_{b_2}$
and
$[\alpha_l^m]_{b_{s^2}}$
, respectively, each defined analogously to
$[\alpha_l^m]_{\unicode{x03BE}}$
. The results for six distinct values of these parameters (
$\unicode{x03BE},\unicode{x03C8},b_1,b_2,b_{s^2}$
) are presented in each column.

Figure 2. Sensitivity of
$B_2^1$
to variations in the second-order growth indices
$\unicode{x03BE}$
and
$\unicode{x03C8}$
. The first and second rows show the ratio
$[\alpha_2^1]_{\unicode{x03BE}}=[B_2^1]_{\unicode{x03BE}}/[B_2^1]_\textrm{GR}$
and
$[\alpha_2^1]_{\unicode{x03C8}}=[B_2^1]_{\unicode{x03C8}}/[B_2^1]_\textrm{GR}$
, respectively, across all possible unique triangle configurations. Here
$[B_2^1]_{\unicode{x03BE}}$
and
$[B_2^1]_{\unicode{x03C8}}$
are the bispectrum multipoles at different values of
$\unicode{x03BE}$
and
$\unicode{x03C8}$
, respectively, while
$[B_2^1]_\textrm{GR}$
corresponds to the values considering the GR. Different columns correspond to specific values of
$\unicode{x03BE}$
and
$\unicode{x03C8}$
, as specified in the respective panels. Each subplot illustrates results in
${\unicode{x03BC}}-t$
space, where each point corresponds to a distinct triangle shape (refer to Figure 2 of Bharadwaj et al. Reference Bharadwaj, Mazumdar and Sarkar2020), with triangle size fixed at
$k_1=0.1$
$\textrm{Mpc}^{-1}$
fixed. The third, fourth, and fifth rows show the
$[\alpha_2^1]_{b_1}$
,
$[\alpha_2^1]_{b_2}$
, and
$[\alpha_2^1]_{b_{s^2}}$
, respectively, depicting the impact of linear and quadratic bias parameters on the bispectrum multipoles. The solid black contours highlight the locations where
$[B_2^1]_\textrm{GR}$
has a negative sign.
Each subplot in Figure 2 shows results across the entire allowed range of
${\unicode{x03BC}}$
and t values, bounded by
$0.5 \leq {\unicode{x03BC}}, t \leq 1$
and
$2{\unicode{x03BC}} t \geq 1$
. This region covers triangles of all possible shapes, with each point
$({\unicode{x03BC}}, t)$
representing a unique triangle configuration. For a detailed discussion on the locations of various triangle shapes in the
${\unicode{x03BC}}-t$
plane, refer to Figure 2 in Bharadwaj et al. (Reference Bharadwaj, Mazumdar and Sarkar2020). Briefly, the upper left corner (
${\unicode{x03BC}}=0.5, t=1$
) represents the equilateral triangles. The upper and lower boundaries, marked by
$t=1$
and
$2{\unicode{x03BC}} t=1$
, respectively, depict L isosceles triangles (the two largest sides are equal;
$t=1$
) and S isosceles triangles (the two smallest sides are equal;
$t=0.5/{\unicode{x03BC}}$
). The right boundary (
${\unicode{x03BC}}=1$
) represents linear triangles, where the three sides are nearly collinear. Furthermore, the upper right (
${\unicode{x03BC}}=1,t=1$
) and lower right (
${\unicode{x03BC}}=1,t=0.5$
) corners correspond to the squeezed and stretched triangles. The relation
${\unicode{x03BC}}=t$
corresponds to the right-angle triangles, which separates the regions with acute triangles (
${\unicode{x03BC}} \lt t$
) from those of obtuse triangles (
${\unicode{x03BC}} \gt t$
).
Considering the first row of Figure 2, we observe that the magnitude of the ratio
$|[\alpha_2^1]_{\unicode{x03BE}}|$
remains below 8 across most of the
${\unicode{x03BC}}-t$
plane for all values of
$\unicode{x03BE}$
. Generally, the ratio stays close to unity at the stretched limit (
${\unicode{x03BC}}=1,t=0.5$
), and it increases as triangles are deformed from obtuse (
${\unicode{x03BC}}\gt t$
) to acute (
${\unicode{x03BC}} \lt t$
). The peak value lies near the right-angle triangles. As the triangle shapes are further deformed toward the equilateral limit,
$|[\alpha_2^1]_{\unicode{x03BE}}|$
begins to decline and approaches a value close to unity. In particular, for the small
$\unicode{x03BE}\in [-0.05,0.05]$
, the ratio
$|[\alpha_2^1]_{\unicode{x03BE}}|$
remains in the order of unity for all triangles, indicating that the minor variations in
$\unicode{x03BE}$
have little effect on the
$B_2^1$
multipole. When
$\unicode{x03BE}=\pm 0.25$
, a few acute triangles near right-angle limit see the ratio rise to around 4. At the extreme values
$\unicode{x03BE}=\pm0.75$
, the
$|[\alpha_{2}^{1}]_{\unicode{x03BE}}|$
for these acute configurations reaches between 4 and 12. Thus, the
$B_2^1$
multipole exhibits the highest sensitivity to
$\unicode{x03BE}$
for these acute triangle configurations, suggesting that such triangles provide an optimal window to probe
$\unicode{x03BE}$
using
$B_2^1$
.
Considering next the
$[\alpha_2^1]_{\unicode{x03C8}}$
shown in the second row, we see that the results for
$\unicode{x03C8}\lt0$
closely resemble those for
$\unicode{x03BE}\gt0$
and vice versa. The overall magnitude of ratio
$[\alpha_2^1]_{\unicode{x03C8}}$
is greater than that of
$[\alpha_2^1]_{\unicode{x03BE}}$
for the identical values of the parameters
$\unicode{x03C8}$
and
$\unicode{x03BE}$
. For most cases, the
$|[\alpha_2^1]_{\unicode{x03C8}}|$
remains below 12. However, for
$\unicode{x03C8}=-0.75$
, the
$|[\alpha_2^1]_{\unicode{x03C8}}|$
has values in the range
$16-20$
for a few acute triangles, and the values decline towards the equilateral limit to
$\sim 2-4$
. The
$[\alpha_2^1]_{\unicode{x03C8}=-0.75}\sim 0-2$
for all obtuse triangles. The
$|[\alpha_2^1]_{\unicode{x03C8}}|$
declines as
$\unicode{x03C8}$
is increased towards
$\unicode{x03C8}_\textrm{GR}$
. At
$\unicode{x03C8}=\pm 0.05$
,
$[\alpha_2^1]_{\unicode{x03C8}}$
remains between
$0-2$
for all triangles. As
$\unicode{x03C8}$
is increased further to
$0.75$
, the
$[\alpha_2^1]_{\unicode{x03C8}}$
begin to surge and reached
$\sim 20$
near right-angle limit.
In addition to the variation in magnitude, the sign of
$B_2^1$
multipole also serves as a crucial indicator of deviation from GR. We find that the ratios
$[\alpha_2^1]_{\unicode{x03BE}}$
and
$[\alpha_2^1]_{\unicode{x03C8}}$
becomes negative for certain acute triangle configurations, implying that
$[B_2^1]_{\unicode{x03BE}}$
and
$[B_2^1]_{\unicode{x03C8}}$
have the opposite signs to that of
$[B_2^1]_\textrm{GR}$
. This observation suggests that if galaxy surveys measure the
$B_2^1$
with a sign contrary to GR predictions, it provides clear evidence of a departure from GR. Moreover, the sign of
$B_2^1$
not only signals a deviation from GR but also distinguishes between the scenarios of
$\unicode{x03BE}\lt0$
(or
$\unicode{x03C8}\gt0$
) and
$\unicode{x03BE}\gt0$
(or
$\unicode{x03C8}\lt0$
). The specific triangle configuration for which the sign of
$B_2^1$
flips provides key insight into the sign of second-order indices
$\unicode{x03BE}$
and
$\unicode{x03C8}$
.
Here, we have identified the triangle shapes that exhibit this sign-flip signature for both
$\unicode{x03BE}\lt0$
(equivalently
$\unicode{x03C8}\gt0$
) and
$\unicode{x03BE}\gt0$
(equivalently
$\unicode{x03C8}\lt0$
). In the plots, we use solid black contours to demarcate the regions where the GR predicts a negative value for the
$B_2^1$
multipole, i.e. where
$[B_2^1]_\textrm{GR}\lt0$
. The contours span all obtuse triangle configurations (
${\unicode{x03BC}} \gt t$
) and extend slightly into the acute regime near the right-angle limit (
${\unicode{x03BC}} \approx t$
). For
$\unicode{x03BE} \lt 0$
(or
$\unicode{x03C8}\gt0$
), we observe that
$[\alpha_2^1]_{\unicode{x03BE}}$
(or
$[\alpha_2^1]_{\unicode{x03C8}}$
) becomes negative for acute triangles close to the equilateral limit. This indicates that
$[B_2^1]_{\unicode{x03BE}}$
(or
$[B_2^1]_{\unicode{x03C8}}$
) is negative for these triangle configurations where GR predicts a positive value
$[B_2^1]_\textrm{GR}\gt0$
, providing a signature of modification to GR with
$\unicode{x03BE} \lt 0$
or
$\unicode{x03C8}\gt0$
. Conversely, for
$\unicode{x03BE} \gt 0$
(or
$\unicode{x03C8}\lt0$
), the ratio
$[\alpha_2^1]_{\unicode{x03BE}}$
(or
$[\alpha_2^1]_{\unicode{x03C8}}$
) becomes negative near the right-angle limit, where
$[B_2^1]_\textrm{GR}$
is negative. This implies that
$[B_2^1]_{\unicode{x03BE}}$
(or
$[B_2^1]_{\unicode{x03C8}}$
) has flipped to a positive value, offering a distinct indicator of deviation from GR associated with
$\unicode{x03BE} \gt 0$
or
$\unicode{x03C8}\lt0$
. Notably, the number of triangle configurations exhibiting this sign flip increases with the strength of the deviation, i.e. with larger values of
$|\unicode{x03BE}|$
and
$|\unicode{x03C8}|$
. Specifically, for
$\unicode{x03C8}=0.75$
, the multipole
$[B_2^1]_{\unicode{x03C8}}$
is negative across all triangle configurations, in contrast to the GR prediction where
$[B_2^1]_\textrm{ GR}$
remains positive for many configurations near the equilateral limit. In a nutshell, the
$B_2^1$
exhibits significant variations with
$\unicode{x03BE}$
and
$\unicode{x03C8}$
across most triangle configurations, and its negative sign for equilateral triangles and positive values near right-angle triangles serves as a remarkable signature of modified gravity.
The bias parameters can also affect the
$B_2^1$
and may potentially produce false signals of modified gravity if the bias of tracers is not accurately known a priori. Further, we analyse how linear bias
$b_1$
and quadratic bias parameters
$b_2, b_{s^2}$
affect the
$B_2^1$
. Third row in Figure 2 presents ratio
$[\alpha_2^1]_{b_1}$
. We see that for small values of
$b_1$
in the range 0.5 to 1, the
$[\alpha_2^1]_{b_1}$
is less than 1 for all triangles. For large values where
$b_1=1.5$
and
$1.75$
, the
$[\alpha_2^1]_{b_1}$
remains within the range
$1-2$
and
$2-4$
respectively. The
$[\alpha_2^1]_{b_1=2}$
is also within the 2–4 range for most triangles, except for a few acute triangles where it increases to approximately
$\sim 4-8$
. Considering the fourth and fifth rows, we see that the
$|[\alpha_2^1]_{b_2}|$
and
$|[\alpha_2^1]_{b_{s^2}}|$
remain below 2 for most of the triangles. The ratio
$|[\alpha_2^1]_{b_{s^2}}|$
becomes negative for a few acute triangle configurations near the right-angle limit, suggesting that a large offset in the parameter
${b_{s^2}}$
can potentially mimic a false signal of modified gravity. However, within the range of all bias parameters considered in this analysis, we do not observe any sign flip in
$B_2^1$
near the equilateral limit – a behaviour that is characteristic of modified gravity.
This concludes our discussion of the
$B_2^1$
multipole. Results for the other quadrupole moments –
$B_2^0$
and
$B_2^2$
– are presented in Appendix A.
4.2.2 Hexadecapole (
$l=4$
) moments
Next, we discuss the hexadecapole (
$l=4$
) moments. The results for the
$l = 4, m = 3$
multipole are shown in Figure 3, while the corresponding results for the other
$l = 4$
multipoles are provided in Appendix A. Considering the first row, we see that the trend of the ratio
$[\alpha_4^3]_{\unicode{x03BE}}$
in the
${\unicode{x03BC}}-t$
plane is somewhat similar to the
$[\alpha_2^1]_{\unicode{x03BE}}$
case. The notable difference is that the magnitude of the ratio
$|[\alpha_4^3]_{\unicode{x03BE}}|$
is largest near the right-angle configurations (
${\unicode{x03BC}} = t$
), reaching values as high as
$\sim 20$
for
$\unicode{x03BE} = \pm 0.75$
. This amplitude gradually decreases as the triangle shape deviates from the right-angle limit, becoming either more acute (
${\unicode{x03BC}} \lt t$
) or more obtuse (
${\unicode{x03BC}} \gt t$
). Nonetheless,
$|[\alpha_4^3]_{\unicode{x03BE}}|$
exhibits a slight increase near the most obtuse configurations – close to the stretched limit – with values rising to approximately
$\sim 1$
–2 for
$\unicode{x03BE} = 0.75$
and
$\sim 2$
–4 for
$\unicode{x03BE} = -0.75$
. For smaller values of
$\unicode{x03BE}$
, the peak values near the right-angle configuration are reduced, with
$|[\alpha_4^3]_{\unicode{x03BE} = \pm 0.25}| \sim 10$
and
$|[\alpha_4^3]_{\unicode{x03BE} = \pm 0.05}| \sim 3$
.

Figure 3. Sensitivity of
$B_4^3$
to variations in the second-order growth indices
$\unicode{x03BE}$
and
$\unicode{x03C8}$
. The rows, from top to bottom, correspond to
$[\alpha_4^3]_{\unicode{x03BE}}$
,
$[\alpha_4^3]_{\unicode{x03C8}}$
,
$[\alpha_4^3]_{b_1}$
,
$[\alpha_4^3]_{b_2}$
, and
$[\alpha_4^3]_{b_{s^2}}$
, respectively. The panels are identical to those in Figure 2.
We now consider
$[\alpha_4^3]_{\unicode{x03C8}}$
, shown in the second row. Much like the quadrupole case, the behaviour observed in the
${\unicode{x03BC}}-t$
space reflects a pattern similar to that when
$\unicode{x03C8} \rightarrow -\unicode{x03BE}$
. The overall trend of
$[\alpha_4^3]_{\unicode{x03C8}}$
closely resembles that of
$[\alpha_4^3]_{\unicode{x03BE}}$
, with the magnitude
$|[\alpha_4^3]_{\unicode{x03C8}}|$
showing a prominent peak near right-angle configurations. However, the magnitude of
$|[\alpha_4^3]_{\unicode{x03C8}}|$
is generally larger than that of
$|[\alpha_4^3]_{\unicode{x03BE}}|$
across the shape space. At the right-angle limit, the peak values of
$|[\alpha_4^3]_{\unicode{x03C8}}|$
reach approximately 20, 16, and 12 for
$\unicode{x03C8} = \pm 0.75$
,
$\pm 0.25$
, and
$\pm 0.05$
, respectively. Beyond this pronounced peak,
$|[\alpha_4^3]_{\unicode{x03C8}}|$
exhibits a secondary rise in amplitude near the stretched limit (
${\unicode{x03BC}} \rightarrow 1$
,
$t \rightarrow 0.5$
), particularly notable for higher values of
$\unicode{x03C8} = \pm 0.75$
, where
$|[\alpha_4^3]_{\unicode{x03C8}}|$
reaches up to 16.
Similar to the quadrupole case, the
$B_4^3$
multipole also exhibits sign reversals across a range of triangle configurations. Notably, this effect is more prominent here, with a greater number of triangles exhibiting sign flips near the right-angle (
${\unicode{x03BC}} = t$
) and stretched limit (
${\unicode{x03BC}} \rightarrow 1$
,
$t \rightarrow 0.5$
). Under GR, we see that the
$[B_4^3]_\textrm{GR}$
is positive for acute triangles and negative for obtuse triangles, except for a few triangles near the stretched limit. However,
$[B_4^3]_{\unicode{x03BE}}$
for
$\unicode{x03BE} \lt 0$
and
$[B_4^3]_{\unicode{x03C8}}$
for
$\unicode{x03C8} \gt 0$
become negative for many acute triangles – regions where GR predicts a positive value. Conversely,
$[B_4^3]_{\unicode{x03BE}}$
for
$\unicode{x03BE} \gt 0$
and
$[B_4^3]_{\unicode{x03C8}}$
for
$\unicode{x03C8} \lt 0$
, turn positive for many obtuse triangles where GR predicts a negative value. This clear contrast in the sign of
$B_4^3$
between GR and the DHOST scenario serves as a strong indicator of modifications in gravity.

Figure 4. Sensitivity of
$B_6^0$
to variations in the second-order growth indices
$\unicode{x03BE}$
and
$\unicode{x03C8}$
. The top, middle, and bottom rows show
$[\alpha_6^0]_{\unicode{x03BE}}$
,
$[\alpha_6^0]_{\unicode{x03C8}}$
, and
$[\alpha_6^0]_{b_1}$
, respectively. The panels are the same as in Figure 2, except for the last two columns, as
$B_6^0$
is unaffected by quadratic bias and second-order tidal bias parameters.
Nonetheless, the sign flip may also be induced by the galaxy bias parameters, as discussed earlier. The third, fourth, and fifth rows in the Figure 3 shows
$[\alpha_4^3]_{b_1}$
,
$[\alpha_4^3]_{b_2}$
and
$[\alpha_4^3]_{b_{s^2}}$
, respectively. We see that the impact of considering different
$b_1$
,
$b_2$
, and
$b_{s^2}$
parameters is severe in the case of
$l=4,m=3$
moment. The
$[\alpha_4^3]_{b_1}$
is small for acute triangles near the equilateral limit, where its magnitude is
$\lt 1$
if
$b_1\lt1.2$
, and it is around
$1-4$
if
$b_1\gt1.2$
. However, the
$[\alpha_4^3]_{b_1}$
is highest (
$\sim 16$
for
$b_1\lt1.2$
and
$\sim 20$
for
$b_1\gt1.2$
) near the right-angle and stretched limit. Furthermore, we observe that
$[B_4^3]_{b_1}$
changes sign for many obtuse triangles when
$b_1 \lt 1.2$
and for a few acute and stretched triangles when
$b_1 \gt 1.2$
. This behaviour closely resembles the sign flip of
$[\alpha_4^3]_{\unicode{x03C8}}$
. Similarly, the influence of second-order bias parameters is most pronounced near the right-angle and stretched configurations, while it is relatively small (
$|[\alpha_4^3]_{b_2}| \approx |[\alpha_4^3]_{b_{s^2}}| \sim 1-2$
) near equilateral and some obtuse triangle configurations. Notably, the trend of
$[\alpha_4^3]_{b_2}$
closely follows that of
$[\alpha_4^3]_{\unicode{x03BE}}$
. In the case of
$b_{s^2}$
, the
$[\alpha_4^3]_{b_{s^2}}$
changes sign for a few obtuse triangles when
$b_{s^2} \lt 0$
and for a few acute triangles when
$b_{s^2} \gt 0$
, similar to
$\unicode{x03C8}$
. However, it does not fully mimic
$\unicode{x03C8}$
, as it does not change the sign above the lower boundary of the black curve for
$b_{s^2} \lt 0$
or near the stretched limit for
$b_{s^2} \gt 0$
. Instead, for
$b_{s^2} \lt 0$
, it changes the sign near the stretched limit, and for
$b_{s^2} \gt 0$
, it changes the sign above the lower boundary of the black curve, which are the features of the
$\unicode{x03BE}$
parameter. Thus,
$b_{s^2}$
exhibits characteristics of both the
$\unicode{x03BE}$
and
$\unicode{x03C8}$
parameters. Consequently, these bias parameters must be accurately accounted for to utilise
$B_4^3$
as a probe for the DHOST theories effectively.
4.2.3 Tetrahexacontapole (
$l=6$
) moments
We now turn to the analysis of the
$l=6$
multipoles. Beginning with the case of
$l=6, m=0$
, illustrated in Figure 4, we find that the
$B_6^0$
multipole exhibits high sensitivity to variations in the parameters
$\unicode{x03BE}$
and
$\unicode{x03C8}$
. This sensitivity is most pronounced for acute configurations
$({\unicode{x03BC}} \lt t)$
, particularly near the equilateral limit, where
$B_6^0$
demonstrates substantial changes even for small
$\unicode{x03BE}$
and
$\unicode{x03C8}$
. Specifically, the maximum value of
$|[\alpha_6^0]_{\unicode{x03BE}}|$
reaches approximately
$\sim 8$
for
$\unicode{x03BE} = \pm 0.05$
, while
$|[\alpha_6^0]_{\unicode{x03C8}}|$
rises steeply to attain values as high as
$\sim 12$
for
$\unicode{x03C8} = \pm 0.05$
near the equilateral limit. However, for other triangle configurations, the magnitudes of the ratios remain relatively small:
$[\alpha_6^{0}]_{\unicode{x03BE} = -0.05} \approx [\alpha_6^{0}]_{\unicode{x03C8} = 0.05} \sim 1$
and
$[\alpha_6^{0}]_{\unicode{x03BE} = 0.05} \approx [\alpha_6^{0}]_{\unicode{x03C8} = -0.05} \sim 2$
. The values of
$|[\alpha_6^0]_{\unicode{x03BE}}|$
and
$|[\alpha_6^0]_{\unicode{x03C8}}|$
increase gradually with higher values of
$\unicode{x03BE}$
and
$\unicode{x03C8}$
, respectively. For equilateral triangles, both
$|[\alpha_6^0]_{\unicode{x03BE}}|$
and
$|[\alpha_6^0]_{\unicode{x03C8}}|$
exceed 20 when
$\unicode{x03BE}$
and
$\unicode{x03C8}$
attain values
$\pm 0.75$
. The number of triangle configurations where both
$|[\alpha_6^0]_{\unicode{x03BE}}|$
and
$|[\alpha_6^0]_{\unicode{x03C8}}|$
exceed 1 increases significantly in the acute triangle region as
$\unicode{x03BE}$
and
$\unicode{x03C8}$
upsurge. In contrast, for the obtuse triangles,
$[\alpha_6^0]_{\unicode{x03BE}}$
and
$[\alpha_6^0]_{\unicode{x03C8}}$
remains largely below 2.

Figure 5. Sensitivity of
$B_6^1$
to variations in the second-order growth indices
$\unicode{x03BE}$
and
$\unicode{x03C8}$
. The top, middle, and bottom rows correspond to
$[\alpha_6^1]_{\unicode{x03BE}}$
,
$[\alpha_6^1]_{\unicode{x03C8}}$
, and
$[\alpha_6^1]_{b_1}$
, respectively. The panels are identical to those in Figure 4.

Figure 6. Sensitivity of
$B_6^2$
to variations in the second-order growth indices
$\unicode{x03BE}$
and
$\unicode{x03C8}$
. The top, middle, and bottom rows correspond to
$[\alpha_6^2]_{\unicode{x03BE}}$
,
$[\alpha_6^2]_{\unicode{x03C8}}$
, and
$[\alpha_6^2]_{b_1}$
, respectively. The panels are identical to those in Figure 4.

Figure 7. Sensitivity of
$B_6^3$
to variations in the second-order growth indices
$\unicode{x03BE}$
and
$\unicode{x03C8}$
. The top, middle, and bottom rows correspond to
$[\alpha_6^3]_{\unicode{x03BE}}$
,
$[\alpha_6^3]_{\unicode{x03C8}}$
, and
$[\alpha_6^3]_{b_1}$
, respectively. The panels are identical to those in Figure 4.

Figure 8. Sensitivity of
$B_6^4$
to variations in the second-order growth indices
$\unicode{x03BE}$
and
$\unicode{x03C8}$
. The top, middle, and bottom rows correspond to
$[\alpha_6^4]_{\unicode{x03BE}}$
,
$[\alpha_6^4]_{\unicode{x03C8}}$
, and
$[\alpha_6^4]_{b_1}$
, respectively. The panels are identical to those in Figure 4.
Similar to the cases of both the quadrupole and hexadecapole moments, the sign of
$B_6^0$
also serves as a valuable indicator of deviations from GR. The
$[B_6^0]_{\text{GR}}$
is negative for a few triangles near the equilateral and squeezed limit. We see that the
$[\alpha_6^0]_{\unicode{x03BE}\lt0}$
and
$[\alpha_6^0]_{\unicode{x03C8}\gt0}$
are negative for many triangles in the region below the black contour near the equilateral limit, while
$[\alpha_6^0]_{\unicode{x03BE}\gt0}$
and
$[\alpha_6^0]_{\unicode{x03C8}\lt0}$
are negative for equilateral triangles. This implies that a measured positive
$B_6^0$
for equilateral triangles or a negative
$B_6^0$
for acute triangles (below black contour) serves as a definitive indicator of deviations from GR.
As noted earlier for the
$l=2$
and
$l=4$
multipoles, the apparent signal of modified gravity in
$B_6^0$
can be induced by the bias parameters. However, the
$l=6$
multipoles offer a distinct advantage over the
$l=2$
and
$l=4$
multipoles as they are unaffected by the quadratic bias
$b_2$
and second-order tidal bias
$b_{s^2}$
parameters. Considering Equation (28), we note that the terms involving
$b_2$
and
$b_{s^2}$
are associated with the
$S_{12}$
term, which is a polynomial in
${\unicode{x03BC}}_i^2$
with the highest power term being the
${\unicode{x03BC}}_i^4$
. So, bispectrum multipoles with
$l\leq 4$
are affected by the higher-order bias terms. Only the linear bias
$b_1$
influences the
$l=6$
multipoles. The third row of the Figure 4 shows
$[\alpha_6^0]_{b_1}$
. We see that, for
$b_1\lt1.2$
, the
$|[\alpha_6^0]_{b_1}|\lt1$
for most of the triangle configurations except at the equilateral limit where it has value
$1-2$
. For
$b_1\gt1.2$
,
$|[\alpha_6^0]_{b_1}|\sim 1-2$
. We see that the small values of
$b_1$
can induce the negative sign in the
$B_6^0$
below the equilateral limit. Overall, the influence of bias parameters on
$B_6^0$
is relatively minor and does not pose a significant concern.
Figure 5 shows the results for
$l=6,m=1$
multipole. We see that the results are very similar to the
$l=6,m=0$
multipole. The
$|[\alpha_6^{1}]_{\unicode{x03BE}}|$
and
$|[\alpha_6^{1}]_{\unicode{x03C8}}|$
are large for acute triangles, especially the equilateral limit, which is sensitive to very small values of
$\unicode{x03BE}$
and
$\unicode{x03C8}$
. In contrast to
$m=0$
case, the
$[B_6^1]_{\text{GR}}$
is negative for all triangle configurations. We see that the
$[\alpha_6^{1}]_{\unicode{x03BE}\gt0}$
and
$[\alpha_6^{1}]_{\unicode{x03C8}\lt0}$
are positive for all triangles, and
$[\alpha_6^{1}]_{\unicode{x03BE}\lt0}$
and
$[\alpha_6^{1}]_{\unicode{x03C8}\gt0}$
are negative for many acute triangles. Thus, the positive values of the observed
$B_6^1$
for acute triangles serve as indicators of modifications in GR. The third row shows that the linear bias parameter has a small impact on
$B_6^1$
and
$[\alpha_6^1]_{b_1}\lt 1$
for all triangles when
$b_1=0.5-1.2$
and
$[\alpha_6^1]_{b_1}\sim 1-2$
for most of the triangles (with the exception of equilateral triangles, where
$[\alpha_6^1]_{b_1}\sim 2-4$
) when
$b_1=1.2-2$
. Similar to the
$m=0$
case, for small values of
$b_1$
,
$[\alpha_6^1]_{b_1}$
is negative near the equilateral limit.
Figure 6 shows that the
$l=6,m=2$
multipole exhibits significant variation with
$\unicode{x03BE}$
and
$\unicode{x03C8}$
. The
$[\alpha_6^{2}]_{\unicode{x03BE}=-0.75}$
and
$[\alpha_6^{2}]_{\unicode{x03C8}=0.75}$
exceed 20 for equilateral triangles. As triangles deformed away from the equilateral triangles,
$[\alpha_6^{2}]_{\unicode{x03BE}=-0.75}$
and
$[\alpha_6^{2}]_{\unicode{x03C8}=0.75}$
becomes
$\sim -20$
. The magnitude of both decreases as we move from acute
$({\unicode{x03BC}} \lt t)$
to obtuse
$({\unicode{x03BC}}\gt t)$
triangles. In the region bounded by black contour, the
$[B_6^2]_\textrm{GR}$
is negative, however,
$[B_6^2]_{\unicode{x03BE}=-0.75}$
and
$[B_6^2]_{\unicode{x03C8}=0.75}$
are positive for all the acute triangles, which indicate the deviation from GR. As
$\unicode{x03BE}$
increases towards
$-0.05$
(or
$\unicode{x03C8}$
is decreased to 0.05), the overall
$[\alpha_6^{2}]_{\unicode{x03BE}}$
(or
$[\alpha_6^{2}]_{\unicode{x03C8}}$
) declines, however its magnitude remains between
$2-8$
at equilateral limit. For positive values of
$\unicode{x03BE}$
(or negative values of
$\unicode{x03C8}$
),
$[\alpha_6^{2}]_{\unicode{x03BE}}$
(or
$[\alpha_6^{2}]_{\unicode{x03C8}}$
) becomes negative for the many triangles near the equilateral limit, and its magnitude remains similar to the case of negative
$\unicode{x03BE}$
values. Considering the third row, we observe that the
$[\alpha_6^2]_{b_1}\lt2$
across all cases in the range
$b_1=0.5-2$
. An additional significant advantage of the
$B_6^2$
multipole is that the sign of
$[\alpha_6^2]_{b_1}$
remains unaffected with
$b_1$
in the region where
$\unicode{x03BE}$
and
$\unicode{x03C8}$
induce the negative sign. Thus, the negative
$B_6^2$
at equilateral and the positive
$B_6^2$
at many acute triangles provide a robust signal for modified gravity.
Further, considering Figure 7, we see that the
$l=6,m=3$
case is very similar to
$m=2$
, with the exception that
$[B_6^3]_\textrm{GR}$
is negative for equilateral triangles and a few linear triangles near the squeezed limit. Similar to
$m=2$
,
$[\alpha_6^{3}]_{\unicode{x03BE}}$
and
$[\alpha_6^{3}]_{\unicode{x03C8}}$
is negative corresponding to a few acute triangles. The magnitude
$|[\alpha_6^{3}]_{\unicode{x03BE}}|$
and
$|[\alpha_6^{3}]_{\unicode{x03C8}}|$
is maximum (
$\sim 20$
) around the black contour near the equilateral limit, and it declines as triangles are distorted from this configuration on either direction in
${\unicode{x03BC}}-t$
space. For the obtuse triangles, the
$[\alpha_6^{3}]_{\unicode{x03BE}}$
is relatively small, remaining below 1 when
$\unicode{x03BE}\lt0$
(or
$\unicode{x03C8}\gt0$
) and ranging between 1 and 4 when
$\unicode{x03BE}\gt0$
(or
$\unicode{x03C8}\lt0$
). The
$|[\alpha_6^{3}]_{b_1}|\lt2$
for all the cases. The positive
$B_6^3$
at equilateral and the negative
$B_6^3$
at a few acute triangles below the black contour signal the modification in the gravity.
Lastly, we discuss the results for
$l=6,m=4$
presented in Figure 8. Considering the first row, we observe that
$[\alpha_6^{4}]_{\unicode{x03BE}=-0.75}$
reaches its highest values for equilateral triangles, ranging between 12 and 16. It gradually declines as the configuration approaches obtuse triangles but increases again near (
${\unicode{x03BC}} = 0.9$
,
$t = 0.6$
). It remains mostly 2 to 4 for the bulk of the triangles near the S-isosceles limit (
$2{\unicode{x03BC}} t=1$
). Notably, there are a few triangles near the L-isosceles limit (
$t=1$
) and at (
${\unicode{x03BC}} = 0.9$
,
$t = 0.6$
), where
$[\alpha_6^{4}]_{\unicode{x03BE}} \sim -12$
. For smaller modification in
$\unicode{x03BE}$
, the magnitude of
$[\alpha_6^{4}]_{\unicode{x03BE}}$
decreases significantly. At
$\unicode{x03BE} = -0.05$
, its values remain between
$0-2$
and it becomes positive for all triangle configurations when
$\unicode{x03BE}$
is small (e.g.
$\unicode{x03BE} = -0.05$
and
$\unicode{x03BE} = 0.05$
). At
$\unicode{x03BE} = 0.25$
,
$[\alpha_6^{4}]_{\unicode{x03BE}}$
becomes
$\sim -2$
for equilateral triangles, increases to
$\sim 8$
for certain L-isosceles triangles, and remains
$\sim 2$
near stretched triangle configurations. For the majority of other triangle shapes,
$|[\alpha_6^{4}]_{\unicode{x03BE}}|$
remains less than 1. As
$\unicode{x03BE}$
approaches
$0.75$
, the magnitude
$|[\alpha_6^{4}]_{\unicode{x03BE}}|$
increases further, reaching a maximum of approximately
$\sim 8$
near the equilateral limit, while consistently exhibiting a negative sign for all S-isosceles triangle configurations.
Considering the second row, we observe that the trend of
$[\alpha_6^{4}]_{\unicode{x03C8}}$
for
$\unicode{x03C8}=-0.75$
is similar to that for
$\unicode{x03BE} = 0.75$
, but the maximum value near the equilateral limit exceeds
$-20$
. The magnitude
$|[\alpha_6^{4}]_{\unicode{x03C8}}| \sim 8$
near the stretched limit and increases to 20 near (
${\unicode{x03BC}} = 0.9$
,
$t = 0.6$
). For the majority of triangle configurations near the S-isosceles limit,
$[\alpha_6^{4}]_{\unicode{x03C8}}$
remains mostly between
$-4$
and
$-8$
. the results for
$\unicode{x03C8} = -0.25$
,
$-0.05$
,
$0.05$
, and
$0.25$
are very similar to those for
$\unicode{x03BE} = 0.75$
,
$0.25$
,
$-0.25$
, and
$-0.75$
, respectively. For
$\unicode{x03C8}=0.75$
, the highest values near equilateral triangles and at (
${\unicode{x03BC}} = 0.9$
,
$t = 0.6$
) range in
$12-20$
.
$[\alpha_6^{4}]_{\unicode{x03C8}} \sim 4-8$
for most of the S-isosceles triangles and
$\sim 1-2$
for linear triangles.
The
$b_1$
majorly affects the
$B_6^4$
corresponding to stretched triangles and a few obtuse triangles. Considering the third row, we see that the
$|[\alpha_6^{4}]_{b_1}|\sim 1-4$
near stretched triangles for
$0.5\leq b_1\leq2$
. The maximum value of
$|[\alpha_6^{4}]_{b_1}|$
is near the black contour around (
${\unicode{x03BC}} = 0.9$
,
$t = 0.6$
), which goes up to 12 for
$b_1=0.5$
and 2 and decreases to 8 and 4 for
$b_1=0.75$
and 1, respectively. For triangles with
$t\gt0.75$
,
$[\alpha_6^{4}]_{b_1}\lt1$
when
$b_1\lt1.2$
and
$[\alpha_6^{4}]_{b_1}\sim 1-2$
when
$b_1\gt1.2$
. The
$b_1$
induces the sign flip near stretched triangles and a few obtuse triangles only. Thus, the positive
$B_6^4$
values observed for equilateral triangles and for triangles with
$0.75 \lt t \lt 0.95$
, along with the negative
$B_6^4$
values for certain L-isosceles triangles (above the black contour), serve as a definitive indicator of modified gravity.
Table 2 summarises the key sign-flip signatures of bispectrum multipoles
$B_l^m$
for testing GR. It encapsulates, for each multipole, the relevant triangle configurations, the contrast between the expected GR sign and the DHOST-induced flip, and their implications for the second-order growth indices.
4.2.4 Higher multipoles
We find that the bispectrum multipoles above
$l=6,m=4$
remain unchanged with respect to
$\unicode{x03BE}$
and
$\unicode{x03C8}$
. Consequently,
$[\alpha_6^5]_{\unicode{x03BE}}=[\alpha_6^6]_{\unicode{x03BE}}=[\alpha_6^5]_{\unicode{x03C8}}=[\alpha_6^6]_{\unicode{x03C8}}=1$
and
$[\alpha_8^m]_{\unicode{x03BE}}=[\alpha_8^m]_{\unicode{x03C8}}=1$
(for all m) across all triangle configurations. In Equation (28), the
$\unicode{x03BE}$
and
$\unicode{x03C8}$
dependence is implicit in the kernels
$G_{12}$
and
$\Delta G_{12}$
which are multiplied by the functions R and
$S_{12}$
, respectively. The R and
$S_{12}$
, which are polynomials of the order
${\unicode{x03BC}}^6_i$
and
${\unicode{x03BC}}^4_i$
, respectively, only contribute to the multipoles
$l\leq 6$
and
$l\leq 4$
, respectively. As a result, the
$\unicode{x03BE}$
dependence is not captured by
$l\gt6$
multipoles.
5. Summary and discussion
In this work, we have explored the potential of the redshift space bispectrum multipole moments of biased tracers as a tool for probing modified gravity theories beyond linear approximation. Here we employ the standard perturbation theory, incorporating second-order perturbations in both the density and velocity fields to capture higher-order properties of gravitational interactions. The analysis is conducted within the framework of the most general scalar-tensor theory of gravity, known as the DHOST theory (Langlois Reference Langlois2019; Kobayashi Reference Kobayashi2019). The influence of modified gravity under the DHOST framework is characterised by three key parameters: a linear growth index,
$\unicode{x03B3}$
, and two second-order growth indices,
$\unicode{x03BE}$
and
$\unicode{x03C8}$
, which encapsulate additional information about structure growth beyond the scope of linear theory (Yamauchi et al. Reference Yamauchi, Yokoyama and Tashiro2017; Yamauchi & Sugiyama Reference Yamauchi and Sugiyama2022; Sugiyama et al. Reference Sugiyama2023b). For the standard GR, these parameters take the values
$\unicode{x03B3}_{\text{GR}} = 6/11$
,
$\unicode{x03BE}_{\text{GR}} = 3/572$
, and
$\unicode{x03C8}_{\text{GR}} = 0$
. Any deviation from these values indicates a departure from GR.
Modifications to gravity affect the clustering properties of biased tracers and the coupling between modes. The bispectrum, being sensitive to the non-linearity and mode coupling introduced by gravitational evolution, serves as a powerful probe for detecting deviations from GR (Peebles Reference Peebles1980; Scoccimarro Reference Scoccimarro2000; Sugiyama et al. Reference Sugiyama2023b). These signatures of modified gravity in the bispectrum become even more pronounced in redshift space, where the anisotropies introduced by peculiar velocities are sensitive solely to the theory of gravity, independent of tracer bias. We decompose the anisotropic redshift space bispectrum in the spherical harmonic basis
$Y_l^m(\hat{\boldsymbol{p}})$
(Scoccimarro et al. Reference Scoccimarro, Couchman and Frieman1999; Bharadwaj et al. Reference Bharadwaj, Mazumdar and Sarkar2020) and study the sensitivity of all possible bispectrum multipole moments
$B_l^m(k_1,{\unicode{x03BC}},t)$
to modifications in gravity across all triangle configurations (
${\unicode{x03BC}},t$
).
We show that the higher-order multipoles, with
$l=2,4,6$
, are comparatively more sensitive to modified gravity than the spherically averaged
$l=0$
moment, which has predominantly been the focus of prior studies. Considering the linear growth, the
$B_0^0$
remains close to the values predicted by the GR for a range of
$\unicode{x03B3}$
. This is consistent with the earlier findings (Borisov & Jain Reference Borisov and Jain2009; Bernardeau & Brax Reference Bernardeau and Brax2011; Tatekawa & Tsujikawa Reference Tatekawa and Tsujikawa2008). The percent level measurements are required to detect deviations in
$\unicode{x03B3}$
using the monopole. Similarly, the
$l=2$
multipoles also show limited sensitivity to the
$\unicode{x03B3}$
. In contrast, the
$l=6$
multipoles, particularly for equilateral and squeezed triangles, are more sensitive to the
$\unicode{x03B3}$
. In addition, the
$B_4^2$
and
$B_4^3$
multipoles for stretched triangles show high sensitivity to the
$\unicode{x03B3}$
. Notably, the highest non-zero multipole is
$(l=8, m=6)$
. However, all
$(l=8)$
multipoles do not provide any independent information on
$\unicode{x03B3}$
beyond what is already captured by the multipoles with
$l \lt 8$
.
The higher-order multipoles exhibit exceptional sensitivity to the non-linear growth indices,
$\unicode{x03BE}$
and
$\unicode{x03C8}$
. Specifically, the multipoles
$B_2^1$
,
$B_4^3$
,
$B_6^0$
,
$B_6^1$
,
$B_6^2$
,
$B_6^3$
, and
$B_6^4$
encapsulate substantial information about departures from GR. The sensitivity of the multipoles varies with triangle configurations. The acute (
$t \gt {\unicode{x03BC}}$
) triangle configurations, in general, demonstrate the strongest sensitivity to
$\unicode{x03BE}$
and
$\unicode{x03C8}$
. For
$B_2^1$
and
$B_4^3$
, the sensitivity peaks near right-angle triangle configurations. In contrast, for the
$l=6$
multipoles, equilateral triangles are the most sensitive, with sensitivity generally decreasing as triangles become more deformed towards squeezed or stretched shapes. An exception is
$B_6^4$
, which shows high sensitivity for stretched triangles as well. We find that the
$(l=6,m\gt4)$
and all
$l=8$
multipoles are unable to capture second-order growth.
Table 2. Summary of the sign-flip features observed in bispectrum multipoles
$B_l^m$
as indicators of deviations from GR. The first column lists the multipoles exhibiting prominent sign-flip behaviour. The second column specifies the triangle configurations associated with the sign flip for each multipole. The third column details the direction of the sign flip, contrasting the expected sign under GR with the altered sign induced in the DHOST scenario. The fourth column presents the physical interpretation of the sign flip and its implications for the second-order growth indices.

We demonstrate that, for specific triangular configurations, the values of these multipoles exhibit opposite signs in modified gravity compared to those predicted by GR. This sign reversal serves as a robust indicator of deviations from GR. The characteristic signatures of modified gravity are as follows: for
$B_2^1$
, negative values for acute triangles and positive values near right-angle triangles; for
$B_4^3$
, negative values for acute and stretched triangles and positive values for obtuse triangles; for
$B_6^0$
, negative values for acute triangles and positive values for equilateral triangles; for
$B_6^1$
, positive values for all acute triangles; for
$B_6^2$
, negative values for equilateral triangles and positive values for other acute triangles; for
$B_6^3$
, positive values near the equilateral limit and negative values for other acute triangles; and for
$B_6^4$
, negative values for many L isosceles and stretched triangles and positive values near the S isosceles limit. These distinct signatures, tabulated in Table 2, underscore the utility of higher-order multipoles in probing deviations from GR.
Although various bispectrum multipoles effectively capture robust signatures of modified gravity, the linear bias
$b_1$
and quadratic bias parameters
$b_2$
and
$b_{s^2}$
have the tendency to replicate these signals if they are not correctly accounted for. While the impact of bias on
$B_2^1$
is minimal, it is significantly more pronounced for the case of
$B_4^3$
. We show that the
$l=6$
multipoles offer a distinct advantage over
$l=2$
and
$l=4$
multipoles as they are not affected by the quadratic bias
$b_2$
and second-order tidal bias
$b_{s^2}$
. Consequently,
$l=6$
multipoles hold considerable potential for probing and constraining theories of modified gravity, emphasising the need to leverage their capabilities in analyses. However, it remains crucial to recognise that similar bispectrum features may arise from either bias or modified gravity. Therefore, a careful disentanglement of these contributions is essential to reliably constrain both the gravitational theory and bias parameters of tracers in the analysis. The ongoing and future surveys such as EuclidFootnote
a
(Laureijs et al. Reference Laureijs2011), DESIFootnote
b
(Levi et al. Reference Levi2013), and SKAFootnote
c
(Weltman et al. Reference Weltman2020) will provide us with an unprecedented volume of data, enabling the detection of modified gravity signatures and constraining gravitational theories with percent-level precision. In future work, we will utilise the fast estimator (Gill & Bharadwaj Reference Gill and Bharadwaj2025) to compute these higher-order bispectrum moments from galaxy survey data and extract critical information about the underlying theory of gravity.

Figure A1. Sensitivity of
$B_0^0$
to variations in the second-order growth indices
$\unicode{x03BE}$
and
$\unicode{x03C8}$
. The panels are identical to those in Figure 2.
Acknowledgements
The author is grateful to Barbie Sangha for helpful discussions and encouragement throughout this project. We thank the anonymous reviewers for their valuable comments.
Data availability statement
The codes and packages for numerical computation involved in this work will be shared on reasonable request to the author.
Appendix A. Results for
$\boldsymbol{l=0}$
and other
$\boldsymbol{l}=\boldsymbol{2,4}$
multipoles
Figures A1, A2, A3, A4, A5, A6, and A7 shows the sensitivity of
$B_0^0$
,
$B_2^0$
,
$B_2^2$
,
$B_4^0$
,
$B_4^1$
,
$B_4^2$
, and
$B_4^4$
, respectively, to the variations in second-order growth indices
$\unicode{x03BE}$
and
$\unicode{x03C8}$
. The results indicate that these multipoles exhibit minimal variation with changes in gravity at the second order. The ratio
$[\alpha_l^m]_{\unicode{x03BE}}$
and
$[\alpha_l^m]_{\unicode{x03C8}}$
remain below 2 in most cases, although they can reach values of up to
$\sim 4$
for a few acute triangle configurations, especially for
$B_2^0,B_2^2,B_4^0,B_4^2,B_4^4$
. Notably, there is no consistent indication of a sign flip signalling the modified gravity in these multipoles, except for
$B_4^2$
which show sign change near triangles near
${\unicode{x03BC}}\approx 0.9, t\approx0.6$
for
$\unicode{x03C8}\gt0$
(or
$\unicode{x03BE}\lt0$
) and near stretched triangles for
$\unicode{x03C8}\lt0$
(or
$\unicode{x03BE}\gt0$
). Additionally, the signal of modified gravity in these multipoles is strongly influenced by the bias parameters
$b_1$
,
$b_2$
, and
$b_{s^2}$
.

Figure A2. Sensitivity of
$B_2^0$
to variations in the second-order growth indices
$\unicode{x03BE}$
and
$\unicode{x03C8}$
. The panels are identical to those in Figure 2.

Figure A3. Sensitivity of
$B_2^2$
to variations in the second-order growth indices
$\unicode{x03BE}$
and
$\unicode{x03C8}$
. The panels are identical to those in Figure 2.

Figure A4. Sensitivity of
$B_4^0$
to variations in the second-order growth indices
$\unicode{x03BE}$
and
$\unicode{x03C8}$
. The panels are identical to those in Figure 2.

Figure A5. Sensitivity of
$B_4^1$
to variations in the second-order growth indices
$\unicode{x03BE}$
and
$\unicode{x03C8}$
. The panels are identical to those in Figure 2.

Figure A6. Sensitivity of
$B_4^2$
to variations in the second-order growth indices
$\unicode{x03BE}$
and
$\unicode{x03C8}$
. The panels are identical to those in Figure 2.

Figure A7. Sensitivity of
$B_4^4$
to variations in the second-order growth indices
$\unicode{x03BE}$
and
$\unicode{x03C8}$
. The panels are identical to those in Figure 2.