The halo bispectrum as a sensitive probe of massive neutrinos and baryon physics

The power spectrum has been a workhorse for cosmological studies of large-scale structure. However, the present-day matter distribution is highly non-Gaussian and significant cosmological information is also contained in higher-order correlation functions. Meanwhile, baryon physics (particularly AGN feedback) has previously been shown to strongly affect the two-point statistics but there has been limited exploration of its effects on higher-order functions to date. Here we use the BAHAMAS suite of cosmological hydrodynamical simulations to explore the effects of baryon physics and massive neutrinos on the halo bispectrum. In contrast to matter clustering which is suppressed by baryon physics, we find that the halo clustering is typically enhanced. The strength of the effect and the scale over which it extends depends on how haloes are selected. On small scales (k>1 $h$ Mpc$^{-1}$, dominated by satellites of groups/clusters), we find that the bispectrum is highly sensitive to the efficiency of star formation and feedback, making it an excellent testing ground for galaxy formation models. We show that the effects of feedback and the effects of massive neutrinos are largely separable (independent of each other) and that massive neutrinos strongly suppress the halo bispectrum on virtually all scales up to the free-streaming length (apart from the smallest scales, where baryon physics dominates). The strong sensitivity of the bispectrum to neutrinos on the largest scales and galaxy formation physics on the smallest scales bodes well for upcoming precision measurements from the next generation of wide-field surveys.


INTRODUCTION
Providing a quantitative explanation for the formation and evolution of large-scale structure (LSS) in the Universe is one of the key tests of modern cosmological theories. The standard model of cosmology, the so-called ΛCDM model, has been remarkably successful in reproducing a wide range of observations, including the observed properties of the cosmic microwave background (e.g. Planck Collaboration et al. 2020) and low-redshift probes such as baryon acoustic oscillations (e.g., Eisenstein et al. 2005) and redshift-space distortions (see e.g. Alam et al. 2017). While ΛCDM has been shown to reproduce most current LSS observations well, the physical nature of both dark matter and dark energy remain elusive. Furthermore, there have been a number of recent mild tensions reported in the best-fit parameter values for certain cosmological parameters, including the Hubble constant (see Verde et al. 2019 for a review) and the LSS parameter 8 ≡ 8 √︁ Ω /0.3 (where 8 is defined as the amplitude of the (linear) power spectrum on the scale of 8 ℎ −1 Mpc and Ω is the present-day matter density parameter) as derived from measurements of cosmic shear, galaxy clustering, and the abundance of massive galaxy clusters (see, e.g., discussion in McCarthy et al. 2018). Whether these tensions are hinting at the presence of (unaccounted for) systematic errors in some of the cosmological analyses, that there is new physics beyond the standard model, or that they represent moderately large statistical fluctuations, is presently unclear.
The drive to test ΛCDM further, and to look for new physics beyond the standard model, requires that we make increasingly precise measurements of 'tried and tested' LSS statistics (e.g., 2-point clustering of galaxy clustering, weak lensing, etc.), but also that we devise new tests of LSS which are capable of probing different aspects of the matter distribution, which in turn could hopefully break important degeneracies between the fitted cosmological parameters. The halo bispectrum (or 3-point function of haloes/galaxies), which we focus on in the present study, is one such example of a LSS test that is becoming increasingly observationally feasible (e.g., Peebles & Groth 1975;Gaztanaga 1994;Scoccimarro et al. 2001;Verde et al. 2002;Gaztañaga et al. 2005 Slepian et al. 2017a,b;Pearson & Samushia 2018;Gualdi et al. 2018Gualdi et al. , 2019Veropalumbo et al. 2021) and for which previous theoretical studies have shown contains a significant source of cosmological information beyond what may be obtained from standard 2-point clustering tests alone (e.g Song et al. 2015;Yankelevich & Porciani 2019;Chudaykin & Ivanov 2019;Barreira 2020;Gualdi & Verde 2020;Moradinezhad Dizgah et al. 2021;Hahn et al. 2020;Heinrich & Doré 2020;Eggemeier et al. 2021;Hahn & Villaescusa-Navarro 2021;Agarwal et al. 2021).
In addition to the extra complexity of following the gravitational evolution of matter when density fluctuations become large, there is also the challenge of accurately accounting for nongravitational physics that comes into play. Specifically, feedback processes associated with the formation of stars and the growth of black holes are well-established, both empirically and via theoretical models and full cosmological hydrodynamical (hereafter hydro) simulations, to have a strong effect on the overall matter distribution on scales of up to a few tens of Mpc (e.g., van Daalen et al. 2011;Schneider & Teyssier 2015;Mummery et al. 2017;Chisari et al. 2019;van Daalen et al. 2020). Previous studies that focused on 2-point statistics have shown that, unless these effects are properly modelled, they will result in biases in the derived cosmological parameters (e.g., Semboloni et al. 2011).
Given that galaxy formation physics has been shown to significantly affect classical 2-point statistics, it is reasonable to expect that it should also affect the 3-point correlation function/bispectrum. But how large are the effects? Is the bispectrum more or less sensitive to these processes than previously studied LSS quantities? And how do we mitigate against these effects (if the goal is cosmological constraints) or exploit them (if the goal is to better understand galaxy formation)? These are the questions we focus on in the present study.
Note that, thus far, there has been remarkably little attention devoted to the effects of baryon physics on the halo bispectrum using cosmological simulations. This is plausibly explained by the fact that very large cosmological volumes (box sizes of at least several hundred Mpc on a side) are required to reliably measure the bispectrum. Combined with the constraint that the simulations must also use full hydrodynamics and a model for galaxy formation in order to self-consistently incorporate the impact of baryons, there are very few simulations that presently meet these criteria. In fact, we are aware of only two suites of hydro simulations of sufficient volume for this purpose: the BAHAMAS simulations (McCarthy et al. , 2018 and the Magneticum simulations Bocquet et al. 2016). However, to date neither have been used to examine the halo/galaxy bispectrum. Semboloni et al. (2013) presented a study of baryonic effects on the matter bispectrum using the OWLS simulations (Schaye et al. 2010). Their study also showed how the combination of the power spectrum and bispectrum can be used to break degeneracies between baryonic and cosmological effects. Foreman et al. (2020), have used BAHAMAS to also explore the effects of baryons on the matter bispectrum, finding that baryon physics does significantly affect the matter bispectrum and that the effects differ in magnitude and scale dependence compared to their effects on the matter power spectrum. Aricò et al. (2021) used the measurements of Foreman et al. (2020) to show that their 'baryonification' method for altering the outputs of collisionless (DM-only) simulations can be used to simultaneously model the effects of baryons on the matter power spectrum and the reduced matter bispectrum. We note here that, while the matter bispectrum is related to the halo bispectrum (and that the former is considerably easier to measure in simulations, given a large number of tracer particles), they are clearly distinct quantities, as haloes/galaxies are biased tracers of the matter distribution. Indeed, we will show later that the impact of baryons on halo clustering is qualitatively different (opposite sign) to that on the matter distribution. Note that, observationally, the halo/galaxy bispectrum may be measured from galaxy redshift surveys, whereas the matter bispectrum can be inferred from a combination of galaxy and weak lensing (cosmic shear) surveys (e.g., Fu et al. 2014).
In this paper, we use the BAHAMAS suite of cosmological hydro simulations to explore the impact of baryon physics (particularly active galactic nucleus (AGN) feedback) on the halo bispectrum. We also examine the role of massive neutrinos, noting that previous theoretical studies have highlighted the sensitivity of the bispectrum to this component (e.g., Hahn et al. 2020;Hahn & Villaescusa-Navarro 2021;Chen et al. 2021), and the potential degeneracy between baryons and neutrinos.
The paper is structured as follows. In Section 2 we describe the cosmological simulations used in this study. In Section 3 we describe our methods for measuring the halo power spectrum and bispectrum and, for the latter, discuss the dependence on triangular configuration. In Section 4 we present our main results, showing the dependence of the bispectrum on baryon physics, halo mass, and neutrino mass, and we explore the degeneracy between baryons and neutrinos. In Section 5 we summarise and discuss our findings.

BAHAMAS
In the present study, we use the BAHAMAS 1 cosmological simulations (McCarthy et al. , 2018 to explore the impact of baryons and neutrinos on the halo bispectrum. BAHAMAS is a suite of largevolume hydro simulations that simultaneously explores variations in baryonic physics and cosmology, including extensions to ΛCDM such as massive neutrinos (Mummery et al. 2017), dynamical dark energy , and a running scalar spectral index ). Here we use 7 runs from McCarthy et al. (2018) which explore variations in AGN feedback and the summed mass of the neutrinos in the context of a WMAP 9-year cosmology (Hinshaw et al. 2013). Table 1 presents the cosmological parameters of BAHAMAS simulations used in the current work.
Each of the hydro runs used in this study has a box size of 400 ℎ −1 Mpc and includes 2 × 1024 3 particles, equally split between baryons and dark matter, which have a mass ratio equal to Ω b /Ω cdm . (As described below, the neutrinos are not followed using particles.) In addition, for each hydro run, there is a corresponding collisionless simulation which, as described in van Daalen et al. (2020), is represented by two collisionless fluids (one for the baryons and one for the dark matter).
The Boltzmann code CAMB 2 (Lewis et al. 2000;Howlett et al. 2012) was used to compute the transfer functions and a modified version of N-GenIC was used to create the initial conditions at a starting redshift of = 127. N-GenIC was modified to include second-order Lagrangian Perturbation Theory corrections and support for massive neutrinos 3 . Separate transfer functions, computed by CAMB, are used for each individual component (baryons, neutrinos, and CDM) when producing the initial conditions. Note also that the same random phases are used for each of the simulations, such that comparisons between different runs are not subject to cosmic variance.
The simulations were carried out with a version of the Lagrangian TreePM-SPH code 3 (last described in Springel 2005), which was modified to include subgrid physics as part of the OWLS project (Schaye et al. 2010). The gravitational softening is fixed to 4 ℎ −1 kpc (in physical coordinates below = 3 and in comoving coordinates at higher redshifts) and the SPH smoothing is done using the nearest 48 neighbours.
To include the effects of massive neutrinos, BAHAMAS uses the semi-linear algorithm developed by Ali-Haïmoud & Bird (2013). The algorithm evaluates neutrino perturbations on the fly at every time step using a linear perturbation integrator, which is sourced from the full non-linear baryons+CDM potential and is added to the total gravitational force. The dynamical responses of the neutrinos to the baryons+CDM and of the baryons+CDM to the neutrinos are mutually and self-consistently included.
The hydro simulations include subgrid prescriptions for metaldependent radiative cooling, star formation and stellar evolution (including chemical enrichment from Type II and Ia supernovae and Asymptotic Giant Branch stars). The simulations include stellar feedback and prescriptions for supermassive black hole growth and AGN feedback. See Schaye et al. (2010) and references therein for a detailed description of the subgrid implementations.
As described in McCarthy et al. (2017), the parameters controlling the efficiencies of the stellar and AGN feedback were adjusted so that the simulations reproduce the present-day galaxy stellar mass function (GSMF) for * > 10 10 M and the amplitude of the gas mass fraction−halo mass relation of groups and clusters, as inferred from high-resolution X-ray observations. These quantities were selected to ensure that the collapsed structures in the simulations have the correct baryon content in a global sense, which van Daalen et al. (2020) have shown is critical for capturing the impact of baryons on the matter power spectrum. We point out that the parameters governing the feedback efficiencies were not recalibrated when varying the cosmological parameters away from the fiducial WMAP 9-yr cosmology (with massless neutrinos). But as shown in McCarthy et al. (2018), the internal properties of collapsed structures (stellar masses, gas masses, etc.) are, to first order, insensitive to the variations in cosmology. Nevertheless, we explore the potential degeneracy between baryon physics and neutrino free-streaming in Section 4.2.

Power spectrum and bispectrum definitions
The power spectrum and the bispectrum are the Fourier transforms of two-and three-point correlation functions, respectively, of the 3 https://github.com/sbird/S-GenIC dimensionless overdensity ( ) parameter, which is defined as: where ( ) is a continuous random field which gives the local density per unit comoving volume in the expanding Universe andī s the mean density with¯= ( ) (the brackets · · · denotes the averages taken over an ideal ensemble of realisations). We define the power spectrum and the bispectrum, respectively, as: where D ( ) is the three dimensional Dirac delta function with 123 = 1 + 2 + 3 , which implies that the bispectrum is defined only for closed triangles of wavevectors.
The halo power spectrum and the halo bispectrum measured from numerical simulation or observational data can be significantly affected by shot noise. Assuming the distribution of haloes derives from Poisson sampling, the measured spectra (which is a sum of the true spectra and shot noise terms) can be defined as (e.g. Oddo et al. 2020, and references therein),: where h is the halo number density (e.g. Matarrese et al. 1997). Note that the assumption of a Poisson distribution for the bispectrum shot noise term is accurate to first-order only, since the clustering of haloes is not strictly Poissonian. We discuss in detail the role of Poisson and non-Poisson shot noise in Section 3.4.

Halo power spectrum and bispectrum measurements
In order to measure halo power spectra and bispectra, we modify the publicly-available 4 code of Foreman et al. (2020). The original code measures the matter power spectrum and matter bispectrum from the BAHAMAS, IllustrisTNG 5 , Illustris 6 and EAGLE 7 simulations.
is built upon the (Hand et al. 2018) simulation analysis package, and uses the standard FFT-based bispectrum estimator (Scoccimarro 2000;Sefusatti et al. 2016;Tomlinson et al. 2019). By default, the BAHAMAS module of reads in the positions and masses of particles in order to compute the matter power spectrum and bispectrum. We modify it for this study to read in the positions (defined here as the centre of potential) and masses of haloes as opposed to particles. When computing the halo power spectra and bispectra with by default we weight each halo equally. In Appendix A we explore the effects of weighting haloes by their mass and show that such a weighting scheme tends to reduce the impact of baryons on the bispectrum, as higher mass Table 1. Cosmological parameter values for the BAHAMAS simulations used here. For all the models a WMAP 9-based cosmology is adopted. The columns are: (1) model name; (2) summed mass of the 3 active neutrino species (we adopt a normal hierarchy for the individual masses); (3) subgrid AGN heating temperature; (4) Hubble's constant; (5) present-day baryon density; (6) present-day dark matter density; (7) present-day neutrino density, computed as Ω = /(93.14 eV ℎ 2 ); (8) spectral index of the initial power spectrum; (9) amplitude of the initial matter power spectrum at a CAMB pivot of 2 × 10 −3 Mpc −1 ; (10) present-day (linearly-evolved) amplitude of the matter power spectrum on a scale of 8 ℎ −1 Mpc (note that we use rather than 8 to compute the power spectrum used for the initial conditions, thus the ICs are 'CMB normalised'). In addition to the cosmological parameters, we also list the following simulation parameters: (11) dark matter particle mass; (12) initial baryon particle mass. (1) ( haloes, which are less affected by ejective AGN feedback, contribute more to the computed spectra. In this work, we analyse all self-gravitating substructures, often referred to as 'subhaloes', regardless of whether they are central or satellite subhaloes. Such a selection is closer to a galaxy-based selection than selecting central subhaloes only since all subhaloes of sufficient mass are expected to host galaxies. Without risk of confusion, going forward we will collectively refer to such selfgravitating substructures as just 'haloes'. For BAHAMAS, a standard friends-of-friends (FoF) algorithm (Davis & Peebles 1983) was run to identify FoF groups and the SUBFIND (Springel et al. 2001;Dolag et al. 2009) algorithm was then used to identify all self-gravitating substructures (haloes).
Note that, because the simulations use the same random phases in the initial conditions, it is possible to identify the same haloes in both the hydro run and its corresponding collisionless run. This allows us to isolate the impact of baryonic and cosmological effects in a straightforward way. We refer to the case where we analyse a common set of haloes between multiple simulations as a 'matched' sample. An 'unmatched' sample refers to the case where we simply select haloes based on some criteria (e.g., total mass) from each simulation without regard for whether they correspond to the same haloes 8 . Observationally, it is often the case that haloes/galaxies are selected by, for example, their total masses, thus it is also useful to analyse the results for the unmatched case.
To construct matched samples of haloes, we use the unique dark matter particle IDs, using the 50 most bound particles in each halo of a reference simulation (e.g., the collisionless run). Whichever halo in the target simulation (e.g., the hydro simulation) contains the largest fraction of these dark matter particles from the reference simulation is considered that haloes match. Matches are performed bĳectively, however, and only haloes that are matched both backwards and forwards are retained in the final matched sample.
The size of the simulation box and its mass resolution limits the range of scales over which the power spectra and bispectra can be measured. f = 2 / box = 0.0157ℎ Mpc −1 corresponds to the so-called the fundamental -mode of the BAHAMAS simulations, and is the longest wavelength mode that can be measured given the finite box size. In terms of the binning strategy between these limits, the 8 In general they will not be the same, since feedback and cosmological effects alter the masses of haloes. choice of bin width, Δ , represents a trade-off between the accuracy of results and the computational cost. Note that the computational cost comes in two forms, which is the evaluation of the spectra from a grid (with finer bins requiring more evaluations) and the construction of the grid itself, which can be the most costly aspect depending on the required resolution. It has been shown before that the bispectrum is sensitive to the size of Δ in various ways (e.g. Yankelevich & Porciani 2019;Oddo et al. 2020). In many theoretical works, a binning of Δ = f is used. However, such narrow bins may not be the best choice for measuring the bispectrum from simulations, as the computational cost scales inversely with the width of the k bins used. In addition to this, narrow -bins also leads to additional numerical noise. The common solution to this problem is to have narrow (wide) -bin for small (large) . For example, Foreman et al. (2020) use Δ = f for < 40 f and Δ = 6 f for 40 f < < max , where max is the maximum wavevector for which the analysis is done. We adopt max = 3ℎ Mpc −1 throughout.
Fortunately, the computational resources to measure the halo bispectrum are lower than for the matter bispectrum in the example above, since the construction of the grid is considerably cheaper for the former. Therefore, when computing the halo bispectrum we use the same bin size of Δ = 2 f over the entire range of -scales. However, note that when plotting power spectrum and bispectrum ratios in the figures described later, we apply some smoothing to the curves to reduce the impact of noise. Specifically, we rebin the spectra by a factor 3 and apply a Savitzky-Golay filter.

Choice of triangular configuration
A potentially important question is what triangular configuration 1 , 2 , 3 should be adopted when computing the bispectrum. Results in the literature are often presented for the equilateral case only. However, there is no particularly strong motivation behind this choice, except perhaps for convenience. We present here a brief analysis of the dependence of (the baryons effects on) the halo bispectrum on the choice of triangular configuration.
In Fig. 1 we show the dependence of the effects of baryonic physics on the halo bispectrum as a function of triangular configuration (note that we discuss the physical origin of these effects in Section 4). We do this by plotting the ratio of the bispectrum measured from the fiducial hydro simulation with respect to the corresponding collisionless (DMO) simulation for the case with massless neutrinos. For this test, we calculate the bispectrum for all haloes with masses above 10 11.5 ℎ −1 M . This choice, which roughly corresponds to a stellar mass of 10 10 M , is dictated primarily by the mass resolution of the simulations. For the unmatched case, this mass cut is applied to both the DMO and hydro simulations separately, and thus somewhat different populations of haloes are selected in the two runs. For the matched case, however, we apply the mass cut to the DMO simulation and then select the corresponding matched haloes from the hydro simulation. For the purposes of this comparison we fix 1 to 1.54 ℎ Mpc −1 , since (i) there is an equilateral triangle with such side length; (ii) there is a sufficient number of independent triangles such that the bispectrum is well measured; and (iii) it falls approximately half way between the full range of scales that we consider in the rest of the paper (0.0157 ℎ Mpc −1 to 3 ℎ Mpc −1 ). In Appendix B we explore how the choice of triangular configuration impacts the results for different choices of 1 .
The dots in Fig. 1 (and in Figs. B1-B2) appear in an irregular way because the bispectrum measurements from the simulations do not contain all the possible triangular configurations for the chosen values of 1 . We adopt the same order of wave-vector length as in which is 1 ≥ 2 ≥ 3 . For theoretical calculations all the triangles which follow triangles' condition on the sides, in our case, 1 ≤ 2 + 3 would exist and will fill all the 2 / 1 − 3 / 1 space. However, not all the configurations exist inside the simulation and therefore there are some gaps in the distribution.
The value of the hydro / DMO in Fig. 1 is colour-coded using two colour bars: one for the unmatched haloes (top row) and another for the matched haloes (bottom row) cases. The arrows pointing to the top right corner of each subplot indicate the ratio for equilateral triangle configuration. The right (left) column shows the results for the case where the shot noise has (not) been subtracted.
We can see from Fig. 1 that the impact of the choice of configuration appears to depend on the selection of haloes adopted (unmatched vs. matched) and on whether the bispectrum has been shot noise subtracted. For the shot noise-subtracted case with un-matched haloes (arguably the most useful comparison, since the spectra should be shot noise subtracted and the unmatched case is closer to an observational selection), the ratio of the bispectrum ranges from 1.2 to 1.6, with the equilateral case yielding close to the maximal ratio. The fact that the ratio depends on the choice of triangular configuration suggests that it may be possible to extract additional information about the impact of baryons by computing the bispectra for a range of different configurations. We plan to explore this possibility further in future work and proceed with the standard equilateral configuration for the remainder of this study, noting that there is some sensitivity to the results on the choice of configuration.

Shot noise and resolution considerations
Here we present an analysis of the impact of shot noise on our results and its dependence on the adopted triangular configuration. We also discuss the possible impact of non-Poisson shot noise as well as possible limitations due to finite mass and force resolution as well as finite volume effects. In short, we conclude that resolution and box size limitations are unlikely to be important for our results and conclusions, while the contribution of shot noise to the halo power spectra and bispectra becomes important (i.e., comparable to the true clustering signal) on scales of ≈ 1 ℎ Mpc −1 . Thus, a degree of caution is warranted when examining the absolute power beyond this scale, although we expect the relative trends in ( ) and ( 1 , 2 , 3 ) (e.g., as we increase feedback or the summed neutrino mass) to be robust. Readers who are mainly interested in the impact of baryon physics and neutrino free-streaming on the halo power spectrum and bispectrum may wish to skip ahead to Section 4.
We first consider the impact of Poisson shot noise. In Fig. 2 we show the absolute power spectrum and bispectrum (both DMO and hydro) with and without shot noise subtraction, which illustrates clearly the relative contribution of shot noise as a function of scale. In the case of the bispectrum, we show the two Poisson terms separately (see equation 5 for details). The plot also shows individual Two vertical lines are added to Fig. 2 to indicate the 50% shot noise contribution for both the power spectrum and bispectrum.
(The values at which Poisson shot noise contributes 1% and 10% are also presented in Table 2.) From this figure we conclude that the contribution from Poisson shot-noise becomes comparable to that of the true clustering signal at scales of ≈ 1 ℎ Mpc −1 . Thus, a degree of caution is warranted when examining the absolute power beyond this scale, although we expect the relative trends (e.g., as we increase feedback or the summed neutrino mass) to be robust.
Note that the shot noise contribution to the bispectrum depends on the adopted triangular configuration. Fig. 3 demonstrates that the shot noise ratio is smaller for the squeezed, elongated and folder triangles, and is generally larger for isosceles and equilateral configurations. The largest shot noise case corresponds to the scalene configuration.
In our fiducial analysis in Section 4, we consider the Poissonian shot noise contribution only. However, on small scales we expect the pure Poissonian assumption to break down since haloes can be strongly clustered (e.g., satellites in a galaxy group). Previous work has shown that the non-Poisson contribution to the shot noise is generally subdominant to that of the Poissonian component. For example, using the halo model, Ginzburg et al. (2017) have calculated that the non-Poisson shot noise is approximately half that of the Poisson shot noise contribution at largest values (smallest scales) considered, a result which they confirmed using N-body simulations. At smaller values (large scales), the contribution from non-Poisson shot noise is generally negligible. Thus, we expect that, at worse, our estimate of the shot noise is biased low by 50% and, therefore, our estimate of the true intrinsic bispectrum signal small scales ( ≈ 1 ℎ Mpc −1 ) may be biased high up to 25%. However, we do not expect our relative findings (i.e., ratios of bispectra) to be as adversely affected, since the absolute spectra will all be similarly affected.
We now turn our attention to resolution considerations. Note that the force resolution limit (taken to be the Plummer equivalent gravitational softening of 4 kpc ℎ −1 ) of the simulations is ≈ 1507 ℎ Mpc −1 , which is well beyond the scales we analyse. Thus, limitations due to force resolution are unimportant for our analyses. However, the force resolution is not the only resolution limitation when considering halo clustering. In particular, there are at least three relevant effects: halo exclusion, finite simulation volume, and the lower mass limit of the halo finder catalog. Halo exclusion implies there will be no clustering below the halo radius because the haloes will be touching, which will also impact the noise on small scales (see Baldauf et al. 2013 for more details). In our analyses, however, we examine the clustering of subhaloes rather than FOF groups. Subhaloes identified with SUBFIND can be in close proximity to each other and can even overlap spatially (particles are assigned to a single subhalo based on an energy unbinding criterion). Thus, we do not anticipate halo exclusion effects to be relevant for our analyses. With regards to simulation volume effects, the finite volume can result in biased halo number densities for very large haloes, since they are rare. For the majority of this study, however, we cut off haloes at mass 10 12.5 ℎ −1 M which is a safe option even in moderate size boxes (see van . Even without this cut, our results are likely to be insensitive to the rarest and most massive haloes (which, in BAHAMAS, are ∼ 10 15 M ), since there are orders of magnitude more haloes at lower masses which dominate the signal (see Fig. 5 below). Finally, with regards to the lower mass limit, we selected a lower cut off of 10 11.5 ℎ −1 M corresponding to approximately 100 particles. These are therefore relatively well-resolved haloes, in terms of being able to robustly estimate their masses and positions. We have also verified that the BAHAMAS halo mass function follows that of other results in the literature (e.g., Tinker et al. 2010) down to this limit.

RESULTS
In this section, we present our main results on the impact of baryon physics and neutrino free-streaming on the halo bispectrum. We first examine the impact of baryons for the massless neutrino case in Section 4.1 and then we explore the effects of neutrinos and their degeneracy with baryonic effects in Section 4.2.

Impact of baryon physics
In Fig. 4 we plot the ratio between the fiducial hydro and the DMO BAHAMAS simulations (both with massless neutrinos) of the halo power spectrum and bispectrum for equilateral triangles. We consider both the raw spectra computed by which includes a shot noise contribution and shot noise-subtracted spectra, where we have assumed the shot noise follows a Poisson distribution as given in equations (4) and (5). We note that when evaluating the shot noise contribution to the bispectrum in equation (5), the power spectra should be shot noise-subtracted (according to eqn. 4). Since the power spectrum differs between the hydro and DMO cases, the bispectrum shot noise term will also, in general, be different for the hydro and DMO cases, even for a matched set of haloes. Furthermore, it is clear that the bispectrum shot noise is not a constant (i.e., it depends on scale), unlike that for the power spectrum. As can be seen in Fig. 4, shot noise does not greatly affect large scales, but starting from > ∼ 0.1ℎ Mpc −1 it begins to have an increasingly important effect. Therefore, we will subtract the shot noise by default for the rest of the paper and we focus on the shot noise-subtracted results below.
Previous studies have shown that including galaxy formation physics in hydro simulations leads to a characteristic suppression in the non-linear matter power spectrum, which is due primarily to the ejection of baryons by AGN feedback (e.g. van Mummery et al. 2017;Chisari et al. 2019;Foreman et al. 2020;Stafford et al. 2020;van Daalen et al. 2020). This effect was also shown to produce a suppression in the matter bispectrum by Foreman et al. (2020) (see also Aricò et al. 2021). However, we show here that, at least at face value, haloes appear to show the opposite effect. In particular, for the unmatched case, the 'hydro' power spectrum and bispectrum are always greater than that of the corresponding 'DMO' simulation, implying that haloes selected to be above a certain mass appear to be more clustered in a hydro simulation than in the corresponding collisionless one. Note that the apparent slight suppression in the bispectrum at large scales (small ) is due to sampling noise from finite box size effects.
In the right panel of Fig. 4 we consider the ratio of the hydro to the DMO case for a given ('matched') set of haloes common to the two simulations, rather than selecting haloes above a given mass limit (as in the left panel). Here we see that on relatively large scales ( < ∼ 0.5 ℎ Mpc −1 ) the ratio now goes to unity, indicating no appreciable difference in the clustering of haloes in the two simulations on large scales. The same conclusion was made by van  when examining the 2-point correlation functions of a matched set of haloes in the OWLS simulations.
As we are examining a common set of haloes in this comparison, it implies that the differences for the unmatched case (left panel) are primarily driven by a difference in the haloes that are selected when using a fixed mass cut in that case. In particular, AGN feedback removes baryons from haloes in the hydro case, reducing their overall masses (e.g., Cui et al. 2014;Cusworth et al. 2014;Velliscig et al. 2014;Schaller et al. 2015). Consequently, some haloes that would (in the absence of feedback) have been just above the mass cut are now below the threshold for selection and are therefore not selected. The net result is that introducing baryon physics (feedback) results in the selection of intrinsically (in the absence of feedback) more massive systems. Since more massive haloes are more strongly biased in their clustering with respect to the mean matter density (e.g., Tinker et al. 2010), the halo clustering is enhanced on large scales in the hydro simulation with respect to the DMO simulation when selecting haloes above a given mass limit. For a common set of haloes (as in the right panel), however, there is no appreciable difference in the clustering on large scales.
On small scales ( > ∼ 1 ℎ Mpc −1 ), however, the clustering is  altered by hydrodynamics even for a common set of haloes. In particular, we find that the power spectrum is slightly suppressed (at the 5% level up to ≈ 3 ℎ Mpc −1 ), whereas the bispectrum is enhanced (at the 10% level up to ≈ 3 ℎ Mpc −1 ). On these scales, the clustering is almost certainly dominated by the so-called '1-halo' term, corresponding to satellites of massive haloes (groups and clusters). Thus, deviations from unity are likely signifying differences in how the satellites are spatially distributed in massive haloes in the hydro case with respect to the DMO case. Physically, such differences can potentially arise due to changes in the efficiency of tidal disruption, with gas ejection likely making satellites more susceptible to stripping at large (satellite-centric) scales, while star formation and adiabatic contraction would likely make disruption of the central regions of satellites more difficult to strip (with respect to the DMO case). However, we cannot easily rule out that differences in the performance of the substructure finder (SUBFIND) between the DMO and hydro cases is also potentially playing a role here. To test this, higher resolution simulations and/or alternative substructure finders (e.g., phase space-based finders) could be employed, but we leave this for future work.

Halo mass dependence
In the analysis presented above we explored the impact of baryon physics on the halo bispectrum at = 0 when selecting all haloes above 10 11.5 ℎ −1 M , corresponding roughly to a stellar mass selection of * > ∼ 10 10 M . Here we attempt to examine the halo mass dependence of the results, noting that the effects of stellar and AGN feedback on the baryon fractions are expected to be halo mass dependent (e.g., Schaller et al. 2015). We, therefore, expect the effects of feedback on the halo bispectrum to be sensitive to the halo mass range examined. However, a limitation of the current study is that very large samples of haloes are required to accurately measure the bispectrum, thus breaking the sample up into multiple Figure 6. Upper panels: The ratio of the power spectrum (dot-dashed lines) and the bispectrum (solid lines) measured from the hydro to DMO BAHAMAS simulations with shot noise subtraction for equilateral triangles for haloes masses in range 10 11.5 ℎ −1 M ≤ < 10 12.5 ℎ −1 M . Different line colours represent variations the efficiency of AGN feedback. Lower panels: The ratio ( ) of the power spectrum and the bispectrum ratios between low/high AGN feedback models with respect to the fiducial (tuned AGN) model. mass bins will yield less statistically robust results. Nevertheless, we proceed here by examining the bispectrum (and power spectrum) in two adjacent mass bins of 10 11.5 ℎ −1 M ≤ < 10 12.5 ℎ −1 M and 10 12.5 ℎ −1 M ≤ < 10 13.5 ℎ −1 M and we compare these with our previous selection of ≥ 10 11.5 ℎ −1 M (Fig. 5). Note that above the mass limit 10 11.5 ℎ −1 M there are 721,459 haloes (or 100%) in the DMO simulation and that the lower mass bin (10 11.5 ℎ −1 M ≤ < 10 12.5 ℎ −1 M ) contains 638,832 of these haloes (or 88.6%) while the higher mass bin contains only 76,016 haloes (or 10.5%).
In Fig. 5 we show the ratio of the halo power spectra and bispectra (hydro to DMO) for the different halo mass selections for both the unmatched (left) and matched (right) cases. The lower mass bin roughly follows the same trend as our total halo sample for the unmatched haloes, though the effects of baryons on both the power spectrum and bispectrum are slightly more pronounced in the low mass bin compared to the total halo selection. The results for the higher mass bin (10 12.5 ℎ −1 M ≤ < 10 13.5 ℎ −1 M ) are, unfortunately, too noisy to make any strong quantitative conclusions. Larger simulation volumes are required to examine the halo bispectrum of high-mass objects. Qualitatively, though, the results are consistent with a reduced impact due to baryon physics in comparison to the lower mass bin. One can also infer this from the comparison of the total selection and lower mass bin comparison, since the latter is more strongly affected by baryons and the only difference in the selections is that the total case also contains high-mass haloes.
For the case of matched haloes (i.e., where haloes are selected according to some mass criterion in the DMO simulations and bĳectively matched to the hydro simulation), shown in the right panel of Fig. 5, we find much the same behaviour as in the right panel of Fig. 4. Specifically, for a given set of haloes, baryon physics does not significantly alter the halo power spectrum or bispectrum on large scales ( < ∼ 0.5 ℎ Mpc −1 ) independent of halo mass. On small scales (the 1-halo regime), however, the spectra are altered and the effects are largest for the lower mass bin (in agreement with the left panel) and we again find that the bispectrum is enhanced due to baryon effects whereas the power spectrum is suppressed (similar to the matter power spectrum).

Feedback variations
In addition to the fiducial calibrated BAHAMAS model, McCarthy et al. (2018) introduced two additional simulations which varied the strength of the AGN feedback: 'low AGN' and 'high AGN'. Specifically, they varied the subgrid AGN heating temperature parameter by ±0.2 dex around the fiducial log 10 [Δ heat ] = 7.8, which approximately brackets the upper and lower envelopes of the observed gas fractions of galaxy groups and clusters, as inferred from resolved X-ray observations (see fig. 3 of McCarthy et al. 2018). Note that Δ heat is the most important parameter in the model for determining the gas fractions of groups/clusters. Higher values of Δ heat correspond to more energetic (and more bursty) AGN feedback episodes and vice-versa.
We show in Fig. 6 the effect of varying Δ heat on the halo power spectrum and bispectrum for a halo mass selection of 10 11.5 ℎ −1 M ≤ < 10 12.5 ℎ −1 M , once again through the ratio of these statistics measured from the hydro simulations computed with respect to the DMO case. For the unmatched case, the power spectrum and bispectrum ratios for all three AGN models show an enhancement due to feedback. Interestingly, it is the 'low AGN' model which shows the largest effect particularly on small scales, whereas for the matter power spectrum and bispectrum the opposite is true (e.g., fig. 13 of Foreman et al. 2020, see also van Daalen et al. 2020). In the case of the matter clustering, the interpretation is straightforward: more energetic feedback leads to increased baryon ejection which reduces the overall matter clustering. For the halo clustering, particularly on small scales which are dominated by the 1-halo term (satellites of groups/clusters), it is plausible that particularly energetic feedback leads to satellites being more susceptible to environmental processing (tidal heating/stripping), reducing the power on small scales compared to a run with less energetic feedback. It is interesting that in the matched halo case (right panel of Fig. 6) the fiducial and low AGN runs show enhanced halo bispectra (and power spectra, for low AGN) relative to the DMO case, whereas the high AGN shows suppressed halo clustering. This suggests there is a relatively fine line between baryon effects enhancing the clustering on small scales relative to DMO (e.g., due to star formation and adiabatic contraction of the satellites subhaloes) and enhanced tidal stripping due to baryon ejection. The ratio of the power spectrum (dot-dashed lines) and the bispectrum (solid lines) measured from the hydro to DMO BAHAMAS simulations with shot noise subtraction for equilateral triangles for haloes masses in range 10 11.5 ℎ −1 M ≤ < 10 12.5 ℎ −1 M . Different line colours represent different summed neutrino masses . Note that the summed neutrino mass is varied simultaneously here for the hydro and DMO runs. Lower panels: As in lower panels of Fig. 6 but ratio ( ) is of the massive neutrino cases with respect to the fiducial massless neutrino case.
The upshot of the sensitivity of the halo clustering due to baryon effects is that the simulations probably cannot be used to robustly predict whether there should be an enhancement or a suppression relative to DMO simulations on small scales (for a common set of haloes), as the answer appears to depend sensitively on the balance between feedback and environmental physics. This is in contrast to the findings of Foreman et al. (2020) (see, e.g., their fig. 13), in which the matter bispectra of the three BAHAMAS hydro simulations consistently show a suppression relative to the DMO case, whereas our measurements of the halo bispectrum show either a suppression or an enhancement depending on the AGN feedback strength. This potentially complicates the interpretation of observational measurements, since a probe such as the cosmic shear bispectrum would observe a suppression while the galaxy 3-point function might be suppressed or enhanced. A joint analysis of the matter and halo bispectra would therefore be useful for simultaneously constraining the effects of feedback and environmental physics on small scales. In this regard, it is noteworthy that the effects are typically much stronger in the bispectrum than in the power spectrum (consistent with the findings of Foreman et al. 2020), suggesting that future measurements of the small-scale bispectrum are a promising tool for studying these effects.

Impact of neutrino free-streaming
Of the various possible extensions to the standard ΛCDM model, massive neutrinos are perhaps the most well motivated, as the results of atmospheric and solar oscillation experiments imply that the three active species of neutrinos have a minimum summed mass, Σ , of at least 0.06 eV (0.1 eV) when adopting a normal (inverted) hierarchy (see Lesgourgues & Pastor 2006 for a review). At present, cosmological studies place the strongest upper limits on the summed neutrino mass. Galaxy clustering studies, in particular, have a rich history of placing important constraints on the summed neutrino mass (e.g., Tegmark et al. 2004;Cole et al. 2005;Sánchez et al. 2006;Beutler et al. 2014;Abbott et al. 2018), with its sensitivity coming from the fact that neutrinos suppress the growth of structure (i.e., reduce the clustering amplitude relative to a massless neutrino cosmology) on scales smaller than the associated free-streaming scale. In the present study, we use the massive neutrinos extension of BAHAMAS (Mummery et al. 2017;McCarthy et al. 2018) to explore the joint effects of baryon physics and neutrinos, and their possible degeneracy, on the halo power spectrum and bispectrum. In addition to the massless neutrino cases explored above, BAHAMAS has four variations of summed neutrino mass, with Σ = 0.06, 0.12, 0.24 or 0.48 eV ( Table 1). As mentioned in Section 2, neutrinos are implemented in the simulations using the semi-linear method of Ali-Haïmoud & Bird (2013), which is sometimes referred to as a 'linear response' method. For further details of the implementation within BAHAMAS, we refer the reader to McCarthy et al. (2018).
In Fig. 7 we show the ratio of the power spectrum and the bispectrum (hydro to DMO) where the summed neutrino mass is systematically varied for both the hydro and DMO runs (simultaneously). In essence, this comparison explores whether the relative impact of baryon physics on the halo clustering is dependent on the summed neutrino mass. If the ratio varies significantly between the simulations, this implies that the impact of baryon physics is dependent on the summed neutrino mass, whereas if the ratio does not vary between the simulations then the effects of baryon physics are independent of (i.e., separable from) those of massive neutrinos. As shown in Fig. 7 (focusing on the bottom subpanels, in particular, which show the ratio with respect to the massless neutrino case) there is no evidence for a strong systematic dependence on the summed neutrino mass, with the results being very similar to the massless case. The noise in the halo bispectrum measurements means we can only place an upper limit on the potential effect but, for the matched case, the impact of baryons on the halo bispectrum agrees amongst all of the neutrino simulations to typically a few per cent accuracy (and to higher accuracy for the halo power spectrum). Thus, the impact of baryons and neutrinos appear to be separable effects for the halo power spectrum and bispectrum, as previously found for the matter clustering in Mummery et al. (2017) and van Daalen et al. (2020).
In Fig. 8 we explore the combined effects of baryon physics and neutrino free-streaming on the halo clustering, by taking the ratio of the hydro simulations of varying neutrino mass to the DMO run with massless neutrinos. Thus, both feedback and neutrinos will impact this metric. In fact, there are four effects present when considering the left panel of Fig. 8 (unmatched haloes); as both feedback and free-streaming affect the halo mass (i.e., alter the halo selection) and both feedback and free-streaming affect the clustering (although we have shown that the clustering is only affected on small scales by baryon physics after we control for the halo mass alteration). Neutrinos alter the halo mass because they represent a form of dark matter that is effectively unable to collapse, whereas all of the dark matter can collapse in the massless neutrino case. Mummery et al. (2017) have shown that it is the largest haloes (massive clusters) that are most strongly affected in terms of halo mass.
In the left panel of Fig. 8 we see that the halo power spectra and bispectra are typically enhanced with respect to the DMO with massless neutrinos case. Increasing the summed neutrino mass tends to reduce this enhancement, implying that the free-streaming reduction in clustering amplitude is a more significant effect than the reduction in halo mass that is due to neutrinos. However, for this range of neutrino masses, the change in halo mass due to baryon ejection is sufficient to overcome the impact of neutrinos, which is why there is an overall enhancement in the clustering rather than a suppression. The one exception to this is the Σ = 0.48 eV case, which shows a slight suppression on quasi-linear scales. Here the neutrino suppression is sufficiently strong to overcome the enhancement due to baryon ejection.
In the right panel of Fig. 8 we do the same comparison but now for a common set of haloes, selected from the DMO massless neutrino run and bĳectively matched to each of the hydro simulations with varying neutrino masses. As explained before, by examining a common set of haloes we are removing the impact of changes in halo mass, due to both feedback and neutrino free-streaming, on the halo selection. Here we see a strong suppression in the clustering which is present on virtually all scales plotted, apart from the very smallest scales where the enhancement due to baryon physics begins to dominate. In fact, the suppression will be present out to approximately the free-streaming scale, which corresponds to fs ≈ 0.05 ℎ Mpc −1 for Σ = 0.06 eV (see eqn. 9 of Ali-Haïmoud & Bird 2013). In common with the impact of feedback, we find that the bispectrum is affected much more significantly than is the power spectrum, although we note that estimates of the bispectrum are more affected by noise (sample variance) than is the power spectrum. In principle, though, upcoming precision measurements of the bispectrum may be able to place strong constraints on both cosmology and mod-els of galaxy formation simultaneously, by using small scales to probe feedback and environmental processing and large scales to help constrain the neutrino mass. Certainly the combination of the power spectrum and bispectrum will allow for much more powerful constraints than may be obtained from the power spectrum alone.
Finally, in the above discussion we have neglected the potential degeneracies between the summed neutrino mass, 8 , and baryons. It is known, for example, that using clustering data alone which is restricted to quasi-linear and non-linear scales (e.g., as in current cosmic shear surveys) that Σ and 8 can be strongly degenerate. Because the BAHAMAS neutrino simulations are 'CMB normalised', varying Σ leads to a variation in 8 as well. To explore the degeneracy in detail would require additional simulations which hold 8 fixed as the neutrino mass is varied (or vice-versa), which is beyond the scope of the present study. We note, however, that with the increasing quality of weak lensing and spectroscopic galaxy clustering surveys, one can exploit the fact that neutrino freestreaming, variations in 8 , and the impact of baryons on the matter and halo clustering all have different scale and redshift dependencies (see, e.g., Mummery et al. 2017 who compared the scale and redshift dependencies of baryons and massive neutrinos on various clustering statistics). We therefore expect upcoming clustering studies to be able to place much stronger constraints on these effects than is currently possible.

SUMMARY AND CONCLUSIONS
The primary goal of this paper was to demonstrate the power of the halo bispectrum as a tool for testing cosmological models and galaxy formation physics (particularly AGN feedback). To this end, we measured the halo power spectrum and the halo bispectrum from the BAHAMAS hydro and dark matter simulation suites. Aside from the fiducial calibrated BAHAMAS hydro simulation (which was calibrated to reproduce the galaxy stellar mass function and the gas fractions of galaxy groups/clusters), the suite contains runs which vary the strength of the AGN feedback and the summed mass of neutrinos, allowing us to explore the sensitivity of the halo clustering to these effects. We modified the publicly available matter clustering code (Foreman et al. 2020) to measure the halo clustering in BAHAMAS.
The main results of our study may be summarised as follows: • The bispectrum measured from the simulations in general depends on the adopted triangular configuration. It is common place to adopt an equilateral configuration mostly for convenience. We explored the dependence of our main results (the relative impact of baryon physics on the halo bispectrum) on the adopted triangular configuration in Fig. 1 (see also Appendix B), finding that the equilateral choice tends to maximise the impact of baryons. We adopted the equilateral choice for the remainder of the analysis but note that the results are somewhat sensitive to the choice of triangular configuration. In the future, it may be possible to exploit this sensitivity as a further test of baryon and cosmological physics.
• When selecting haloes above a given mass, baryon physics (particularly AGN feedback) tends to enhance the halo power spectrum and bispectrum relative to a collisionless (dark matter-only, DMO) case (see left panel of Fig. 4). This enhancement is driven by the change in halo mass due to AGN feedback. Specifically, haloes just above the selection threshold in the DMO case have their masses reduced due to gas ejection and therefore not selected in the hydro simulation case. Consequently, the mean bias of haloes selected in the hydro case is larger than that of the DMO case, giving rise to an apparent enhancement in the clustering. Comparing the clustering of a common set of haloes in the hydro and DMO cases (see right panel of Fig. 4) demonstrates that baryon physics does not appreciably affect the clustering on large scales of < ∼ 0.5 ℎ Mpc −1 . On small scales, dominated by the 1-halo term (satellites of massive groups and clusters), the power spectrum and bispectrum can be significantly affected (≈ 10% level) and, interestingly, we find for the fiducial calibrated BAHAMAS model that the bispectrum is enhanced with respect to the DMO case, while the power spectrum is suppressed.
• We find that the effects described above are a function of halo mass, with lower masses (∼ 10 12 M ) being more significantly affected by baryon physics (see Fig. 5).
• We explored the sensitivity of the results to variations in galaxy formation physics by using two variations from the BAHAMAS suite ('low AGN', 'high AGN'). Interestingly, we find that reducing the efficiency of AGN feedback results in an increased enhancement in the halo bispectrum on small scales (see right panel of Fig. 6), which is plausibly due to the fact satellites are more resilient against tidal disruption in this case. Increasing the efficiency of AGN feedback can actually change the sign of the effect, such that the halo bispectrum is suppressed with respect to the DMO case. This suggests that there is a very fine balance between star formation/adiabatic contraction (which tend to make haloes less susceptible to tidal heating/stripping) and gas ejection (which make them more susceptible). Small-scale measurements of the halo bispectrum, therefore, offer a means of simultaneously probing feedback and environmental physics.
• Using the massive neutrinos extension of BAHAMAS we explored the joint effects of neutrino free-streaming and baryon physics (and their potential degeneracy) on the halo clustering. We showed that the effects of baryon physics are largely independent of (separable from) the effects of neutrinos (see Fig. 7). For a common set of haloes, we find neutrinos strongly suppress the halo bispectrum on all but the smallest scales probed in the present study (the smallest scales are dominated by baryon physics which can either enhance or suppress the bispectrum, as discussed above). See the right panel of Fig. 8.
Our results have demonstrated that the halo power spectrum and particularly the halo bispectrum are highly sensitive to neutrino effects on the largest scales and baryon physics (AGN feedback and star formation) on the smallest scales. This bodes well for upcoming precision measurements of the bispectrum for surveys like Euclid, LSST, and DESI. To capitalise on these measurements and the implied sensitivity what is required now is a large suite of very large volume cosmological simulations that systematically and simultaneously explore variations in cosmological parameters and galaxy formation physics. This is the aim of ongoing work.

APPENDIX A: HALO MASS-WEIGHTED SPECTRA
In the main analysis, when computing the halo power spectra and bispectra we treat all haloes equally (same weighting). Here we explore an alternative weighting scheme, which is to weight each halo by its mass. It was shown previously by Seljak et al. (2009) that the effects of shot noise can be significantly reduced in the halo power spectrum by adopting such a weighting scheme. In there is an option to include particle weight into analysis when computing the grids. We simply substitute in the halo mass. For the matched case we used DMO halo masses for both DMO and hydro grids.
We demonstrate our findings in Fig. A1. While we see no obvious advantage in terms of the shot noise for the bispectrum, it is interesting to note that the effects of baryon physics are reduced significantly when weighting haloes by their mass instead of equally weighting each halo. This is likely because higher mass haloes are less significantly affected by ejective feedback. Thus, if the aim is to mitigate against the effects of baryons, employing a halo mass weighting scheme goes some way towards achieving this objective.

APPENDIX B: TRIANGULAR CONFIGURATION DEPENDENCE AS A FUNCTION OF SCALE
Here we extend the analysis presented in Section 3.3 for different values of 1 . In Section 3.3 we adopted 1 = 1.54 ℎ Mpc −1 , which is approximately half way between between the minimum and maximum scales we can sample in the simulation box. Here we Figure A1. The same as at Fig. 4 but including cases where the haloes are weighted by their masses when computing the power spectra and bispectra. study the dependence of the bispectrum ratio for different triangular configurations on large scales ( 1 = 1.01 ℎ Mpc −1 ; see Fig. B1) and for small scales ( 1 = 2.92 ℎ Mpc −1 ; see Fig. B2). Again the values of 1 are chosen in a way that the equilateral configuration 1 = 2 = 3 must exist. For small values of 1 , the number of possible triangles in the box is low. We therefore plot in Fig. B1 the first wave vector value where the number of triangles is representative.
Overall, we see similar trends to that presented in Fig. 1, indicating that the dependence on the adopted triangular configuration is similar for different choices of 1 . Figure B1. The same as in Fig. 1 but for a lower value of 1 = 1.01 ℎ Mpc −1 . Figure B2. The same as in Fig. 1 and Fig. B1 but for a higher value of 1 = 2.92 ℎ Mpc −1 .