case, so this result might be fortuitous and should be interpreted carefully.
Moreover, the fact that the validity ranges are consistently lower than for
the one-loop model, despite having a larger parameter space, suggests that the
higher-derivative terms mainly act as a convenient way of absorbing the bias
loop corrections. This conclusion is also supported by various higher-
derivative parameters attaining a significant dependence on
$k_{\mathrm{max}}$.

####  V.2.2 Improvements over the power spectrum

After finding the validity ranges for the joint power spectrum and bispectrum
analysis, we can turn to the question of how much extra information on $A_{s}$
we gain compared to the power spectrum alone. For that purpose we show in Fig.
13 the $68\,\%$ and $95\,\%$ confidence limits on $A_{s}$ and $b_{1}$
(rescaled back to the original linear bias parameter) obtained for these two
cases, where we have used the one-loop bispectrum model with fixed NLE bias
operators and similarly assumed coevolution for $\gamma_{21}$ (the only
contributing NLE operator) in the individual power spectrum analysis. The
constraints are respectively plotted at the maximum valid $k_{\mathrm{max}}$
value of the joint fit as determined from the data shown in Fig. 12. For each
of the samples the power spectrum contours (gray) tightly follow the expected
degeneracy direction, $b_{1}\,\sqrt{A_{s}}=\mathrm{const.}$, which prevents
any better constraint on $A_{s}$. Inclusion of an increasing number of
bispectrum triangle configurations (larger $k_{\mathrm{max}}$), however,
breaks this degeneracy more and more efficiently, which leads to substantial
reductions in the uncertainties. As indicated for $A_{s}$ by the green arrows
in the plot, these range between factors of 2.2 for HALO1 and 5.7 for HALO3.
Moreover, we have checked that enlarging the cutoff scale for the power
spectrum fit, e.g. to $k_{\mathrm{max}}=0.3\,h/\mathrm{Mpc}$, has little to no
impact on this outcome, and $\Delta A_{s}$ generally shrinks by at least a
factor of 3 (except for HALO1).

Nevertheless, as mentioned at the beginning of this section, these results are
not fully representative in the sense that this analysis was performed in real
opposed to redshift space, where the combination of power spectrum multipoles
already helps in disentangling $b_{1}$ from $A_{s}$ (assuming $\Lambda$CDM,
otherwise from a combination of $A_{s}$ and the growth rate of structure). In
addition, correlations between different triangle configurations that we have
ignored in the covariance matrix of the bispectrum, could further diminish the
improvement factors quoted in Fig. 13. On the other hand, similar correlations
in the power spectrum part of the covariance matrix would act in the opposite
direction and the cross-covariance between both statistics, which was also not
included in our analysis, has been shown to have a strong _positive_ effect on
constraining the amplitude of fluctuations Sefusatti  _et al._ (2006); Byun
_et al._ (2017). Our results are therefore at the very least indicative of the
large potential of the bispectrum, and in particular we have shown here for
the first time that this potential is not heavily diminished by the increased
parameter space that is required for describing the bispectrum on mildly
nonlinear scales.

##  VI Conclusions

We have analyzed real-space measurements of the power spectrum and bispectrum
from a large set of galaxy and halo mock catalogs, and performed fits of their
joint likelihood using models build from perturbation theory. In particular,
we have accounted here for the first time for all next-to-leading (one-loop)
corrections due to galaxy bias in the bispectrum, which puts its description
on equal footing with models routinely used in the analysis of the power
spectrum Gil-Marín  _et al._ (2015b); Beutler  _et al._ (2017); Grieb  _et
al._ (2017); Ivanov  _et al._ (2019); D’Amico  _et al._ (2019); Tröster  _et
al._ (2020), and allowed us to include more nonlinear triangle shapes than in
previous studies Feldman  _et al._ (2001); Saito  _et al._ (2014); Gil-Marín
_et al._ (2017b, a); Slepian  _et al._ (2017a); Pearson and Samushia (2018);
Sugiyama  _et al._ (2021). Since modeling the nonlinear corrections requires
the introduction of additional bias parameters, and thus a larger number of
nuisance parameters that need to be marginalized over, our analysis pursued
two main goals: 1) identifying whether the information from nonlinear scales
in the bispectrum offsets the penalty of the extra parameters, and 2) testing
the applicability and universality of coevolution and peak-background split
(PBS) relations in order to reduce the bias parameter space.

###  VI.1 At fixed cosmology

Holding initially all cosmological parameters fixed at their fiducial values,
we have found that the one-loop galaxy bias corrections to the bispectrum
become relevant at scales $k\gtrsim 0.17\,h/\mathrm{Mpc}$. This is largely
consistent across our various samples, for each of which we adopted
statistical uncertainties corresponding to an effective survey volume of
$6\,(\mathrm{Gpc}/h)^{3}$, and also agrees qualitatively with the scale at
which analogous corrections start to significantly affect the power spectrum.
The posteriors of both statistics remain fully compatible and based on the
goodness-of-fit, as well as the accurate recovery of the linear bias parameter
$b_{1}$, we have quantitatively ascertained that joint fits of the power
spectrum and bispectrum can be applied over the entire tested range of scales,
i.e. $k_{\mathrm{max}}=0.3\,h/\mathrm{Mpc}$, for most of the samples. Notably,
we gain a factor of 1.5 to 2 when contrasting the constraining power on
$b_{1}$ at the maximum scale at which our metrics regard the model as valid
with that of a model excluding loop corrections. This indicates a clear
benefit from the information on nonlinear scales, despite the larger nuisance
parameter space.

From the same analysis we have carefully extracted measurements of all bias
parameters up to fourth order, and for the parameters shared between the power
spectrum and bispectrum — the nonlinear and tidal biases $b_{2}$ and
$\gamma_{2}$, as well as the nonlocal parameter $\gamma_{21}$ — these range
among the most stringent that have been reported so far. Comparing these
measurements with the commonly applied local Lagrangian approximations for
$\gamma_{2}$ and $\gamma_{21}$, we have found clear deviations that are in
agreement with the results of other, independent analyses Lazeyras and Schmidt
(2018); Abidi and Baldauf (2018). This suggests the presence of Lagrangian
tidal bias and, especially, the importance of including its effect in joint
analyses of the power spectrum and bispectrum (see also Oddo  _et al._
(2020)). At the high-mass end of our halo samples (large $b_{1}$), we have
also detected a departure from the excursion-set prediction of Sheth  _et al._
(2013) for $\gamma_{2}$, which was working well when considering only data
from the power spectrum Eggemeier  _et al._ (2020). On the other hand, for all
of our samples $\gamma_{21}$ is in excellent agreement with the more general
coevolution relation, which takes the non-zero Lagrangian tidal bias into
account. The PBS relation of Lazeyras  _et al._ (2016) for $b_{2}$ (calibrated
for dark matter halos) matches very well with the results obtained from the
halo samples, but disagrees strongly with those from all three galaxy
catalogs, which we have attributed to the nonlinearity of the $b_{2}(b_{1})$
relation.

Constraints on higher-order bias parameters are generally much weaker and
individually often consistent with zero, since they only enter as nonlinear
corrections to the bispectrum. However, an analysis of the parameter
correlation matrix revealed that there are two very well constrained
combinations of third-order parameters — one of them in good agreement with
the PBS relation for $b_{3}$ — which implies that they cannot simply be
dropped from the model.

Taking these results as guidance in order to find suitable reductions of the
parameter space, we have shown that the application of coevolution relations
to all galaxy bias contributions that are generated by nonlocal gravitational
evolution (NLE bias operators corresponding to $\gamma_{21}$ and all four
fourth-order parameters) presents an exceptionally good approximation for our
data. The estimated validity ranges are not diminished, while the
uncertainties on $b_{1}$ shrink by $15$ \- $20\,\%$. Further improvements can
be achieved by fixing the scale-dependent stochastic component in the power
spectrum using the empirical relation provided in Eggemeier  _et al._ (2020).
Hence, with these approximations, modeling the bispectrum at one-loop order
only takes three more parameters in excess of those already required by the
one-loop power spectrum.

###  VI.2 Varying $A_{s}$

In a second part of our analysis we have allowed the amplitude of
fluctuations, $A_{s}$, to vary alongside the bias parameters. Judging the
model performance in that case by goodness-of-fit and an accurate estimation
of $A_{s}$, we have found slightly reduced validity ranges for both the tree-
level and one-loop bias descriptions compared to the results at fixed
cosmology, although the latter remains reliable up to at least
$k_{\mathrm{max}}=0.25\,h/\mathrm{Mpc}$ for the majority of samples. We have
further shown that the penalty on $A_{s}$ from marginalizing over bias
parameters is more severe than on $b_{1}$, leading for some samples to
comparable or even worse constraints on $A_{s}$ from the one-loop model than
what we obtain from the tree-level bias description.

However, we have also demonstrated that the coevolution relations for NLE
operators continue to hold and their application turns the situation clearly
in favor of the one-loop model. In particular, it improves the statistical
uncertainties on $A_{s}$ by factors of 1.5 to 2.4 compared to an analysis
constrained to larger scales using the tree-level model. In comparison with
the power spectrum alone the joint analysis shrinks the constraints on $A_{s}$
by factors of $\sim 4$ \- 6 for most samples, and evidently highlights its
value for breaking the prominent $b_{1}$ \- $A_{s}$ degeneracy. Even though
this degeneracy is already partly broken by the power spectrum monopole and
quadrupole in a redshift-space analysis, likely leading to less pronounced
improvements than those quoted above, the bispectrum will still be
instrumental in disentangling the effects of $A_{s}$ and the growth rate of
structures, which we will consider in future work.

###  VI.3 On scale-dependent stochastic and higher-derivative contributions

We have also analyzed the effect of including scale-dependent stochastic
(three extra parameters) and higher-derivative (five extra parameters)
contributions, both to the tree-level and one-loop galaxy bias model of the
bispectrum. While the performance of the tree-level model is not significantly
enhanced when augmented by scale-dependent stochasticity, we have seen that it
leads to slightly improved results and less $k_{\mathrm{max}}$ dependent
fitting parameters for the one-loop model when applied to the more heavily
biased samples. This is in accordance with the same finding for the power
spectrum presented in Eggemeier  _et al._ (2020).

At fixed cosmology higher-derivative terms can extend the validity of the
tree-level bias description far into the nonlinear regime for multiple
samples, which implies that they are able to efficiently absorb the missing
bias loop corrections. However, this comes at the expense of a significantly
increased uncertainty on $b_{1}$, such that in this case there is little
benefit in going beyond the large-scale analysis. Interestingly, the same
effect is not as severe with varying $A_{s}$, giving minimum measurement
errors on $A_{s}$ that are larger than those obtained by the one-loop model by
only a few per-cent for the galaxy samples. On the contrary, for the halo
samples we did not find the same improvements from the inclusion of higher-
derivative contributions, so we caution that their ability to absorb and
marginalize over loop corrections may well depend on the particular sample
under consideration.

###### Acknowledgements.

The authors thank Cristiano Porciani and Emiliano Sefusatti for useful
discussions. AE acknowledges support from the European Research Council (grant
number ERC-StG-716532-PUNCA), RES from the STFC (grant number ST/P000525/1,
ST/T000473/1), MC and AP from the Spanish Ministry of Science MINECO under
grant PGC2018-102021, and AGS from the Excellence Cluster ORIGINS, which is
funded by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) under Germany’s Excellence Strategy - EXC-2094 - 390783311. This
research made use of matplotlib, a Python library for publication quality
graphics Hunter (2007).

##  APPENDIX A FOURTH-ORDER BIAS OPERATORS

The one-loop bispectrum receives contributions from four fourth-order bias
operators, whose expressions are collected in the following. As explained in
Sec. II.1, each of them can be written in terms of Galileons and higher-order
LPT potentials:

$\displaystyle\mathrm{i)}\quad$ $\displaystyle\delta\,{\cal
G}_{2}(\varphi_{2},\varphi_{1})\,,$ (92) $\displaystyle\mathrm{ii)}\quad$
$\displaystyle{\cal G}_{3}(\varphi_{2},\varphi_{1},\varphi_{1})\,,$ (93)
$\displaystyle\mathrm{iii)}\quad$ $\displaystyle{\cal
G}_{2}(\varphi_{2},\varphi_{2})\,,$ (94) $\displaystyle\mathrm{iv)}\quad$
$\displaystyle{\cal G}_{2}(\varphi_{3},\varphi_{1})=\frac{1}{18}{\cal
G}_{2}(\varphi_{3}^{(a)},\varphi_{1})+\frac{5}{42}{\cal
G}_{2}(\varphi_{3}^{(b)},\varphi_{1})$ $\displaystyle\hskip
56.00014pt+\frac{1}{14}\nabla_{i}\left(\bm{\nabla}\times\bm{A}_{3}\right)_{j}\nabla_{ij}\varphi_{1}\,,$
(95)

where the last operator is a combination of the two scalar third-order
potentials $\varphi_{3}^{(a)}$ and $\varphi_{3}^{(b)}$, and the vector
potential $\bm{A}_{3}$ Eggemeier  _et al._ (2019). They are defined as

$\displaystyle\nabla^{2}\,\varphi_{3}^{(a)}$ $\displaystyle=-{\cal
G}_{3}(\varphi_{1})\,,$ (96) $\displaystyle\nabla^{2}\,\varphi_{3}^{(b)}$
$\displaystyle=-{\cal G}_{2}(\varphi_{2},\varphi_{1})\,,$ (97)
$\displaystyle\nabla^{2}\bm{A}_{3}$
$\displaystyle=-\hat{e}_{i}\,\epsilon_{ijk}\,\left(\nabla_{jl}\,\varphi_{1}\right)\,\left(\nabla_{kl}\,\varphi_{2}\right)\,,$
(98)

where $\epsilon_{ijk}$ is the fully anti-symmetric Levi-Civita symbol and
$\hat{e}_{i}$ denotes the unit vector in $i$-direction.

In Fourier space we can express any of these operators as the integral

$\begin{split}{\cal
O}^{(4)}(\bm{k})=\;&(2\pi)^{3}\int_{\bm{k}_{1},\ldots,\bm{k}_{4}}\delta_{D}(\bm{k}-\bm{k}_{1234})\,\\\
&\times\,{\cal K}^{(4)}_{{\cal
O}}(\bm{k}_{1},\ldots,\bm{k}_{4})\prod_{i=1}^{4}\delta_{L}(\bm{k}_{i})\,,\end{split}$
(99)

and the kernel functions corresponding to the four operators are given by
Eggemeier  _et al._ (2019)

$\displaystyle{\cal K}^{(4)}_{\delta{\cal
G}_{2}(\varphi_{2},\varphi_{1})}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4})$
$\displaystyle=\frac{1}{12}\Big{[}K(\bm{k}_{1},\bm{k}_{23})\,K(\bm{k}_{2},\bm{k}_{3})+\text{sym.}(12)\Big{]}\,,$
(100) $\displaystyle{\cal K}^{(4)}_{{\cal
G}_{3}(\varphi_{2},\varphi_{1},\varphi_{1})}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4})$
$\displaystyle=\frac{1}{6}\Big{[}L(\bm{k}_{1},\bm{k}_{2},\bm{k}_{34})\,K(\bm{k}_{3},\bm{k}_{4})+\text{sym.}(6)\Big{]}\,,$
(101) $\displaystyle{\cal K}^{(4)}_{{\cal
G}_{2}(\varphi_{2},\varphi_{2})}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4})$
$\displaystyle=\frac{1}{3}\Big{[}K(\bm{k}_{12},\bm{k}_{34})\,K(\bm{k}_{1},\bm{k}_{2})\,K(\bm{k}_{3},\bm{k}_{4})+\text{sym.}(3)\Big{]}\,,$
(102) $\displaystyle{\cal K}^{(4)}_{{\cal
G}_{2}(\varphi_{3},\varphi_{1})}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4})$
$\displaystyle=\frac{1}{12}\left[\frac{1}{18}K(\bm{k}_{1},\bm{k}_{234})\left(\frac{15}{7}\,K(\bm{k}_{23},\bm{k}_{4})\,K(\bm{k}_{2},\bm{k}_{3})-L(\bm{k}_{2},\bm{k}_{3},\bm{k}_{4})\right)\right.$
$\displaystyle\hskip
11.99998pt\left.+\frac{1}{14}\Big{(}M(\bm{k}_{1},\bm{k}_{23},\bm{k}_{4},\bm{k}_{234})-M(\bm{k}_{1},\bm{k}_{234},\bm{k}_{23},\bm{k}_{4})\Big{)}K(\bm{k}_{2},\bm{k}_{3})+\text{sym.}(12)\right]\,,$
(103)

where $\mathrm{sym.}(n)$ stands for the $n-1$ terms that need to be added in
order to symmetrize the expressions over the four participating wave vectors.
In Eq. (103) we have further made use of the kernel function

$M(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4})\equiv\frac{(\bm{k}_{1}\cdot\bm{k}_{2})\,(\bm{k}_{2}\cdot\bm{k}_{3})\,(\bm{k}_{3}\cdot\bm{k}_{4})\,(\bm{k}_{4}\cdot\bm{k}_{1})}{(k_{1}\,k_{2}\,k_{3}\,k_{4})^{2}}\,.$
(104)

Finally, when computing the loop corrections for the two-point propagator (see
Eq. 9), we need to make sure to only include the “finite” part of the kernels
and subtract any potential sensitivities to the nonlinear regime in the large-
scale limit (which are absorbed by the bias parameters). This only occurs for
the first operator (Eq. 92) and one can show that

$\begin{split}{\cal K}_{\delta{\cal
G}_{2}(\varphi_{2},\varphi_{1})}^{(4,\mathrm{lim})}(\bm{k_{1}},\bm{k}_{2},\bm{q})&\equiv\lim_{\bm{q}\to
0}{\cal K}_{\delta{\cal
G}_{2}(\varphi_{2},\varphi_{1})}^{(4)}(\bm{k}_{1},\bm{k}_{2},\bm{q},-\bm{q})\\\
&=\frac{1}{3}K(\bm{k}_{1},\bm{q})\,K(\bm{k}_{2},\bm{q})\\\ &\hskip
10.00002pt+\frac{1}{6}K(\bm{k}_{1},\bm{k}_{2})\,K(\bm{k}_{1}+\bm{k}_{2},\bm{q})\,,\end{split}$
(105)

and so we define

${\cal K}_{\delta{\cal G}_{2}(\varphi_{2},\varphi_{1})}^{(4,F)}\equiv{\cal
K}_{\delta{\cal G}_{2}(\varphi_{2},\varphi_{1})}^{(4)}-{\cal K}_{\delta{\cal
G}_{2}(\varphi_{2},\varphi_{1})}^{(4,\mathrm{lim})}\,,$ (106)

where we have suppressed the argument
$(\bm{k}_{1},\bm{k}_{2},\bm{q},-\bm{q})$.

Table 5: Constraints on third- and fourth-order bias parameters from a joint analysis of the power spectrum and bispectrum with cutoff scales as described in Sec. IV.3. Sample | $b_{3}$ | $\gamma_{3}$ | $\gamma_{2}^{\times}$ | $\gamma_{211}$ | $\gamma_{22}$ | $\gamma_{31}$ |  $\gamma_{21}^{\times}$
---|---|---|---|---|---|---|---
MGS | $-1.5\pm 4$ | $-0.91\pm 2$ | $-0.76\pm 2$ | $0.0098\pm 1$ | $0.039\pm 0.05$ | $0.39\pm 0.8$ |  $-1.4\pm 2$
LOWZ | $-9.5\pm 9$ | $-0.54\pm 3$ | $-1.6\pm 3$ | $-0.48\pm 0.9$ | $-0.03\pm 0.2$ | $4.4\pm 2$ | $-1\pm 4$
CMASS | $-5.5\pm 10$ | $-1\pm 3$ | $-1.4\pm 3$ | $-0.52\pm 1$ | $0.086\pm 0.1$ | $1.6\pm 2$ | $-2.1\pm 5$
HALO1 | $-3.4\pm 10$ | $-2.6\pm 5$ | $-2\pm 5$ | $2.7\pm 5$ | $0.12\pm 0.1$ | $-0.19\pm 2$ | $-2\pm 5$
HALO2 | $-11\pm 20$ | $3\pm 4$ | $2\pm 6$ | $0.43\pm 2$ | $0.51\pm 0.5$ | $-0.53\pm 5$ | $6.4\pm 8$
HALO3 | $-11\pm 30$ | $-0.57\pm 8$ | $-1.9\pm 9$ | $-0.1\pm 3$ | $0.27\pm 0.5$ | $1.7\pm 7$ | $-0.79\pm 10$
HALO4 | $14\pm 20$ | $11\pm 6$ | $8.5\pm 8$ | $3.4\pm 7$ | $0.017\pm 1$ | $6.9\pm 10$ |  $21\pm 20$

##  APPENDIX B EVOLUTION OF CONSERVED TRACERS

To complete the expressions given in Sec. II.2, we here briefly review the
general solution of the coupled equations for matter and tracer densities,
$\delta$ and $\delta_{g}$, as well as matter velocity $\theta$, where we
assume the conservation of tracers and the absence of velocity bias
($\theta_{g}=\theta$). For further information, see Scoccimarro (2001); Chan
_et al._ (2012); Eggemeier  _et al._ (2019).

We start by introducing the triplet

$\Psi_{a}(\bm{k},\tau)\equiv\Big{[}\delta(\bm{k},\tau),\,\theta(\bm{k},\tau)/f{\cal
H},\,\delta_{g}(\bm{k},\tau)\Big{]}\,,$ (115)

where $\tau$ is the conformal time, ${\cal H}$ the conformal Hubble rate, and
$f$ the growth rate of structures. Using this triplet and changing time
variable to the logarithm of the linear growth factor,
$\eta\equiv\ln{D(\tau)}$, the continuity and Euler equations can be written in
the following compact form

$\begin{split}&\frac{\partial\Psi_{a}(\bm{k},\eta)}{\partial\eta}+\Omega_{ab}\,\Psi_{b}(\bm{k},\eta)\\\
&\hskip
1.99997pt=(2\pi)^{3}\int_{\bm{k}_{1},\bm{k}_{2}}\delta_{D}(\bm{k}-\bm{k}_{12})\,\gamma_{abc}(\bm{k}_{1},\bm{k}_{2})\\\
&\hskip
3.00003pt\times\,\Psi_{b}(\bm{k}_{1},\eta)\,\Psi_{c}(\bm{k}_{2},\eta)\,.\end{split}$
(116)

The couplings between densities and velocities at a linear level are described
by the matrix

$\Omega_{ab}\equiv\frac{1}{2}\left[\begin{array}[]{ccc}0&-2&0\\\ -3&1&0\\\
0&-2&0\end{array}\right]\,,$ (120)

whereas nonlinear interactions are encoded in
$\gamma_{abc}(\bm{k}_{1},\bm{k}_{2})$, whose only non-zero components are
given by

$\displaystyle\gamma_{121}=\gamma_{323}$
$\displaystyle=\alpha(\bm{k}_{1},\bm{k}_{2})\equiv\frac{\bm{k}_{12}\cdot\bm{k}_{1}}{k_{1}^{2}}\,,$
(121) $\displaystyle\gamma_{222}$
$\displaystyle=\beta(\bm{k}_{1},\bm{k}_{2})\equiv\frac{k_{12}^{2}\,\left(\bm{k}_{1}\cdot\bm{k}_{2}\right)}{2k_{1}^{2}\,k_{2}^{2}}\,,$
(122)

and $\gamma_{112}(\bm{k}_{1},\bm{k}_{2})=\gamma_{121}(\bm{k}_{2},\bm{k}_{1})$.
Letting $\phi_{a}(\bm{k})$ denote the initial conditions for the three fields
at $\eta=0$, it can be shown via Laplace transformation that Eq. (116) has an
integral solution,

$\begin{split}\Psi_{a}(\bm{k},\eta)=\,&g_{ab}(\eta)\,\phi_{b}(\bm{k})+(2\pi)^{3}\int_{0}^{\eta}\text{d}\eta^{\prime}\,g_{ab}(\eta-\eta^{\prime})\\\
&\times\int_{\bm{k}_{1},\bm{k}_{2}}\delta(\bm{k}-\bm{k}_{12})\,\gamma_{bcd}(\bm{k}_{1},\bm{k}_{2})\,\Psi_{c}(\bm{k}_{1},\eta^{\prime})\,\Psi_{d}(\bm{k}_{2},\eta^{\prime})\,,\end{split}$
(123)

with $g_{ab}(\eta)$ denoting the linear (time) propagator, which is given by a
combination of growing and decaying modes:

$\displaystyle g_{ab}(\eta)=\;$
$\displaystyle\frac{\mathrm{e}^{\eta}}{5}\left[\begin{array}[]{ccc}3&2&0\\\
3&2&0\\\
3&2&0\end{array}\right]-\frac{\mathrm{e}^{-3\eta/2}}{5}\left[\begin{array}[]{ccc}-2&2&0\\\
3&-3&0\\\ -2&2&0\\\ \end{array}\right]$ (130)
$\displaystyle+\left[\begin{array}[]{ccc}0&0&0\\\ 0&0&0\\\
-1&0&1\end{array}\right]\,.$ (134)

Plugging in the multipoint propagator expansion for $\Psi_{a}$ and $\phi_{a}$
into Eq. (123) one can derive a set of recursion relations at arbitrary order
of perturbation theory, which lead to the expressions in Eq. (13) to (15).

##  APPENDIX C CONSTRAINTS ON HIGHER-ORDER BIAS PARAMETERS

In Table 5 we have collected the remaining constraints on the third- and
fourth-order bias parameters for all of our samples. These values derive from
the same joint power spectrum and bispectrum analysis that was described in
Sec. IV.3 and thus complete Table 3.

## References

  * Levi  _et al._ (2013) M. Levi, C. Bebek, T. Beers, R. Blum, R. Cahn, D. Eisenstein, B. Flaugher, K. Honscheid, R. Kron, O. Lahav, P. McDonald, N. Roe, D. Schlegel, and representing the DESI collaboration, arXiv e-prints , arXiv:1308.0847 (2013), arXiv:1308.0847 [astro-ph.CO] .
  * Laureijs  _et al._ (2011) R. Laureijs, J. Amiaux, S. Arduini, J. L. Auguères, J. Brinchmann, R. Cole, M. Cropper, C. Dabin, L. Duvet, A. Ealet, B. Garilli, P. Gondoin, L. Guzzo, J. Hoar, H. Hoekstra, and Others, arXiv e-prints , arXiv:1110.3193 (2011), arXiv:1110.3193 [astro-ph.CO] .
  * Maartens  _et al._ (2015) R. Maartens, F. B. Abdalla, M. Jarvis, and M. G. Santos, arXiv e-prints , arXiv:1501.04076 (2015), arXiv:1501.04076 [astro-ph.CO] .
  * Sefusatti  _et al._ (2006) E. Sefusatti, M. Crocce, S. Pueblas, and R. Scoccimarro, Phys. Rev. D 74, 023522 (2006), arXiv:astro-ph/0604505 .
  * Song  _et al._ (2015) Y.-S. Song, A. Taruya, and A. Oka, J. Cosmology Astropart. Phys 2015, 007 (2015), arXiv:1502.03099 [astro-ph.CO] .
  * Gagrani and Samushia (2017) P. Gagrani and L. Samushia, MNRAS 467, 928 (2017), arXiv:1610.03488 [astro-ph.CO] .
  * Byun  _et al._ (2017) J. Byun, A. Eggemeier, D. Regan, D. Seery, and R. E. Smith, MNRAS 471, 1581 (2017), arXiv:1705.04392 [astro-ph.CO] .
  * Yankelevich and Porciani (2019) V. Yankelevich and C. Porciani, MNRAS 483, 2078 (2019), arXiv:1807.07076 [astro-ph.CO] .
  * Agarwal  _et al._ (2020) N. Agarwal, V. Desjacques, D. Jeong, and F. Schmidt, arXiv e-prints , arXiv:2007.04340 (2020), arXiv:2007.04340 [astro-ph.CO] .
  * Chudaykin and Ivanov (2019) A. Chudaykin and M. M. Ivanov, J. Cosmology Astropart. Phys 2019, 034 (2019), arXiv:1907.06666 [astro-ph.CO] .
  * Hahn  _et al._ (2020) C. Hahn, F. Villaescusa-Navarro, E. Castorina, and R. Scoccimarro, J. Cosmology Astropart. Phys 2020, 040 (2020), arXiv:1909.11107 [astro-ph.CO] .
  * Kamalinejad and Slepian (2020) F. Kamalinejad and Z. Slepian, arXiv e-prints , arXiv:2011.00899 (2020), arXiv:2011.00899 [astro-ph.CO] .
  * Hahn and Villaescusa-Navarro (2020) C. Hahn and F. Villaescusa-Navarro, arXiv e-prints , arXiv:2012.02200 (2020), arXiv:2012.02200 [astro-ph.CO] .
  * Scoccimarro  _et al._ (2012) R. Scoccimarro, L. Hui, M. Manera, and K. C. Chan, Phys. Rev. D 85, 083002 (2012), arXiv:1108.5512 [astro-ph.CO] .
  * Sefusatti  _et al._ (2012) E. Sefusatti, M. Crocce, and V. Desjacques, MNRAS 425, 2903 (2012), arXiv:1111.6966 [astro-ph.CO] .
  * Tellarini  _et al._ (2016) M. Tellarini, A. J. Ross, G. Tasinato, and D. Wands, J. Cosmology Astropart. Phys 2016, 014 (2016), arXiv:1603.06814 [astro-ph.CO] .
  * Karagiannis  _et al._ (2018) D. Karagiannis, A. Lazanu, M. Liguori, A. Raccanelli, N. Bartolo, and L. Verde, MNRAS 478, 1341 (2018), arXiv:1801.09280 [astro-ph.CO] .
  * Karagiannis  _et al._ (2020) D. Karagiannis, J. Fonseca, R. Maartens, and S. Camera, arXiv e-prints , arXiv:2010.07034 (2020), arXiv:2010.07034 [astro-ph.CO] .
  * Moradinezhad Dizgah  _et al._ (2020) A. Moradinezhad Dizgah, M. Biagetti, E. Sefusatti, V. Desjacques, and J. Noreña, arXiv e-prints , arXiv:2010.14523 (2020), arXiv:2010.14523 [astro-ph.CO] .
  * Gil-Marín  _et al._ (2015a) H. Gil-Marín, J. Noreña, L. Verde, W. J. Percival, C. Wagner, M. Manera, and D. P. Schneider, MNRAS 451, 539 (2015a), arXiv:1407.5668 [astro-ph.CO] .
  * Gil-Marín  _et al._ (2017a) H. Gil-Marín, W. J. Percival, L. Verde, J. R. Brownstein, C.-H. Chuang, F.-S. Kitaura, S. A. Rodríguez-Torres, and M. D. Olmstead, MNRAS 465, 1757 (2017a), arXiv:1606.00439 [astro-ph.CO] .
  * Slepian  _et al._ (2017a) Z. Slepian, D. J. Eisenstein, F. Beutler, C.-H. Chuang, A. J. Cuesta, J. Ge, H. Gil-Marín, S. Ho, F.-S. Kitaura, C. K. McBride, R. C. Nichol, W. J. Percival, S. Rodríguez-Torres, A. J. Ross, R. Scoccimarro, H.-J. Seo, J. Tinker, R. Tojeiro, and M. Vargas-Magaña, MNRAS 468, 1070 (2017a).
  * Slepian  _et al._ (2017b) Z. Slepian, D. J. Eisenstein, J. R. Brownstein, C.-H. Chuang, H. Gil-Marín, S. Ho, F.-S. Kitaura, W. J. Percival, A. J. Ross, G. Rossi, H.-J. Seo, A. Slosar, and M. Vargas-Magaña, MNRAS 469, 1738 (2017b), arXiv:1607.06097 [astro-ph.CO] .
  * Pearson and Samushia (2018) D. W. Pearson and L. Samushia, MNRAS 478, 4500 (2018), arXiv:1712.04970 [astro-ph.CO] .
  * Smith  _et al._ (2008) R. E. Smith, R. K. Sheth, and R. Scoccimarro, Phys. Rev. D 78, 023523 (2008), arXiv:0712.0017 [astro-ph] .
  * Rampf and Wong (2012) C. Rampf and Y. Y. Y. Wong, J. Cosmology Astropart. Phys 2012, 018 (2012), arXiv:1203.4261 [astro-ph.CO] .
  * Bernardeau  _et al._ (2012) F. Bernardeau, M. Crocce, and R. Scoccimarro, Phys. Rev. D 85, 123519 (2012), arXiv:1112.3895 [astro-ph.CO] .
  * Assassi  _et al._ (2014) V. Assassi, D. Baumann, D. Green, and M. Zaldarriaga, J. Cosmology Astropart. Phys 8, 056 (2014), arXiv:1402.5916 .
  * Baldauf  _et al._ (2015) T. Baldauf, L. Mercolli, M. Mirbabayi, and E. Pajer, J. Cosmology Astropart. Phys 2015, 007 (2015), arXiv:1406.4135 [astro-ph.CO] .
  * Angulo  _et al._ (2015a) R. E. Angulo, S. Foreman, M. Schmittfull, and L. Senatore, Journal of Cosmology and Astro-Particle Physics 2015, 039 (2015a), arXiv:1406.4143 [astro-ph.CO] .
  * Lazanu  _et al._ (2016) A. Lazanu, T. Giannantonio, M. Schmittfull, and E. P. S. Shellard, Phys. Rev. D 93, 083517 (2016), arXiv:1510.04075 [astro-ph.CO] .
  * Hashimoto  _et al._ (2017) I. Hashimoto, Y. Rasera, and A. Taruya, Phys. Rev. D 96, 043526 (2017), arXiv:1705.02574 [astro-ph.CO] .
  * Desjacques  _et al._ (2018a) V. Desjacques, D. Jeong, and F. Schmidt, J. Cosmology Astropart. Phys 2018, 035 (2018a), arXiv:1806.04015 [astro-ph.CO] .
  * Eggemeier  _et al._ (2019) A. Eggemeier, R. Scoccimarro, and R. E. Smith, Phys. Rev. D 99, 123514 (2019), arXiv:1812.03208 [astro-ph.CO] .
  * Desjacques  _et al._ (2018b) V. Desjacques, D. Jeong, and F. Schmidt, Phys. Rep. 733, 1 (2018b).
  * Kaiser (1984) N. Kaiser, ApJ 284, L9 (1984).
  * Coles (1993) P. Coles, MNRAS 262, 1065 (1993).
  * Fry and Gaztanaga (1993) J. N. Fry and E. Gaztanaga, ApJ 413, 447 (1993), arXiv:astro-ph/9302009 .
  * Catelan  _et al._ (1998) P. Catelan, F. Lucchin, S. Matarrese, and C. Porciani, MNRAS 297, 692 (1998), arXiv:astro-ph/9708067 .
  * Catelan  _et al._ (2000) P. Catelan, C. Porciani, and M. Kamionkowski, MNRAS 318, L39 (2000), astro-ph/0005544 .
  * McDonald and Roy (2009) P. McDonald and A. Roy, J. Cosmology Astropart. Phys 8, 020 (2009), arXiv:0902.0991 [astro-ph.CO] .
  * Matsubara (2011) T. Matsubara, Phys. Rev. D 83, 083518 (2011), arXiv:1102.4619 [astro-ph.CO] .
  * Guo and Jing (2009) H. Guo and Y. P. Jing, ApJ 702, 425 (2009), arXiv:0907.0282 [astro-ph.CO] .
  * Pollack  _et al._ (2012) J. E. Pollack, R. E. Smith, and C. Porciani, MNRAS 420, 3469 (2012), arXiv:1109.3458 .
  * Chan  _et al._ (2012) K. C. Chan, R. Scoccimarro, and R. K. Sheth, Phys. Rev. D 85, 083509 (2012), arXiv:1201.3614 [astro-ph.CO] .
  * Baldauf  _et al._ (2012) T. Baldauf, U. Seljak, V. Desjacques, and P. McDonald, Phys. Rev. D 86, 083540 (2012), arXiv:1201.4827 [astro-ph.CO] .
  * Sheth  _et al._ (2013) R. K. Sheth, K. C. Chan, and R. Scoccimarro, Phys. Rev. D 87, 083002 (2013), arXiv:1207.7117 [astro-ph.CO] .
  * Pollack  _et al._ (2014) J. E. Pollack, R. E. Smith, and C. Porciani, MNRAS 440, 555 (2014), arXiv:1309.0504 [astro-ph.CO] .
  * Saito  _et al._ (2014) S. Saito, T. Baldauf, Z. Vlah, U. Seljak, T. Okumura, and P. McDonald, Phys. Rev. D 90, 123522 (2014), arXiv:1405.1447 [astro-ph.CO] .
  * Biagetti  _et al._ (2014) M. Biagetti, V. Desjacques, A. Kehagias, and A. Riotto, Phys. Rev. D 90, 103529 (2014), arXiv:1408.0293 [astro-ph.CO] .
  * Angulo  _et al._ (2015b) R. Angulo, M. Fasiello, L. Senatore, and Z. Vlah, Journal of Cosmology and Astro-Particle Physics 2015, 029 (2015b), arXiv:1503.08826 [astro-ph.CO] .
  * Bel  _et al._ (2015) J. Bel, K. Hoffmann, and E. Gaztañaga, MNRAS 453, 259 (2015), arXiv:1504.02074 [astro-ph.CO] .
  * Lazeyras and Schmidt (2018) T. Lazeyras and F. Schmidt, Journal of Cosmology and Astro-Particle Physics 2018, 008 (2018), arXiv:1712.07531 [astro-ph.CO] .
  * Abidi and Baldauf (2018) M. M. Abidi and T. Baldauf, J. Cosmology Astropart. Phys 7, 029 (2018), arXiv:1802.07622 .
  * Oddo  _et al._ (2020) A. Oddo, E. Sefusatti, C. Porciani, P. Monaco, and A. G. Sánchez, J. Cosmology Astropart. Phys 2020, 056 (2020), arXiv:1908.01774 [astro-ph.CO] .
  * Chiang  _et al._ (2015) C.-T. Chiang, C. Wagner, A. G. Sánchez, F. Schmidt, and E. Komatsu, J. Cosmology Astropart. Phys 2015, 028 (2015), arXiv:1504.03322 [astro-ph.CO] .
  * Fry (1996) J. N. Fry, ApJ 461, L65+ (1996).
  * Sánchez  _et al._ (2017) A. G. Sánchez, R. Scoccimarro, M. Crocce, J. N. Grieb, S. Salazar-Albornoz, C. Dalla Vecchia, M. Lippich, F. Beutler, J. R. Brownstein, C.-H. Chuang, D. J. Eisenstein, F.-S. Kitaura, M. D. Olmstead, W. J. Percival, F. Prada, S. Rodríguez-Torres, A. J. Ross, L. Samushia, H.-J. Seo, J. Tinker, R. Tojeiro, M. Vargas-Magaña, Y. Wang, and G.-B. Zhao, MNRAS 464, 1640 (2017), arXiv:1607.03147 .
  * Gil-Marín  _et al._ (2017b) H. Gil-Marín, W. J. Percival, L. Verde, J. R. Brownstein, C.-H. Chuang, F.-S. Kitaura, S. A. Rodríguez-Torres, and M. D. Olmstead, MNRAS 465, 1757 (2017b), arXiv:1606.00439 [astro-ph.CO] .
  * Beutler  _et al._ (2017) F. Beutler, H.-J. Seo, S. Saito, C.-H. Chuang, A. J. Cuesta, D. J. Eisenstein, H. Gil-Marín, J. N. Grieb, N. Hand, F.-S. Kitaura, C. Modi, R. C. Nichol, M. D. Olmstead, W. J. Percival, F. Prada, A. G. Sánchez, S. Rodriguez-Torres, A. J. Ross, N. P. Ross, D. P. Schneider, J. Tinker, R. Tojeiro, and M. Vargas-Magaña, MNRAS 466, 2242 (2017), arXiv:1607.03150 [astro-ph.CO] .
  * Grieb  _et al._ (2017) J. N. Grieb, A. G. Sánchez, S. Salazar-Albornoz, R. Scoccimarro, M. Crocce, C. Dalla Vecchia, F. Montesano, H. Gil-Marín, A. J. Ross, F. Beutler, S. Rodríguez-Torres, C.-H. Chuang, F. Prada, F.-S. Kitaura, A. J. Cuesta, D. J. Eisenstein, W. J. Percival, M. Vargas-Magaña, J. L. Tinker, R. Tojeiro, J. R. Brownstein, C. Maraston, R. C. Nichol, M. D. Olmstead, L. Samushia, H.-J. Seo, A. Streblyanska, and G.-b. Zhao, MNRAS 467, 2085 (2017), arXiv:1607.03143 .
  * Eggemeier  _et al._ (2020) A. Eggemeier, R. Scoccimarro, M. Crocce, A. Pezzotta, and A. G. Sánchez, arXiv e-prints , arXiv:2006.09729 (2020), arXiv:2006.09729 [astro-ph.CO] .
  * (63) A. Pezzotta et al., in prep. .
  * Senatore (2015) L. Senatore, Journal of Cosmology and Astro-Particle Physics 2015, 007 (2015), arXiv:1406.7843 [astro-ph.CO] .
  * Mirbabayi  _et al._ (2015) M. Mirbabayi, F. Schmidt, and M. Zaldarriaga, J. Cosmology Astropart. Phys 7, 030 (2015), arXiv:1412.5169 .
  * Kofman and Pogosyan (1995) L. Kofman and D. Pogosyan, ApJ 442, 30 (1995), astro-ph/9403029 .
  * Bertschinger (1995) E. Bertschinger, “Cosmological dynamics,” NASA STI/Recon Technical Report N (1995), arXiv:astro-ph/9503125 [astro-ph] .
  * McDonald (2006) P. McDonald, Phys. Rev. D 74, 103512 (2006), astro-ph/0609413 .
  * Crocce and Scoccimarro (2006a) M. Crocce and R. Scoccimarro, Phys. Rev. D 73, 063519 (2006a), arXiv:astro-ph/0509418 .
  * Crocce and Scoccimarro (2006b) M. Crocce and R. Scoccimarro, Phys. Rev. D 73, 063520 (2006b), arXiv:astro-ph/0509419 .
  * Mo  _et al._ (1997) H. J. Mo, Y. P. Jing, and S. D. M. White, MNRAS 284, 189 (1997).
  * Paranjape and Sheth (2012) A. Paranjape and R. K. Sheth, MNRAS 426, 2789 (2012), arXiv:1206.3506 [astro-ph.CO] .
  * Wagner  _et al._ (2015) C. Wagner, F. Schmidt, C. T. Chiang, and E. Komatsu, MNRAS 448, L11 (2015), arXiv:1409.6294 [astro-ph.CO] .
  * Lazeyras  _et al._ (2016) T. Lazeyras, C. Wagner, T. Baldauf, and F. Schmidt, Journal of Cosmology and Astro-Particle Physics 2016, 018 (2016), arXiv:1511.01096 [astro-ph.CO] .
  * Matsubara (1995) T. Matsubara, ApJS 101, 1 (1995), astro-ph/9501056 .
  * Bernardeau  _et al._ (2008) F. Bernardeau, M. Crocce, and R. Scoccimarro, Phys. Rev. D 78, 103521 (2008), arXiv:0806.2334 .
  * Scoccimarro (1997) R. Scoccimarro, ApJ 487, 1 (1997).
  * Desjacques (2008) V. Desjacques, Phys. Rev. D 78, 103503 (2008), arXiv:0806.0007 [astro-ph] .
  * Desjacques  _et al._ (2010) V. Desjacques, M. Crocce, R. Scoccimarro, and R. K. Sheth, Phys. Rev. D 82, 103529 (2010), arXiv:1009.3449 [astro-ph.CO] .
  * Bernardeau  _et al._ (2002) F. Bernardeau, S. Colombi, E. Gaztañaga, and R. Scoccimarro, Phys. Rep. 367, 1 (2002), astro-ph/0112551 .
  * Lazeyras and Schmidt (2019) T. Lazeyras and F. Schmidt, J. Cosmology Astropart. Phys 2019, 041 (2019), arXiv:1904.11294 [astro-ph.CO] .
  * Pueblas and Scoccimarro (2009) S. Pueblas and R. Scoccimarro, Phys. Rev. D 80, 043504 (2009), arXiv:0809.4606 .
  * Dekel and Lahav (1999) A. Dekel and O. Lahav, ApJ 520, 24 (1999), arXiv:astro-ph/9806193 [astro-ph] .
  * Taruya and Soda (1999) A. Taruya and J. Soda, ApJ 522, 46 (1999), arXiv:astro-ph/9809204 [astro-ph] .
  * Matsubara (1999) T. Matsubara, ApJ 525, 543 (1999), arXiv:astro-ph/9906029 [astro-ph] .
  * Mo and White (1996) H. J. Mo and S. D. M. White, MNRAS 282, 347 (1996), arXiv:astro-ph/9512127 .
  * Sheth and Lemson (1999) R. Sheth and G. Lemson, MNRAS 304, 767 (1999).
  * Smith  _et al._ (2007) R. E. Smith, R. Scoccimarro, and R. K. Sheth, Phys. Rev. D 75, 063512 (2007), arXiv:astro-ph/0609547 .
  * Baldauf  _et al._ (2013) T. Baldauf, U. Seljak, R. E. Smith, N. Hamaus, and V. Desjacques, Phys. Rev. D 88, 083507 (2013), arXiv:1305.2917 [astro-ph.CO] .
  * Schmidt (2016) F. Schmidt, Phys. Rev. D 93, 063512 (2016), arXiv:1511.02231 [astro-ph.CO] .
  * Gil-Marín  _et al._ (2015b) H. Gil-Marín, J. Noreña, L. Verde, W. J. Percival, C. Wagner, M. Manera, and D. P. Schneider, MNRAS 451, 539 (2015b), arXiv:1407.5668 [astro-ph.CO] .
  * Scoccimarro (2015) R. Scoccimarro, Phys. Rev. D 92, 083532 (2015), arXiv:1506.02729 [astro-ph.CO] .
  * Sefusatti  _et al._ (2016) E. Sefusatti, M. Crocce, R. Scoccimarro, and H. M. P. Couchman, MNRAS 460, 3624 (2016), arXiv:1512.07295 [astro-ph.CO] .
  * Peebles (1980) P. Peebles,  _The large-scale structure of the universe_ (Princeton University Press, 1980).
  * Scoccimarro  _et al._ (1998) R. Scoccimarro, S. Colombi, J. N. Fry, J. A. Frieman, E. Hivon, and A. Melott, ApJ 496, 586 (1998), arXiv:astro-ph/9704075 .
  * Scoccimarro  _et al._ (1999) R. Scoccimarro, H. M. P. Couchman, and J. A. Frieman, ApJ 517, 531 (1999), arXiv:astro-ph/9808305 [astro-ph] .
  * Feldman  _et al._ (1994) H. A. Feldman, N. Kaiser, and J. A. Peacock, ApJ 426, 23 (1994), arXiv:astro-ph/9304022 .
  * Tegmark (1997) M. Tegmark, Phys. Rev. Lett. 79, 3806 (1997), arXiv:astro-ph/9706198 [astro-ph] .
  * Sefusatti  _et al._ (2010) E. Sefusatti, M. Crocce, and V. Desjacques, MNRAS 406, 1014 (2010), arXiv:1003.0007 [astro-ph.CO] .
  * Gelman and Rubin (1992) A. Gelman and D. B. Rubin, Statistical Science 7, 457 (1992).
  * Lewis (2019) A. Lewis,  (2019), arXiv:1910.13970 [astro-ph.IM] .
  * Scoccimarro and Frieman (1999) R. Scoccimarro and J. A. Frieman, ApJ 520, 35 (1999), arXiv:astro-ph/9811184 .
  * Sefusatti and Scoccimarro (2005) E. Sefusatti and R. Scoccimarro, Phys. Rev. D 71, 063001 (2005), arXiv:astro-ph/0412626 [astro-ph] .
  * Bernardeau  _et al._ (2014) F. Bernardeau, A. Taruya, and T. Nishimichi, Phys. Rev. D 89, 023502 (2014).
  * Nishimichi  _et al._ (2016) T. Nishimichi, F. Bernardeau, and A. Taruya, Physics Letters B 762, 247 (2016), arXiv:1411.2970 [astro-ph.CO] .
  * Nishimichi  _et al._ (2017) T. Nishimichi, F. Bernardeau, and A. Taruya, Phys. Rev. D 96, 123515 (2017), arXiv:1708.08946 [astro-ph.CO] .
  * Ivanov  _et al._ (2019) M. M. Ivanov, M. Simonović, and M. Zaldarriaga, arXiv e-prints , arXiv:1909.05277 (2019), arXiv:1909.05277 [astro-ph.CO] .
  * D’Amico  _et al._ (2019) G. D’Amico, J. Gleyzes, N. Kokron, D. Markovic, L. Senatore, P. Zhang, F. Beutler, and H. Gil-Marín, arXiv e-prints , arXiv:1909.05271 (2019), arXiv:1909.05271 [astro-ph.CO] .
  * Tröster  _et al._ (2020) T. Tröster, A. G. Sánchez, M. Asgari, C. Blake, M. Crocce, C. Heymans, H. Hildebrandt, B. Joachimi, S. Joudaki, A. Kannawadi, C.-A. Lin, and A. Wright, A&A 633, L10 (2020), arXiv:1909.11006 [astro-ph.CO] .
  * Feldman  _et al._ (2001) H. A. Feldman, J. A. Frieman, J. N. Fry, and R. Scoccimarro, Physical Review Letters 86, 1434 (2001), arXiv:astro-ph/0010205 .
  * Sugiyama  _et al._ (2021) N. S. Sugiyama, S. Saito, F. Beutler, and H.-J. Seo, MNRAS 501, 2862 (2021), arXiv:2010.06179 [astro-ph.CO] .
  * Hunter (2007) J. D. Hunter, Computing in Science & Engineering 9, 90 (2007).
  * Scoccimarro (2001) R. Scoccimarro, in  _The Onset of Nonlinearity in Cosmology_, Vol. 927, edited by J. N. Fry, J. R. Buchler, and H. Kandrup (2001) pp. 13–23, arXiv:astro-ph/0008277 [astro-ph] .


# Verifying the Australian MWA EoR pipeline I: sky model and correlated
measurement density

J. L. B. Line  International Centre for Radio Astronomy Research, Curtin
University, Perth, WA 6102, Australia  [ C. M. Trott  International Centre for
Radio Astronomy Research, Curtin University, Perth, WA 6102, Australia  J. H.
Cook  International Centre for Radio Astronomy Research, Curtin University,
Perth, WA 6102, Australia  B. Greig  School of Physics, University of
Melbourne, Parkville, VIC 3010, Australia  N. Barry  International Centre for
Radio Astronomy Research, Curtin University, Perth, WA 6102, Australia  C. H.
Jordan  International Centre for Radio Astronomy Research, Curtin University,
Perth, WA 6102, Australia

###### Abstract

We present the first of two papers dedicated to verifying the Australian Epoch
of Reionisation pipeline (AusEoRPipe) through simulation. The AusEoRPipe aims
to disentangle radiation emitted by gas surrounding the very first stars from
contaminating foreground astrophysical sources, and has been in development
for close to a decade. In this paper, we build an accurate sky model that can
be used by the WODEN simulation software to create visibilities containing a
predictable signal. We verify that the power spectrum estimator CHIPS can
recover this signal in the absence of foregrounds. We also investigate how
measurements in Fourier-space are correlated, and how their gridded density
affects the power spectrum. We measure and fit for this effect using Gaussian-
noise simulations of the MWA phase I layout. We find a gridding density
correction factor of 2.651 appropriate for integrations equal to or greater
than 30 minutes of data, which contain observations with multiple primary beam
pointings and LSTs. Paper II of this series will use the results of this paper
to test the AusEoRPipe in the presence of foregrounds and instrumental
effects.

###### keywords:

Astronomy data analysis, Reionisation, GPU computing

ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO-3D)
J. L. B<EMAIL_ADDRESS>ARC Centre of Excellence for All Sky
Astrophysics in 3 Dimensions (ASTRO-3D) ARC Centre of Excellence for All Sky
Astrophysics in 3 Dimensions (ASTRO-3D) ARC Centre of Excellence for All Sky
Astrophysics in 3 Dimensions (ASTRO-3D) ARC Centre of Excellence for All Sky
Astrophysics in 3 Dimensions (ASTRO-3D) ARC Centre of Excellence for All Sky
Astrophysics in 3 Dimensions (ASTRO-3D)

##  1 Introduction

Constraining when and how the first stars were formed has been a goal of
astronomers for decades. By measuring cosmically-redshifted radiation from the
hydrogen gas surrounding those first stars, one can conceivably map the
effects of the ionising radiation coming from them, to infer their properties
as a function of redshift. Unfortunately, the detection of radiation from the
Epoch of Reionisation (EoR) (redshifted to the 50-250 MHz frequency range), is
impeded by a myriad of astrophysical sources, including Active Galactic
Nuclei, Radio Galaxies, Supernova Remnants, and the diffuse synchrotron
radiation emitted from the Milky Way. These foregrounds drown out the signal
Furlanetto2006b (16), and must either be subtracted from the data or avoided.
One approach that somewhat naturally separates foregrounds from the signal is
via measuring the power spectrum (PS). The PS can be used to garner spatial
information of the effects of reionisation by the first stars, and infer their
properties. Furthermore, Fourier transforming over frequency to obtain a PS
results in the foregrounds and signal manifesting in different areas of the PS
given their differing spectral behaviours. Recent upper limits on the EoR have
come from: Hydrogen Epoch of Reionisation Array (DeBoer2016, 13, HERA,); Low-
Frequency Array (VanHaarlem2013, 19, LOFAR,); The Murchison Widefield Array
(Tingay2013, 40, 44, MWA,); New Extension in Nançay Upgrading LOFAR
(Zarka2012, 47, NenuFar,).

There are a number of challenges which must be overcome to make the
statistical detection of the signal. As well as nullifying the aforementioned
foregrounds (Cook2022, 12, 1, e.g.), one must deal with ionospheric refraction
(Elder2021, 14, 10, e.g.), terrestrial interference (Offringa2019, 36, 46,
e.g.), and a myriad of instrumental calibration effects (Kern2020, 23, 11, 31,
24, e.g.), all of which cause limiting systematics. These challenges have
forced a decade of not only hardware but also software development, to
calibrate and treat the data with extreme precision to uncover the underlying
signal. Recent best-effort PS limits include Mertens2020 (28, 41, 38, 20, 34).

Given that measuring the PS is as much a software as it is an observational
experiment, great care must be taken to understand what biases and systematics
the software pipeline imparts upon the data. Given the limiting systematics
listed above, a natural route is via simulation, where one can control which
systematics are injected into a data set, with a ground “truth” to be
recovered. In short, if you inject a known signal into a data set, process
through your pipeline and recover what you injected, you gain confidence in
your pipeline. This approach has been taken by a number of authors, including
Barry2019_pipeline (6) using MWA data and HERA_verify (2) with HERA data. To
date, the Australian EoR pipeline (AusEoRPipe), which is designed to process
MWA data, has not been tested from end to end. This is the goal of this series
of papers. This paper aims to test the gridding and PS estimation aspects of
the pipeline; the following paper in this series will focus on calibration.

The native output of an interferometer is a visibility, a sampling of the
Fourier transform of the visible sky. There are two main ways to simulate a
signal in visibility space. One can either make an analytic model that
produces similar statistics in a PS, to directly generate in visibility space,
or can create an image based model and derive visibilities from that. The
former is far less computationally expensive, but the latter locks in
projection and geometrical effects inherent to true interferometric
observations. Having an image-based model locked to celestial coordinates
allows for simulated observations at various hour angles, to probe the effects
of changing instrumental primary beam patterns and visibility sampling.

Simulating an image-based signal for the MWA is non-trivial however, given the
sky coverage of the instrument. The full-width half-maximum of the main lobe
of the primary beam of the MWA can be greater than $20\times 20$ square
degrees, and $50\times 50$ square degrees down to the 1% beam level.
Simulating a realistic EoR volume capable of covering this sky-area is
computationally demanding. Fortunately, Greig2022 (18) generated a simulated
EoR volume large enough, which we use in this work and detail in Section 3.
Given the resolution of the instrument is of order an arcminute, the number of
pixels required to represent lightcones from this EoR volume is
$\mathcal{O}(10^{7})$. Down-sampling this model can affect the resultant
statistics of the signal. For this reason, we use the WODEN simulation package
Line2022 (26), capable of ingesting sky models with millions of components,
and producing realistic simulated MWA visibilities. As WODEN is GPU-
accelerated, these simulations can be run in a reasonable time frame (see
Section 4.1 for details). Importantly, WODEN is also capable of simulating
both the signal and the foregrounds at the expected power levels, which can be
many orders of magnitude apart. This allows us to test signal recovery in the
presence of foregrounds without boosting the amplitude of the signal.

This paper is organised as follows. In Section 2 we briefly overview the
AusEoRPipe and outline what simulated data are necessary to test them. In
Section 3 we detail the construction of a sky model compatible with WODEN. In
Section 4 we verify the steps taken to project the model into celestial
coordinates, by applying them to a sky model of Gaussian noise. We also
investigate the effects of gridding visibilities when generating power
spectra, and fit for the effect this has on estimating the resultant power. In
Section 5 we test the model and how well the AusEoRPipe can recover the
expected PS. We discuss our results in Section 6, and conclude in Section 7.

##  2 Pipeline overview

The AusEoRPipe is designed to take raw data from the MWA correlator, calibrate
for instrumental and ionospheric effects, subtract foregrounds, and then
estimate a PS. Here we list the different software packages that make up the
AusEoRPipe, and suggest simulations needed to test them. This Section is
intended to give the reader a broad overview of the AusEoRPipe and the overall
goals of this series of papers, to give this paper context. At the end of this
Section, we detail which parts of the work are covered in this paper.

Birli 111https://github.com/MWATelescope/Birli is designed to take raw MWA
data and preprocess it, including geometric correction, averaging, and radio
frequency interference (RFI) flagging (via AOFlagger Offringa2010 (35) or
SINSS Wilensky2019 (45)). Birli has been written to replace Cotter
Offringa2015 (37), to have one unifying preprocessing package that works with
both legacy and new MWA data (Morrison2023, 33, see). Given Birli has been
tested to reproduce results out of Cotter, which has been tested and used
within a number of MWA pipelines, we choose not to test Birli via simulation.
To do so however, one would need to simulate raw MWA data, which for the
original correlator would mean producing visibilities with no phase tracking,
adding a frequency dependent bandpass from the MWA polyphase band filter, and
add in realistic RFI.

hyperdrive 222https://github.com/MWATelescope/mwa_hyperdrive (Jordan et al.
submitted), takes either raw or preprocessed data, calibrates it, and then
subtracts foregrounds. Built to replace the RTS Mitchell2008 (32), hyperdrive
also creates model visibilities to calibrate against via an image-based sky
model. It initially performs a direction independent (DI) calibration step,
that derives a gain for each MWA receiving element (tile). All frequencies are
calibrated independently. After DI calibration, foregrounds can be subtracted.
The Fourier transform along frequency to create the PS is expected to separate
the signal from the foregrounds, as foreground emission is supposed to be
spectrally smooth, whereas the signal has spectral structure. Incomplete
sampling of the $u,v$ plane however causes mode mixing between the two, as
well as other frequency-related instrumental effects. A detection therefore
relies on the foregrounds being subtracted without injecting false spectral
structure, as to preserve the underlying signal. To test this functionality, a
simulation needs to contain both a signal and foregrounds, as well as
frequency dependent instrumental gain errors. One can also add frequency-
dependent effects like cable reflections (Ewall-Wice2016, 15, e.g.) and the
bandpass to investigate how calibrating each frequency channel independently
performs. These effects are direction-independent and so can be added to
visibilities post simulation (functionality which is in development in WODEN).
Further simulations including RFI and ionospheric refraction would test how
robust the calibration is to environmental influences. These effects are
direction dependent and so must be adding during the calculations of the
visibilities, which WODEN is currently incapable of.

CHIPS Trott2016 (42) takes the calibrated and subtracted visibilities and
grids them into a spectral cube using an optimised Blackman-harris kernel. It
then Fourier-transforms along frequency, and then either cylindrically
averages to create a two dimensional PS (2D PS) or spherically averages to
create a one dimensional PS (1D PS). Integrations of tens to hundreds of hours
of data into a 1D PS have set the limits that have been released by the
AusEoRPipe to date. The minimal test for CHIPS is to simulate a set of
visibilities with a known signal, and see if CHIPS can recover it. Beyond
that, taking in calibrated visibilities out of hyperdrive, derived from
simulations with realistic instrumental effects and foregrounds, allows us to
probe the effects of calibration on the resultant PS.

WODEN 333https://github.com/JLBLine/WODEN, as previously mentioned, simulates
MWA visibilities from a sky based model. It can use the same primary beam
model as hyperdrive, the FEE beam Sokolowski2017 (39) (via the
mwa_hyperbeam444https://github.com/MWATelescope/mwa_hyperbeam package), and
creates visibilities that can be fed directly into either hyperdrive or CHIPS.
WODEN calculates the interferometric measurement equation for every component
in the sky model. Each component can either be a point source (a dirac-delta
function upon the sky), an elliptical Gaussian, or a Shapelet model (Line2020,
25, see). The measurement equation encodes baseline and sky projection effects
inherently, making the simulator accurate across the large field of view of
the MWA. The caveat being that the sky model must be broken into the small
components listed here, meaning representations of large-scale structure such
as the galactic diffuse emission need millions of components. Aside from the
primary beam model, the WODEN code base is entirely independent to hyperdrive.
Although they use the same methodology to generate visibilities, this
redundancy in code is by design; any internal bug that may cancel out when
generating visibilities and calibrating them all with the same code base is
avoided.

Along with the aforementioned packages, extensive quality metrics are
generated at various stages (Nunhokee et al. submitted). These metrics are
used to cull any observations that are deemed unusable due to insurmountable
instrumental effects, ionospheric conditions, RFI events, contamination by
bright sources in the primary beam sidelobes, and other effects. The
simulations detailed above could be used to check whether these metrics are
able to catch poor data, and check whether the limits set for each metric
still allow for a detection of the signal.

In this paper, we ignore calibration entirely, and focus on generating a
realistic sky model to simulate through WODEN and directly input into CHIPS.
As visibilities are additive, once we are able to generate accurate
simulations, we can add a variety of foregrounds and instrumental effects to
test the calibration and subtraction step, without the need to rerun the more
expensive simulations. We leave testing hyperdrive to the second paper in this
series.

###  2.1 MWA EoR observing

The MWA is a low-frequency radio interferometer, with receiving elements
consisting of 4$\times$4 grids of dual linear-polarisation bow-tie antennas.
These ‘tiles’ are electronically steered through beamforming. To estimate the
PS, the MWA EoR collaboration have identified a number of fields with lower
sky temperatures, which have now been observed for hundreds of hours
(Trott2020, 41, e.g.), in an attempt to average over thermal noise. Given the
electronic beamforming, this has resulted in a drift-and-shift observational
campaign, where the target field drifts through a number of pointings, with
scheduling keeping the field centre as close to the primary beam centre as
possible (Jacobs2016, 21, see).

In this paper we focus on the EoR0 field, centred at RA, DEC
$=0^{h},-27^{\circ}$, observed between frequencies of 167-198 MHz (known as
the high-band). These parameters have consistently yielded the best MWA limits
(Barry2019, 7, 41, 10, e.g.). We also focus on five pointings, labelled -2
through +2 (Beardsley2016, 8, c.f.), where -2 means two pointings before the
zenith 0 pointing, and +2 meaning two pointings after 0 pointing. There is an
approximate 15∘ shift between pointings, all along the meridian. Examples are
shown in Figure 1. MWA EoR data are typically taken in two minute chunks known
as an observation.

##  3 sky-model

We use the lightcone detailed in Greig2022 (18), which is derived from a
simulated EoR cube with sides of length $7.5\,$Gpc. This cube was generated
using a simplified version of 21CMFACST Mesinger2007 (29, 30), and was
spectrally sampled at 80 kHz to match typical MWA EoR analysis parameters.
This methodology generates a volume that projects to a sky coverage of $\sim
50\times 50$ square degrees, which is essential given the footprint of the MWA
primary beam. A number of steps are necessary to translate the lightcone box
into a WODEN sky model. The box is a collection of 2D projected 21 cm
intensity distributions, each at a unique redshift. These cartesian
projections are taken by slicing a full 3D simulation volume at regular
redshift intervals, as the simulation is evolved with time. This leaves us
with a number of transformations needed to allow WODEN to ingest the model:

  * •

Each 2D slice has a pixel resolution $\Delta x$ in cMpc; each slice therefore
has a different angular resolution, as this quantity is redshift dependent. To
run efficiently, WODEN needs a grid of pixels constantly sampled in RA/DEC,
necessitating a redshift dependent angular interpolation.

  * •

The box is regularly sampled in redshift ($\Delta z$), which results in
irregular sampling in frequency ($\Delta\nu$). Interferometric data is sampled
regularly in frequency, necessitating a second interpolation over frequency.

  * •

The box is in units of mK. WODEN ingests units of Jansky (Jy), which is an
integrated flux density. The conversion from mK to Jy/sterrad is
straightforward, however, the physical volume of each pixel changes as a
function of frequency. Along with the angular extent, the pixel volume is
determined by $\Delta z$. Given the interpolation over angle and frequency,
this volume change must be taken into account, as it effects the resultant
variance of the map. This effect is noted in the script
make_flat_spectrum_eor.py555https://github.com/RadioAstronomySoftwareGroup/pyradiosky/blob/main/scripts/make_flat_spectrum_eor.py
in the Python package pyradiosky.

In the following, we use the cosmology assumed by Greig2022 (18), a
$\Lambda$CDM cosmology with:
$H_{0}=68;\Omega_{M}=0.31;\Omega_{\Lambda}=0.69;\Omega_{b}=0.048$. We use
astropy astropy:2013 (5, 3, 4) for all cosmological calculations.

We choose to interpolate all slices to the finest angular resolution in the
box to retain as much angular structure as possible. We call this the
reference redshift, $z_{\mathrm{ref}}=7.5693$. Given $\Delta\theta=\Delta
x/\mathcal{D}(z)$, where $\mathcal{D}$ is the comoving distance, we
interpolate to an angular resolution of $\sim 27$ arcseconds. We found
interpolating the Cartesian slice directly into a TAN FITS projection
(FITS2002, 9, a gnomic projection; see) centred at RA, DEC
$=0^{h},-27^{\circ}$ returned the expected PS. We experimented with bilinear
and bicubic interpolation, but found significant signal loss in the resultant
PS. It’s possible some variant of Gaussian Process Regression may be more
effective, however for the purposes of this work we found a simple nearest-
neighbour approximation was sufficient. Similarly, we applied a nearest-
neighbour interpolation along the frequency axis, as again we saw signal loss
in the final PS when applying linear or cubic interpolation. Along lines of
sight, the signal rapidly fluctuates between positive and zero, and so these
kinds of interpolation tend to create false signal. Once these interpolations
have been applied, we then scale for the change in volume of the pixel
($\mathrm{V}_{\mathrm{pix}}(z)$) and its effect on the variance. We do this
for each redshift by calculating via the differential comoving volume and
multiplying by $\Delta z$. We then scale each redshift slice by a factor
$C_{\mathrm{pix}}$ given by

$C_{\mathrm{pix}}=\sqrt{\frac{\mathrm{V}_{\mathrm{pix}}(z)}{\mathrm{V}_{\mathrm{pix}}(z_{\mathrm{ref}})}},$
(1)

with the square root ensuring the variance is scaled by the pixel volume.

Once we have interpolated, scaled, and transformed each slice in Jy/sterrad,
we convert each pixel into a point source with units of Jy by multiplying by
the pixel solid angle. Ideally, we would the find some way to tile these maps
to give an all-sky sky model. However, given we cannot interpolate the box
without signal loss, it is computationally infeasible. At the lowest
frequencies, the number of point sources necessary would approach 200 million.
However, given the extreme volume of the simulated box, without tiling, the
sky model already covers the main lobe of the MWA primary beam down to the
$~{}1$% power level for pointings of interest (see Figure 1).

Figure 1: Demonstration of the MWA primary beam and its interaction with the
sky model, via the: a) -2 pointing; b) -1 pointing; c) Zenith pointing. Each
plot is locked to the observer in HA/Dec. The solid black lines represent the
instrumental Stokes I primary beam pattern contoured at 1%, 10%, 50%, and 90%
power levels. The coloured images show the apparent sky model after
attenuation by the primary beam at three different LSTs, with the
corresponding colour lines marking the edges of the full sky model. Both the
beam and sky model are shown at 167 MHz, where the primary beam is largest for
a high-band observation. Note that the +1, +2, pointings as described in
Section 5 are simply westward reflections of the -1, -2 pointings, so aren’t
reproduced here for brevity.

To investigate whether the projection steps detailed in this Section induce
any form of signal loss or bias, we apply them to a purely Gaussian noise
simulation in Section 4. We report of the results of the sky model in Section
5.

##  4 White Gaussian noise sky simulation

The power spectrum $P(k)$ of a White Gaussian noise distribution $\mathcal{N}$
of mean $\mu=0$ and variance $\sigma^{2}$ is proportionate to the variance,
i.e. $P(k)\propto\sigma^{2}$. By creating a sky model of purely White Gaussian
noise, we can therefore predict an output $P(k)$, and check the normalisation
of the AusEoRPipe. Given we have outlined a method to project a box from units
of cMpc in Section 3, we can simply repeat the entire process, starting from a
box of purely White Gaussian noise (from hereon in referred to as the Gaussian
noise simulation). We choose to set $\sigma^{2}=16\textrm{mK}^{2}$ as this
yields a power comparable to those expected from the signal. We typically
report the power spectrum in units of mK2 Mpc3 $h^{-3}$, where $h$ is the
reduced Hubble constant; we use $h=0.68$ in this work. We can therefore
predict our expected noise power spectrum $P_{N}$ via

$P_{N}(k)=\sigma^{2}\mathrm{V}_{\mathrm{pix}}(z_{\mathrm{ref}})h^{3}=8.923\,\mathrm{mK}^{2}\mathrm{Mpc}^{3}h^{-3}$
(2)

We use WODEN 666For this paper we use the docker v2.0 container of WODEN,
installed via docker pull docker://jlbline/woden-2.0 to simulate a single
zenith pointing observation with an 8 s time cadence and 80 kHz frequency
resolution, using a frequency-interpolated version of the FEE beam. We use the
convention that Stokes I = (XX + YY) / 2, where XX and YY are the two linear
polarisations. Given we are using simulated data with isotropic and
unpolarised sky models, there is little difference between XX and YY, and
therefore the power spectra shown are effectively Stokes I. All PS shown in
this paper are from the north-south aligned polarisation. We run CHIPS on the
simulated visibilities, to produce the 2D PS as shown in Figure 2. All CHIPS
plots were generated using the Python package
CHIPS_wrappers777https://github.com/JLBLine/CHIPS_wrappers.

Figure 2 demonstrates a known effect of estimating the PS from gridded
visibilities, where the power inferred is lower where the density of gridded
visibilities is higher, as observed by Barry2019_pipeline (6). Note we have
not corrected for this in Figure 2, which is normally done by multiplying by a
factor two (as found by Barry2019_pipeline (6)), to explicitly show the
effect.

Figure 2: Data from a single zenith observation, without any correction for
gridding density, where: a) shows the 2D PS; b) the ratio to the expected
value; c) the CHIPS gridding weights. $k_{\perp}$ modes are derived from the
instrument sampling in the visibility $u,v$ plane, averaged down to one
dimension, and $k_{\parallel}$ from the fourier transform of the visibilities
along frequency. The CHIPS gridding weights therefore show the $u,v$ gridding
density averaged into one dimension.

Visibilities are gridding using a kernel which is designed to minimise
aliasing and match the area of the MWA primary beam Trott2016 (42). Given the
array layout, and the need to match the kernel to the primary beam, the kernel
has an unavoidable footprint large enough to overlap with neighbouring gridded
visibilities. Depending on the covariance between neighbouring visibilities,
and the density of the gridded visibilities, the power estimated will
therefore vary as a function of $u,v$ gridding location. Using simulations,
Barry2019_pipeline (6) found a factor of two correction was sufficient as a
normalisation factor for modes of interest. This estimate was made using the
FHD/$\varepsilon$ppsilon pipeline and a model of the MWA primary beam as a
gridding kernel, whereas CHIPS uses a Blackman-Harris kernel. This estimate
was also made using 2 s resolution data, rather than 8 s. Given these
differences, we investigate the gridding density correction factor in Section
4.1 using CHIPS.

###  4.1 Gridding density correction factor

The gridding density is directly affected by the $u,v$ coverage, which in turn
is dictated by the array layout and phase centre. The covariance of
neighbouring $u,v$ points will therefore be affected by LST (as the phase
centre is always set to EoR0 field centre). The primary beam pointing will
also have an effect as this changes the amplitude and spectral behaviour of
the visibilities. To thoroughly investigate the gridding density we therefore
simulate a number of observations spanning realistic ranges of LST and
pointings. To do so, we take the first (subset LST1), central (subset LST2),
and final observation (subset LST3) across five different pointings (-2 to +2)
from three real nights of MWA phase I EoR observing (see Figure 3). For the
LST2 subset, we select the observation where the primary beam pointing is
closest to the EoR0 field centre, to maximise the beam coverage over the sky
model coverage. As the visibilities are normally averaged to $8\,$s, $80\,$kHz
we simulate at this cadence to save compute time. Each observation takes a
total of $\sim 32\,$GPU hours. For 45 simulated observations this is a total
of 60 GPU days. We use the Pawsey Garrawarla
cluster888https://pawsey.org.au/systems/garrawarla/, and split each simulation
across 24 GPUs, meaning these simulations take less than 3 days real-time.

Figure 3: Observations simulated in this paper, based on two real nights of
MWA EoR0 observing. The blue diamonds show the simulated LST1 subset, green
hexagons the LST2 subset, and orange circles the LST3 subset. Each hollow
square shows a different two minute snapshot which was not simulated but
exists in the real data set. Dividing lines and labels show the changes from
pointings -2 through +2. Note observational constraints mean there are less +2
pointings.

By running CHIPS on various combinations of these observations, we can
investigate both primary beam pointing and LST effects. All integrations over
multiple observations are done coherently, i.e. all observations are gridded
to the same $u,v$ grid, which is then used to estimate the PS. We take the
median of the ratio of the observed PS to the expected value as a function of
$k_{\perp}$ (i.e. take the median along the y-axis of the middle panel in
Figure 2) to see the effects of gridding, as shown in Figure 4. We take the
median, rather than the mean, as the distribution of ratios along each
$k_{\perp}$ bin display significant skew.

Figure 4: The median ratio of the recovered noise PS to expected value, as a
function of $k_{\perp}$, for: a) the five different pointings, each for a
single observation; b) the three different LST subsets integrated over 5
observations; c) the three different LST subsets integrated over 15
observations; d) all three LST subsets integrated together over 15, 30, and 45
observations. The horizontal dashed line shows the mean of the median ratios
for $k_{\perp}<0.1$ as measured from the bottom right plot. This dashed line
Any dataset which is a multiple of five observations has an even split across
the five pointings.

Figure 4a shows that for a single observation, the estimated power is
consistent across pointings, with lower power estimated a low $k_{\perp}$
(higher gridding density). This shows the beam volume correction applied
internally in CHIPS returns a consistent power level across the five
pointings. Figure 4b shows that the power recovered from the three LST subsets
over 5 observations is roughly consistent, but the power recovered does vary
somewhat. This recovered power depends on the exact $u,v$ coverage, and the
way the primary beam interacts with the sky model. Figure 4c shows that when
integrating over 15 observations for each LST subset, the recovered power at
low $k_{\perp}$ is consistent with the power recovered when integrating over 5
observations. However, the power recovered at higher $k_{\perp}$ is actually
higher than at when integrating over 5 observations. This highlights the fact
that the power estimated comes from a combination of not just the amount of
data gridded, but the covariance between the neighbouring gridded
visibilities. This is further highlighted by showing that integrating the
three LST subsets together across 15 observations yields significantly lower
power at high $k_{\perp}$. High $k_{\perp}$ is derived from longer baselines,
where the drift of a baseline $u,v$ coordinate over time is greatest. The LST1
and LST3 groups are also closely separated in time, yielding close gridding
locations, but overlap as the beam pointing changes. The amplitude and the
spectral behaviour of the primary beam therefore sharply changes on the sky,
reducing the covariance between neighbouring gridded visibilities. This
results in the lower power recovered. Figure 4d shows that when integrating
more and more data with the same combination of LST coverage and pointings,
there is no change in power at low $k_{\perp}$.

Given the stability of the recovered power when combining all LST subsets with
an integration of 15 observations or more, one could fit a functional form to
derive a correction. We investigate this approach in A, and find little to no
difference in using a fit when compared to just applying a scalar
normalisation factor. PS upper limits are derived from integrating hours of
data to overcome thermal noise, and modes where
$k_{\perp}>0.1\,h\,\textrm{Mpc}^{-1}$ are typically cut when using the
AusEoRPipe. To obtain a normalisation factor for the gridding density we
therefore simply take the mean of the median observed ratios for the
integrations over all LSTs and pointings where
$k_{\perp}<0.1\,h\,\textrm{Mpc}^{-1}$. This value is plotted as the horizontal
grey dashed line in Figure 4. Inverting this gives a normalisation factor of
2.651. We apply this normalisation and create a number of 1D PS for various
integrations, shown in Figure 5, and compare the outcomes to the previous
correction factor of two. This shows that the new normalisation factor causes
too much power to be recovered for a single observation, and when only using
the LST2 subset, but recovers the correct power to within 10% for most $k$
modes when integrating over all LSTs and pointings. The factor two correction
(coincidentally) does well for a single observation, but does not recover the
correct power for higher integrations.

Figure 5: 1D PS for various integrations of the Gaussian noise simulation,
where a) a factor of two was used to normalise for gridding density and b) a
factor of 2.5845 was used. The vertical dashed line shows the maximum
$k_{\perp}$-mode.

##  5 Recovered signal

To test the model detailed in Section 3, we simulate 30 observations, matching
the first two nights detailed in Figure 3, with the same parameters as used
for the Gaussian noise simulation (detailed in Section 4). We leave out
simulating the third night purely to save on computational resources. To
qualitatively assess the accuracy of the model through WODEN and into CHIPS,
we compare a CHIPS 2D PS to one derived directly from the lightcone box, shown
in Figure 6. We produce the 2D PS directly from the lightcone using the OSIRIS
package999https://github.com/JaidenCook/OSIRIS Cook2022 (12), where no primary
beam or instrument sampling was applied, meaning the entire $u,v$ plane was
sampled. We see that the CHIPS PS is broadly consistent with the OSIRIS PS.
Noticeably, it seems there is a potential signal loss at high $k_{\parallel}$
in the CHIPS PS (the darker blue region at the top of Figure 6b). There is
also excess power seen along high $k_{\perp}$, around the dashed black line,
which is expected as the instrument baseline sampling causes mode mixing,
moving power up from lower to higher $k_{\parallel}$.

To test the gridding normalisation, we apply it to various integrations and
produce 1D PS, shown in Figure 7. We compare them to an expected signal, again
derived directly from the lightcone box. Excellent agreement is shown, and
similarly to the Gaussian noise simulation, the gridding density normalisation
estimates too much power for short integrations, but renders consistent
results for longer integrations with the same mix of LSTs and pointings.

In Figure 8 we show the ratio of the final recovered signal to the expected
value, for a 60 min integration. While good agreement is shown, there does
seem to be signal loss around $0.1<k<1$, which was not seen for the Gaussian
noise simulation. This is likely due to the frequency interpolation of the
original lightcone as detailed in Section 3. For reference, Figure 8 shows the
maximum $k_{\perp}$ mode that went into the average from 2D to 1D. $k$ modes
above this are binned from only high $k_{\parallel}$ modes, which is where
signals with the most spectral structure manifest. As this area is where the
signal loss is found, it is likely the frequency interpolation is causing
this. As the Gaussian noise sky model has no correlation over frequency, it
should not suffer from this interpolation, which is why it doesn’t display the
same signal loss.

Figure 8 also shows the effect of performing the ‘wedge-cut’, where modes
normally dominated by foreground sources are removed. While broadly
consistent, the signal loss seems to worsen after the cuts. This is likely due
to the reduced number of bins going into the average from 2D to 1D. As can be
seen from Figure 6c, the horizon cut (marked by the solid black line) removes
more samples for smaller $k_{\parallel}$ modes. However, lower sampling is not
the only possible reason. The manifestation of the primary beam, as well as
the exact $u,v$ distribution, could both cause a systematic bias at specific
spots in the 2D PS, which could contribute here.

Figure 6: 2D PS from the simulation where: a) produced using OSIRIS directly
from the lightcone box from Greig2022 (18); b) a CHIPS 2D PS made from
integrating over 30 simulated observations; c) the same PS from (b) but with
standard cuts made to remove foreground contamination. The solid black line
indicates an estimate of the horizon, and the dashed black line the full-width
half-max of the MWA primary beam. Figure 7: Recovered signal from various
integration lengths of simulated observations, when normalised with a factor
2.5845. The panel bottom right is included to illustrate the normalisation
effects between $0.1<k<1.0\,h\,\textrm{Mpc}^{-1}$. The two minute data set is
for a zenith pointing; the 10 minute data set is for the LST2 subset; all
other data sets are an even split between the five pointings and LST subsets
as described in Section 5. Figure 8: Ratios of the recovered signal to the
expected value, for a 60 min integration. The orange with crosses line shows
the ratio when using the factor two gridding density correction, and the blue
with circles when using the fitted correction. The vertical dashed line shows
the maximum $k_{\perp}$-mode, meaning any $k$-mode above this was derived
purely from the Fourier transform of the visibilities along frequency.

##  6 Discussion

The Gaussian noise simulation results in Section 4 confirm that gridding
visibilities without accounting for their full covariance from overlapping
footprints requires a normalisation factor. They reveal however that this
normalisation factor depends on an interplay between the array layout, the LST
range of the integrated observations, as well as the primary beam pointings,
as these all effect the covariance between the visibilities. To truly correct
for this effect, one would need to calculate this covariance, and propagate it
through the gridding and PS estimation steps (Liu2014, 27, see). We instead
simply measure and fit for the effect, given we can estimate the expected
power from the Gaussian noise simulation. It could be argued that simply
taking the median of the ratios where $k_{\perp}<0.1\,h\,\textrm{Mpc}^{-1}$
from Figure 4d is too simplistic, given we can fit a functional form as shown
in Appendix A. We choose to err on the side of simplicity however, given the
single scaling factor we derive recovers the correct Gaussian-noise power to
within 10% for most $k$-modes. This normalisation factor comes with a set of
caveats, which require further investigation to be understood:

  * •

We have only tested the MWA phase I layout. Given the covariance between $u,v$
may well be different for different array layouts, we cannot be sure this
correction factor would work with the MWA phase II compact layout, for example

  * •

Upon initial investigation, 2 s resolution data was seen to require a smaller
normalisation factor for the same integration length. It is possible that
gridding four more times more data will eventually reach the same plateau as 8
s data, but it might change the level of covariance in the gridded
visibilities. Testing this will require four times the computational resources

  * •

It is possible that the gridding normalisation will change for a different
frequency range, as the rate of change of the $u,v$ points with time will be
different. To test this with the method outlined in this paper, another
lightcone box will need to be generated to match the required frequency
coverage

  * •

We attribute the signal loss seen in Figure 8 to the frequency interpolation
of the original lightcone box. However, there could be other sources of signal
loss, which include the primary beam shape altering the signal distribution,
and a dependence on details of the gridding within a $k_{\bot}$ bin. To remove
frequency interpolation and check this signal loss issue, the original
lightcone box would need to be regenerated with constant frequency sampling.
This is non-trivial, as it would require 21CMFAST simulations with adaptable
voxel sizes. Alternatively, one could run the noise sky simulation test using
a type of noise with a spectral structure (e.g. Brownian noise). As long as
the shape of the noise could be easily predicted in PS space, one could test
for signal loss

It should be noted that this method depends on not only the Gaussian noise
model being correctly scaled, but that our prediction for the power level is
correct. However, given the normalisation also yields the expected power level
for the model, with that prediction made directly from the lightcone, this
gives us confidence that our prediction and scaling are correct. We further
note that we ran the WODEN simulations in so-called float mode, which sets
some of the internal precision to single instead of double precision, with an
appreciable speed up (Line2022, 26, see). We ran a number of comparisons to
fully double simulations all the way through to the PS and saw no change in
the results.

It is worth considering that gridding density is not necessarily the only
effect that could be causing a loss in signal. The model itself, the WODEN
simulator, and CHIPS are potentially sources of signal loss. However, given
that the Gaussian noise sky model is as simple a model as possible, WODEN has
been heavily unit tested101010See
https://woden.readthedocs.io/en/latest/testing/cmake_testing.html for details
on testing, and we are fitting a CHIPS-specific normalisation factor, we are
confident the normalisation derived here is applicable to real data put
through the AusEoRPipe.

##  7 Conclusion

We have projected a lightcone into an celestial coordinate-locked sky model
for use with the WODEN simulator. The model is stored in a FITS table format
easily adaptable to other simulation packages, and is available from the PASA
Datastore. In testing the validity of this model, we generated an equivalent
White Gaussian noise model. We used this to test the effects of gridding
visibilities using a kernel to produce power spectra, through simulated
observations. We confirmed the known effect that power estimated from gridded
visibilities is lower where the gridding density is higher. We found the
effect to be dependent on a combination of array layout, primary beam
pointing, and LST of integrated observations. We also observed that the
reduction in power estimated plateaus after integrating more than 30 minutes
of data. We found a single scalar normalisation factor of 2.651 was sufficient
to recover the expected power to within 10% for most $k$-modes in a 1D PS. We
then generated simulated observations of the sky model, and verified that
CHIPS was able to recover the expected signal. We found that frequency
interpolation of the original lightcone box causes some signal loss at higher
$k$, but with the 2.651 normalisation factor, the recovered power is again
close to within 10% of the expect signal. We conclude that this sky model is
sufficient to accurately test the AusEoRPipe, and can be used in conjunction
with foreground models and simulated instrumental effects to test calibration.

This scientific work uses data obtained from Inyarrimanha Ilgari Bundara / the
Murchison Radio-astronomy Observatory. We acknowledge the Wajarri Yamaji
People as the Traditional Owners and native title holders of the Observatory
site. Establishment of CSIRO’s Murchison Radio-astronomy Observatory is an
initiative of the Australian Government, with support from the Government of
Western Australia and the Science and Industry Endowment Fund. Support for the
operation of the MWA is provided by the Australian Government (NCRIS), under a
contract to Curtin University administered by Astronomy Australia Limited.
This work was supported by resources provided by the Pawsey Supercomputing
Research Centre with funding from the Australian Government and the Government
of Western Australia.

During this work we made extension use of the kvis Gooch1996 (17) and DS9
Joye2003 (22) FITS file image viewers.

#### Funding Statement

This research was supported by the Australian Research Council Centre of
Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through
project number CE170100013. CMT is supported by an ARC Future Fellowship under
grant FT180100321.

#### Data Availability Statement

The and White Gaussian noise sky models are available in the PASA Datastore in
a FITS table format. The CHIPS codebase is available upon request to the
authors; all other code bases are linked throughout the paper.

##  Appendix A Fitting for gridding density correction factor

As noted in Section 4.1, the power estimated from the Gaussian noise
simulations is consistent for integrations of 15, 30, and 45 observations
(when using the same combination of LSTs and pointings). We can therefore
correct for gridding density by fitting the observed ratio of measured to
expected power as a function of $k_{\perp}$. We fit a broken power-law (see
Figure 9) as defined by111111See
https://gist.github.com/cgobat/12595d4e242576d4d84b1b682476317d for more
information on the functional form

$A\frac{k_{\perp}}{k_{\mathrm{break}}}^{-\alpha_{1}}\left[0.5\left(1+\left(\frac{k_{\perp}}{k_{\mathrm{break}}}\right)^{\frac{1}{\Delta}}\right)\right]^{(\alpha_{1}-\alpha_{2})\Delta},$
(3)

where $A$ is an amplitude factor, $k_{\mathrm{break}}$ is the $k_{\perp}$
value where the power-law breaks, $\alpha_{1}$ is the power-law index for
$k_{\perp}<k_{\mathrm{break}}$, $\alpha_{2}$ is the power-law index for
$k_{\perp}>k_{\mathrm{break}}$, and $\Delta$ is the sharpness of the break. We
set $\alpha_{1}=0$ as the gridding density plateaus at low $k_{\perp}$, and
set $A=0.37728$, the median of the ratios where
$k_{\perp}<0.1\,h\,\textrm{Mpc}^{-1}$ from Figure 4d. We obtain $\alpha_{2}$
by performing a least squares fit with a single power-law model
$Ck_{\perp}^{\alpha_{2}}$, where $C$ is some constant, on ratios where
$k_{\perp}>0.07\,h\,\textrm{Mpc}^{-1}$. We can then find $k_{\mathrm{break}}$
by finding where the power-law is equal to the amplitude $A$:

$\displaystyle A=Ck_{\mathrm{break}}^{\alpha_{2}},$ (4) $\displaystyle
k_{\mathrm{break}}=10^{\frac{\log_{10}(A)-C}{\alpha_{2}}}.$ (5)

We fit this function to the median value of the ratios for each $k_{\perp}$
bin, as we found the values in each bin to have non-zero skew. Taking the
median should therefore be a better representation of the data rather than the
mean. We summarise the fitted parameters in Table 1. To create a normalisation
function, we simply take the inverse of this broken power-law.

Figure 9: Fitting the median recovered ratio of the noise simulation as a function of $k_{\perp}$. Circles, squares, and triangles show the median ratio for integrations of 15, 30, and 45 observations respectively. The shaded regions are bound by the median absolute deviation. The dashed black line shows the fitted broken power-law. Table 1: Parameters for broken power-law fit Parameter |  Value
---|---
$a_{1}$ |  0.0
$\Delta$ |  0.01
$A$ |  0.377280447218
$a_{2}$ |  -0.161193559625
$k_{\mathrm{break}}$ |  0.064354113566

We compare how well this fit predicts the observed ratios in Figure 10. Here
we plot the ratios as a function of the gridding weights, rather than
$k_{\perp}$. We use the seaborn Waskom2021 (43)
kdeplot121212https://seaborn.pydata.org/generated/seaborn.kdeplot.html
function to generate two dimensional kernel density estimates (KDE) to show
the distribution of observed ratios in blue. We also plot the KDE of the
ratios predicted by the fitted broken power-law in orange. We see that the
fitted function predicts values close to the median value at high gridding
weights. As data with a high gridding weight is naturally up-weighted, these
data points contribute more to the estimated power in the 1D PS. The net
result of this is shown in Figure 11, where we compare the recovered 1D PS
from the Gaussian noise simulation when correcting with the single scalar
normalisation factor, and the fitted functional form. No difference is seen at
low $k$, with minimal difference seen at high $k$. We conclude that as the
functional form provides little difference to any limits, the scalar
normalisation factor is sufficient for correcting for gridding density.

Figure 10: Kernel density estimates of the recovered ratio of the noise
simulation as a function of CHIPS gridding weights (blue), and those predicted
by fitted broken power-law (orange) for: a) the 15 observation integration of
all LSTS; b) 30 observations; c) 45 observations.

Figure 11: The recovered 1D PS from the Gaussian noise simulation when
correcting with the single scalar normalisation factor, and the fitted
functional form as shown in Figure 9

.

## References

  * (1) Anshuman Acharya et al.  “21-cm Signal from the Epoch of Reionization: A Machine Learning upgrade to Foreground Removal with Gaussian Process Regression”  In _MNRAS_ , 2023  DOI: 10.1093/mnras/stad3701
  * (2) James E. Aguirre et al.  “Validation of the HERA Phase I Epoch of Reionization 21 cm Power Spectrum Software Pipeline”  In _ApJ_ 924.2, 2022, pp. 85  DOI: 10.3847/1538-4357/ac32cd
  * (3) Astropy Collaboration et al.  “The Astropy Project: Building an Open-science Project and Status of the v2.0 Core Package”  In _AJ_ 156.3, 2018, pp. 123  DOI: 10.3847/1538-3881/aabc4f
  * (4) Astropy Collaboration et al.  “The Astropy Project: Sustaining and Growing a Community-oriented Open-source Project and the Latest Major Release (v5.0) of the Core Package”  In _apj_ 935.2, 2022, pp. 167  DOI: 10.3847/1538-4357/ac7c74
  * (5) Astropy Collaboration et al.  “Astropy: A community Python package for astronomy”  In _A &A_ 558, 2013, pp. A33  DOI: 10.1051/0004-6361/201322068
  * (6) N. Barry et al.  “The FHD/Eppsilon Epoch of Reionisation power spectrum pipeline”  In _PASA_ 36, 2019, pp. e026  DOI: 10.1017/pasa.2019.21
  * (7) N. Barry et al.  “Improving the Epoch of Reionization Power Spectrum Results from Murchison Widefield Array Season 1 Observations”  In _ApJ_ 884.1, 2019, pp. 1  DOI: 10.3847/1538-4357/ab40a8
  * (8) A.. Beardsley et al.  “First Season MWA EoR Power spectrum Results at Redshift 7”  In _ApJ_ 833.1, 2016, pp. 102  DOI: 10.3847/1538-4357/833/1/102
  * (9) M.. Calabretta and E.. Greisen  “Representations of celestial coordinates in FITS”  In _A &A_ 395, 2002, pp. 1077–1122  DOI: 10.1051/0004-6361:20021327
  * (10) J. Chege et al.  “Optimising MWA EoR data processing for improved 21-cm power spectrum measurements—fine-tuning ionospheric corrections”  In _PASA_ 39, 2022, pp. e047  DOI: 10.1017/pasa.2022.34
  * (11) A. Chokshi et al.  “Dual polarization measurements of MWA beampatterns at 137 MHz”  In _MNRAS_ 502.2, 2021, pp. 1990–2004  DOI: 10.1093/mnras/stab156
  * (12) J.. Cook, C.. Trott and J… Line  “Investigating the contribution of extended radio sources to the Epoch of Reionization power spectrum”  In _MNRAS_ 514.1, 2022, pp. 790–805  DOI: 10.1093/mnras/stac1330
  * (13) David R DeBoer et al.  “Hydrogen Epoch of Reionization Array (HERA)”  In _PASP_ 129.974, 2017, pp. 045001  DOI: 10.1088/1538-3873/129/974/045001
  * (14) H.. Edler, F. de Gasperin and D. Rafferty  “Investigating ionospheric calibration for LOFAR 2.0 with simulated observations”  In _A &A_ 652, 2021, pp. A37  DOI: 10.1051/0004-6361/202140465
  * (15) A. Ewall-Wice et al.  “First limits on the 21 cm power spectrum during the Epoch of X-ray heating”  In _MNRAS_ 460.4, 2016, pp. 4320–4347  DOI: 10.1093/mnras/stw1022
  * (16) Steven R. Furlanetto, S. Peng Oh and Frank H. Briggs  “Cosmology at low frequencies: The 21cm transition and the high-redshift Universe”  In _Physics Reports_ 433.4, 2006, pp. 181–301  DOI: https://doi.org/10.1016/j.physrep.2006.08.002
  * (17) Richard Gooch  “Karma: a Visualization Test-Bed”  In _Astronomical Data Analysis Software and Systems V_ 101, Astronomical Society of the Pacific Conference Series, 1996, pp. 80
  * (18) Bradley Greig et al.  “Generating extremely large-volume reionization simulations”  In _MNRAS_ 516.4, 2022, pp. 5588–5600  DOI: 10.1093/mnras/stac2506
  * (19) M.. Haarlem et al.  “LOFAR: The LOw-Frequency ARray”  In _A &A_ 556, 2013, pp. A2  DOI: 10.1051/0004-6361/201220873
  * (20) HERA Collaboration et al.  “Improved Constraints on the 21 cm EoR Power Spectrum and the X-Ray Heating of the IGM with HERA Phase I Observations”  In _ApJ_ 945.2, 2023, pp. 124  DOI: 10.3847/1538-4357/acaf50
  * (21) Daniel C. Jacobs et al.  “THE MURCHISON WIDEFIELD ARRAY 21 cm POWER SPECTRUM ANALYSIS METHODOLOGY”  In _ApJ_ 825.2, 2016, pp. 114  DOI: 10.3847/0004-637X/825/2/114
  * (22) W.. Joye and E. Mandel  “New Features of SAOImage DS9”  In _Astronomical Data Analysis Software and Systems XII_ 295, Astronomical Society of the Pacific Conference Series, 2003, pp. 489
  * (23) Nicholas S. Kern et al.  “Mitigating Internal Instrument Coupling for 21 cm Cosmology. II. A Method Demonstration with the Hydrogen Epoch of Reionization Array”  In _ApJ_ 888.2, 2020, pp. 70  DOI: 10.3847/1538-4357/ab5e8a
  * (24) Matthew Kolopanis, Jonathan C. Pober, Daniel C. Jacobs and Samantha McGraw  “New EoR power spectrum limits from MWA Phase II using the delay spectrum method and novel systematic rejection”  In _MNRAS_ 521.4, 2023, pp. 5120–5138  DOI: 10.1093/mnras/stad845
  * (25) J… Line et al.  “Modelling and peeling extended sources with shapelets: A Fornax A case study”  In _PASA_ 37, 2020, pp. e027  DOI: 10.1017/pasa.2020.18
  * (26) Jack Line  “WODEN: A CUDA-enabled package to simulate low-frequency radio interferometric data”  In _The Journal of Open Source Software_ 7.69, 2022, pp. 3676  DOI: 10.21105/joss.03676
  * (27) Adrian Liu, Aaron R. Parsons and Cathryn M. Trott  “Epoch of reionization window. II. Statistical methods for foreground wedge reduction”  In _Phys. Rev. D_ 90.2, 2014, pp. 023019  DOI: 10.1103/PhysRevD.90.023019
  * (28) F.. Mertens et al.  “Improved upper limits on the 21 cm signal power spectrum of neutral hydrogen at z $\approx$ 9.1 from LOFAR”  In _MNRAS_ 493.2, 2020, pp. 1662–1685  DOI: 10.1093/mnras/staa327
  * (29) Andrei Mesinger and Steven Furlanetto  “Efficient Simulations of Early Structure Formation and Reionization”  In _ApJ_ 669.2, 2007, pp. 663–675  DOI: 10.1086/521806
  * (30) Andrei Mesinger, Steven Furlanetto and Renyue Cen  “21CMFAST: a fast, seminumerical simulation of the high-redshift 21-cm signal”  In _MNRAS_ 411.2, 2011, pp. 955–972  DOI: 10.1111/j.1365-2966.2010.17731.x
  * (31) M. Mevius et al.  “A numerical study of 21-cm signal suppression and noise increase in direction-dependent calibration of LOFAR data”  In _MNRAS_ 509.3, 2022, pp. 3693–3702  DOI: 10.1093/mnras/stab3233
  * (32) Daniel a. Mitchell et al.  “Real-Time Calibration of the Murchison Widefield Array”  In _ISTSP_ 2.5, 2008, pp. 707–717  DOI: 10.1109/JSTSP.2008.2005327
  * (33) I.. Morrison et al.  “MWAX: A new correlator for the Murchison Widefield Array”  In _PASA_ 40, 2023, pp. e019  DOI: 10.1017/pasa.2023.15
  * (34) S. Munshi et al.  “First upper limits on the 21 cm signal power spectrum from cosmic dawn from one night of observations with NenuFAR”  In _A &A_ 681, 2024, pp. A62  DOI: 10.1051/0004-6361/202348329
  * (35) A.. Offringa  “AOFlagger: RFI Software”, Astrophysics Source Code Library, record ascl:1010.017, 2010, pp. ascl:1010.017  ASCL: 1010.017
  * (36) A.. Offringa, F. Mertens and L… Koopmans  “The impact of interference excision on 21-cm epoch of reionization power spectrum analyses”  In _MNRAS_ 484.2, 2019, pp. 2866–2875  DOI: 10.1093/mnras/stz175
  * (37) A.. Offringa et al.  “The Low-Frequency Environment of the Murchison Widefield Array: Radio-Frequency Interference Analysis and Mitigation”  In _PASA_ 32, 2015, pp. e008  DOI: 10.1017/pasa.2015.7
  * (38) M. Rahimi et al.  “Epoch of reionization power spectrum limits from Murchison Widefield Array data targeted at EoR1 field”  In _MNRAS_ 508.4, 2021, pp. 5954–5971  DOI: 10.1093/mnras/stab2918
  * (39) M. Sokolowski et al.  “Calibration and Stokes Imaging with Full Embedded Element Primary Beam Model for the Murchison Widefield Array”  In _PASA_ 34, 2017, pp. e062  DOI: 10.1017/pasa.2017.54
  * (40) S… Tingay et al.  “The Murchison Widefield Array: The Square Kilometre Array Precursor at Low Radio Frequencies”  In _PASA_ 30, 2013, pp. e007  DOI: 10.1017/pasa.2012.007
  * (41) Cathryn M. Trott et al.  “Deep multiredshift limits on Epoch of Reionization 21 cm power spectra from four seasons of Murchison Widefield Array observations”  In _MNRAS_ 493.4, 2020, pp. 4711–4727  DOI: 10.1093/mnras/staa414
  * (42) Cathryn M. Trott et al.  “CHIPS: The Cosmological HI Power Spectrum Estimator”  In _ApJ_ 752.2, 2016, pp. 34  DOI: 10.3847/0004-637X/818/2/139
  * (43) Michael L. Waskom  “seaborn: statistical data visualization”  In _Journal of Open Source Software_ 6.60 The Open Journal, 2021, pp. 3021  DOI: 10.21105/joss.03021
  * (44) Randall B. Wayth et al.  “The Phase II Murchison Widefield Array: Design overview”  In _PASA_ 35, 2018, pp. e033  DOI: 10.1017/pasa.2018.37
  * (45) Michael J. Wilensky et al.  “Absolving the SSINS of Precision Interferometric Radio Data: A New Technique for Mitigating Faint Radio Frequency Interference”  In _PASP_ 131.1005, 2019, pp. 114507  DOI: 10.1088/1538-3873/ab3cad
  * (46) Michael J. Wilensky et al.  “Evidence of Ultrafaint Radio Frequency Interference in Deep 21 cm Epoch of Reionization Power Spectra with the Murchison Wide-field Array”  In _ApJ_ 957.2, 2023, pp. 78  DOI: 10.3847/1538-4357/acffbd
  * (47) P. Zarka, J.. Girard, M. Tagger and L. Denis  “LSS/NenuFAR: The LOFAR Super Station project in Nançay”  In _SF2A-2012: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics_ , 2012, pp. 687–694