# Interacting dark energy from the joint analysis of the power spectrum and
bispectrum multipoles with the EFTofLSS

Maria Tsedrik,1 Chiara Moretti,1,2,5 Pedro Carrilho,1 Federico Rizzo,2,3,4
Alkistis Pourtsidou1,5
1Institute for Astronomy, The University of Edinburgh, Royal Observatory,
Edinburgh EH9 3HJ, UK
2Istituto Nazionale di Astrofisica, Osservatorio Astronomico di Trieste, via
Tiepolo 11, 34143 Trieste, Italy
3Istituto Nazionale di Fisica Nucleare, Sezione di Trieste, via Valerio 2,
34127 Trieste, Italy
4Institute for Fundamental Physics of the Universe, Via Beirut 2, 34151
Trieste, Italy
5Higgs Centre for Theoretical Physics, School of Physics and Astronomy, The
University of Edinburgh, Edinburgh EH9 3FD, UK E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Interacting dark energy models have been suggested as alternatives to the
standard cosmological model, $\Lambda$CDM. We focus on a phenomenologically
interesting class of dark scattering models that is characterised by pure
momentum exchange between dark energy and dark matter. This model extends the
parameter space with respect to $\Lambda$CDM by two parameters, $w$ and $A$,
which define the dark energy equation of state and the strength of the
coupling between dark energy and dark matter, respectively. In order to test
non-standard cosmologies with Stage-IV galaxy clustering surveys, it is
crucial to model mildly nonlinear scales and perform precision vs accuracy
tests. We use the Effective Field Theory of Large-Scale Structure, and we
perform validation tests by means of an MCMC analysis using a large set of
N-body simulations. We find that adding the bispectrum monopole to the power
spectrum multipoles improves the constraints on the dark energy parameters by
$\sim 30\%$ for $k_{\mathrm{max},B}^{l=0}=0.11$ $h$ Mpc-1 without introducing
biases in the parameter estimation. We also find that the same improvement can
be achieved with more moderate scale cuts and the use of bias relations, or
with the addition of the bispectrum quadrupole. Finally, we forecast
constraints assuming Stage-IV surveys specifications at $z=1$ in the
interacting dark energy scenario: $\sigma_{w}=0.08$ and $\sigma_{A}=2.51$ b
GeV-1, as well as for $w$CDM cosmology: $\sigma_{w}=0.1$.

###### keywords:

cosmology: large-scale structure of Universe – cosmology: cosmological
parameters – cosmology: dark energy

††pubyear: 2022††pagerange: Interacting dark energy from the joint analysis of
the power spectrum and bispectrum multipoles with the EFTofLSS–10

##  1 Introduction

In the next few years, Stage-IV galaxy surveys such as the Dark Energy Survey
(DES)111https://www.darkenergysurvey.org/ (Abbott et al., 2021), the Dark
Energy Spectroscopic Instrument (DESI)222https://www.desi.lbl.gov/ (Aghamousa
et al., 2016), Euclid333http://euclid-ec.org (Laureijs et al., 2011), the
Nancy Grace Roman Space Telescope444https://www.nasa.gov/roman (Spergel et
al., 2015), and the Vera C. Rubin Observatory’s Legacy Survey of Space and
Time (LSST)555https://www.lsst.org/ (Ivezić et al., 2019) promise to provide
high precision measurements, which will allow us to study the large-scale
structure of the Universe with unprecedented accuracy. A key goal of these
surveys is to test for deviations from the standard cosmological model,
$\Lambda$CDM. The $\Lambda$CDM model postulates two exotic components: cold
dark matter (CDM), and dark energy in the form of a cosmological constant
($\Lambda$). The nature of these dark sector components is unknown. Moreover,
in recent years issues related to tensions in the determination of the $H_{0}$
and $\sigma_{8}$ cosmological parameters between early and late Universe
probes have emerged. The search for alternatives to $\Lambda$CDM that could
also explain the apparent tensions has led to a plethora of exotic dark
energy, dark matter, and modified gravity theories (for reviews see Copeland
et al., 2006; Bertone et al., 2005; Clifton et al., 2012; Perivolaropoulos &
Skara, 2021). In preparation for the stringent tests of these theories with
forthcoming data, validation tests of the exotic models against simulations
have to be performed. One of these non-standard cosmological models, called
Interacting Dark Energy, is the focus of this study.

The term Interacting Dark Energy (IDE) describes a broad class of models, in
which a non-gravitational coupling between dark energy and dark matter is
allowed. A large suite of such models can be constructed by making different
choices for the form of the coupling function (for instance, Amendola, 2000;
Farrar & Peebles, 2004; Clemson et al., 2012; Pourtsidou et al., 2013).
Interacting dark energy models have also been considered as candidate models
to alleviate the $\sigma_{8}$ and $H_{0}$ tension, with various degrees of
success (Pourtsidou & Tram, 2016; Baldi & Simpson, 2017; Barros et al., 2019;
Di Valentino et al., 2020a, b; Kase & Tsujikawa, 2020; Amendola & Tsujikawa,
2020; Beltrán Jiménez et al., 2021). The vast majority of models exhibits
background energy transfer. As a result, they are severely constrained by
Cosmic Microwave Background (CMB), Baryonic Acoustic Oscillations (BAO) and
Supernovae (SNeIa) measurements (for instance, Xia, 2009; Di Valentino et al.,
2017; van de Bruck & Thomas, 2019; Pan et al., 2019; Gómez-Valent et al.,
2020; Yang et al., 2021; Nunes et al., 2022). In this work we consider a
popular class of IDE models, which is uncoupled in the background and exhibits
only momentum exchange between dark energy and dark matter, at the level of
the linear perturbations. This feature makes the model capable of fitting both
large-scale structure (LSS) and CMB measurements extremely well for a wide
range of couplings (Pourtsidou & Tram, 2016; Baldi & Simpson, 2017; Spurio
Mancini & Pourtsidou, 2022).

In contrast to modified gravity (MG) theories, IDE models do not influence the
speed of gravitational waves or the gravitational potential of planets, hence
they do not need screening mechanisms in order to surpass Solar System tests.
Furthermore, instead of an enhancement in the growth of structure at nonlinear
scales (as in the case of popular MG models, e.g., Sotiriou & Faraoni, 2010),
some interacting models predict a suppression of growth and hence lead to an
alleviation of the $\sigma_{8}$ tension.

Non-standard cosmologies generally have a larger parameter space (more free
parameters) than $\Lambda$CDM, and these may result to observable effects at
linear and nonlinear scales. For example, Baldi & Simpson (2015) showed that
nonlinear scales are much more informative in terms of constraints for dark
matter - dark energy scattering than the linear regime. This highlights the
importance of accurately modelling nonlinear scales for extended models. The
specific IDE model under consideration in this work was studied by Carrilho et
al. (2021). The authors used a set of approximate (COLA) simulations run with
standard cosmology to perform null tests of three perturbative approaches to
compute the halo power spectrum multipoles. The authors conclude that a) a
false $\sim 3\sigma$ detection of exotic dark energy would occur should the
nonlinear modelling be incorrect b) precise and accurate constraints on IDE
parameters are possible when appropriate scale cuts are considered, and c) the
Effective Field Theory of Large-Scale Structure (EFTofLSS) approach (Baumann
et al., 2012; Carrasco et al., 2012; d’Amico et al., 2020; Rizzo et al., 2022)
provides a model accurate up to a wave-number of $k_{\mathrm{max}}\sim 0.3$
$h$ Mpc-1. In this study we apply the same model and extend upon the previous
work by adding the bispectrum monopole and quadrupole to the analysis, and by
using better simulations.

The joint analysis of the power spectrum and bispectrum has been considered
for Stage-IV surveys analyses that aim to constrain $\Lambda$CDM. The addition
of the bispectrum is useful for breaking bias degeneracies and tightly
constraining bias parameters. An example is the significant increase in the
constraining power on the mass fluctuation amplitude $A_{s}$ when the real
space analysis includes the bispectrum (Eggemeier et al., 2021): indeed, the
bispectrum is able to break the strong degeneracy between $A_{s}$ and the
linear bias, $b_{1}$. This is however no longer the case in redshift space,
since the aforementioned degeneracy is already broken at the level of the
power spectrum multipoles. There are various examples of the joint analysis in
redshift space in the literature. To name a few, a Fisher Matrix forecast
analysis was performed in Yankelevich & Porciani (2019) with tree-level
Standard Perturbation Theory (SPT) (Bernardeau et al., 2002); a full
likelihood analysis was done in Gualdi & Verde (2020) and Gualdi et al. (2021)
with the TNS modelling (Taruya et al., 2010) and the Quijote suite of N-body
simulations (Villaescusa-Navarro et al., 2020), as well as in d’Amico et al.
(2020); Philcox & Ivanov (2022); Ivanov et al. (2022), with data from the BOSS
survey (Beutler et al., 2017) and the EFTofLSS model. The above forecasts and
data analyses concentrated on the $\Lambda$CDM model. In this paper we analyse
the informational content of the bispectrum in addition to the power spectrum
in redshift space for a non-minimal, interacting dark energy model. To our
knowledge, this is the first analysis of this kind for exotic dark energy.

This paper is organized as follows: section 2 is dedicated to the specific
interacting dark energy model we use in our work, its main features and impact
on structure formation; in section 3 we describe the EFTofLSS model for the
power spectrum and bispectrum with a brief overview on the nuisance parameters
and their co-dependencies. In section 4 we describe our likelihood pipeline
and performance metrics. Our results are presented in section 5. We conclude
in section 6.

##  2 Interacting Dark Energy Modelling

In this work we focus on the subclass of IDE models which exhibits no coupling
at the background level and introduces a pure momentum transfer at the level
of linear perturbations. In general, non-gravitational interaction between
dark matter and dark energy can be introduced either at the level of the
action (Pourtsidou et al., 2013; Tamanini, 2015; Kase & Tsujikawa, 2020) or
more phenomenologically at the level of the fluid equations (e.g., Väliviita
et al., 2008; Simpson, 2010). This leads to two different approaches for our
model of interest.

The first approach corresponds to the “Type 3” class of models from the
action-based derivation in Pourtsidou et al. (2013). There, the pure exchange
of momentum is generated by a coupling between the dark matter velocity field
and the covariant derivative of the dark energy scalar field. In this case,
the two parameters describing the extension with respect to the standard
cosmological model are given by the coupling strength $\beta$ and the
exponential factor of the quintessence potential $\lambda$. These parameters
are included in the action as follows:

$S_{\phi}=\int\differential
t~{}\differential^{3}x~{}a^{3}\left[\frac{1}{2}(1-2\beta)\dot{\phi}^{2}-\frac{1}{2}|\nabla\phi|^{2}-V_{0}\mathrm{e}^{-\lambda\phi}\right]\,,$
(1)

where $\phi$ is the scalar field (quintessence), the last term is the
quintessence potential of the single exponential form, the dot represents a
derivative with respect to cosmic time and $a$ is the scale factor. This model
is explicitly studied in Pourtsidou & Tram (2016), Spurio Mancini & Pourtsidou
(2022) and has been implemented into the CLASS (Lesgourgues, 2011; Blas et
al., 2011) and CAMB (Lewis & Challinor, 2011) Boltzmann codes. These analyses
vary $\beta$ and $\lambda$, while the quintessence potential normalisation,
$V_{0}$, is varied automatically by the Boltzmann solver in order to match a
given $\Omega_{\phi}$ today.

The second approach is derived in the elastic scattering formalism by Simpson
(2010). The main idea is that a dark matter particle moves through the
isotropic dark energy fluid and experiences a drag-force proportional to the
scattering cross-section between the particle and the fluid. Since a dark
matter particle is more massive than the energy exchange with the dark energy
fluid, they interact via elastic scattering, in a way that resembles Thomson
scattering between non-relativistic electrons and photons in the radiation-
dominated era.

Such interaction involves no energy transfer and hence it leaves the energy
conservation equation unchanged with respect to the standard cosmology.
Moreover, we assume that the dark energy density and velocity fields are
homogeneous (i.e. $\delta_{\mathrm{DE}}=\theta_{\mathrm{DE}}=0$), since dark
energy perturbations are damped within the cosmic horizon. This follows from
the assumption of the dark energy speed of sound being equal to the speed of
light, and was confirmed numerically in Baldi & Simpson (2015). While these
fluctuations become non-negligible at super-horizon scales, even the most
extreme modifications do not contradict the CMB measurements due to the
presence of cosmic variance (Pourtsidou & Tram, 2016). In presence of momentum
exchange between the two components of the dark sector, and with the
approximation of homogeneity in dark energy, this IDE model modifies only the
Euler equation to

$a~{}\partial_{a}\Theta+\left(2+(1+w)~{}\xi~{}\frac{\rho_{\mathrm{DE}}}{H}+\frac{a~{}\partial_{a}H}{H}\right)\Theta+\frac{\nabla^{2}\Phi}{a^{2}H^{2}}=0\,,$
(2)

where $\xi=\sigma_{\mathrm{D}}/m_{\mathrm{c}}$ with $\sigma_{\mathrm{D}}$ is
the cross-section of the dark energy - dark matter interaction and
$m_{\mathrm{c}}$ the mass of dark matter particle, and
$\Theta=\theta_{\mathrm{c}}/aH$ with $\theta_{\mathrm{c}}$ the dark matter
velocity divergence. $H=\dot{a}/a$ is the Hubble rate and it describes the
expansion history in the flat universe composed only of dark matter and dark
energy by

$H^{2}=H_{0}^{2}\left(\Omega_{\mathrm{c},0}~{}a^{-3}+\Omega_{\mathrm{DE},0}~{}a^{-3(1+w)}\right)\,,$
(3)

where as above $w$ is the dark energy equation of state parameter. $H_{0}$,
$\Omega_{\mathrm{c},0}$, $\Omega_{\rm{DE},0}$ are the values of Hubble
parameter, dark matter density and dark energy density parameters today.
Additionally, the gravitational potential stays unchanged and is given by the
following Poisson equation:

$\nabla^{2}\Phi=\frac{3}{2}a^{2}H^{2}\Omega_{\mathrm{c}}\delta_{\mathrm{c}}\,,$
(4)

where $\Omega_{\mathrm{c}}$ is the background dark matter density parameter at
the scale factor $a$ and $\delta_{\mathrm{c}}$ is the corresponding density
contrast.

Figure 1: Top: the ratio of the linear growth factor to the $\Lambda$CDM case
$D_{\Lambda\mathrm{CDM}}=0.475$ for the different values of IDE parameters.
Bottom: the ratio of the logarithmic growth rate to the $\Lambda$CDM case
$f_{\Lambda\mathrm{CDM}}=0.861$ for the different values of IDE parameters.
Both growth parameters are computed at redshift $z=1$. Parameter $A$ is given
in the units of b GeV-1. The red patches in the upper left and the lower right
corners denote the forbidden area of the IDE parameter space due to the
condition $A/(1+w)\geq 0$.

Although there is no direct mapping between the elastic scattering and the
“Type 3” models, a particular subclass of “Type 3” models can reproduce the
drag-like behaviour (Skordis et al., 2015). In Baldi & Simpson (2017) the
authors show how these models can be matched approximately. In this work we
use the phenomenological approach with parameters $w$ and $A=\xi(1+w)$, which
characterise the equation of state for dark energy and the coupling strength,
respectively. We assume that $w$ is constant, it affects the expansion history
in Equation 3 as well as structure formation in Equation 2. At the same time,
$A$ impacts only the perturbation growth and is chosen in this form (i.e. a
combination of $\xi$ and $w$) for the purposes of parameter inference. This is
particularly important for $w\approx-1$, where we can constrain $A$ while
$\xi$ could take an arbitrary large value.

Parameters $w$ and $A$ modify the linear growth factor $D(a)$ and the
logarithmic growth rate $f=\differential\ln{D}/\differential\ln{a}$. For a
constant value of $w$ in the late Universe with redshift $z\leq 10$, in the
canonical (phantom) case of $w_{\mathrm{DE}}>-1\,(w_{\mathrm{DE}}<-1)$ the
growth factor is decreased (increased) with respect to $\Lambda$CDM. The
effect gets more prominent as the redshift decreases, and for larger values of
$\xi$. This leads to a scale-independent suppression (enhancement) of the
power spectrum at large scales. On smaller scales, the modification comes from
the drag term in Equation 2: for $w_{\mathrm{DE}}>-1\,(w_{\mathrm{DE}}<-1)$
there is a dissipation (injection) of kinetic energy within collapsed
structures. In other words, for $w_{\mathrm{DE}}>-1\,(w_{\mathrm{DE}}<-1)$
particles inside collapsed objects experience additional friction (drag),
hence it becomes harder (easier) for them to stream away, with a strong and
scale-dependent enhancement (suppression) of the growth of structure in the
nonlinear regime. This behaviour, anticipated from the theory, has been shown
to hold true in N-body simulations by Baldi & Simpson (2015) and Baldi &
Simpson (2017).

In this work we consider only one redshift $z=1$, at which the impact of the
IDE parameters is observable but not as prominent as at lower redshifts. In
Figure 1 we show the impact of different values of $w$ and $A$ on the growth
parameters at the redshift of our interest. The ratios are taken with respect
to the $\Lambda$CDM values: $D_{\Lambda\mathrm{CDM}}=0.475$ and
$f_{\Lambda\mathrm{CDM}}=0.861$ for our choice of cosmological parameters. The
red patches in the upper left and the lower right corners denote the forbidden
area of the IDE parameter space: the sign of $A$ and $(1+w)$ has to be the
same, since $\xi=A/(1+w)$ from Equation 2 has to be non-negative. For $w<-1$
dark energy starts to dominate over matter at later times than in
$\Lambda$CDM, hence with decreasing $w$ the growth rate quickly becomes $1$
(as in the matter-dominated epoch) nearly independent from the value of $A$.
The linear growth factor is also independent of $A$ and is slightly larger
than the $\Lambda$CDM case for $w<-1$. For $w>-1$ the expected power
suppression at linear scales is present, with a stronger effect for larger
values of $A$.

##  3 Power Spectrum and Bispectrum Modelling

Halos and galaxies are biased tracers of the matter density field (Desjacques
et al., 2018). Their clustering properties up to mildly nonlinear scales can
be described using the perturbative framework. The EFTofLSS formalism is based
on perturbation theory, with the addition of counterterms to account for the
impact of unknown small-scale physics, such as galaxy formation, on large-
scale modes. The functional form of the counterterms is specified by the
symmetries of the density field, while the nuisance parameters are determined
by the data. We note that application to data would not be possible without
huge effort in the development on the theoretical side (see for example
Bernardeau et al. 2002; Crocce & Scoccimarro 2008; Taruya et al. 2010;
Senatore & Zaldarriaga 2014, 2015; Vlah et al. 2015; Baldauf et al. 2015b;
Perko et al. 2016; Vlah & White 2019, and for a more comprehensive list for
the development of the EFTofLSS see the footnote in D’Amico et al. 2022).

###  3.1 Bias expansion

Both dark matter halos (the focus of this study) and galaxies are biased
tracers of the underlying matter distribution. The halo overdensity can be
expanded as $\delta_{h}=\sum_{n}b_{n}/n!\delta^{n}$, where $\delta^{n}$ are
higher powers of the matter density contrast and $b_{n}$ are the corresponding
local bias terms (Fry & Gaztanaga, 1993). Halos are formed by the
gravitational collapse of matter from a spatially finite region. Anisotropies
created by this process impact the local distribution of objects, and thus the
bias parameters, with tidal effects (Chan et al., 2012; Baldauf et al., 2012)
and introduce scale dependency by breaking the local assumption (McDonald &
Roy, 2009). Effects from an ellipsoidal gravitational collapse are encoded in
the non-local operators, dubbed $\mathcal{G}_{2}$ and $\Gamma_{3}$: the first
one represents the tidal stress tensor produced by the gravitational
potential, while the second one denotes the difference between the tidal
stress tensors from the gravitational and velocity potentials. Effects due to
the finite size of the collapsing region are taken into account by the higher-
order derivative of the matter field $\nabla^{2}\delta$. Additionally, the
relation between matter and galaxies is also affected by stochasticity, whose
impact can be modeled as an additional contribution to Poisson shot noise on
large scales (Seljak et al., 2009; Paech et al., 2017), while on small scales
it can exhibit a scale-dependence due to the halo-exclusion effect (Baldauf et
al., 2013a; Eggemeier et al., 2020). We consider the bias expansion given by
(Desjacques et al., 2018; Assassi et al., 2014):

$\delta_{h}=b_{1}\delta+\frac{b_{2}}{2}\delta^{2}+b_{\mathcal{G}_{2}}\mathcal{G}_{2}+b_{\Gamma_{3}}\Gamma_{3}+b_{\nabla^{2}\delta}\nabla^{2}\delta+\epsilon+\epsilon_{\delta}\delta\,,$
(5)

where $b_{1}$ and $b_{2}$ are the linear and quadratic bias, and the last two
terms correspond to stochasticity contributions. Note that the higher
derivative term gives rise to a $k^{2}P_{L}(k)$ correction, which is
degenerate with one of the counterterms. Therefore, we omit
$b_{\nabla^{2}\delta}$ in our analysis, as previously done in Oddo et al.
(2021). Moreover, we omit from Equation 5 all operators that do not contribute
to our chosen theoretical models for the power spectrum and bispectrum.

###  3.2 Power spectrum and bispectrum

The mildly nonlinear power spectrum and bispectrum in real space using the
EFTofLSS can be modeled with three ingredients: a) the leading-order
contribution in SPT at one-loop for the power spectrum with the addition of
the EFTofLSS counterterms, and at tree-level for the bispectrum, b) a bias
expansion, and c) an infrared-resummation routine to take the damping of the
oscillatory features into account (Baldauf et al., 2015a; Blas et al., 2016).
In addition to that, to model redshift space quantities we should take the
distortions due to peculiar velocities into consideration. This introduces an
additional dependency of the observables on the direction of separation
between the objects with respect to the line-of-sight, $\boldsymbol{\hat{s}}$.
For the power spectrum this dependency is characterized by the cosine of the
angle formed by the separation wave-number with the line-of-sight,
$\mu=(\boldsymbol{k}\cdot\boldsymbol{\hat{s}})/k$. For the bispectrum we need
two cosines, $\mu_{1}$ and $\mu_{2}$, of $\boldsymbol{k}_{1}$ and
$\boldsymbol{k}_{2}$ with the line of sight $\boldsymbol{\hat{s}}$, with the
condition
$\boldsymbol{k}_{1}+\boldsymbol{k}_{2}+\boldsymbol{k}_{3}=\boldsymbol{0}$.
However, it is more convenient to describe triangle configurations with three
wave-numbers ($k_{1}$, $k_{2}$, $k_{3}$), and two angles,
$\theta=\arccos{\mu_{1}}$ and $\xi$ being the azimuthal rotational angle of
$\boldsymbol{k}_{2}$ around $\boldsymbol{k}_{1}$ (Scoccimarro et al., 1999a).
In total, the redshift space bispectrum depends on five variables, three of
which determine the shape of the triangle while the remaining two define its
orientation with respect to the line-of-sight.

The complete set of equations we use in this work will be detailed in a
forthcoming paper by Moretti et al. (2022), where we also give a detailed
description of the code used to perform the analysis. Here we briefly
highlight the main formulae and techniques used. These formulae are the same
as in $\Lambda$CDM, as the interaction under consideration only modifies the
linear growth factor, $D$ (appearing below via $P_{L}$) and the linear growth
rate of structure, $f$.

The expression for the halo power spectrum in redshift space that we adopt in
our analysis is given by:

$\displaystyle P_{hh}(k,\mu)$ $\displaystyle=Z_{1}(\boldsymbol{k})P_{L}(k)$
$\displaystyle+2\int\differential^{3}\boldsymbol{q}[Z_{2}(\boldsymbol{q},\boldsymbol{k}-\boldsymbol{q})]^{2}P_{L}(q)P_{L}(|\boldsymbol{k}-\boldsymbol{q}|)$
$\displaystyle+6Z_{1}(\boldsymbol{k})P_{L}(k)\int\differential^{3}\boldsymbol{q}Z_{3}(\boldsymbol{k},\boldsymbol{q},-\boldsymbol{q})P_{L}(q)$
$\displaystyle+P_{\mathrm{ctr}}(\boldsymbol{k})+P_{\mathrm{stoch}}(k)\,,$ (6)

with $Z_{1}(\boldsymbol{k})$, $Z_{2}(\boldsymbol{k}_{1},\boldsymbol{k}_{2})$
and $Z_{3}(\boldsymbol{k}_{1},\boldsymbol{k}_{2},\boldsymbol{k}_{3})$ the
redshift space kernels for the SPT loop corrections (Scoccimarro et al.,
1999b), and $P_{L}$ being the IR-resummed linear power spectrum (Baldauf et
al., 2015a; Blas et al., 2016). The stochastic power spectrum is given by
(Eggemeier et al., 2020)

$P_{\mathrm{stoch}}(k)=(1+\alpha_{P}+\epsilon_{k^{2}}k^{2})\bar{n}^{-1}\,,$
(7)

with two free parameters, $\alpha_{P}$ and $\epsilon_{k^{2}}$, describing a
constant deviation from purely Poisson shot noise $\bar{n}^{-1}$, and higher-
order scale dependent corrections generated to account for the short-range
non-locality, respectively. The physical meaning of $\alpha_{P}$ can be
explained with the halo exclusion effect (Smith et al., 2007; Baldauf et al.,
2013b); its value can be either positive (e.g., in galaxy populations with
high satellite fractions) or negative (e.g., in central galaxies of massive
halos) (Baldauf et al., 2013b). The term $(1+\alpha_{P})\,\bar{n}^{-1}$
corresponds to the noise parameter $N$ usually used in the literature, which
denotes a scale-independent noise power spectrum. In this work we include
scale-dependency into the shot noise, a feature which is supported by
numerical simulations (Ginzburg et al., 2017). Furthermore, as we show below,
the data prefers a strong scale-dependence in the shot noise (see subsection
4.1 and Figure 2).

The redshift space loop corrections lead to 28 independent integrals, which we
compute using the Fast-PT algorithm (McEwen et al., 2016; Fang et al., 2017).
The EFTofLSS counterterms for this model are

$\displaystyle P_{\mathrm{ctr}}(k,\mu)$
$\displaystyle=-2\tilde{c}_{0}k^{2}P_{L}(k)-2\tilde{c}_{2}k^{2}f\mu^{2}P_{L}(k)$
$\displaystyle-2\tilde{c}_{4}k^{2}f^{2}\mu^{4}P_{L}(k)+P_{\mathrm{ctr},\nabla^{4}\delta}(k,\mu)\,,$
(8)

with an additional counterterm proportional to $\mu^{4}k^{4}P_{L}(k)$ to
include higher-order contributions and model to some extent the Finger of God
effect (FoG):

$P_{\mathrm{ctr},\nabla^{4}\delta}(k,\mu)=c_{\nabla^{4}\delta}k^{4}f^{4}\mu^{4}(b_{1}+f\mu^{2})^{2}P_{L}(k)\,.$
(9)

This power spectrum model includes in total 10 parameters: {$b_{1}$, $b_{2}$,
$b_{\mathcal{G}_{2}}$, $b_{\Gamma_{3}}$, $\tilde{c}_{0}$, $\tilde{c}_{2}$,
$\tilde{c}_{4}$, $\alpha_{P}$, $\epsilon_{k^{2}}$, $c_{\nabla^{4}\delta}$ }.

For the tree-level bispectrum the expression at tree-level is given by

$\displaystyle
B_{g}(\boldsymbol{k}_{1},\boldsymbol{k}_{2},\boldsymbol{\hat{s}})=B_{\mathrm{stoch}}(\boldsymbol{k}_{1},\boldsymbol{k}_{2},\boldsymbol{\hat{s}})$
$\displaystyle+\sum_{\boldsymbol{k}_{1}\leq\boldsymbol{k}_{i}\leq\boldsymbol{k}_{j}\leq\boldsymbol{k}_{3}}2Z_{1}(\boldsymbol{k}_{i})Z_{1}(\boldsymbol{k}_{j})Z_{2}(\boldsymbol{k}_{i},\boldsymbol{k}_{j})P_{L}(k_{i})P_{L}(k_{j})\,,$
(10)

with the stochastic contribution

$B_{\mathrm{stoch}}(\boldsymbol{k}_{1},\boldsymbol{k}_{2},\boldsymbol{\hat{s}})=(1+\alpha_{B})\bar{n}^{-1}\sum_{i=1}^{3}Z_{1}^{2}(\boldsymbol{k}_{i})P_{L}(k_{i})\,,$
(11)

where $\alpha_{B}$ describes a deviation from the Poisson limit. In contrast
to Ivanov et al. (2022), we do not include the FoG modelling and its
corresponding parameter in our tree-level model, since such a contribution is
expected only at 1-loop level. This leaves the halo bispectrum to be described
by 4 parameters: {$b_{1}$, $b_{2}$, $b_{\mathcal{G}_{2}}$, $\alpha_{B}$ }.

Overall, we express our theoretical prediction in terms of multipoles, in
order to eliminate the $\mu$-dependence. For the power spectrum we compute

$P_{l}(k)=\frac{2l+1}{2}\int_{-1}^{1}\differential\mu
P_{hh}(k,\mu)\mathcal{P}_{l}(\mu)$ (12)

where $\mathcal{P}_{l}(\mu)$ are the Legendre polynomials of order $l$ and
$P_{hh}(k,\mu)$ is the redshift space power spectrum given by subsection 3.2.
In this analysis we use the monopole ($l=0$), quadrupole ($l=2$), and
hexadecapole ($l=4$). Studies have shown that these three multipoles contain
most of the cosmological information (e.g., Taruya et al., 2010; Beutler et
al., 2014; Markovic et al., 2019).

An analogous expansion can be applied to the bispectrum, by expanding in terms
of spherical harmonics as

$B_{h}(k_{1},k_{2},k_{3},\theta,\xi)=\sum_{l}\sum_{m=-l}^{l}B_{lm}(k_{1},k_{2},k_{3})Y_{l}^{m}(\theta,\xi)$
(13)

with the same convention for triangle description as above (Scoccimarro et
al., 1999a). The expansion coefficients are given by

$B_{lm}(k_{1},k_{2},k_{3})=\int_{-1}^{1}\differential\cos{\theta}\int_{0}^{2\pi}\differential\xi
B_{h}(k_{1},k_{2},k_{3},\theta,\xi)Y_{l}^{m}(\theta,\xi)\,.$ (14)

As is typical in the literature, when we describe the redshift space
multipoles of the bispectrum, we consider only the case with $m=0$, as it has
been shown that such seemingly strong compression does not lead to a
significant loss of information (Gagrani & Samushia, 2017). This allows us to
write the bispectrum version of Equation 12, since the spherical harmonics
reduce to Legendre polynomials in the $m=0$ case:

$\displaystyle B_{l}(k_{1},k_{2},k_{3})=\frac{2l+1}{2}$
$\displaystyle\crossproduct\int_{-1}^{1}\differential\cos{\theta}\left[\frac{1}{2\pi}\int_{0}^{2\pi}\differential\xi
B_{h}(k_{1},k_{2},k_{3},\theta,\xi)\right]\mathcal{P}_{l}(\cos{\theta})\,.$
(15)

Let us briefly list the known degeneracies, which follow from previous
studies, our tests and theory directly:

  1. 1.

$f-b_{1}$ are anti-correlated: both parameters control the amplitude on large
scales, they occur combined with various powers in the model, hence if one
gets larger the other should decrease.

  2. 2.

The counterterms $\tilde{c}_{0}$, $\tilde{c}_{2}$, $\tilde{c}_{4}$ and
$c_{\nabla^{4}\delta}$ are all slightly degenerate with one another, but
especially strong is the anti-correlation between $\tilde{c}_{4}$ and
$\tilde{c}_{2}$. This becomes clear after the projection of subsection 3.2
into multipoles: for the monopole the first three terms in this equation yield
$-1/2\,(2\,\tilde{c}_{0}+0.7\,\tilde{c}_{2}f+0.4\,\tilde{c}_{4}f^{2})\,k^{2}P_{L}$,
while for the quadrupole we obtain
$-5/2\,(0.3\,\tilde{c}_{2}f+0.2\,\tilde{c}_{4}f^{2})\,k^{2}fP_{L}$. For
instance, for $f\approx 1$ the counterterm parameters $\tilde{c}_{4}$ and
$\tilde{c}_{2}$ cancel each other out if $\tilde{c}_{4}\approx-\tilde{c}_{2}$.
In these ad-hoc calculations, we have neglected the angle-dependence in the
linear power spectrum due to IR-resummation, which however will affect the
counterterms in the same way and does not impact our argument. Additionally,
$\tilde{c}_{4}$ cannot be constrained without the power spectrum hexadecapole
data. This can be seen directly from subsection 3.2 where there are three
terms with coefficients $\tilde{c}_{0}$, $\tilde{c}_{2}$, $\tilde{c}_{4}$ that
all of them have a form of $k^{2}P_{L}(k)$. From this it follows that we need
all three multipoles to constrain all three parameters.

  3. 3.

The nonlinear bias parameters $b_{2}$, $b_{\mathcal{G}_{2}}$ and
$b_{\Gamma_{3}}$ are all degenerate: $b_{2}$ is correlated with
$b_{\mathcal{G}_{2}}$, while $b_{\mathcal{G}_{2}}$ and $b_{\Gamma_{3}}$ are
strongly anti-correlated. The power spectrum model includes overall 30
independent contributions (28 for the loop-corrections and 2 for the shot
noise), in which the bias parameters come as products with various powers and
various combinations: for instance, $b_{2}$ appears as $b_{1}b_{2}$,
$b_{2}^{2}$, $b_{2}f$, so that it can be constrained with power spectrum data
only, while $b_{\Gamma_{3}}$ appears only twice and only linearly in
combination with $b_{1}$ and $f$, which makes it harder to constrain. However,
the inclusion of the tree-level bispectrum solves this issue, since the
bispectrum model does not include $b_{\Gamma_{3}}$. Additionally, due to the
absence of $b_{\Gamma_{3}}$, $b_{2}$ and $b_{\mathcal{G}_{2}}$ are better
constrained since the strong degeneracy is not present anymore. Therefore, we
expect measurements of $b_{2}$ and $b_{\mathcal{G}_{2}}$ to be more precise
from bispectrum measurements only, although $b_{1}$ is better constrained by
power spectrum only measurements.

  4. 4.

The noise parameters are anti-correlated with one another since they both
control the deviation from Poisson shot noise (constant and scale-dependent).
They dominate on small scales together with other effects like FoG and the
monopole counterterm, hence we expect correlations between $\alpha_{P}$,
$\epsilon_{k^{2}}$ and $\tilde{c}_{0}$, $c_{\nabla^{4}\delta}$.

##  4 MCMC analysis

###  4.1 Data

In this work we analyse halo power spectrum and bispectrum measurements in
redshift space from the large set of 298 Minerva N-body simulations (Grieb et
al., 2016; Lippich et al., 2019). These follow the evolution of $1000^{3}$
dark matter particles in a cubic box of side $L=1500$ Mpc $h^{-1}$,
corresponding to a total effective volume of $\sim 1,000$ Gpc3 $h^{-3}$. The
mass cut in the halo catalogues is set to $1.12\cdot
10^{13}\,M_{\odot}\,h^{-1}$, with a resulting halo mean number density of
$\bar{n}=2.13\cdot 10^{-4}$ $h^{3}$ Mpc-3. We limit this study to a single
redshift $z=1$. The binning scheme we apply is given by:
$k_{i}=(2+(i-1))k_{f}$, with the fundamental frequency $k_{f}=2\pi L^{-1}$,
the bin size is $\Delta k=k_{f}$, and the total number of bins is $N_{k}=128$
for the power spectrum and $N_{k}=29$ for the bispectrum (with a total number
of triangles $N_{t}=2766$). For a detailed description of the measurements,
the halo catalog construction and the estimators, we refer the reader to Oddo
et al. (2021, 2020); Rizzo et al. (2022).

The fiducial cosmology of the simulations is the $\Lambda$CDM best-fit of the
combined analysis of WMAP and BOSS DR9 results (Sánchez et al., 2013) and is
given by: $h=0.695$, $\Omega_{\mathrm{m}}=0.285$, $\Omega_{\mathrm{b}}=0.046$,
$n_{s}=0.9632$, $\sigma_{8}=0.828$. This means that the fiducial IDE
parameters we should recover are $w=-1$ and $A=0$.

Figure 2: Left: average measurements over the 298 Minerva simulations of the
halo power spectrum multipoles – monopole (blue crosses), quadrupole (purple
dots) and hexadecapole (orange triangles), the solid lines correspond to the
posterior-averaged joint model with $k^{l=0,2}_{\mathrm{max},P}=0.225$ $h$
Mpc-1, $k^{l=4}_{\mathrm{max},P}=0.1$ $h$ Mpc-1 and
$k^{l=0}_{\mathrm{max},B}=0.111$ $h$ Mpc-1. The large scattering in the
residuals of the hexadecapole is a consequence of the effective binning scheme
we adopt in the theoretical prediction. Right: the halo bispectrum monopole
average (green crosses). The average is taken over the 298 N-body simulations
and the Poisson shot noise contribution is not subtracted. Shaded areas denote
the uncertainties on the mean. The bottom panels show the relative error on
the residuals of the chain-averaged joint model and the mean of the data. The
dashed vertical lines correspond to the value of $k_{\mathrm{max}}$ used in
the fit.

###  4.2 Covariance

The covariance matrix is constructed from 10,000 mock halo catalogues produced
with the Pinocchio code (Monaco et al., 2002; Munari et al., 2017), 298 of
which have initial conditions chosen to match the initial conditions of the
298 Minerva simulations, while the rest have different initial seeds. The
total halo power spectrum of the mock catalogs is matched to the one of the
N-body simulations at large scales by adjusting the mass threshold (Oddo et
al., 2020). The large number of mocks allows us to suppress the noise in the
off-diagonal components, giving a more accurate covariance, which allows us to
better assess the goodness of the fit of the theoretical models we study.

The uncertainties associated to a measurement inferred from this covariance
correspond to the volume of a single realization. However, in our analysis we
fit the 298 realizations simultaneously by summing up their log-likelihoods
(see subsection 4.3), thus we effectively fit for the full simulation volume
of $\sim 1000$ Gpc3 $h^{-3}$. This volume is unrealistically large when
compared to a typical redshift survey, but can be used to study the
theoretical systematics of the model. To perform a more realistic analysis,
the volume can be decreased in two ways:

  1. 1.

Fitting fewer measurements: it is sufficient to decrease the number of
measurements (e.g., in our case 4 realizations correspond to $13.5$ Gpc3
$h^{-3}$), while using the same covariance for a single measurement. However,
during various tests we found that the resulting parameters and the goodness-
of-fit vary substantially depending on the particular set of simulations
selected. We interpret this as a consequence of the fact that the volume of a
single simulation is large, since we obtain large amount of low wave-number
modes, which are heavily affected by cosmic variance.

  2. 2.

Rescaling the covariance: this mimics the effect of analysing a subset of
measurements, but with the same suppression in cosmic variance as if we were
fitting all of them. Unfortunately, due to the artificial rescaling of the
covariance we cannot use the traditional measure of goodness-of-fit, the
$\chi^{2}$-statistic, because the increased error bars do not match the
scatter in the measurements.

In order to mimic the error budget of Stage-IV spectroscopic galaxy surveys,
we rescale the measured covariance $C_{ij}$ by a factor $\eta$:

$C_{ij}^{\mathrm{resc}}=\frac{\eta}{N_{\mathrm{R}}+N_{\mathrm{C}}}\sum_{\alpha=1}^{N_{\mathrm{R}}+N_{\mathrm{C}}}(X^{\alpha}_{i}-\overline{X_{i}})(Y^{\alpha}_{j}-\overline{Y_{j}})\,,$
(16)

where $N_{\mathrm{R}}$ and $N_{\mathrm{C}}$ are the number of measurements
from the simulations and catalogs, $X,Y\in\\{P_{l_{1}},B_{l_{2}}\\}$ with $l$
being the multipole order, $\alpha$ is the realization index, $i$ and $j$
denote indices of the wave-number or triangle configuration bin. This factor
$\eta$ is a ratio between the effective volumes of all simulation realisations
and a typical Stage-IV survey with $V\sim 8$ Gpc${}^{3}\,h^{-3}$ at $z=1$,
giving $\eta\approx 126$. The effective volume is computed as (Tegmark, 1997)

$V_{\mathrm{eff}}=\left[\frac{\bar{n}P_{hh}(k^{*})}{1+\bar{n}P_{hh}(k^{*})}\right]^{2}V\,,$
(17)

where $\bar{n}$ is the mean number density, $V$ is the volume of the measured
simulation and $k^{*}=0.1\,h$ Mpc-1 is the reference scale. This allows us to
roughly match our simulation data to that expected from future surveys, giving
similar constraining power and signal-to-noise ratio. It also allows us to
extend the range of validity of the theoretical model up to higher values of
$k_{\mathrm{max}}$ with respect to the cuts we would have to adopt in the case
of the non-rescaled covariance and the full volume of the simulations.

It is important to notice that shot noise is not subtracted from our
measurements (both the Minerva and the Pinocchio ones). This is particularly
prominent in the power spectrum monopole, as can be seen at small scales in
Figure 2. In this figure we show the average of the measurements with the
corresponding error bars from the rescaled covariance, the chain averaged
model from the joint Bayesian analysis (described below) and their residuals.
Notice that the very small scatter in the residuals of the monopole and
quadrupole (middle-left panel) is an artifact of our inflated errorbars, while
for the hexadecapole (bottom-left panel) we have some residual scatter due to
the fact that we evaluate the model on the effective modes rather than
performing the full average binning.

###  4.3 Likelihood evaluation

In a Bayesian framework (Bayes, 1763; MacKay, 2003), the probability of a
model given some data is proportional to the product between the probability
of obtaining the data under the assumption that the model is correct and the
probability of the model parameters to have the corresponding values. In other
words, the posterior distribution is proportional to the likelihood function
multiplied by the prior distribution.

In our analysis we choose the priors on nuisance parameters to be non-
informative: they are uniform and very broad. We consider as fiducial values
for these parameters the best fit values obtained in previous studies (Oddo et
al., 2021; Rizzo et al., 2022). Tests were performed to verify that the prior
choice does not influence the parameter inference process. However, the priors
on the IDE parameters $w$ and $A$ are informative and depend on each other.
The equation of state parameter $w$ has a flat prior with size $[-2,-0.5]$.
For $w<-2$ dark energy starts to dominate over matter at a later epoch than
our redshift of interest $z=1$, so that the growth rate converges quickly to 1
independent from $\xi$ as in the matter-dominated epoch. As for the upper
bound, we have observational evidence of the late time acceleration, which
implies $w<-1/3$. The parameter $A$ also has a flat prior with size
$[-20,20]$, but with the additional condition that $A$ must have the same sign
as $1+w$. The upper limit for $A$ is given under the consideration of the
upper bound for $w$ and the most extreme value of $\xi=50$ b GeV-1 tested out
in simulations (Baldi & Simpson, 2015), while the lower limit is chosen to
have a symmetric prior around the fiducial value of $A=0$. The additional
condition on the sign is explained by the fact that the parameter
$\xi=A/(1+w)$ has to be non-negative, since it represents the ratio of two
positive quantities: the cross section of dark sector interactions and the
mass of dark energy particles (see discussion at the end of section 2).

Given $N_{\mathrm{R}}=298$ independent realizations, we evaluate the total
likelihood as the product of individual likelihoods, under the assumption that
each is well described by a multivariate Gaussian. As already mentioned,
instead of fitting the mean of all realizations with the total volume
covariance, we fit each individual simulation and compute

$\displaystyle-2\ln{\mathcal{L}_{\mathrm{tot}}}=-2\sum_{\alpha}^{N_{\mathrm{R}}}\ln{\mathcal{L}_{\alpha}}=-\sum_{\alpha}^{N_{\mathrm{R}}}\chi^{2}_{\alpha}=$
$\displaystyle\sum_{\alpha=1}^{N_{\mathrm{R}}}\sum_{i,j=1}^{N_{\mathrm{b}}}\left(X_{i}^{\alpha}-X^{\mathrm{theo}}_{i}\right)C_{X,ij}^{-1}\left(X_{j}^{\alpha}-X^{\mathrm{theo}}_{j}\right)\,,$
(18)

where $X^{\alpha}_{i}$ is the power spectrum or bispectrum multipole
measurement from the $\alpha$-realization of the $i$-th Fourier bin from the
total number of bins $N_{\mathrm{b}}$, $X^{\mathrm{theo}}_{i}$ is the
theoretical prediction and $C_{ij}$ is the rescaled covariance matrix from
equation Equation 16. We also test applying the correction to the likelihood
function suggested in Sellentin & Heavens (2016), which should take into
account the finite number of mocks used in the construction of the covariance
matrix:

$-2\ln{\mathcal{L}_{\alpha}}=N_{\mathrm{M}}\ln{\left(1+\frac{\chi^{2}_{\alpha}}{N_{\mathrm{M}}-1}\right)}\,,$
(19)

where $N_{\mathrm{M}}$ is the number of mock catalogs. We found no impact on
the posterior distributions, which confirms our expectation of a negligible
effect given the large number of mocks catalogues used.

Our likelihood pipeline uses the emcee package (Foreman-Mackey et al., 2013)
to sample the posterior distribution. We run the MCMC with 300 walkers, and
assume them to have reached convergence after 30000 steps in the analysis
which only includes the power spectrum multipoles, and after 15000 steps when
the bispectrum multipole(s) is (are) included.

###  4.4 Model evaluation

The evaluation of the theoretical model should in principle be performed on
the same grid that is used for the measurements, and the model should be
binned in the same way. This can be understood as the resolution of the
following discrepancy: the theoretical prediction assumes an infinite universe
with an infinite number of $k$-values, while in our case the measurements are
performed in boxes with periodic boundary conditions, hence over a restricted
number of $k$-bins. In general, there are four approaches to resolve this:

  1. 1.

Full binning: compute an exact average of the theoretical predictions for the
observables over each bin in Fourier space. In other words, calculate the
theoretical model at each $q_{i}\in k_{i}$, with $k_{i}$ being a Fourier bin
and $q_{i}$ denoting all discrete wave-numbers in a bin of size $\Delta
k_{i}$, and then find a mean value. In our work the bin size equals the
fundamental frequency $\Delta k=2\pi/L$. Advantage: this is the most
consistent way; disadvantage: it is computationally demanding especially when
re-computing the model at each step of the likelihood evaluation.

  2. 2.

Effective binning: evaluate the theoretical predictions only at one
“effective” wave-number, which is computed as the average over each bin
centered at $k$ for the power spectrum multipoles, while for the bispectrum
the “effective” triplet is computed in a similar but hierarchical way (see
equation B.2 in Rizzo et al., 2022). Advantage: less computationally
expensive; disadvantage: it might introduce systematic errors comparable to
the statistical uncertainties for the full volume of our data set (Oddo et
al., 2020).

  3. 3.

Expansion approach (Oddo et al., 2021): this is an extension of the effective
approach which includes a Taylor expansion of the theoretical model around the
effective wave-number up to second order. Advantage: high accuracy, especially
for the hexadecapole of the power spectrum; disadvantage: it can be
computationally demanding if a large number of terms is involved.

  4. 4.

Integration approach: similar to the averaging approach, but the integral is
computed instead of building an average as a sum over the multipoles, divided
by the number of all discrete vectors in a bin. This approach is popular in
the literature (e.g., Ivanov et al., 2022; Eggemeier et al., 2021) but
similarly to the effective binning it can introduce non-negligible systematic
errors, especially for folded triangles. However, this issue can be bypassed
by introducing “discreetness weights” as corrections (Ivanov et al., 2022).

In our study we choose the effective binning for power spectrum and bispectrum
multipoles. Our tests show that this type of binning does not have an impact
on the inferred parameters due to the enhanced size of our uncertainties after
rescaling the covariance, while being the most computationally efficient.

We also note that in the bispectrum case we include all triangular
configurations with $k_{1}\geq k_{2}\geq k_{3}$ in our analysis, and do not
separate them into squeezed/equilateral/isosceles configurations as is common
in the literature. Previous studies (e.g., Oddo et al., 2020; Eggemeier et
al., 2021; Smith et al., 2008) show that the signal-to-noise ratio and the
accuracy of a model increase for high $k_{\mathrm{max}}$ if nearly equilateral
triangle bins are excluded, while still including elongated triangle
configurations. Whether or not this selection improves our constraints of dark
energy models when adopting realistic error bars is a subject of our future
studies. However, we expect the effect to be minimal.

Figure 3: Left: performance metrics as a function of the maximum wave-number
of the monopole and quadrupole (blue crosses). Middle: performance metrics
with the inclusion of the hexadecapole (orange crosses) and with the best
scale cut $k^{l=0,2}_{\mathrm{max},P}=0.225$ $h$ Mpc-1 from the left panel.
Right: joint models — the power spectrum monopole and quadrupole with the
bispectrum monopole (light-green crosses), the power spectrum monopole and
quadrupole with the bispectrum monopole and one counterterm parameter set to
zero (dark-green dots), all power spectrum multipoles with
$k^{l=4}_{\mathrm{max},P}=0.1$ $h$ Mpc-1 and the bispectrum monopole (dark-red
triangles). The dashed grey lines correspond to the values of
$k_{\mathrm{max}}$. The confidence intervals for FoB are calculated by direct
integration of a two-dimensional Gaussian over an ellipse between $\pm$ FoB
and equating the integral to the 68% and 95% percentile thresholds, which
results in FoB equal to 1.52 and 2.49, respectively.

###  4.5 Model selection

Our goal is to determine the appropriate scale cuts and forecast the
constraining power of future surveys for the IDE parameters, $w$ and $A$. It
is crucial to determine the largest range of scales which can be used to
extract information from the data in an unbiased manner and tighten the
cosmological constraints. This can be also achieved by removing or alleviating
degeneracies between the parameters of the theoretical model. For instance,
bias relations (i.e., the relations between bias parameters based on their
connection due to the underlying physics of gravitational collapse) are able
to reduce the parameter space without introducing bias in the estimated
parameters. In order to evaluate whether a certain scale cut or a reduction of
degeneracies in parameter space provide a better description of the
observables we need a way to quantitatively compare the posterior
distributions from our Markov chains. Since the $\chi^{2}$-statistic is not
valid anymore as a measurement of the goodness-of-fit, due to the rescaling of
the covariance matrix (see discussion in subsection 4.2), we make use of two
performance metrics to select the best scale cut and/or reduced model (see a
similar approach in Eggemeier et al., 2020; Osato et al., 2019):

  1. 1.

the Figure of Bias (FoB) quantifies the relative separation of the measured
parameters from their fiducial values in terms of the variance of the
posterior distribution. We calculate the 68% - 95% percentile thresholds for
the FoB by assuming the posterior distribution of the parameters of interest
to be Gaussian. Hence, for one parameter the 68% threshold equals 1, while for
two parameters the integration of two-dimensional Gaussian over an ellipse
results in a FoB of 1.52. It is defined as

$\mathrm{FoB}=\sqrt{(\theta_{\mathrm{fid}}-\bar{\theta})S^{-1}(\theta_{\mathrm{fid}}-\bar{\theta})}\,,$
(20)

where $\bar{\theta}$ and $\theta_{\mathrm{fid}}$ are vectors labeling,
respectively, the posterior averages and fiducial values of the parameters we
want to measure, and $S=\mathrm{cov}(\theta)$ is the covariance matrix of the
parameters calculated from the chains. In our case, $\theta=[w,A]$ and
$\theta_{\mathrm{fid}}=[-1,0]$;

  2. 2.

the Figure of Merit (FoM) shows the merit of the considered model with respect
to the parameters varied in the fit. It is the inverse of the volume of the
68% contours of the parameters, effectively giving a global measure of the
inverse size of the uncertainties obtained on the parameters (Wang, 2008).
This is given by

$\mathrm{FoM}=\frac{1}{\sqrt{\mathrm{det}(S)}}\,.$ (21)

##  5 Results

###  5.1 Joint base model

In total our base model contains 13 parameters: {$w$, $A$, $b_{1}$, $b_{2}$,
$b_{\mathcal{G}_{2}}$, $b_{\Gamma_{3}}$, $\tilde{c}_{0}$, $\tilde{c}_{2}$,
$\tilde{c}_{4}$, $\alpha_{P}$, $\epsilon_{k^{2}}$, $c_{\nabla^{4}\delta}$,
$\alpha_{B}$ }. Firstly, we aim to find the highest $k_{\mathrm{max}}$ at
which the model parameters are not biased while the constraints on the dark
energy parameters are the most informative (i.e., FoM reaches its maximum).
For this purpose we run our likelihood pipeline, calculate the performance
metrics and check the posterior distributions for the nuisance parameters. The
latter is done in order to avoid any deviation larger than two standard
deviations from the fiducial values for $b_{1}$, $b_{2}$,
$b_{\mathcal{G}_{2}}$, $b_{\Gamma_{3}}$ $\alpha_{P}$, $\epsilon_{k^{2}}$ and
$\alpha_{B}$ as determined in the previous analyses of Oddo et al. (2020,
2021); Rizzo et al. (2022). Note that the counterterms in real space and in
redshift space are not the same.

We vary the value of $k_{\mathrm{max}}$ jointly for the power spectrum
monopole and quadrupole, as is typically done in analyses of observational
data (Beutler et al., 2014). The scale cut is varied in the range of
$k^{l=0,2}_{\mathrm{max},\,P}\in[0.1,0.3]$ $h$ Mpc-1 and is found to maximize
the FoM while keeping the FoB within one standard deviation at $0.225$ $h$
Mpc-1. This can be seen in the left panel of Figure 3, where blue crosses
correspond to the joint analysis of the lower order power spectrum multipoles.
Keeping this value fixed, we then add the hexadecapole with
$k^{l=4}_{\mathrm{max},\,P}\in[0,0.2]$ $h$ Mpc-1. We expect our model to
describe it less efficiently than the lower order multipoles, therefore the
range of scales included for $P_{4}$ is smaller than for $P_{0}$ and $P_{2}$.
Unsurprisingly the addition of the hexadecapole allows us to constrain
$\tilde{c}_{4}$ and hence improve significantly the constraints on
$\tilde{c}_{2}$, as well as on the other counterterms, though less
significantly. However, we see in the middle panel of Figure 3 no improvement
in the constraints of the IDE parameters, the hexadecapole only slightly
decreases the FoM and increases the FoB. This finding is in agreement with the
conclusions of Carrilho et al. (2021).

Next we include the bispectrum monopole, $B_{0}$. This is shown in the right
panel of Figure 3. We perform the analysis by adding it to $P_{0}+P_{2}$,
since the power spectrum hexadecapole does not contribute to the improvement
of the constraints in the IDE parameters. The results are denoted by the light
green crosses connected by a solid line. We observe an increase in the FoB and
no improvement in the FoM for scales up to $k^{l=0}_{\mathrm{max},\,B}=0.08$
$h$ Mpc-1 for $P_{0}+P_{2}+B_{0}$. However, for
$k^{l=0}_{\mathrm{max},\,B}>0.08$ $h$ Mpc-1 we notice a steep growth in the
FoM as we add more $k$-bins. Clearly, inclusion of the small-scale information
contained in the higher $k$-bins provides better constraining power,
especially since the number of triangles increases significantly with each
bin. The results do not show a bias in any of the parameters of the model, at
least up to the 27-th $k$-bin (corresponding to
$k^{l=0}_{\mathrm{max},\,B}=1.12$ $h$ Mpc-1) from the total 29 bins of our
measurements.

Additionally, we perform a joint analysis with the hexadecapole of the power
spectrum included up to $0.1$ $h$ Mpc-1, which is rather conservative and
common in the literature. This model with all power spectrum multipoles and
$B_{0}$ shows a mild but continuous improvement of the FoM in comparison to
the model with only lower order power spectrum multipoles (see the dark red
triangles in the right panel). Noticeably, adding the power spectrum
hexadecapole in the joint analysis with the bispectrum monopole always results
in a smaller FoB (and thus bias) of the IDE parameters. We check whether the
same effect could be achieved by setting the otherwise unconstrained
$\tilde{c}_{4}$ to zero in the model with the monopole and quadrupole only. As
expected, the constraints of the $P_{0}+P_{2}+B_{0}$ model are not sensitive
to this change. This is shown by the dark-green dots in the right panel, which
overlap with the light-green crosses of the base model with $\tilde{c}_{4}\neq
0$. We therefore conclude that the inclusion of the hexadecapole of the power
spectrum, in combination with the bispectrum monopole, can help (especially at
lower scale cuts) by breaking degeneracies between the counterterms and bias
parameters.

Figure 4: Marginal posterior distributions for the IDE parameters with the
base model and scale cuts as specified in the plot. The power spectrum
monopole and quadrupole analysis is denoted by the dotted blue line,
constraints from all power spectrum multipoles are given by the solid orange
line, the joint analysis of the power spectrum monopole plus quadrupole and
bispectrum monopole is presented with the dotted-dashed light-green line, and
the full joint analysis is shown with the dark-red line. The fact that the
contour shows non-vanishing values for the forbidden regions $A/(1+w)<0$ is an
artifact of the smoothing used for plotting. Thin grey lines correspond to the
fiducial values. Figure 5: Joint analysis of all power spectrum multipoles and
the bispectrum monopole. Left: performance metrics as a function of the
maximum wave-number of the power spectrum monopole and quadrupole with
$k^{l=0}_{\mathrm{max},B}=0.08$ $h$ Mpc-1. Right: performance metrics as a
function of the maximum wave-number of the bispectrum monopole with
$k^{l=0,2}_{\mathrm{max},P}=0.225$ $h$ Mpc-1 and
$k^{l=4}_{\mathrm{max},P}=0.1$ $h$ Mpc-1. Dark-red triangles correspond to the
base model with 13 parameters, olive crosses denote the tidal bias relation,
purple stars represent the $b_{2}$-relation, the dark-cyan squares denote the
combined relation of $b_{2}$ and $b_{\mathcal{G}_{2}}$. The confidence
intervals for FoB are calculated as in Figure 3.

The full posterior plot is shown in Figure 7 of Appendix A. Additionally, we
present in Figure 4 the triangle plot for the marginalized posterior
distribution of the IDE parameters. The characteristic “butterfly” pattern, in
which the posterior distribution is present only in two quadrants – upper-
right and lower-left, is the consequence of the priors on the IDE parameters
with the condition $A/(1+w)\geq 0$, as described in subsection 4.3. The $w-A$
contour shows a strong degeneracy: as argued in section 2, a very negative
value of $w$ leads to a later start of the dark energy dominated epoch, hence
to a longer matter domination which is associated with increased growth of
structures. A negative value of $A$ leads to extra drag in Equation 2, and
thus to an enhancement of the growth of structures in the linear regime (the
opposite is true for the nonlinear regime, as discussed in section 2). On
large scales the drag term acts in the same direction as the gravitational
acceleration, since the velocity field is always aligned with the spatial
gradient of the gravitational potential in the linear regime. The opposite
happens for $w>-1$ and $A>0$. This is in line with what can already be seen in
Figure 1 and its diagonal (from upper left to lower right) pattern in the
growth rate $f$, and is particularly prominent in the upper right corner (for
positive values of $A$ and $(1+w)$). There we observe the same level of
suppression for either a larger $w$ value with $A$ fixed, or for a larger $A$
value with $w$ fixed. It is therefore apparent how the two effects are
strongly degenerate, which results in larger uncertainties on them, especially
in biased cases. Such a conclusion is in agreement with the findings of
Carrilho et al. (2021).

In addition to the base model with 13 parameters, we run the joint analysis
adopting some assumptions that are common in the literature:
$b_{\Gamma_{3}}=0$ and / or $\epsilon_{k^{2}}=0$. Assuming that noise is
scale-independent (as opposed to leaving the $\epsilon_{k^{2}}$ parameter
free), results in no improvement and no bias in the constraints. This is a
consequence of the rescaled error bars, since we know from the real space
analysis with the full volume of the simulations by Oddo et al. (2021) that
the scale-dependent term in shot noise is relevant at larger Fourier modes due
to accounting for additional corrections beyond the assumed one-loop model.
The assumption of $b_{\Gamma_{3}}=0$ does not yield improved constraints, but
it biases the noise parameters $\alpha_{P}$, $\epsilon_{k^{2}}$ and
$\alpha_{B}$ by more than $2\sigma$ with respect to the fiducial values. This
is in line with the joint analysis performed on the same data set in real
space, which measures $b_{\Gamma_{3}}=0$ to be inconsistent with zero at more
than $2\sigma$ (Oddo et al., 2021). The assumption of both conditions results
in bias in the noise parameters and $b_{\mathcal{G}_{2}}$ at more than
$2\sigma$ level, while $\alpha_{B}$ gets biased at more than $3\sigma$ level.
From this we conclude that, when including the bispectrum monopole,
$b_{\Gamma_{3}}$ should not be set to zero.

###  5.2 Bias relations

Although the strong degeneracy between $b_{\mathcal{G}_{2}}$ and
$b_{\Gamma_{3}}$ is broken by the inclusion of the bispectrum monopole in the
analysis, we are still interested in the impact of bias relations for two
reasons: firstly, to reduce the dimensionality of the parameter space;
secondly, we know that both $A$ and $w$ are correlated with $b_{1}$, hence a
tighter constraint on the latter will lead to tighter constraints on the IDE
parameters.

Figure 6: Joint analysis of all power spectrum multipoles with the bispectrum
monopole (dark-red triangles) and quadrupole (purple stars). The performance
metrics, FoM and FoB, are presented as a function of the maximum wave-number
on both multipoles of the bispectrum. The confidence intervals for the FoB are
calculated as in Figure 3.

We investigate the full joint analysis of $P_{0}+P_{2}+P_{4}+B_{0}$ for the
following bias relations:

  1. 1.

the fit from the separate universe simulations presented in Lazeyras et al.
(2016):

$b_{2}(b_{1},b_{\mathcal{G}_{2}})=0.412-2.143b_{1}+0.929b_{1}^{2}+0.008b_{1}^{3}+\frac{4}{3}b_{\mathcal{G}_{2}}\,;$
(22)

  2. 2.

the fit from the excursion set approach of Eggemeier et al. (2020):

$b_{\mathcal{G}_{2}}(b_{1})=0.524-0.547b_{1}+0.046b_{1}^{2}\,;$ (23)

  3. 3.

the combination of both relations: $b_{2}(b_{1},b_{\mathcal{G}_{2}}(b_{1}))$.

These relations have been shown to improve the fits of the power spectrum and
bispectrum in Oddo et al. (2020, 2021); Rizzo et al. (2022). Results for the
individual relations as well as for the combined one are shown in Figure 5 in
terms of the FoM and FoB metrics. We notice that the
$b_{2}(b_{1},b_{\mathcal{G}_{2}})$ relation, denoted by purple stars, does not
impact our constraints both when varying the maximum wave-number for the power
spectrum (left panel) and for the bispectrum (right panel). This can be
explained by a combination of two factors. Firstly, the fact that in our base
model the $b_{2}$-relation always stays consistent with the fiducial value of
$b_{2}$, (that measured from the real space analysis of Oddo et al., 2020,
2021; Rizzo et al., 2022), which implies that we are not introducing any new
information with the bias relation. Then, while a different behaviour has been
seen in the analysis of Oddo et al. (2021) in Fig. 8 with very small error
bars, the fact that ours are “realistic” (larger) for Stage-IV surveys means
that we cannot detect this effect. On the other hand, the tidal bias relation
(olive crosses) and its combination with the $b_{2}$-relation (dark-cyan
squares) lead to a reduction in the constraints on $w$, $A$, $b_{1}$ and
$b_{2}$. This improvement is particularly prominent in the right panel, where
the number of modes in the bispectrum monopole is varied while the scale-range
of power spectrum multipoles is fixed. At $k^{l=0}_{\mathrm{max},B}=0.1$ $h$
Mpc-1 we see the $b_{\mathcal{G}_{2}(b_{1})}$ has a higher value of FoM, which
equals $\simeq 7.5$ and is larger than FoM$\simeq 4$ of the base model with
the same largest Fourier mode. However, for $k^{l=0}_{\mathrm{max},B}>0.1$ $h$
Mpc-1 cases (ii) and (iii) cause a continuous shift of the IDE parameters to
the positive values of $A$ and $(1+w)$. Nonetheless, even for
$k^{l=0}_{\mathrm{max},B}=0.111$ $h$ Mpc-1 all nuisance and IDE parameters
stay within 2 standard deviations of their fiducial values (see the full
posterior distribution in Figure 8), FoM reaches $\simeq 11$ for both
relations. From the lower right panel we see that the combined bias relation
results in the less biased dark energy parameters, in contrast to the tidal-
relation only.

We note that the bias relations applied here are derived in the context of the
$\Lambda$CDM cosmology. From the simulations performed in Baldi & Simpson
(2015), we know that at the redshift of our interest deviations of the halo
mass function from the standard cosmology are negligible. Moreover, the
measurements in our analysis do not include modifications in the gravitational
interaction. Therefore, the application of the bias relations is justified.
Nevertheless, we aim to study the impact of the interactions between dark
matter and dark matter on the bias relations in our future work.

###  5.3 Bispectrum quadrupole

Recent progress in measuring the bispectrum multipoles from simulations (Rizzo
et al., 2022) showed that, while the monopole constrains mostly the bias
parameters $b_{1}$, $b_{2}$ and $b_{\mathcal{G}_{2}}$, the quadrupole can
break the $b_{1}-f$ degeneracy and can constrain the growth factor even in the
case of a bispectrum only analysis. These constraints on the growth rate from
the bispectrum multipoles are weaker than the ones that can be obtained from
the power spectrum multipoles. The argument for the inclusion of higher order
multipoles of the bispectrum is also discussed in Gualdi & Verde (2020);
Gualdi et al. (2021).

However, it raises the question: if including the bispectrum quadrupole
improves the constraints on the growth rate, does it also do so for the IDE
parameters $w$ and $A$? In particular, it is interesting to investigate if
this holds also in the case of Stage-IV-like error bars, such as the ones we
employ here. We run our likelihood pipeline for the base model and find that
indeed, for $k^{l=0,2}_{\mathrm{max},B}>0.09$ $h$ Mpc-1, the addition of the
quadrupole improves the constraints on the dark energy parameters. Our results
are shown in Figure 6, where dark-red triangles denote the power spectrum
multipoles plus the bispectrum monopole, while the addition of $B_{2}$ is
denoted by purple stars. We note that the FoB and FoM as functions of
$k_{\mathrm{max}}$ resemble the trend obtained in the previous section with
the bias relations, with the exception that in the present case the slope of
the FoM in the upper panel is steeper at large wave-numbers. For instance, for
$k^{l=0,2}_{\mathrm{max},B}=0.1$ $h$ Mpc-1, we get FoM $\simeq 9$, while with
the tidal bias relation from the previous section and without the bispectrum
quadrupole we obtained FoM $\simeq 7.5$. The full posterior distribution for
the $k^{l=0,2}_{\mathrm{max},B}=0.111$ $h$ Mpc-1 can be seen in Figure 9. The
increase in the FoM is the result of an improvement in the constraints on the
IDE parameters and the linear bias. Moreover, because of the degeneracy
between the linear bias and the other bias parameters, a tighter constraint on
the former leads to the latter being more constrained too. This effect is more
significant than in the case of the bias relations, and avoids resorting to
such relations, which might not be valid in extended cosmological models. This
is a hint that the bispectrum quadrupole might play a crucial role in
constraining extended cosmological models with spectroscopic galaxy clustering
data.

###  5.4 Forecast for $w$CDM

In addition to IDE, we repeat our analysis for the simplest extension to
$\Lambda$CDM model in which $w=w_{0}$ is constant, but can assume values
different from $-1$. Based on what we have learned with the IDE model we
perform the following analyses:

  1. 1.

base model with $k^{l=0}_{\mathrm{max},B}=0.111$ $h$ Mpc-1;

  2. 2.

combined bias relations with $k^{l=0}_{\mathrm{max},B}=0.095$ $h$ Mpc-1;

  3. 3.

bispectrum monopole and quadrupole with $k^{l=0,2}_{\mathrm{max},B}=0.09$ $h$
Mpc-1.

We obtain the following constraints on the dark energy equation of state
parameter for the three runs respectively: $w=-0.99\pm 0.10$, $-0.98\pm 0.10$,
$-0.99\pm 0.09$. The full posterior is shown in Figure 10. All three cases
result in similar level of uncertainty on the inferred parameters, except on
the shot noise parameter of the bispectrum for which both the bias relations
and the inclusion of $B_{2}$ yield tighter constraints. However, the
improvement with respect to the analysis with the power spectrum only is
apparent, roughly $38\%$, demonstrating once again that the inclusion of the
bispectrum does increase the cosmological information that can be extracted
from the data.

##  6 Conclusions

In this paper we have performed a detailed analysis of the role of the
bispectrum multipoles in constraining an interacting dark energy model. This
is relevant for Stage-IV spectroscopic surveys, as a validation test and a
comprehensive case study for the joint power spectrum and bispectrum analyses
of extended cosmological models.

We modelled the power spectrum and bispectrum multipoles with the EFTofLSS at
one-loop and tree-level, respectively. We computed the linear growth factor
and the logarithmic growth rate by solving the growth equation for the IDE
model. As input data (the fiducial data vector in our likelihood pipeline) we
used a large set of simulations at a single redshift $z=1$, complemented with
a rescaled numerical covariance constructed from $10,000$ mock catalogues. Our
fits used an effective volume that mimics the error budget of Stage-IV
spectroscopic galaxy surveys, and we run our MCMC analyses with a Gaussian
likelihood function.

Using this setup, we first added more observables in sequence and studied them
in terms of their ability to constrain the IDE parameters $w$ and $A$ without
introducing a bias in the dark energy and nuisance parameters of the model.
For the latter, we adopted as fiducial values the best fit values from
previous studies performed in the context of standard $\Lambda$CDM (Oddo et
al., 2020, 2021; Rizzo et al., 2022; Moretti et al., 2022). We found that the
tree-level model for the bispectrum monopole is still valid up to $k\simeq
0.11~{}h~{}{\rm Mpc}^{-1}$, close to the largest value available in our
measurements. In contrast to adding the hexadecapole to the lower-order
multipoles of the power spectrum, the addition of the bispectrum monopole
improved the constraints on the IDE parameters by $\sim 30\%$ for the equation
of state parameter $w$, and $\sim 26\%$ for the coupling parameter $A$. We
also found that the $P_{0}+P_{2}+P_{4}+B_{0}$ combination yields better
constraints and less bias in the inferred parameters than $P_{0}+P_{2}+B_{0}$.

In general, the large number of triangular configurations over which the
bispectrum model has to be evaluated when including such nonlinear scales
results in a rather computationally expensive analysis. Hence, we studied two
approaches, which could allow us to get the same level of constraints but at
smaller scales. The first approach aims at reducing the parameter space, which
in our base case included 13 parameters, by means of bias relations. We tested
the validity of two bias relations, and found that the tidal bias relation and
its combination with the $b_{2}$ relation are the most suitable for this goal.
In the second approach we added the bispectrum quadrupole to our Bayesian
analysis. We found that either applying the bias relations or including the
bispectrum quadrupole allows us to efficiently constrain the IDE parameter at
more moderate values of $k_{\mathrm{max}}$ than in the base model
$P_{0}+P_{2}+P_{4}+B_{0}$. For instance, to achieve the same level of
constraints with FoM $\sim 6$ one can use either of the following: a) use only
the bispectrum monopole in the base model, but include it up to high Fourier
mode of $k_{\mathrm{max},B}^{l=0}=0.11$ $h$ Mpc-1, b) apply the tidal bias
relation in the scale range with the highest mode
$k_{\mathrm{max},B}^{l=0}=0.095$ $h$ Mpc-1, c) add the quadrupole measurements
and evaluate the model up to $k_{\mathrm{max},B}^{l=0,2}=0.09$ $h$ Mpc-1.

We then forecasted the constraining power of Stage-IV spectroscopic galaxy
surveys on the IDE parameters and the dark energy equation of state parameter
in wCDM. At $z=1$ we found $\sigma_{w}=0.08$ and $\sigma_{A}=2.51$ b GeV-1 for
the IDE case, while for the $w$CDM case, we found s $\sigma_{w}=0.1$. The
first values are in a good agreement with the results of Carrilho et al.
(2021), where $\sigma_{w}=0.06$ and $\sigma_{A}=2$ b GeV-1 were found when
only power spectrum multipoles were used in the MCMC analysis. However, in
Carrilho et al. (2021) the simulated data contained more than one order of
magnitude less shot noise in comparison to our case, while the volume of the
measurements was smaller with $V\sim 4$ Gpc3 $h^{-3}$. Additionally, the
covariance in the previous work was computed analytically up to the linear
order, while the one used in this work is based on mocks, which is more
realistic and contributes to larger uncertainties on the inferred parameters.
For the $w$CDM cosmology, our forecasted constraint with the bispectrum and
only one redshift bin outperforms the full-shape analysis (FS) with the
EFTofLSS from the BOSS survey of D’Amico et al. (2021). We stress that the
analysis of BOSS data takes advantage of two redshift bins and three different
sky cuts, with the addition of BAO measurements, while the volume of $6.5$
Gpc3 $h^{-3}$ and number density of $2.2\times 10^{-4}$ $h^{3}$ Mpc-3 are very
comparable with our values. On the other hand, we did not include effects such
as the impact of the window function and Alcock-Paczinsky in our analysis. In
D’Amico et al. (2021), the authors mention that FS analysis with the power
spectrum alone is not sufficient to constraint $w$, due to strong
degeneracies. They measure $w=-1.101^{+0.14}_{-0.11}$, which is still a
stronger limit than the one obtained from Planck2018 data with lensing (Planck
Collaboration et al., 2020): $w=-1.57^{+0.16}_{-0.33}$. These values of $w$
from the power spectrum analysis are in agreement with our forecast, when the
bispectrum is not included (see Figure 10).

A number of interesting questions remain to be answered. The first natural
step would be to perform this analysis with the data from the BOSS survey.
Additionally, it would be interesting to study the time-dependent dark energy
equation of state, for example, in the CPL parametrization (Linder, 2003;
Chevallier & Polarski, 2001), which would allow for different effects of the
interaction. This study requires data at different redshift bins, in order to
constrain $w_{0}$ and $w_{a}$. In the case of IDE, the time-dependent equation
of state promises to have weaker impact on the nonlinear scales than the
constant-$w$ models Baldi & Simpson (2017). Another potential direction of our
study is the inclusion of one-loop effects in the bispectrum model, as was
recently done for $\Lambda$CDM cosmology in D’Amico et al. (2022) and Philcox
et al. (2022). In the latter the authors show that the one-loop corrections
for the bispectrum start to be significant already at
$k_{\mathrm{max},B}^{l=0}>0.1$ $h$ Mpc-1. This is also expected to extend the
reach of the model and therefore its constraining power.

## Acknowledgements

We are grateful to Claudio Dalla Vecchia and Ariel Sánchez for running and
making available the Minerva simulations, performed on the Hydra and Euclid
clusters at the Max Planck Computing and Data Facility in Garching, and to
Pierluigi Monaco, for producing the Pinocchio mocks, run on the GALILEO
cluster at CINECA thanks to an agreement with the University of Trieste. We
thank Andrea Oddo and Emiliano Sefusatti for the effort in building the
initial likelihood code, and for useful conversations. We acknowledge use of
the Cuillin computing cluster, Royal Observatory, University of Edinburgh. We
acknowledge use of open source software: Python (Van Rossum & Drake Jr, 1995;
Hunter, 2007), numpy (Van Der Walt et al., 2011), scipy (Virtanen et al.,
2020), astropy (Astropy Collaboration et al., 2018), corner (Foreman-Mackey,
2016), GetDist (Lewis, 2019). MT’s research is supported by a doctoral
studentship in the School of Physics and Astronomy, University of Edinburgh.
AP is a UK Research and Innovation Future Leaders Fellow [grant MR/S016066/1].
PC and CM’s research is supported by a UK Research and Innovation Future
Leaders Fellowship [grant MR/S016066/1]. For the purpose of open access, the
author has applied a Creative Commons Attribution (CC BY) licence to any
Author Accepted Manuscript version arising from this submission.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Abbott et al. (2021) Abbott T. M. C., et al., 2021, Dark Energy Survey Year 3 Results: Cosmological Constraints from Galaxy Clustering and Weak Lensing (arXiv:2105.13549)
  * Aghamousa et al. (2016) Aghamousa A., et al., 2016, The DESI Experiment Part I: Science,Targeting, and Survey Design (arXiv:1611.00036)
  * Amendola (2000) Amendola L., 2000, Phys. Rev. D, 62, 043511
  * Amendola & Tsujikawa (2020) Amendola L., Tsujikawa S., 2020, J. Cosmology Astropart. Phys., 2020, 020
  * Assassi et al. (2014) Assassi V., Baumann D., Green D., Zaldarriaga M., 2014, J. Cosmology Astropart. Phys., 2014, 056
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Baldauf et al. (2012) Baldauf T., Seljak U., Desjacques V., McDonald P., 2012, Phys. Rev. D, 86, 083540
  * Baldauf et al. (2013a) Baldauf T., Seljak U., Smith R. E., Hamaus N., Desjacques V., 2013a, Phys. Rev. D, 88, 083507
  * Baldauf et al. (2013b) Baldauf T., Seljak U., Smith R. E., Hamaus N., Desjacques V., 2013b, Phys. Rev. D, 88, 083507
  * Baldauf et al. (2015a) Baldauf T., Mirbabayi M., Simonović M., Zaldarriaga M., 2015a, Phys. Rev. D, 92, 043514
  * Baldauf et al. (2015b) Baldauf T., Mercolli L., Mirbabayi M., Pajer E., 2015b, J. Cosmology Astropart. Phys., 2015, 007
  * Baldi & Simpson (2015) Baldi M., Simpson F., 2015, MNRAS, 449, 2239
  * Baldi & Simpson (2017) Baldi M., Simpson F., 2017, MNRAS, 465, 653
  * Barros et al. (2019) Barros B. J., Amendola L., Barreiro T., Nunes N. J., 2019, J. Cosmology Astropart. Phys., 2019, 007
  * Baumann et al. (2012) Baumann D., Nicolis A., Senatore L., Zaldarriaga M., 2012, J. Cosmology Astropart. Phys., 2012, 051
  * Bayes (1763) Bayes T., 1763, Phil. Trans. of the Royal Soc. of London, 53, 370
  * Beltrán Jiménez et al. (2021) Beltrán Jiménez J., Bettoni D., Figueruelo D., Teppa Pannia F. A., Tsujikawa S., 2021, Phys. Rev. D, 104, 103503
  * Bernardeau et al. (2002) Bernardeau F., Colombi S., Gaztañaga E., Scoccimarro R., 2002, Phys. Rep., 367, 1
  * Bertone et al. (2005) Bertone G., Hooper D., Silk J., 2005, Physics Reports, 405, 279
  * Beutler et al. (2014) Beutler F., et al., 2014, MNRAS, 443, 1065
  * Beutler et al. (2017) Beutler F., et al., 2017, MNRAS, 466, 2242
  * Blas et al. (2011) Blas D., Lesgourgues J., Tram T., 2011, JCAP, 07, 034
  * Blas et al. (2016) Blas D., Garny M., Ivanov M. M., Sibiryakov S., 2016, J. Cosmology Astropart. Phys., 2016, 052
  * Carrasco et al. (2012) Carrasco J. J. M., Hertzberg M. P., Senatore L., 2012, Journal of High Energy Physics, 2012, 82
  * Carrilho et al. (2021) Carrilho P., Moretti C., Bose B., Markovič K., Pourtsidou A., 2021, J. Cosmology Astropart. Phys., 2021, 004
  * Chan et al. (2012) Chan K. C., Scoccimarro R., Sheth R. K., 2012, Phys. Rev. D, 85, 083509
  * Chevallier & Polarski (2001) Chevallier M., Polarski D., 2001, International Journal of Modern Physics D, 10, 213
  * Clemson et al. (2012) Clemson T., Koyama K., Zhao G.-B., Maartens R., Valiviita J., 2012, Phys. Rev. D, 85, 043007
  * Clifton et al. (2012) Clifton T., Ferreira P. G., Padilla A., Skordis C., 2012, Phys. Rep., 513, 1
  * Copeland et al. (2006) Copeland E. J., Sami M., Tsujikawa S., 2006, International Journal of Modern Physics D, 15, 1753
  * Crocce & Scoccimarro (2008) Crocce M., Scoccimarro R., 2008, Phys. Rev. D, 77, 023533
  * D’Amico et al. (2021) D’Amico G., Senatore L., Zhang P., 2021, J. Cosmology Astropart. Phys., 2021, 006
  * D’Amico et al. (2022) D’Amico G., Donath Y., Lewandowski M., Senatore L., Zhang P., 2022, arXiv e-prints, p. arXiv:2206.08327
  * Desjacques et al. (2018) Desjacques V., Jeong D., Schmidt F., 2018, Phys. Rep., 733, 1
  * Di Valentino et al. (2017) Di Valentino E., Melchiorri A., Mena O., 2017, Phys. Rev. D, 96, 043503
  * Di Valentino et al. (2020a) Di Valentino E., Melchiorri A., Mena O., Vagnozzi S., 2020a, Physics of the Dark Universe, 30, 100666
  * Di Valentino et al. (2020b) Di Valentino E., Melchiorri A., Mena O., Vagnozzi S., 2020b, Phys. Rev. D, 101, 063502
  * Eggemeier et al. (2020) Eggemeier A., Scoccimarro R., Crocce M., Pezzotta A., Sánchez A. G., 2020, Physical Review D, 102, 103530
  * Eggemeier et al. (2021) Eggemeier A., Scoccimarro R., Smith R. E., Crocce M., Pezzotta A., Sánchez A. G., 2021, Phys. Rev. D, 103, 123550
  * Fang et al. (2017) Fang X., Blazek J. A., McEwen J. E., Hirata C. M., 2017, J. Cosmology Astropart. Phys., 2017, 030
  * Farrar & Peebles (2004) Farrar G. R., Peebles P. J. E., 2004, ApJ, 604, 1
  * Foreman-Mackey (2016) Foreman-Mackey D., 2016, The Journal of Open Source Software, 1, 24
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Fry & Gaztanaga (1993) Fry J. N., Gaztanaga E., 1993, ApJ, 413, 447
  * Gagrani & Samushia (2017) Gagrani P., Samushia L., 2017, MNRAS, 467, 928
  * Ginzburg et al. (2017) Ginzburg D., Desjacques V., Chan K. C., 2017, Phys. Rev. D, 96, 083528
  * Gómez-Valent et al. (2020) Gómez-Valent A., Pettorino V., Amendola L., 2020, Phys. Rev. D, 101, 123513
  * Grieb et al. (2016) Grieb J. N., Sánchez A. G., Salazar-Albornoz S., Dalla Vecchia C., 2016, MNRAS, 457, 1577
  * Gualdi & Verde (2020) Gualdi D., Verde L., 2020, J. Cosmology Astropart. Phys., 2020, 041
  * Gualdi et al. (2021) Gualdi D., Gil-Marín H., Verde L., 2021, J. Cosmology Astropart. Phys., 2021, 008
  * Hunter (2007) Hunter J. D., 2007, Computing In Science & Engineering, 9, 90
  * Ivanov et al. (2022) Ivanov M. M., Philcox O. H. E., Nishimichi T., Simonović M., Takada M., Zaldarriaga M., 2022, Phys. Rev. D, 105, 063512
  * Ivezić et al. (2019) Ivezić v., et al., 2019, Astrophys. J., 873, 111
  * Kase & Tsujikawa (2020) Kase R., Tsujikawa S., 2020, Physics Letters B, 804, 135400
  * Laureijs et al. (2011) Laureijs R., et al., 2011, Euclid Definition Study Report (arXiv:1110.3193)
  * Lazeyras et al. (2016) Lazeyras T., Wagner C., Baldauf T., Schmidt F., 2016, J. Cosmology Astropart. Phys., 2016, 018
  * Lesgourgues (2011) Lesgourgues J., 2011, arXiv e-prints, p. arXiv:1104.2932
  * Lewis (2019) Lewis A., 2019, GetDist: a Python package for analysing Monte Carlo samples (arXiv:1910.13970)
  * Lewis & Challinor (2011) Lewis A., Challinor A., 2011, CAMB: Code for Anisotropies in the Microwave Background, Astrophysics Source Code Library, record ascl:1102.026 (ascl:1102.026)
  * Linder (2003) Linder E. V., 2003, Phys. Rev. Lett., 90, 091301
  * Lippich et al. (2019) Lippich M., et al., 2019, MNRAS, 482, 1786
  * MacKay (2003) MacKay D. J. C., 2003, Information Theory, Inference, and Learning Algorithms. Copyright Cambridge University Press
  * Markovic et al. (2019) Markovic K., Bose B., Pourtsidou A., 2019, Open J. Astrophys., 2, 13
  * McDonald & Roy (2009) McDonald P., Roy A., 2009, J. Cosmology Astropart. Phys., 2009, 020
  * McEwen et al. (2016) McEwen J. E., Fang X., Hirata C. M., Blazek J. A., 2016, J. Cosmology Astropart. Phys., 2016, 015
  * Monaco et al. (2002) Monaco P., Theuns T., Taffoni G., 2002, MNRAS, 331, 587
  * Moretti et al. (2022) Moretti C., Rizzo F., Pardede K., Oddo A., Sefusatti E., Porciani C., Monaco P., 2022, in preparation
  * Munari et al. (2017) Munari E., Monaco P., Sefusatti E., Castorina E., Mohammad F. G., Anselmi S., Borgani S., 2017, MNRAS, 465, 4658
  * Nunes et al. (2022) Nunes R. C., Vagnozzi S., Kumar S., Di Valentino E., Mena O., 2022, Phys. Rev. D, 105, 123506
  * Oddo et al. (2020) Oddo A., Sefusatti E., Porciani C., Monaco P., Sánchez A. G., 2020, J. Cosmology Astropart. Phys., 2020, 056
  * Oddo et al. (2021) Oddo A., Rizzo F., Sefusatti E., Porciani C., Monaco P., 2021, J. Cosmology Astropart. Phys., 2021, 038
  * Osato et al. (2019) Osato K., Nishimichi T., Bernardeau F., Taruya A., 2019, Phys. Rev. D, 99, 063530
  * Paech et al. (2017) Paech K., Hamaus N., Hoyle B., Costanzi M., Giannantonio T., Hagstotz S., Sauerwein G., Weller J., 2017, MNRAS, 470, 2566
  * Pan et al. (2019) Pan S., Yang W., Di Valentino E., Saridakis E. N., Chakraborty S., 2019, Phys. Rev. D, 100, 103520
  * Perivolaropoulos & Skara (2021) Perivolaropoulos L., Skara F., 2021, arXiv e-prints, p. arXiv:2105.05208
  * Perko et al. (2016) Perko A., Senatore L., Jennings E., Wechsler R. H., 2016, arXiv e-prints, p. arXiv:1610.09321
  * Philcox & Ivanov (2022) Philcox O. H. E., Ivanov M. M., 2022, Phys. Rev. D, 105, 043517
  * Philcox et al. (2022) Philcox O. H. E., Ivanov M. M., Cabass G., Simonović M., Zaldarriaga M., Nishimichi T., 2022, arXiv e-prints, p. arXiv:2206.02800
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Pourtsidou & Tram (2016) Pourtsidou A., Tram T., 2016, Phys. Rev. D, 94, 043518
  * Pourtsidou et al. (2013) Pourtsidou A., Skordis C., Copeland E. J., 2013, Phys. Rev. D, 88, 083505
  * Rizzo et al. (2022) Rizzo F., Moretti C., Pardede K., Eggemeier A., Oddo A., Sefusatti E., Porciani C., Monaco P., 2022, arXiv e-prints, p. arXiv:2204.13628
  * Sánchez et al. (2013) Sánchez A. G., et al., 2013, MNRAS, 433, 1202
  * Scoccimarro et al. (1999a) Scoccimarro R., Couchman H. M. P., Frieman J. A., 1999a, ApJ, 517, 531
  * Scoccimarro et al. (1999b) Scoccimarro R., Couchman H. M. P., Frieman J. A., 1999b, ApJ, 517, 531
  * Seljak et al. (2009) Seljak U., Hamaus N., Desjacques V., 2009, Phys. Rev. Lett., 103, 091303
  * Sellentin & Heavens (2016) Sellentin E., Heavens A. F., 2016, MNRAS, 456, L132
  * Senatore & Zaldarriaga (2014) Senatore L., Zaldarriaga M., 2014, arXiv e-prints, p. arXiv:1409.1225
  * Senatore & Zaldarriaga (2015) Senatore L., Zaldarriaga M., 2015, J. Cosmology Astropart. Phys., 2015, 013
  * Simpson (2010) Simpson F., 2010, Phys. Rev. D, 82, 083505
  * Skordis et al. (2015) Skordis C., Pourtsidou A., Copeland E. J., 2015, Phys. Rev. D, 91, 083537
  * Smith et al. (2007) Smith R. E., Scoccimarro R., Sheth R. K., 2007, Phys. Rev. D, 75, 063512
  * Smith et al. (2008) Smith R. E., Sheth R. K., Scoccimarro R., 2008, Phys. Rev. D, 78, 023523
  * Sotiriou & Faraoni (2010) Sotiriou T. P., Faraoni V., 2010, Reviews of Modern Physics, 82, 451
  * Spergel et al. (2015) Spergel D., et al., 2015, arXiv preprint arXiv:1503.03757
  * Spurio Mancini & Pourtsidou (2022) Spurio Mancini A., Pourtsidou A., 2022, MNRAS, 512, L44
  * Tamanini (2015) Tamanini N., 2015, Phys. Rev. D, 92, 043524
  * Taruya et al. (2010) Taruya A., Nishimichi T., Saito S., 2010, Phys. Rev. D, 82, 063522
  * Tegmark (1997) Tegmark M., 1997, Phys. Rev. Lett., 79, 3806
  * Väliviita et al. (2008) Väliviita J., Majerotto E., Maartens R., 2008, J. Cosmology Astropart. Phys., 2008, 020
  * Van Der Walt et al. (2011) Van Der Walt S., Colbert S. C., Varoquaux G., 2011, preprint,  (arXiv:1102.1523)
  * Van Rossum & Drake Jr (1995) Van Rossum G., Drake Jr F. L., 1995, Python reference manual. Centrum voor Wiskunde en Informatica Amsterdam
  * Villaescusa-Navarro et al. (2020) Villaescusa-Navarro F., et al., 2020, ApJS, 250, 2
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Vlah & White (2019) Vlah Z., White M., 2019, JCAP, 03, 007
  * Vlah et al. (2015) Vlah Z., White M., Aviles A., 2015, JCAP, 09, 014
  * Wang (2008) Wang Y., 2008, Phys. Rev. D, 77, 123525
  * Xia (2009) Xia J.-Q., 2009, Phys. Rev. D, 80, 103514
  * Yang et al. (2021) Yang W., Pan S., Saló L. A., de Haro J., 2021, Phys. Rev. D, 103, 083520
  * Yankelevich & Porciani (2019) Yankelevich V., Porciani C., 2019, MNRAS, 483, 2078
  * d’Amico et al. (2020) d’Amico G., Gleyzes J., Kokron N., Markovic K., Senatore L., Zhang P., Beutler F., Gil-Marín H., 2020, J. Cosmology Astropart. Phys., 2020, 005
  * van de Bruck & Thomas (2019) van de Bruck C., Thomas C. C., 2019, Phys. Rev. D, 100, 023515

##  Appendix A Full posterior distributions

We collect here the full posterior distributions for the different runs
described in the text.

Figure 7: Posterior distributions for the base model with the scale cuts given
in the triangle plot. The power spectrum monopole and quadrupole analysis is
denoted by the dotted blue line, constraints from all power spectrum
multipoles are given by the solid orange line, joint analysis of power
spectrum monopole plus quadrupole and bispectrum monopole is presented by the
dotted-dashed light-green line, and the full joint analysis is shown by the
dark-red contours. Thin grey lines correspond to the fiducial values known
from the fiducial cosmology and previous analysis of this data set in Oddo et
al. (2020, 2021); Rizzo et al. (2022) . Figure 8: Posterior distributions for
the $P_{0}+P_{2}+P_{4}+B_{0}$ models with bias relations at the scale cuts
given on the triangle plot. The joint analysis with the base model is shown by
the dark-red dotted-dashed lines, olive dotted lines denote the tidal bias
relation, purple solid lines represent the $b_{2}$-relation, the dark-cyan
contours denote the combined relation of $b_{2}$ and $b_{\mathcal{G}_{2}}$.
Thin grey lines correspond to the fiducial values known from the fiducial
cosmology and previous analysis of this data set in Oddo et al. (2020, 2021);
Rizzo et al. (2022). Figure 9: Posterior distributions for the base model with
and without the bispectrum quadrupole at the scale cuts given on the triangle
plot. The joint analysis with the bispectrum monopole is shown by the dark-red
dotted-dashed lines, the joint analysis with the bispectrum monopole and
quadrupole is denoted by the purple contours. Thin grey lines correspond to
the fiducial values known from the fiducial cosmology and previous analysis of
this data set in Oddo et al. (2020, 2021); Rizzo et al. (2022). Figure 10:
Posterior distributions for the $w$CDM model for four models at the scale cuts
given on the triangle plot. The joint analysis with the bispectrum monopole is
shown by the solid dark-red lines, the joint analysis with the bispectrum
monopole with the combined bias relation is given by the dotted-dashed dark-
cyan lines, the joint analysis with the bispectrum monopole and quadrupole is
denoted by the purple contours, constraints from all power spectrum multipoles
are shown by the dotted orange line. Thin grey lines correspond to the
fiducial values known from the fiducial cosmology and previous analysis of
this data set in Oddo et al. (2020, 2021); Rizzo et al. (2022).