Modelling large-scale halo bias using the bispectrum
aa r X i v : . [ a s t r o - ph . C O ] D ec Mon. Not. R. Astron. Soc. , 000–000 (0000) Printed 3 October 2018 (MN L A TEX style file v2.2)
Modelling large-scale halo bias using the bispectrum
Jennifer E. Pollack ⋆ , Robert E. Smith , † , & Cristiano Porciani ‡ Argelander Institut f¨ur Astronomie der Universit¨at Bonn, Auf dem H¨ugel 71, D-53121 Bonn, Germany Institute for Theoretical Physics, University of Zurich, Zurich, CH 8037
ABSTRACT
We study the relation between the density distribution of tracers for large-scale struc-ture and the underlying matter distribution – commonly termed bias – in the ΛCDMframework. In particular, we examine the validity of the local model of biasing atquadratic order in the matter density. This model is characterized by parameters b and b . Using an ensemble of N -body simulations, we apply several statistical meth-ods to estimate the parameters. We measure halo and matter fluctuations smoothedon various scales. We find that, whilst the fits are reasonably good, the parametersvary with smoothing scale. We argue that, for real-space measurements, owing to themixing of wavemodes, no smoothing scale can be found for which the parameters areindependent of smoothing. However, this is not the case in Fourier space. We measurehalo and halo-mass power spectra and from these construct estimates of the effectivelarge-scale bias as a guide for b . We measure the configuration dependence of the halobispectra B hhh and reduced bispectra Q hhh for very large-scale k -space triangles. Fromthis data we constrain b and b , taking into account the full bispectrum covariancematrix. Using the lowest-order perturbation theory, we find that for B hhh the best-fitparameters are in reasonable agreement with one another as the triangle scale is var-ied; although, the fits become poor as smaller scales are included. The same is true for Q hhh . The best-fit values were found to depend on the discreteness correction. This ledus to consider halo-mass cross-bispectra. The results from these statistics supportedour earlier findings. We then developed a test to explore whether the inconsistency inthe recovered bias parameters could be attributed to missing higher-order correctionsin the models. We prove that low-order expansions are not sufficiently accurate tomodel the data, even on scales k ∼ . h Mpc − . If robust inferences concerningbias are to be drawn from future galaxy surveys, then accurate models for the fullnonlinear bispectrum and trispectrum will be essential. Key words: cosmology: theory, large-scale structure
The accurate estimation and modelling of higher-orderclustering statistics in current and future galaxy redshiftsurveys has the potential to act as a powerful probefor cosmological physics. The higher-order connected mo-ments, beginning at lowest order with the three-pointcorrelation function and its Fourier analogue, the bis-pectrum, when interpreted within the gravitational insta-bility paradigm, encode important information regardingthe growth of large-scale structure (Matarrese et al. 1997;Scoccimarro et al. 1998). Their measurements also provideinsight into the statistical nature of the primordial fluc- ⋆ E-mail: [email protected] † [email protected] ‡ [email protected] tuations (Fry & Scherrer 1994; Sefusatti & Komatsu 2007;Nishimichi et al. 2010; Baldauf et al. 2011) and the cosmo-logical parameters (Sefusatti et al. 2006). Another attributeof three-point statistics, and the focus of our study, is theircapability to probe the manner in which an observabletracer population of objects, such as galaxies, is related tothe unobservable matter distribution – termed the ‘bias’(Kaiser 1984; Dekel & Rees 1987; Fry & Gaztanaga 1993;Dekel & Lahav 1999; Catelan et al. 2000).If the primordial fluctuations were Gaussian, as appearsto be the case (Komatsu et al. 2010), then the statisticalproperties of the initial fields are fully characterized by thepower spectrum, with all higher-order connected correlatorsvanishing. However, gravitational instability leads to thecoupling of Fourier modes and this generates a hierarchyof non-vanishing connected correlators, each of which hasa precise characteristic mathematical structure. The matter c (cid:13) J. E. Pollack, R. E. Smith & C. Porciani bispectrum is thus an inherently nonlinear quantity, whosesignal depends on closed triangles in Fourier space. In the-ory, this should vanish at early times and on scales largeenough where linear theory is valid. If galaxy bias is localand linear, then the bispectrum of the observable tracersis proportional to the matter bispectrum. If on the otherhand, bias is local and nonlinear then the triangle configu-ration dependence of the signal is modified, and this happensin a very precise and calculable way. Thus the bispectrumcan be used to constrain the bias (Fry & Gaztanaga 1993;Matarrese et al. 1997; Scoccimarro et al. 1999).There is a long and rich history of measurements ofthree-point statistics from galaxy surveys, going all theway back to Peebles & Groth (1975). However, attemptsto constrain the nonlinearity of galaxy bias from galaxyredshift surveys have only been performed over the lastdecade. Feldman et al. (2001) and Scoccimarro et al. (2001)both analyzed the IRAS survey using the bispectrum, andfound a negative quadratic bias; although, due to smallsample size the constraints were rather weak. Verde et al.(2002) analyzed the 2dFGRS survey, also using the bispec-trum approach, and claimed that the flux-limited samplewas an unbiased tracer of the dark matter. A subsequentanalysis of the final 2dFGRS data set by Gazta˜naga et al.(2005), using the 3-point correlation function, contradictedthis: using information from weakly non-linear scales ( R ∼ − h − Mpc) the unbiased case ( b = 1 and b = 0)was excluded at the order of 9 σ . More recently a num-ber of authors have analyzed various data releases of theSDSS (Nishimichi et al. 2007 – DR3, McBride et al. 2011– DR6 and Marin 2010 – DR7). These all claim a non-zero quadratic bias term for most of the samples within thedataset. Obviously, these variations of the results with sur-vey and statistical method require an explanation.Whilst the local bias model can be used to test whetherthe bias is linear or nonlinear, a significant detection of non-zero nonlinear bias does not imply that we have understoodthe bias. In order to believe that these measurements aremeaningful, we need to be sure that the local model is indeedthe correct model for interpreting data. This is currently anopen question. Attempting to shed light on this subject isone of the aims of this paper. Over the past few years, thelocal model of galaxy bias has been scrutinized by a numberof authors (Heavens et al. 1998; Gazta˜naga & Scoccimarro2005; Smith et al. 2007; Guo & Jing 2009; Manera et al.2010; Roth & Porciani 2011; Manera & Gazta˜naga 2011).However, no firm conclusions have yet been reached.In this paper, we use a large ensemble of 40 mid-resolution, large volume, pure dark matter N -body simu-lations, to test the validity of the local bias model. In thisstudy we compare a selection of different methods for deter-mining the bias. We first present a point-wise comparisonof the halo and matter density fields smoothed on certainscales. We also utilize the power spectrum to estimate alarge-scale effective bias parameter. Then we expend mostof our efforts on using the bispectrum and reduced bispec-trum approach for constraining the bias. Besides the auto-bispectra, we also present, for the first time, measurementsof the halo-matter cross-bispectra: B hhm and B hmm , and Q hhm and Q hmm . The value of these new statistics becomesapparent when correcting for shot-noise effects. Finally, weperform a numerical test that allows us to sharply illumi- nate the importance of terms in the theory that are beyondthe tree-level expansions typically used.The paper is divided up as follows: in § §
3, we provide details of the numerical simulations usedin this work. In §
4, we present the results for the bias pa-rameters from the various commonly used simple estimators.In § § § § § In the fluid approximation, the gravitational collapse of col-lisionless cold matter structures in the expanding Universe,can be fully characterized by specifying the evolution of thedensity δρ ( x ) and the peculiar velocity δ v ( x ) perturbations(Bernardeau et al. 2002). Focusing primarily on the densityfield, we work with models of the matter density contrast: δ ( x , t ) ≡ ρ ( x , t ) − ρ ( t ) ρ ( t ) , (1)where ¯ ρ ( t ) is the mean matter density of the Universe. InFourier space, we define its corresponding Fourier represen-tation, δ ( k ), accordingly, as δ ( x ) = Z d k (2 π ) δ ( k ) e − i k · x ⇔ δ ( k ) = Z d x δ ( x ) e i k · x . (2)It can be shown that the nonlinear equations of motionfor δ ( k ) can be solved exactly by perturbative expansionsof the type (Juszkiewicz 1981; Vishniac 1983; Goroff et al.1986): δ ( k ) = ∞ X n =1 a n ( t ) δ n ( k ) , (3)and, where δ n ( k ) is given by, δ n ( k ) = Z d q (2 π ) . . . Z d q n (2 π ) (2 π ) δ D ( k − q − · · · − q n ) × F n ( q , . . . , q n ) δ ( q ) · · · δ ( q n ) . (4)The density kernel F n is the dimensionless, homoge-neous, mode coupling function that couples together theamplitudes and phases of n initial Fourier wavemodes { δ ( q ) , . . . , δ ( q n ) } . As was shown by (Goroff et al. 1986;Makino et al. 1992; Jain & Bertschinger 1994), the n thorder kernel may be constructed recursively from thelower-order solutions. Linear theory is thus represented by F ( q ) = 1, and the first nonlinear correction by F , where F ( k , k ) = 57 + k · k k k (cid:18) k k + k k (cid:19) + 27 ( k · k ) k k . (5)The above approach defines the standard perturbation the-ory (hereafter SPT). Before moving on, we note that theabove statements are only exactly true for the Einstein-de Sitter model. However, it has been shown that the F kernel is almost independent of cosmology (Fry 1994;Bouchet et al. 1995; Hivon et al. 1995). We therefore adoptEq. (5) when dealing with the density at second order. c (cid:13) , 000–000 odelling halo bias using the bispectrum Owing to the stochastic nature of the density field, we arenot interested in reproducing a specific density field per se ,but instead in characterizing its statistical properties. In thiswork we focus on 2- and 3-point correlation functions inFourier space. These we may write as: h δ ( k ) δ ( k ) i ≡ (2 π ) δ D ( k ) P mm ( k ) ; (6) h δ ( k ) δ ( k ) δ ( k ) i ≡ (2 π ) δ D ( k ) B mmm ( k , k , k ) , (7)where P mm ( k ) and B mmm ( k , k , k ) constitute definitionsof the power and bispectrum. For the Dirac delta functionswe used the short-hand notation δ D ( k ...n ) ≡ δ D ( k + . . . + k n ) and these guarantee that P and B are translationally in-variant. This is an important property for estimation, sinceit means that we should consider only closed pairs and tri-angles in Fourier space: P k i = 0.The perturbative expansion of the density field de-scribed in the previous section implies that P mm and B mmm may also be described in a perturbative fashion. Hence, h δ ( k ) δ ( k ) i = h [ δ ( k ) + δ ( k ) + . . . ] × [ δ ( k ) + δ ( k ) + . . . ] i ; (8) h δ ( k ) . . . δ ( k ) i = h [ δ ( k ) + δ ( k ) + . . . ] . . . × [ δ ( k ) + δ ( k ) + . . . ] i . (9)Since we are assuming that the initial Fourier modes areGaussianly distributed, i.e. the phase of each initial modeis uniformly random φ ∈ [0 , π ], modes must cancel inpairs. Hence, Wick’s theorem applies, and so odd prod-ucts of initial Fourier modes must vanish: h δ ( k ) δ ( k ) i = h δ ( k ) δ ( k ) δ ( k ) i = 0 . This leads us to write the per-turbative expansions for P mm and B mmm as: P mm ( k ) = P (0)mm ( k ) + P (1)mm ( k ) + . . . ; (10) B mmm ( k , k ) = B (0)mmm ( k , k ) + B (1)mmm ( k , k ) + . . . (11)We shall refer to the lowest order terms in the expansions as‘tree-level’ terms. For P , P (0) is simply the linear spectrum,while for B the tree-level term can be written: B (0)mmm ( k , k ) = 2 P (0)mm ( k ) P (0)mm ( k ) F ( k , k ) + 2 cyc . (12)In this work we shall mainly be dealing with tree-level quan-tities; we now set: P (0)mm = P mm and B (0)mmm = B mmm , unlessotherwise stated.Another statistical quantity commonly used toexplore galaxy clustering is the reduced bispectrum(Peebles & Groth 1975; Scoccimarro et al. 1998), which canbe defined: Q mmm ( k , k , k ) ≡ B mmm ( k , k , k ) P mm ( k ) P mm ( k ) + 2 cyc . (13)As will be made clear below, the importance of this statisticbecomes apparent when one considers non-Gaussian termsthat are generated by simple quadratic products of Gaussianfields. In this case Q mmm simply scales as a constant. In this study we investigate the relation between the clus-tering of dark matter haloes and total matter. If galaxies are only formed in dark matter haloes, as is the usual assump-tion for all models of galaxy formation (White & Rees 1978),then understanding the clustering of haloes is an essentialcomponent of any theory of galaxy biasing (Smith et al.2007). In the local model of halo biasing, the number den-sity of dark matter haloes of mass scale M , smoothed overa scale R , can be expressed as a function of the local matterdensity, smoothed on the scale R . This function may thenbe Taylor expanded to give (Fry & Gaztanaga 1993; Coles1993; Mo et al. 1997; Smith et al. 2007): δ h ( x | M, R ) = ∞ X j =0 b j ( M ) j ! [ δ ( x | R )] j , (14)where we defined the smoothed halo over-density to be, δ h ( x | M, R ) ≡ [ n h ( x | M, R ) − n h ( M )] /n h ( M ). Owing to thefact that h δ h i = 0, the constant coefficient b ( M ) = − P ∞ j =2 b j ( M ) h δ j i /j ! (Fry & Gaztanaga 1993). Note thaton Fourier transforming δ h ( x | M, R ) the constant b onlycontributes to the unmeasurable k = 0 mode. The terms b ( M ) and b ( M ) represent the linear and first nonlinearbias parameters, respectively.In Fourier space, Eq. (14) can be written as: δ h ( k | M, R ) = b ( M ) δ ( k | R )+ b ( M )2 Z d q (2 π ) δ ( q | R ) δ ( k − q | R ) + . . . , (15)where δ i ( q j | R ) ≡ W ( | q j | R ) δ i ( q j ). If one inserts the SPTexpansions for the density into the local model, then, up tosecond order in the density and bias, one finds: δ h ( k | M, R ) = b ( M ) [ δ ( k | R ) + δ ( k | R )]+ b ( M )2 Z d q (2 π ) δ ( q | R ) δ ( k − q | R ) . (16)Using this approach one may then find a perturbative ex-pansion for the halo power and bispectra: P hh ( M ) = P (0)hh ( M ) + P (1)hh ( M ) + . . . (17) B hhh ( M ) = B (0)hhh ( M ) + B (1)hhh ( M ) + . . . . (18)Again, we refer to the lowest order terms in these expansionsas tree-level terms, and for these we have: e P (0)hh ( k | M ) = b ( M ) e P mm ( k ) ; (19) e B (0)hhh ( k , k | M ) = b ( M ) e B mmm ( k , k ) + b ( M ) b ( M ) × h e P mm ( k ) e P mm ( k ) + 2 cyc i , (20)where in the above expressions we have derived the spec-tra of the smoothed fields: e P ≡ W ( kR ) P ( k ), and e B ≡ W ( k R ) W ( k R ) W ( k R ) B . However, when we estimate thebispectrum from data we do not smooth the fields apartfrom the CIC assignment scheme used to obtain the densitycontrast field. As pointed out by Smith et al. (2007, 2008)and Sefusatti (2009), one way to overcome this is to adoptthe ansatz: P (0)hh ( k | M ) = e P (0)hh ( k | M, R ) W ( kR ) ; (21) B (0)hhh ( k , k , k | M ) = e B (0)hhh ( k , k , k | M, R ) W ( k R ) W ( k R ) W ( k R ) . (22) c (cid:13)000
4, we present the results for the bias pa-rameters from the various commonly used simple estimators.In § § § § § In the fluid approximation, the gravitational collapse of col-lisionless cold matter structures in the expanding Universe,can be fully characterized by specifying the evolution of thedensity δρ ( x ) and the peculiar velocity δ v ( x ) perturbations(Bernardeau et al. 2002). Focusing primarily on the densityfield, we work with models of the matter density contrast: δ ( x , t ) ≡ ρ ( x , t ) − ρ ( t ) ρ ( t ) , (1)where ¯ ρ ( t ) is the mean matter density of the Universe. InFourier space, we define its corresponding Fourier represen-tation, δ ( k ), accordingly, as δ ( x ) = Z d k (2 π ) δ ( k ) e − i k · x ⇔ δ ( k ) = Z d x δ ( x ) e i k · x . (2)It can be shown that the nonlinear equations of motionfor δ ( k ) can be solved exactly by perturbative expansionsof the type (Juszkiewicz 1981; Vishniac 1983; Goroff et al.1986): δ ( k ) = ∞ X n =1 a n ( t ) δ n ( k ) , (3)and, where δ n ( k ) is given by, δ n ( k ) = Z d q (2 π ) . . . Z d q n (2 π ) (2 π ) δ D ( k − q − · · · − q n ) × F n ( q , . . . , q n ) δ ( q ) · · · δ ( q n ) . (4)The density kernel F n is the dimensionless, homoge-neous, mode coupling function that couples together theamplitudes and phases of n initial Fourier wavemodes { δ ( q ) , . . . , δ ( q n ) } . As was shown by (Goroff et al. 1986;Makino et al. 1992; Jain & Bertschinger 1994), the n thorder kernel may be constructed recursively from thelower-order solutions. Linear theory is thus represented by F ( q ) = 1, and the first nonlinear correction by F , where F ( k , k ) = 57 + k · k k k (cid:18) k k + k k (cid:19) + 27 ( k · k ) k k . (5)The above approach defines the standard perturbation the-ory (hereafter SPT). Before moving on, we note that theabove statements are only exactly true for the Einstein-de Sitter model. However, it has been shown that the F kernel is almost independent of cosmology (Fry 1994;Bouchet et al. 1995; Hivon et al. 1995). We therefore adoptEq. (5) when dealing with the density at second order. c (cid:13) , 000–000 odelling halo bias using the bispectrum Owing to the stochastic nature of the density field, we arenot interested in reproducing a specific density field per se ,but instead in characterizing its statistical properties. In thiswork we focus on 2- and 3-point correlation functions inFourier space. These we may write as: h δ ( k ) δ ( k ) i ≡ (2 π ) δ D ( k ) P mm ( k ) ; (6) h δ ( k ) δ ( k ) δ ( k ) i ≡ (2 π ) δ D ( k ) B mmm ( k , k , k ) , (7)where P mm ( k ) and B mmm ( k , k , k ) constitute definitionsof the power and bispectrum. For the Dirac delta functionswe used the short-hand notation δ D ( k ...n ) ≡ δ D ( k + . . . + k n ) and these guarantee that P and B are translationally in-variant. This is an important property for estimation, sinceit means that we should consider only closed pairs and tri-angles in Fourier space: P k i = 0.The perturbative expansion of the density field de-scribed in the previous section implies that P mm and B mmm may also be described in a perturbative fashion. Hence, h δ ( k ) δ ( k ) i = h [ δ ( k ) + δ ( k ) + . . . ] × [ δ ( k ) + δ ( k ) + . . . ] i ; (8) h δ ( k ) . . . δ ( k ) i = h [ δ ( k ) + δ ( k ) + . . . ] . . . × [ δ ( k ) + δ ( k ) + . . . ] i . (9)Since we are assuming that the initial Fourier modes areGaussianly distributed, i.e. the phase of each initial modeis uniformly random φ ∈ [0 , π ], modes must cancel inpairs. Hence, Wick’s theorem applies, and so odd prod-ucts of initial Fourier modes must vanish: h δ ( k ) δ ( k ) i = h δ ( k ) δ ( k ) δ ( k ) i = 0 . This leads us to write the per-turbative expansions for P mm and B mmm as: P mm ( k ) = P (0)mm ( k ) + P (1)mm ( k ) + . . . ; (10) B mmm ( k , k ) = B (0)mmm ( k , k ) + B (1)mmm ( k , k ) + . . . (11)We shall refer to the lowest order terms in the expansions as‘tree-level’ terms. For P , P (0) is simply the linear spectrum,while for B the tree-level term can be written: B (0)mmm ( k , k ) = 2 P (0)mm ( k ) P (0)mm ( k ) F ( k , k ) + 2 cyc . (12)In this work we shall mainly be dealing with tree-level quan-tities; we now set: P (0)mm = P mm and B (0)mmm = B mmm , unlessotherwise stated.Another statistical quantity commonly used toexplore galaxy clustering is the reduced bispectrum(Peebles & Groth 1975; Scoccimarro et al. 1998), which canbe defined: Q mmm ( k , k , k ) ≡ B mmm ( k , k , k ) P mm ( k ) P mm ( k ) + 2 cyc . (13)As will be made clear below, the importance of this statisticbecomes apparent when one considers non-Gaussian termsthat are generated by simple quadratic products of Gaussianfields. In this case Q mmm simply scales as a constant. In this study we investigate the relation between the clus-tering of dark matter haloes and total matter. If galaxies are only formed in dark matter haloes, as is the usual assump-tion for all models of galaxy formation (White & Rees 1978),then understanding the clustering of haloes is an essentialcomponent of any theory of galaxy biasing (Smith et al.2007). In the local model of halo biasing, the number den-sity of dark matter haloes of mass scale M , smoothed overa scale R , can be expressed as a function of the local matterdensity, smoothed on the scale R . This function may thenbe Taylor expanded to give (Fry & Gaztanaga 1993; Coles1993; Mo et al. 1997; Smith et al. 2007): δ h ( x | M, R ) = ∞ X j =0 b j ( M ) j ! [ δ ( x | R )] j , (14)where we defined the smoothed halo over-density to be, δ h ( x | M, R ) ≡ [ n h ( x | M, R ) − n h ( M )] /n h ( M ). Owing to thefact that h δ h i = 0, the constant coefficient b ( M ) = − P ∞ j =2 b j ( M ) h δ j i /j ! (Fry & Gaztanaga 1993). Note thaton Fourier transforming δ h ( x | M, R ) the constant b onlycontributes to the unmeasurable k = 0 mode. The terms b ( M ) and b ( M ) represent the linear and first nonlinearbias parameters, respectively.In Fourier space, Eq. (14) can be written as: δ h ( k | M, R ) = b ( M ) δ ( k | R )+ b ( M )2 Z d q (2 π ) δ ( q | R ) δ ( k − q | R ) + . . . , (15)where δ i ( q j | R ) ≡ W ( | q j | R ) δ i ( q j ). If one inserts the SPTexpansions for the density into the local model, then, up tosecond order in the density and bias, one finds: δ h ( k | M, R ) = b ( M ) [ δ ( k | R ) + δ ( k | R )]+ b ( M )2 Z d q (2 π ) δ ( q | R ) δ ( k − q | R ) . (16)Using this approach one may then find a perturbative ex-pansion for the halo power and bispectra: P hh ( M ) = P (0)hh ( M ) + P (1)hh ( M ) + . . . (17) B hhh ( M ) = B (0)hhh ( M ) + B (1)hhh ( M ) + . . . . (18)Again, we refer to the lowest order terms in these expansionsas tree-level terms, and for these we have: e P (0)hh ( k | M ) = b ( M ) e P mm ( k ) ; (19) e B (0)hhh ( k , k | M ) = b ( M ) e B mmm ( k , k ) + b ( M ) b ( M ) × h e P mm ( k ) e P mm ( k ) + 2 cyc i , (20)where in the above expressions we have derived the spec-tra of the smoothed fields: e P ≡ W ( kR ) P ( k ), and e B ≡ W ( k R ) W ( k R ) W ( k R ) B . However, when we estimate thebispectrum from data we do not smooth the fields apartfrom the CIC assignment scheme used to obtain the densitycontrast field. As pointed out by Smith et al. (2007, 2008)and Sefusatti (2009), one way to overcome this is to adoptthe ansatz: P (0)hh ( k | M ) = e P (0)hh ( k | M, R ) W ( kR ) ; (21) B (0)hhh ( k , k , k | M ) = e B (0)hhh ( k , k , k | M, R ) W ( k R ) W ( k R ) W ( k R ) . (22) c (cid:13)000 , 000–000 J. E. Pollack, R. E. Smith & C. Porciani
On applying this ‘de-smoothing’ operation to Eqs (19)and (20), one finds: P (0)hh ( k | M ) = b ( M ) P mm ( k ) ; (23) B (0)hhh ( k , k | M ) = b ( M ) B mmm ( k , k ) + b ( M ) b ( M ) × f W k , k P mm ( k ) P mm ( k ) + 2 cyc . (24)where we have defined the function (Sefusatti 2009): f W k , k ≡ W ( | k | R ) W ( | k | R ) W ( | k + k | R ) . (25)Note that in the limit of very large scales or arbitrarily smallsmoothing scales, k i R → i ∈ { , , } , Eq. (24) approx-imates to: B (0)hhh ( k , k | M ) ≈ b ( M ) B mmm ( k , k ) + b ( M ) b ( M ) × [ P mm ( k ) P mm ( k ) + 2 cyc ] . (26)Again, since in this paper we are only considering tree-levelexpressions we shall take P (0)hh → P hh and B (0)hhh → B hhh .Considering now the reduced halo bispectrum, it maybe defined in a similar fashion to Eq. (13): Q hhh ( k , k , k | M ) ≡ B hhh ( k , k | M ) P hh ( k | M ) P hh ( k | M ) + 2 cyc . (27)On inserting our tree-level expressions from Eqs (23)and (24), we find that Q hhh ( M ) = Q mmm b ( M ) + b ( M ) b ( M ) α ( k , k , k ) , (28)where we have defined, α ( k , k , k ) ≡ f W k , k P ( k ) P ( k ) + 2 cyc P ( k ) P ( k ) + 2 cyc . (29)Again, in the limit of very large scales or arbitrarily smallsmoothing scales and α ( k , k , k ) →
1, the above expres-sion approximates to: Q hhh ( k , k , k | M ) ≈ Q mmm ( k , k , k ) b ( M ) + b ( M ) b ( M ) . (30)We now see the utility of the reduced bispectrum: if one con-structs halo/galaxy density fields from a local transforma-tions of the matter density, then the lowest order nonlinearcorrections will lead to a function that is a scaled version ofthe matter Q mmm , plus a constant offset. Moreover, if thedensity field were simply Gaussian, then estimates of Q hhh on large scales would directly measure b /b . N -BODY SIMULATIONS For our investigations of the bias, we use an ensemble of40 large N -body simulations, executed on the zBOX-2 andzBOX-3 supercomputers at the University of Z¨urich. We useonly the z = 0 outputs from the simulations. Each simula-tion was performed using the publicly available Gadget-2 code (Springel 2005), and followed the nonlinear evolutionunder gravity of N = 750 equal-mass particles in a comov-ing cube of length L sim = 1500 h − Mpc.The cosmological model that we simulated was analo-gous to the basic vanilla ΛCDM model determined by theWMAP experiment (Komatsu et al. 2009): matter density Ω m = 0 .
25, vacuum density Ω Λ = 0 .
75, power spectrum nor-malization σ = 0 .
8, power spectral index n = 1, and dimen-sionless Hubble parameter h = 0 .
7. The transfer functionfor the simulations was generated using the publicly avail-able cmbfast code (Seljak & Zaldarriaga 1996; Seljak et al.2003), with high sampling of the spatial frequencies on largescales. Initial conditions were set at redshift z = 49 us-ing the serial version of the publicly available code(Scoccimarro 1998; Crocce et al. 2006).Dark matter halo catalogues were generated for eachsimulation using the Friends-of-Friends (FoF) algorithm(Davis et al. 1985), with the linking-length parameter b =0 .
2, where b is the fraction of the inter-particle spacing. Forthis we employed the fast parallel B-FoF code, provided tous by V. Springel. The minimum number of particles forwhich an object is considered to be a bound halo was setat 20 particles. This gave a minimum host halo mass of M min = 1 . × h − M ⊙ . For our analysis of the bias,we use the full sample of haloes and this corresponded toroughly N h ≈ . × haloes per simulation. Further de-tails of the simulations may be found in Smith (2009). Before we examine halo bias in the context of the bispec-trum, we explore two alternative methods for studying thebias. We first evaluate the second-order local biasing modeldirectly, by comparing in a point-wise fashion the halo andmatter density fields, smoothed over a range of scales. Then,we use power spectra to determine an effective large-scalebias.
One obvious way to examine the local model of biasing isto simply construct a scatter plot of the local density ofdark matter haloes against the local density of dark matterin the simulations (see for example Sheth & Lemson 1999;Dekel & Lahav 1999). As was discussed in § R .We generate the smoothed density fields as follows: weassign particles/haloes to a Fourier grid using the CIC algo-rithm (c.f. § W ( kR ) ≡ exp (cid:2) − ( kR ) / (cid:3) . (31)Finally, on taking the inverse Fourier transform, we obtainthe smoothed δ ( x | R ) and δ h ( x | R ). We perform the aboveprocedure for 28 of the ensemble of simulations and considerthe filter scales: R = { , , } h − Mpc.In Figure 1 we present the bin averaged scatter plotsof δ h ( x | R ) vs. δ ( x | R ), averaged over the realizations. Thecolour contours are shaded by the normalized populationdensity of that pixel, e.g. the central red region indicates thatmost of the points in the simulation are regions of densityclose to average. We also see that as the smoothing scale isdecreased (panels going from left-to-right), that the scatterincreases and that there are more points that have higher c (cid:13) , 000–000 odelling halo bias using the bispectrum Figure 1.
Scatter plots of δ h ( x ) versus δ ( x ) smoothed with a Gaussian filter of various scales averaged over the realizations. From leftto right, the panels correspond to the smoothing scales R = { , , } h − Mpc. The color coding denotes the log of the populationdensity, i.e. the red region corresponds to the largest concentration of points and the white background to null values. The dot-dashedline in each panel denotes the local halo bias model up to second-order with the best-fitting bias parameters averaged over 28 realizations.
R b ± σ b b ± σ b b ± σ b [ h − Mpc] × −
50 1.3 ± ± ± ± ± ± ± ± ± Table 1.
Average of the mean bias parameters and the root-meansquare errors for the local halo bias model up to second-orderaveraged over 28 realizations determined from fitting the scatterplots of the halo and matter density fields smoothed on scales k s = { . , . , . } h Mpc − . and lower density. Conversely, as the filter scale is increasedthe relation becomes tighter and more linear. One obviousconclusion that may be drawn from this behaviour is thatthe bias relation is certainly not deterministic.In order to obtain a more quantitative understanding,we next consider fitting for the parameters of the local biasmodel at second-order. From Eq. (14) we have: δ h ( x | M, R ) = b ( M )+ b ( M ) δ ( x | R )+ b ( M )2 [ δ ( x | R )] . (32)We perform a least-squares analysis on each realization, andthen average over the resulting set of bias parameters to ob-tain the mean parameters: b ( M ), b ( M ), and b ( M ). The1 σ -errors are then estimated in the usual way, as quadraticdeviations from the sample mean. In Fig. 1 we plot the re-sultant best-fit local model as the dot-dashed line in each ofthe three panels.The information on the parameters is summarized inTable 1 as a function of the filter scale R . This clearlyshows that the estimates of b increase as the smooth-ing scale is decreased, whereas those for b appear to beparabolical. Naively, one might expect that the nonlinearbias terms should approach zero as the amount of smooth-ing is increased and nonlinearities are washed out, however,at R = 50 h − Mpc even with σ ( x | R ) <
1, the fluctuations are still significant enough to yield a non-zero b . Note alsothat in all cases b = 0.The local model, as written in Eq. (14), asserts that theparameters b i are independent of the smoothing scale R , andwe, therefore, consider the implications as follows. Supposethat nonlinear bias is exactly as described by Eq. (32), butthat the coefficients are not independent of the smoothingscale. Let us now consider the results that would be obtainedfrom measurements for two smoothing scales R a and R b .From Eq. (32) we would have: δ h ( x | R a ) = b a + b a δ ( x | R a ) + b a δ ( x | R a )] ; (33) δ h ( x | R b ) = b b + b b δ ( x | R b ) + b b δ ( x | R b )] . (34)Supposing now that we desmoothed each of the fields, byFourier transforming and dividing out the appropriate win-dow function. We would then have: δ h ( k ) = b a δ ( k ) + b a Z d q (2 π ) δ ( q ) δ ( | k − q | ) f W q , k − q ( R a ) δ h ( k ) = b b δ ( k ) + b b Z d q (2 π ) δ ( q ) δ ( | k − q | ) f W q , k − q ( R b ) . In order for the above equations to be equivalent, then wemust have b a = b b (35) b a = b b (cid:20) W ( qR b ) W ( qR a ) W ( | k − q | R a ) W ( | k − q | R b ) W ( kR a ) W ( kR b ) (cid:21) . (36)The last of the two equations may only be satisfied if andonly if R a = R b or { kR, qR, | k − q | R } ≪
1. Since the δ h ( x | R )vs. δ ( x | R ) method is inherently a real space measure it in-volves contributions from all Fourier modes. It is thereforedifficult to ensure that b a = b b .We conclude that the above method will not be a safeway to recover bias parameters independent of the smooth-ing scale. We now turn to Fourier space methods. c (cid:13)000
1. Since the δ h ( x | R )vs. δ ( x | R ) method is inherently a real space measure it in-volves contributions from all Fourier modes. It is thereforedifficult to ensure that b a = b b .We conclude that the above method will not be a safeway to recover bias parameters independent of the smooth-ing scale. We now turn to Fourier space methods. c (cid:13)000 , 000–000 J. E. Pollack, R. E. Smith & C. Porciani
Figure 2.
Scale dependence of effective bias parameters b NLhh , b Lhh , b NLhm and b Lhm (c.f. Eqs (38) and (39)), estimated from the auto- andcross-power spectrum as a function of wavemode. For the left and central panels: solid blue and open red symbols denote the bias when P hh is not and is shot noise corrected, respectively. The first panel shows b hh when the nonlinear matter power spectrum is used; thesecond panel shows the same but when the linear matter power spectrum is used; the third panel shows b hm , where the red stars andblue points denote the case where the nonlinear and linear matter power spectra are used, respectively. k [ h Mpc − ] b NLhh b NL , SChh b Lhh b L , SChh ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± Table 2.
Weighted average estimates of the effective bias, b hh (see Eq. (38) for a defintions). We now use various halo power spectra to derive estimatesfor an effective large-scale halo bias.In order to do this, we first measure the Fourier trans-form of the matter and halo density fields as described inAppendix A. The halo-halo, halo-mass and mass-mass powerspectra, { P hh , P hm , P mm } , are then estimated from the databy performing the following sums: b P µν ( k l ) = V µ N ( k ) N ( k ) X m =1 δ µ ( k l ) δ ∗ ν ( k l ) , (37)where { µ, ν } ∈ { m , h } , V µ is the sample volume (which inour case is the simulation volume), and where N ( k ) are thenumber of Fourier modes in a shell of thickness ∆ k .Following Smith et al. (2007), we next construct the es-timators: b b NLhh = 1 N s N s X i =1 s b P hh ( k i ) b P mm ( k i ) ; b b Lhh = 1 N s N s X i =1 s b P hh ( k i ) P Lmm ( k i ) ; (38) b b NLhm = 1 N s N s X i =1 b P hm ( k i ) b P mm ( k i ) ; b b Lhm = 1 N s N s X i =1 b P hm ( k i ) P Lmm ( k i ) , (39)where N s is the number of simulations and P Lmm is the linearmatter power spectrum. Note that in the case of b hh we alsoconsider shot-noise corrected versions of these two estima-tors, i.e. we correct P hh using Eq. (A11). We denote thesebias estimates by b NL , SChh and b L , SChh , respectively. Finally, wedetermine the 1 σ errors by evaluating the variance of eachrealization against the mean. The first panel of Fig. 2 shows b NLhh (solid blue datapoints) and b NL , SChh (open red points). The bias for the shot-noise corrected terms remains roughly constant at ∼ . k ∼ . h Mpc − and with very small er-rors, indicating that the result is highly constrained by thedata. On scales smaller than this the bias is a decreasingfunction of k . Without shot-noise correction, we find thatthe bias is strongly scale dependent, and the bias rapidlyincreases with increasing k .The second panel of Fig. 2 shows the results obtainedfrom using b Lhh (solid blue points) and b L , SChh (open hexago-nal points). The results are similar to those for b NL , SChh , butwith increased cosmic variance on large scales. An oscillationstructure is also present, this can be understood as explainedin Guzik et al. (2007). Nevertheless, comparing the two pro-vides a clear indication of the validity of the tree-level powerspectrum up to k = 0 . h Mpc − .The third panel of Fig. 2 shows the bias results b NLhm (solid red points) and b Lhm (solid blue points). The valueof b NLhm stays roughly constant for the whole scale rangeconsidered in the estimate, while b Lhm is not smooth andclearly shows the imprint of the oscillation structure. Nev-ertheless, we find b Lhm ∼ .
48, to within the errors for k < . h Mpc − . On comparing the results for b NL , SChh and b NLhm , we see that for scales k . . h Mpc − , these esti-mates are compatible and that the effective large-scale biasis roughly b ∼ .
49. Interestingly, these findings are consis-tent with real space measures of the effective large scale biasfrom cell variances (Smith & Marian 2011).In Table 2, we report the weighted average and cor- c (cid:13) , 000–000 odelling halo bias using the bispectrum responding 1 σ error on the effective bias, b hh , computedover the same k -modes corresponding to the magnitude ofthe third wavevector k for each triangle configuration. Thetabulated results for the analysis of the uncorrected dataconfirms the results shown in Figure 2, that bias is indeedscale-dependent. Applying the shot-noise correction yields amore constant effective bias, even for the range of k -modesentering into our bispectrum estimation. Interestingly, thevalue b = 1 . ± .
002 found for k ∈ [0 . , . b obtained from fitting thedensity fields smoothed on scales R = 50 h − Mpc. There-fore, if we opt to infer that the effective bias is equivalentto b over these scales, then the bispectrum (reduced bis-pectrum) should also yield this value for b when fitted overthe same scale ranges (c.f. Table 2). That is, if the local biasmodel is correct and the tree-level bispectrum (reduced bis-pectrum) is a sufficient description of the nonlinearities onthese scales.Before moving on, we point out that one can also usethe halo power spectra to define an effective b (Smith et al.2009), however we shall not explore this here. In this section, we present our main results from the analysisof the halo bispectra.
The computational code used to estimate the matter andhalo power and bispectra is a modified version of the codedeveloped in Smith et al. (2008), which itself is based onthe algorithm of Scoccimarro et al. (1998). The major mod-ification to that code, which we have implemented, is thatno random subsampling of the Fourier modes is performedto estimate the bispectrum. Instead, all modes that con-tribute to a particular triangle configuration are used. Inthis work we use a FFT grid of size N g = 512 to esti-mate the power and bispectra. We only evaluate trianglesthat have k = 2 k , but consider the variation of B withthe angular separation of the two vectors. The largest scaleat which we estimate the bispectrum is k = 0 . h Mpc − ,and this is ≈ . k f , where k f = 2 π/L ≈ . h Mpc − . Fur-ther details of the bispectrum estimation procedure may befound in Appendix A.Figure 3 shows the ensemble-averaged shot-noise cor-rected results for the halo bispectra B hhh (open red squares)and matter bispectra B mmm (solid blue diamonds), mea-sured from the ensemble of N -body simulations. The fourpanels show the results obtained for the scales: k = { . , . , . , . } h Mpc − . The error bars are the 1- σ errors on the mean, derived from the ensemble to ensemblevariation. The solid red line represents the tree-level pre-diction for B mmm as given by Eq. (12). We see that thisappears to be a good description of the B mmm estimates forthe scales that we have considered. We notice that, for thecase k = 0 . h Mpc − , the theory systematically under-predicts the measurements for θ/π ∼ . § Q hhh and Q mmm . Our method for estimating the bias parameters follows anapproach similar to that presented by Scoccimarro (2000)and Porciani & Giavalisco (2002). To start, we take a χ function that is a quadratic form of the type: χ ( b , b ) = N θ X i =1 N θ X j =1 ∆ i ( b , b ) d r − ij ∆ j ( b , b ) , (40)where N θ is the number of angular bins considered and∆ i ≡ b B hhh ( k , k , θ i ) − B modhhh ( k , k , θ i | b , b ) σ hhh ( k , k , θ i ) . (41)Note that in dividing the difference between the estimate ofthe ensemble average bispectrum ( b B ) and the model predic-tion ( B mod ) by the standard deviation ( σ hhh ), d r − ij is in factthe inverse correlation matrix. Recall that the correlationand covariance matrices are related by: r ij = C ij / p C ii C jj .In order to minimize this χ function and so recover thebest-fit bias parameters, we need an estimate of d r − ij , andwe do this using the standard unbiased estimator: b r ij = 1( N sim − N sim X k =1 e ∆ ki e ∆ kj , (42)where N sim is the number of simulations and where e ∆ i ≡ b B ( k )hhh ( k , k , θ i ) − b B hhh ( k , k , θ i ) σ hhh ( k , k , θ i ) , (43)where in the above b B ( k )hhh is the k th estimate of the bispec-trum and where b B hhh is the mean of the ensemble.Next, we use singular-value-decomposition (SVD) to in-vert the b r ij . For the estimate of the inverse correlation ma-trix, we utilize principal component analysis (PCA) to re-move some of the noisy eigenvectors. We select the fractionof principal components that account for 95% of the vari-ance. According to this selection criteria we typically retain15 out of 20 of the most ‘dominant’ eigenmodes. Note thatas pointed out in Hartlap et al. (2007), d C − = b C − . How-ever, since we are using PCA, this should be a subdominantcorrection. Thus, we may approximate Eq. (40) as: χ = N θ X i =1 N θ X j =1 ∆ i ( b , b )[ R T Λ R ] − ij ∆ j ( b , b ) , = N θ X i =1 N θ X j =1 ∆ i ( b , b ) N θ X l =1 R Til Λ − ll R lj ∆ j ( b , b ) , ≈ N θ X l =1 Λ − ll Y l Θ ll , (44)where the correlation matrix r was diagonalized by rota-tion into its eigenbasis, i.e. r = R T Λ R , with Λ repre-senting a diagonal matrix of eigenvalues. We also defined Y l ≡ P N θ i =1 R li ∆ i ( b , b ). Note that in the final approximateexpression we include a matrix Θ ll , this a diagonal matrixwith entries either 1 or 0, depending on whether the eigen-vector is to be retained or cut from the PCA reconstruction.Finally, the χ ( b , b ) function was minimized using theLevenberg-Marquardt routine for non-linear least squaresfitting. c (cid:13)000
The computational code used to estimate the matter andhalo power and bispectra is a modified version of the codedeveloped in Smith et al. (2008), which itself is based onthe algorithm of Scoccimarro et al. (1998). The major mod-ification to that code, which we have implemented, is thatno random subsampling of the Fourier modes is performedto estimate the bispectrum. Instead, all modes that con-tribute to a particular triangle configuration are used. Inthis work we use a FFT grid of size N g = 512 to esti-mate the power and bispectra. We only evaluate trianglesthat have k = 2 k , but consider the variation of B withthe angular separation of the two vectors. The largest scaleat which we estimate the bispectrum is k = 0 . h Mpc − ,and this is ≈ . k f , where k f = 2 π/L ≈ . h Mpc − . Fur-ther details of the bispectrum estimation procedure may befound in Appendix A.Figure 3 shows the ensemble-averaged shot-noise cor-rected results for the halo bispectra B hhh (open red squares)and matter bispectra B mmm (solid blue diamonds), mea-sured from the ensemble of N -body simulations. The fourpanels show the results obtained for the scales: k = { . , . , . , . } h Mpc − . The error bars are the 1- σ errors on the mean, derived from the ensemble to ensemblevariation. The solid red line represents the tree-level pre-diction for B mmm as given by Eq. (12). We see that thisappears to be a good description of the B mmm estimates forthe scales that we have considered. We notice that, for thecase k = 0 . h Mpc − , the theory systematically under-predicts the measurements for θ/π ∼ . § Q hhh and Q mmm . Our method for estimating the bias parameters follows anapproach similar to that presented by Scoccimarro (2000)and Porciani & Giavalisco (2002). To start, we take a χ function that is a quadratic form of the type: χ ( b , b ) = N θ X i =1 N θ X j =1 ∆ i ( b , b ) d r − ij ∆ j ( b , b ) , (40)where N θ is the number of angular bins considered and∆ i ≡ b B hhh ( k , k , θ i ) − B modhhh ( k , k , θ i | b , b ) σ hhh ( k , k , θ i ) . (41)Note that in dividing the difference between the estimate ofthe ensemble average bispectrum ( b B ) and the model predic-tion ( B mod ) by the standard deviation ( σ hhh ), d r − ij is in factthe inverse correlation matrix. Recall that the correlationand covariance matrices are related by: r ij = C ij / p C ii C jj .In order to minimize this χ function and so recover thebest-fit bias parameters, we need an estimate of d r − ij , andwe do this using the standard unbiased estimator: b r ij = 1( N sim − N sim X k =1 e ∆ ki e ∆ kj , (42)where N sim is the number of simulations and where e ∆ i ≡ b B ( k )hhh ( k , k , θ i ) − b B hhh ( k , k , θ i ) σ hhh ( k , k , θ i ) , (43)where in the above b B ( k )hhh is the k th estimate of the bispec-trum and where b B hhh is the mean of the ensemble.Next, we use singular-value-decomposition (SVD) to in-vert the b r ij . For the estimate of the inverse correlation ma-trix, we utilize principal component analysis (PCA) to re-move some of the noisy eigenvectors. We select the fractionof principal components that account for 95% of the vari-ance. According to this selection criteria we typically retain15 out of 20 of the most ‘dominant’ eigenmodes. Note thatas pointed out in Hartlap et al. (2007), d C − = b C − . How-ever, since we are using PCA, this should be a subdominantcorrection. Thus, we may approximate Eq. (40) as: χ = N θ X i =1 N θ X j =1 ∆ i ( b , b )[ R T Λ R ] − ij ∆ j ( b , b ) , = N θ X i =1 N θ X j =1 ∆ i ( b , b ) N θ X l =1 R Til Λ − ll R lj ∆ j ( b , b ) , ≈ N θ X l =1 Λ − ll Y l Θ ll , (44)where the correlation matrix r was diagonalized by rota-tion into its eigenbasis, i.e. r = R T Λ R , with Λ repre-senting a diagonal matrix of eigenvalues. We also defined Y l ≡ P N θ i =1 R li ∆ i ( b , b ). Note that in the final approximateexpression we include a matrix Θ ll , this a diagonal matrixwith entries either 1 or 0, depending on whether the eigen-vector is to be retained or cut from the PCA reconstruction.Finally, the χ ( b , b ) function was minimized using theLevenberg-Marquardt routine for non-linear least squaresfitting. c (cid:13)000 , 000–000 J. E. Pollack, R. E. Smith & C. Porciani
Figure 3.
Ensemble-averaged matter and halo bispectrum measurements for 40 LCDM N -body simulations in real-space in comparisonwith the PT models at tree-level. Each panel shows the shot-noise corrected bispectrum measurements as a function of angle for a varietyof triangle configurations at different scale ranges set by k = { .
03 0 .
04 0 .
05 0 . } h Mpc − and k = 2 k . The blue solid symbolsrepresent the matter bispectrum, whereas the open squares denote the halo bispectrum. The tree-level bispectrum is represented by thesolid orange line and the local halo bias model with the best-fitting parameters listed in Table 4 is the dashed violet line. To the best-fit parameters ( b , b ), we assign both systematicand statistical errors.In our context, the systematic errors correspond to theerrors induced in the best-fit parameters from fitting thedata with a noisy inverse covariance matrix (or correlationmatrix). Owing to the relatively low number of simulations( N sim = 40), we expect that Eq. (42) provides a noisy es-timate of r − ij . In order to estimate the errors this has onthe best-fit parameters we employ the jackknife subsamplingmethod (see for example Norberg et al. 2009). This involvesslicing the total data set into N sub subsamples. Then a re-sampling of the data is obtained by excluding one of the sub-samples from the set. From this resampling we then estimatethe mean statistic of interest and the inverse correlation ma-trix as described in the previous section. The resampled dataset is then used to determine a new estimate of the best-fitbias parameters. This procedure is then repeated for all of the possible N sub resamplings of the data. In our particu-lar case we treat the measurements from each simulation asthe regions to be included or excluded, and this gives us 40jackknife estimates of the bias parameters ( b , b ). The pa-rameter covariance matrix for the systematic errors can becomputed as (Norberg et al. 2009): b C JK [ b i , b j ] = N sub − N sub N sub X k =1 ( b i,k − ˆ¯ b j )( b j,k − ˆ¯ b j ) , (45)where b i,k is the estimate of b i from the k th resampling ofthe data and ˆ¯ b i is the estimate of the mean b i obtained fromall of the resamplings.The statistical error is obtained directly from the non-linear least-squares analysis. The routine mrqmin providesan approximation to the errors on the parameters that cor-responds to a ∆ χ ≈ c (cid:13) , 000–000 odelling halo bias using the bispectrum Figure 4.
Ensemble-averaged matter and halo reduced bispectrum measurements of the 40 LCDM N -body simulations in real-space incomparison with the PT models at tree-level. Point and line styles as in Figure 3. correspond to either ∆ χ = (2 . , . ∼ σ , ∼ σ ) errors for a two-parameter model.Given that we consider the two forms of error: system-atic and statistical), and that one is never consistently largerthan the other, in all forthcoming tables, we choose to re-port only the total error. This is obtained simply from thetwo errors added in quadrature. B Before we report the estimates of the halo bias parameters,we first present a test of the validity of the tree-level modelfor the matter bispectrum. We do this by applying the χ test described above, to the B mmm and Q mmm data, andso fit for b and b . Note that since the total number ofprincipal components retained equals 15, then for a two-parameter model the number of degrees-of-freedom equals13. If the tree-level expressions in the large-scale limit asgiven by Eqs (26) and (30) are correct, then we should expectto find b = 1 and b = 0.Table 3 presents the best-fit nonlinear bias parameters for the four different bispectrum scale ranges discussed ear-lier. In the analysis we fit the shot-noise corrected bispectra.The χ values (end column of the table), confirms that thetree-level expressions B and Q , provide good fits forthe triangle configurations with k = { . , . } h Mpc − .However, for k = { . , . } h Mpc − the fits are poorgiven the χ estimates, and we see that, for both B and Q ,they yield non-zero values for b at 1 σ . The results also im-ply that the failure of the tree-level model on these scalesis more severe for Q than for B . This can be understoodby noting that b from Q mmm shows a prominent departurefrom unity, whereas B mmm does not (although the deviationstill exceeds 2 σ ). We thus conclude that it is likely that thetree-level expressions for the halo bispectra will only be validfor k . h Mpc − , for our chosen bispectra configura-tions. b and b from halo bispectra Table 4 presents the best-fit nonlinear bias parameters andtheir respective 1 σ errors in quadrature, obtained from the c (cid:13) , 000–000 J. E. Pollack, R. E. Smith & C. Porciani k [ h Mpc − ] b ± σ b b ± σ b χ B mmm ± ± Q mmm ± ± B mmm ± ± Q mmm ± ± B mmm ± ± Q mmm ± ± B mmm ± ± Q mmm ± ± Table 3.
Assessment of the validity of the tree-level modellingby fitting the matter bispectra and reduced bispectra. Column 1:bispectrum triangle scale; column 2: statistic, where B mmm and Q mmm are shot-noise corrected; columns 3 and 4: best-fit b and b along with 1 σ errors; column 5: χ . χ analysis of B hhh and Q hhh . Note that we present theresults for both the uncorrected and shot-noise correctedmeasurements, indicated in the table by superscript ‘SC’.Table 4 also shows the χ value of these best-fit parametersas an indication of the goodness-of-fit.In Figures 3 and 4 we also show the tree-level theo-retical models for B hhh and Q hhh (dashed lines), where thebest-fit bias parameters from Table 4 have been used. Thesefigures demonstrate that, at least by-eye, the tree-level mod-els provide a reasonable description of the data. However amore detailed inspection of Table 4 reveals some importantdiscrepancies.For the case of fitting B , the shot-noise correction is lessimportant, as we see that the estimates of b for all bispectraconfigurations with and without shot-noise corrections areconsistent to within the errors, and have b ∼ .
4. However, b shows systematic differences, being more negative if thecorrection is made, and for this we find that b ∼ − .
25. Onthe other hand, for the case of Q , the results clearly showthat the shot-noise subtraction has an important effect onthe recovered values for the bias parameters. If the shot-noise is not corrected, then we see that the estimates for b increase systematically as we go from triangle configurationswith k = 0 . h Mpc − to k = 0 . h Mpc − . Whereas ifit is corrected, then we find b ∼ . b ∼ − . B and Q we see that, whilst the values for b disagree significantly,surprisingly those for b remain consistent at the 1 σ level.The χ function of Eq. (41) may be interpreted asa Gaussian likelihood if we make the transformation, L ( { B hhh }| b , b ) ∝ exp[ − χ / p ( b , b |{ B hhh } ).Figure 5 shows the 1 σ likelihood confidence contoursin the posterior probability for the nonlinear bias parame-ters for the four scales considered according to our methodof analysis described above. The solid lines denote the sizeof the confidence regions at the 68% and 95% level (i.e.∆ χ ≈ . , .
17) when we construct a correlation matrixfrom the 40 realizations without regard for the systematicuncertainty. The dashed-lines demonstrate the magnitude k [ h Mpc − ] b ± σ b b ± σ b χ B hhh ± ± B SChhh ± ± Q hhh ± ± Q SChhh ± ± B hhh ± ± B SChhh ± ± Q hhh ± ± Q SChhh ± ± B hhh ± ± B SChhh ± ± Q hhh ± ± Q SChhh ± ± B hhh ± ± B SChhh ± ± Q hhh ± ± Q SChhh ± ± Table 4.
Best-fit bias parameters from fitting the halo-halo-halobispectra and reduced bispectra. Column 1: bispectrum trianglescale; column 2: statistic, where B hhh and Q hhh are raw, andwhere B SChhh and Q SChhh are shot-noise corrected; columns 3 and 4:best-fit b and b along with 1 σ errors; column 5: χ . at which the 68% and 95% confidence regions expand fol-lowing our generation of a set of jackknife subsamples tomonitor the effect due to the implicit error associated withthe estimated correlation matrix. Hence, we clearly see therelevance of accounting for the uncertainty of the correlationmatrix when obtaining the bias parameter constraints. Thediscrepancy between the resulting jackknife error ellipses for B and Q is less severe than the likelihood contours obtainedfrom the complete sample where the level of agreement im-proves progressing to large scales, yet this might be due tothe fact that the statistical error is more prominent at largerscales. Interestingly, the overlap of the two likelihood regionsat 2 σ for k =0.03, 0.04 and 0.05 h Mpc − occurs with therectangular region or strip denoting the effective bias mea-sure, b NL , SChh , at 1 σ . These set of panels in Figure 5 conveypictorially the information obtained from the results in Ta-ble 4 that the likelihood contours from analysis of Q showan evolution with decreasing scale toward larger and larger b , whereas the constraints on b remain consistent. Lastly,the constraints obtained from analyzing the Q -amplitudesare much weaker than those coming from the bispectrum. In § Q . In thissection, we attempt to develop the use of cross-bispectra, asmeasures of the bias that are less susceptible to discretenesseffects. The use of cross-correlations in large-scale structurework has long been known as a way of reducing discretenesscorrections (Peebles 1980). However, it is only relatively re-cent that it has been applied to study bias (Smith et al.2007; Dalal et al. 2008; Smith 2009; Padmanabhan et al.2009; Desjacques et al. 2009; Pillepich et al. 2010). c (cid:13) , 000–000 odelling halo bias using the bispectrum Figure 5.
Evolution of the likelihood contours for the bias parameters b and b , estimated from B hhh and Q hhh , with scale. Thesolid lines denote the 68% and 95% confidence intervals, obtained from a full exploration of the likelihood surface around the bestfit values; the dashed lines denote the same, but where the jackknife parameter covariance matrix from Eq. (45) has been used todetermine the error contours. The top left, top right, bottom left and bottom right panels show the results for triangle configurationswith k = { . , . , . , . } h Mpc − , respectively. The vertical black lines denote the effective bias parameter b NL , SChh , using thesame wavemodes that enter into the bispectrum estimates.
We may define the halo cross-bispectra as follows: h δ h ( k ) δ h ( k ) δ ( k ) i = (2 π ) δ D ( k ) B hhm ( k , k , k ) ;(46) h δ h ( k ) δ ( k ) δ ( k ) i = (2 π ) δ D ( k ) B hmm ( k , k , k ) . (47)We then symmeterize these quantities by the operations: B (sym)hhm = [ B hhm + B hmh + B mhh ] / B (sym)hmm = [ B hmm + B mhm + B mmh ] / . (49)For ease of notation we shall now simply take B (sym)hhm ≡ B hhm and B (sym)hmm ≡ B hmm , unless otherwise indicated. We maynow also define the cross-reduced bispectra as: Q hhm ≡ B hhm /P P hhm ; (50) Q hmm ≡ B hmm /P P hmm , (51) where we have for the denominators (again symmetrized): P P hhm = 23 [ P hh ( k ) P hm ( k ) + 2 cyc ]+ 13 [ P hm ( k ) P hm ( k ) + 2 cyc ] ; (52) P P hmm = 23 [ P hm ( k ) P mm ( k ) + 2 cyc ]+ 13 [ P hm ( k ) P hm ( k ) + 2 cyc ] . (53)The relations for P P hhm and
P P hmm can easily be con-structed using a graphical approach. Let us consider threenodes two of which are the same and the third is different(we shall think of the nodes as the density fields). Label thesenodes 1, 2, and 3. Now consider all possible ways to connectthe three nodes together by two edges. When two nodes,which are the same, connect together this gives us an auto-power spectrum with a delta function, and when two nodesthat are different connect together this gives us a cross- c (cid:13)000
P P hmm can easily be con-structed using a graphical approach. Let us consider threenodes two of which are the same and the third is different(we shall think of the nodes as the density fields). Label thesenodes 1, 2, and 3. Now consider all possible ways to connectthe three nodes together by two edges. When two nodes,which are the same, connect together this gives us an auto-power spectrum with a delta function, and when two nodesthat are different connect together this gives us a cross- c (cid:13)000 , 000–000 J. E. Pollack, R. E. Smith & C. Porciani power spectrum and delta function. One may then sym-metrize the results by considering all possible relabellingsof the nodes and dividing by three.In Appendix B we calculate the tree-level cross-bispectra, B hmm and B hhm , in the local model of halo bi-asing. The main results are: B (0)hmm ( k , k , k ) ≈ b ( M ) B (0)mmm ( k , k , k )+ b ( M )3 hf W k , k P (0)mm ( k ) P (0)mm ( k ) + 2 cyc i ; (54) B (0)hhm ( k , k , k ) ≈ b ( M ) B (0)mmm ( k , k , k ) + 13 b ( M ) × b ( M ) hf W k , k P (0)mm ( k ) P (0)mm ( k ) + 2 cyc i . (55)In the limit of large scales, and or small smoothing scales,the filter functions f W k , k → B (0)hmm ≈ b B (0)mmm + b h P (0)mm ( k ) P (0)mm ( k ) + 2 cyc i ; (56) B (0)hhm ≈ b B (0)mmm + 13 b b h P (0)mm ( k ) P (0)mm ( k ) + 2 cyc i . (57)At second order in the nonlinear bias and PT, the cross-reduced bispectra are: Q (0)hmm ≈ Q mmm b ( M ) + b ( M )2 b ( M ) + b ( M ) α ( k , k , k ) ;(58) Q (0)hhm ≈ Q mmm b ( M ) + 1 + 2 b ( M )2 b ( M ) + b ( M ) α ( k , k , k ) . (59)In the large-scale limit, α →
1, these expressions become: Q (0)hmm ≈ Q mmm b ( M ) + b ( M )2 b ( M ) + b ( M ) ; (60) Q (0)hhm ≈ Q mmm b ( M ) + 1 + 2 b ( M )2 b ( M ) + b ( M ) . (61) The cross-bispectra B hhm and B hmm can be estimated fol-lowing the algorithm described in § B hhm we have, b B dhhm ( k , k , θ ) = 13 V µ N tri N tri X ( n , n ) × n R e [ δ h ( k n ) δ h ( k n ) δ m ( k n )] + 2 cyc o , (62)and with a similar relation for b B dhmm . The reduced bispectraare estimated by dividing the above bispectrum estimatesby estimates for P P hhm and
P P hmm from Eqs (52) and (53),respectively.One further complication is constructing the correctionsfor shot-noise. This may be performed following the counts-in-cells approach of Peebles (1980) (and see also Smith2009). We find that the symmetrized corrections for B hhm and B hmm can be written as: b ¯ B hhm , shot ≡ n h h P dhm ( k ) + 2 cyc i ; (63) b ¯ B hmm , shot ≡ n m h P dhm ( k ) + 2 cyc i , (64) k [ h Mpc − ] b ± σ b b ± σ b χ B hmm ± ± B SChmm ± ± Q hmm ± ± Q SChmm ± ± B hmm ± ± B SChmm ± ± Q hmm ± ± Q SChmm ± ± B hmm ± ± B SChmm ± ± Q hmm ± ± Q SChmm ± ± B hmm ± ± B SChmm ± ± Q hmm ± ± Q SChmm ± ± Table 5.
Best-fit bias parameters from halo-mass-mass bispec-tra and reduced bispectra. Column 1: bispectrum triangle scale;column 2: statistic, where B hmm and Q hmm are raw and where B SChmm and Q SChmm are shot-noise corrected; columns 3 and 4: best-fit b and b along with 1 σ errors; column 5: χ . k [ h Mpc − ] b ± σ b b ± σ b χ B hhm ± ± B SChhm ± ± Q hhm ± ± Q SChhm ± ± B hhm ± ± B SChhm ± ± Q hhm ± ± Q SChhm ± ± B hhm ± ± B SChhm ± ± Q hhm ± ± Q SChhm ± ± B hhm ± ± B SChhm ± ± Q hhm ± ± Q SChhm ± ± Table 6.
Best-fit bias parameters from halo-halo-mass bispectraand reduced bispectra. Column 1: bispectrum triangle scale; col-umn 2: statistic, where B hhm and Q hhm are raw and where B SChhm and Q SChhm are shot-noise corrected; columns 3 and 4: best-fit b and b along with 1 σ errors; column 5: χ . where n m = N/V µ and n h = N h /V µ , are the number den-sity of matter particles and haloes, respectively. For the re-duced bispectra we must correct the estimates of P P hhm and
P P hmm , which are written in the following form: Q denomhhm , shot = 23 n h h b P dhm ( k ) + 2 cyc i ; (65) Q denomhmm , shot = 23 n m h b P dhm ( k ) + 2 cyc i . (66)As it is the case that n m ≫ n h , we expect that the shot-noise corrections to B hmm will be significantly smaller than c (cid:13) , 000–000 odelling halo bias using the bispectrum Figure 6.
Evolution of the 95% likelihood contours for b and b obtained from the halo auto- and cross-bispectra and reduced bispectraas a function of scale. In each panel, the solid red lines of increasing thickness denote { B hmm , B hhm , B hhh } and the dashed blue linesof increasing thickness denote { Q hmm , Q hhm , Q hhh } . The alphabetical labels { a , b , c , d } correspond to the triangle configurations with k = { . , . , . , . } h Mpc − , respectively. The vertical black lines denote the effective bias parameter b NL , SChh , using the samewavemodes that enter into the bispectrum estimates. for B hhm . Hence, we shall think of this as being an almostperfect measure independent of discreteness.Using these estimators, we compute the ensemble av-erage and ensemble-to-ensemble variations of the halo-masscross-bispectra. We do this for the same bispectra configu-rations that were considered in § We estimate the nonlinear bias parameters and their errorsfrom the cross-bispectra using the same method employedfor the auto-bispectra in § § b and b , that are obtained from fitting theshot-noise corrected bispectra { B SChhh , B
SChhm , B
SChmm } and re-duced bispectra { Q SChhh , Q
SChhm , Q
SChmm } . The four panels showthe results obtained from fitting triangle configurations, k ∈ { . , . , . , . } h Mpc − , with k /k = 2, andthese correspond to the top-left, top-right, bottom left andbottom right panels, respectively. For comparative purposes,the vertical band in each panel, denotes the 1 σ constraint on b NL , SChh , obtained from the shot-noise corrected halo andnonlinear matter power spectra (c.f. § { B SChmm , B
SChhm , B
SChhh , } (solid red lines of increasing thickness), from the figureand the tables, we see that all of the results are reason-ably consistent with one another over the various scaleranges considered. However, when smaller scales are used(i.e. k > . h Mpc − ), the consistency weakens and thebest-fit parameters, obtained from B SChhh and B SChmm , differby & . σ .Evaluating the results for the reduced bispectra { Q SChmm , Q
SChhm , Q
SChhh } (dashed blue lines of increasing thick-ness), the four panels show a strong evolution of the er-ror ellipses with scale. We also note that the level of agree-ment between the different estimators also evolves strongly,becoming weaker and weaker as smaller scales are consid-ered. At the largest scale where k = 0 . h Mpc − , all of the2 σ likelihood contour regions overlap. However, this consis-tency is broken for the next scale range, k = 0 . h Mpc − ,where Q hmm and Q hhm are shifted downwards and furtherto the right favoring a more negative b and higher b . Thetrend continues in this same direction heading to smallerand smaller scales. c (cid:13) , 000–000 J. E. Pollack, R. E. Smith & C. Porciani
Comparing the results from both B and Q together, wesee that only on the largest scales is there any degree of over-all consistency. One way to interpret the results up to now,is that, if we believe b ≈ b NL , SChh , then the agent driving theinconsistency between the parameter estimates, is the break-down of the local bias model at tree-level. Furthermore, thebreakdown of the local tree-level model is more severe forthe reduced bispectrum than for the bispectrum.
This final set of analysis consists of a simple proof of methodtest, and we determine whether, when the underlying biasmodel is known, the ‘true’ bias parameters of the model areindeed recoverable with our approach.
For these tests, and for simplicity, we shall assume that thelocal model of biasing at quadratic order is the correct un-derlying bias model. Nonlinear biased density fields of thistype may be obtained through the following procedure.For each of the z = 0 outputs of the 40 simulations,we assign the nonlinear density field of matter to a cubicalFourier grid using the CIC algorithm. This is then Fouriertransformed. Each Fourier mode is then smoothed using aGaussian filter of scale R . We then inverse Fourier trans-form this field and obtain the smoothed, nonlinear matterdistribution in real space. Using this we next form the sum, δ b ( x | R ) = b b1 δ ( x | R ) + b b2 [ δ ( x | R )] / , (67)where b b1 and b b2 are the artificial bias parameters. Finally,this is Fourier transformed to give us δ b ( k | R ). Thus, given δ ( k | R ) and δ b ( k | R ), we can now use our standard bispec-trum estimators to estimate B bbb , B bbm , B bmm and B mmm .We refer to this procedure as the ‘biasing-by-hand’ test.The major benefits of these tests are that we are able tobetter gauge the effects to which nonlinearities beyond tree-level order influence the measured bispectra and reducedbispectra. We also note that shot-noise plays no rˆole here,since the biased field is created from the matter density fieldwhich is densely sampled. In order to interpret the results from such a constructionwe may use the results presented in Appendix B, with thesmall modification that we do not de-smooth the results. Ifwe define the smoothed bispectra as B ( k , k , k ) ≡ W ( k R ) W ( k R ) W ( k R ) B ( k , k , k ) , (68)then for {B bmm , B bbm , B bbb } , we have: B bmm = b B mmm + b P , m ; (69) B bbm = b B mmm + b b P , m + b P , m ; (70) B bbb = b B mmm + b b P , m + b b P , m + b P , m , (71) where for ease of notation we take b b i = b i and in the abovewe have suppressed the dependence of B , P , m , P , m and P , m , on ( k , k , − k − k ). We have also introduced theauxiliary functions: P n, m ≡ W ( k R ) W ( k R ) W ( k R ) P n, m ; (72) P , m ≡ Z d q (2 π ) f W q , k − q × T ( q , k − q , k , k ) + 2 cyc ; (73) P , m ≡ Z d q (2 π ) d q (2 π ) f W q , k − q f W q , k − q × P , m ( q , k − q , q , k − q , k )+2 cyc ; (74) P , m ≡ Z d q (2 π ) . . . d q (2 π ) f W q , k − q . . . f W q , k − q × P ( q , k − q , q , k − q , q , k − q ) . (75)The attractive aspect of this test can now be understood: ifwe move the terms in Eqns (69), (70) and (71), which areproportional to B mmm from the right to the left-hand-side,then we may rewrite this system as the matrix equation: Y bmm Y bbm Y bbb = b / b b / b /
12 0 b b / b b / b / P , m P , m P , m , (76)where we defined Y bmm ≡ B bmm − b B mmm , etc. This equa-tion may be inverted to give, P , m P , m P , m = 1 b b − b b b b − b Y bmm Y bbm Y bbb . (77)Hence, if we specify b , b and measure the four bispectra B mmm , B bmm , B bbm and B bbb , then we can determine ex-actly P , m , P , m and P , m . Thus we have complete knowl-edge of all components of the nonlinear model at all orders inthe theory. The lowest order perturbation theory expansionsof these statistics are (c.f. Appendix B): B (0)bmm ≈ b B (0)mmm + b h P (0)mm ( k ) P (0)mm ( k ) + 2 cyc i ; (78) B (0)bbm ≈ b B (0)mmm + b b h P (0)mm ( k ) P (0)mm ( k ) + 2 cyc i ;(79) B (0)bbb ≈ b B (0)mmm + b b h P (0)mm ( k ) P (0)mm ( k ) + 2 cyc i , (80)where in the above, we defined P mm ( k ) ≡ W ( k | R ) P mm ( k ). Following the algorithm described in § R = { , , . } h − Mpc. In all cases we apply the same non-linear bias: b = 1 .
63 and b = − .
53. Whilst these valuesare somewhat arbitrary, they were selected to coincide withthe best-fit values to the scatter plot of δ h ( x | R ) vs. δ ( x | R ),smoothed at R ∼ h − Mpc, that we recorded in § B mmm , B bmm , B bbm and B bbb for triangle configurationswith k = 0 . h Mpc − , k /k = 2, over 20 angular bins.From these we use the method described above, to recoverthe higher-order terms: P , m , P , m and P , m .We now define three modelling cases of interest: c (cid:13) , 000–000 odelling halo bias using the bispectrum R [ h − Mpc] b ± σ b b ± σ b χ B bbb ± ± B bbm ± ± B bmm ± ± B bbb ± ± B bbm ± ± B bmm ± ± B bbb ± ± B bbm ± ± B bmm ± ± Table 7.
Constraints on b and b obtained using the ExactTrispectrum model described in the text. The actual inputbias parameters were b = 1 .
63 and b = − .
53. Column 1: thesmoothing scale of the biased density field; Column 2: measuredquantity; Column 3 and 4 best-fit values for b and b along with1 σ errors; Column 5: the median χ . R [ h − Mpc] b ± σ b b ± σ b χ B bbb ± ± B bbm ± ± B bmm ± ± B bbb ± ± B bbm ± ± B bmm ± ± B bbb ± ± B bbm ± ± B bmm ± ± Table 8.
Same as Table 7, but this time the χ analysis is forthe Tree-level model described in the text. • Case 1:
All Order : Eqns (69), (70) and (71) are usedto interpret the data. • Case 2:
Exact Trispectrum : Eqns (69), (70) and (71)are exact up to P , m . All higher-order terms ( P , m , P , m )are dropped from the modelling. • Case 3:
Tree-level : lowest order expansions given byEqns (78), (79) and (80) are used to interpret the data.For each of the models described above, we then apply thesame χ –fitting analysis, as described in § b and b parameters.We begin by first examining the All Order expansionmodel. We confirm that for this case, the true bias parame-ters b = 1 .
63 and b = − .
53 are recovered exactly, albeitwith some uncertainty, however, with a χ = 0, and for allthe smoothing lengths considered. This null test is impor-tant, because it gives us confidence that any departures ofthe fits from the true bias values, can be attributed solelyto a break down of the theoretical modelling.Next we focus on the Exact Trispectrum modelwhere B mmm and P , m are measured from the simulations. InTable 7 we report the best-fitting bias parameters with the1 σ errors expressed in quadrature for the auto- and cross-bispectrum and for the four smoothing scales examined. Forthe case B bmm , a quick inspection of Eq. (69) tells us thatthe modelling should be exact, and indeed we see that thebias parameters are correctly recovered. However, for thecases B bbm and B bbb we see that the absence of the higher-order terms ( P , m , P , m ), induce biases in the parameters. For b the deviation from the true value is relatively small,with the value of the parameter only slightly decreasing insize. For b the deviations are larger, and this parameter be-comes more positive. We also note that the deviations fromthe true values appear to increase as the smoothing scale isdecreased.Finally, we focus on the Tree-level model. Table 8presents the best-fit results for b and b . We see that innearly all cases, there are systematic biases in the recoveryof the nonlinear bias parameters for all of the measured bis-pectra. In particular, for the case of B bmm , the results aremost deviant and poorly constrained. Whereas for B bbb , onlywhen the data has been smoothed on scales R = 20 h − Mpcare the recovered parameters close to the true values.The comparison of the results from this analysis leadsus to conclude that the recovered bias parameters are verysensitive to the inclusion of beyond leading order correctionsin the modelling. Furthermore, accurate nonlinear modellingof, at the very least, the matter bispectrum and trispectrumwill be essential, if we are to safely recover the nonlinearbias parameters from this approach.
We have evaluated the local halo bias model at second-orderusing three different probes: smoothed density fields; powerspectra; and bispectra and reduced bispectra. A summary ofour results for the best-fitting bias parameters determinedfrom shot-noise corrected spectra is shown in Figure 7.In the figure we also compare our estimates for b and b with the analytical predictions for the halo bias parame-ters obtained from the peak background split (PBS) ansatz(Bardeen et al. 1986; Mo & White 1996). The average the-ory bias parameters are obtained through computing theexpression: b i = 1 n Z ∞ M min dM n ( M ) b i ( M ) ; n ≡ Z ∞ M min dM n ( M ) , (81)where n ( M ) is the halo mass function, and M min is set equal to the value of the minimum halo massidentified in the simulations (see §
3) . We evaluate theabove integral using three different fits to N-body simula-tions by Sheth & Tormen (1999), Warren et al. (2006) andPillepich et al. (2010). The corresponding expressions forthe bias parameters as a function of halo mass are presentedin Scoccimarro et al. (2001) and Manera et al. (2010).Considering, the results for b (bottom panel), we seethat when the reduced bispectra, Q hhh , Q hhm and Q hmm are used, the recovered parameters are poorly constrainedand appear incompatible with respect to the other estimatesand are only weakly consistent with one another. On theother hand, the estimates from the bispectra B hhh , B hhm and B hmm are in much better agreement with each other.They are also in close agreement with the predictions fromWarren et al. (2006) and Pillepich et al. (2010), which bothprovided an estimate of b = 1 .
39. However, they slightlyundershoot the values from the effective bias estimates, b NLhh and b Lhh , likewise the smoothed density fields, and theSheth-Tormen prediction. The analytical predictions fromthe Sheth & Tormen (1999) mass-function yielded b = 1 . c (cid:13)000
39. However, they slightlyundershoot the values from the effective bias estimates, b NLhh and b Lhh , likewise the smoothed density fields, and theSheth-Tormen prediction. The analytical predictions fromthe Sheth & Tormen (1999) mass-function yielded b = 1 . c (cid:13)000 , 000–000 J. E. Pollack, R. E. Smith & C. Porciani
Figure 7.
Plot summary of bias measurements on b and b for the second-order local bias model from different estimators:,shot-noise corrected B , Q , P , and finally smoothed, δ R , in com-parison with analytical predictions applying the peak-backgroundsplit ansatz with the Sheth & Tormen (1999) mass function de-noted by the dotted line, as well as the Warren et al. (2006) andPillepich et al. (2010) mass-functions, which are both representedby the single dashed-line. results smoothed on a scale R ∼ h − Mpc. The recoveredvalues of b from the effective bias in the power spectraand the smoothed density fields collectively are in broadagreement, but the latter increase with decreasing smooth-ing scale.In the top panel of the figure, we see that the con-straints on b from the different estimators used in thesimulations are reasonably consistent with one another, al-beit with significant error bars. These estimates also agreewell with the PBS prediction from the Warren et al. (2006)and Pillepich et al. (2010) mass functions, which give anaverage b = − .
24. However, the prediction from theSheth & Tormen (1999) mass function gives b = − . δ h vs. δ . Perhaps, this is a consequence ofnon-locality. In this case we need a more advanced theoret-ical approach to understand the halo clustering. One pos-sibility may be that the bias is local in Lagrangian space(Catelan et al. 2000; Matsubara 2011).Secondly, our simple biasing-by-hand test has enabledus to discern that, for the current set of tests that we haveperformed, the most likely explanation at this point is thatthe tree-level expansions for B hhh and Q hhh are not suffi-ciently accurate enough. Higher-order nonlinear corrections in the modelling must be included, and if possible all orderexpansions for B mmm and P , m would be invaluable.Thirdly, as we have argued, the local bias model onlymakes sense in the context of smoothing. The bias param-eters one recovers from fitting, depend sensitively on thesmoothing scale R . For the biasing-by-hand tests, the exactsmoothing scale was known beforehand. However, in realdata we do not know this a priori . In all of the cases, whenrecovering bias parameters from the bispectra, we have as-sumed that we are on sufficiently large scales such that W ( k i R ) →
1. However, in general R should be a free pa-rameter and as such, marginalized over in the analysis.Manera & Gazta˜naga (2011) performed a similar studyof nonlinear halo bias with the three-point correlation func-tion in configuration space. In contrast to our analysis, theymeasured the bias parameters for different halo mass bins.They found inconsistencies between the predictions of thedifferent estimators considered. They evaluated the scatterplots of δ h vs. δ as function of smoothing scale, and foundthat stability in the local bias parameters ( b , b ) occurredfor smoothing scales, R > − h − Mpc, albeit with largererrors. They also found that the bias predictions derivedfrom δ h vs. δ for R = 60 h − Mpc were in good agreement,to within the errors, with the linear bias measured from eval-uating the two-point correlation function on large and inter-mediate scales. As in the case of our findings, they found thelinear bias measured from evaluating the three-point corre-lation function, expressed in terms of the Q -amplitudes, wasnot consistent with that of the two-point correlation func-tion for the lower mass bins: M < h − M ⊙ . They wereunable to formulate solid conclusions for larger mass ranges.Guo & Jing (2009) explored the differences between es-timates of bias from Q and P . They also found that b basedon analysis of the mock galaxy catalogues was larger forgalaxy reduced bispectra and power spectra Q g than for P g . While Guo & Jing (2009) also noted that this might bedue to the failure of SPT at tree-level, they also reportedthat agreement could be found between estimators if Q mmm measured directly from the simulations was used in place ofthe tree-level expression. However, when we performed thesame test with our data we found no dramatic reconcilia-tion of the two bias estimates. The investigation performedby Guo & Jing (2009) was carried out using only 4 largevolume and 3 smaller volume runs. As a result of having toofew realizations, they assumed the Gaussian approximationfor the covariance matrix in order to perform their study atlarge-scales. In this paper we have used a sample of 40 large volume N -body simulations, with total volume V ∼
135 Gpc h − ,to test the local model of halo biasing, and the extent towhich nonlinearities impact the modelling. We used threedifferent methods for exploring the bias: smoothed densityfields; power spectra; and bispectra and reduced bispectra.We focused mainly on the results from the bispectra. All ofthe reported results were scaled to a single realization of oursimulations, and so are directly relevant for galaxy surveyswith a total volume V ∼ h − Gpc .In § c (cid:13) , 000–000 odelling halo bias using the bispectrum ory and how they connect to density statistics. We thenreviewed the local model of halo biasing, drawing special at-tention to the rˆole that smoothing plays in the theory. Theimportant result being that even at tree-level, the smooth-ing explicitly enters the theory. The expressions for B hhh and Q hhh , which are typically used in all past and currentanalysis, make the assumption that smoothing is unimpor-tant. In subsequent sections we argued that this assumptionis not safe.In § N -body simulations, andthe halo catalogues used in this study.In § δ h ( x | R ) and δ ( x | R ) smoothed on the set of scales R = { , , } h − Mpc. To this data we fitted the localmodel of halo biasing up to second order, including b , b and b . We found that the fits were reasonably good, however thebest-fit parameters showed a running with the filter scale R .We then demonstrated, theoretically, why the nonlinear biasparameters from this approach, could not be made indepen-dent of smoothing scale. We then turned to Fourier spacestatistics, and used the halo auto- and cross-power spec-tra to obtain an effective large-scale bias. We found thatthe effective bias estimators b NLhh and b NLhm were reasonablyscale independent for k < . h Mpc − . However, on scalessmaller than this, b NLhh decreased with increasing wavenum-ber, whereas b NLhm remained surprisingly flat.In § k = { . , . , . , . } h Mpc − and with k /k = 2and θ ∈ [0 , π ]. These triangles all lay in the weakly nonlin-ear regime k = 0 . − . h Mpc − . We modeled these es-timates using tree-level perturbation theory expressions forthe matter bispectrum and nonlinear bias at second order,and assumed smoothing to be unimportant. Our method forestimation of the bias parameters followed a standard min-imum χ approach. We estimated the covariance matrix forthe full ensemble applying principal component analysis tominimize the intrinsic noise. We also performed a jackknifesubsampling routine to propagate the error of the estimatedcovariance matrix onto the errors of the bias parameters.We tested how well the measurements of the matterbispectra B mmm and reduced bispectra Q mmm could be de-scribed by such modelling. The results obtained for the biasparameters b and b showed that the tree-level expressionswere a good description of the data for configurations, k = { . , . } h Mpc − , for which b = 1 and b = 0. However,for smaller scale triangles, k = { . , . } h Mpc − , signif-icant deviations were apparent, and these were manifest as b = 1 and b = 0 at high significance.We then applied the χ test to the halo bispectra B hhh and reduced bispectra Q hhh . We found, for the shot-noisecorrected B hhh , that the estimated values for b ∼ .
40 and b ∼ − .
25 were reasonably consistent with one another.However, the fits became progressively poorer as smallerscales were added, yet the reduced– χ remained . k = { . } h Mpc − . For the shot-noise corrected Q hhh ,we found that the values of b were significantly larger b ∼ .
85, with large errors, and the values evolved withtriangle configuration scale. However, b ∼ − .
3, appearedto be more stable, although again with large errors. For tri-angle configurations with k > . h Mpc − , the fits from B hhh and Q hhh were inconsistent with each other at the ∼ σ level. For both B hhh and Q hhh shot-noise corrections signif-icantly influenced the recovered bias parameters.In § B hhm and B hmm , and reduced bispectra Q hhm and Q hmm .We calculated the tree-level expressions for these quantitiessymmetrized in all of their arguments. We then developedestimators for them. We showed that for B hmm and Q hmm ,provided the matter distribution was densely sampled, theshot-noise corrections were small.We applied the χ analysis from § b and b . We foundthat for B hhm the shot-noise corrected data were all rea-sonably consistent with one another, giving b ∼ .
39 and b ∼ − .
3. For B hmm we found a similar pattern, exceptthat for k > . h Mpc − where we found b ∼ .
0, butwith large errors. The results for Q hhm and Q hmm appearedto vary significantly.Finally in § b and b to some fiducial values, and then measuredthe smoothed matter and halo bispectra and their cross-bispectra, then the higher order matter correlators P ,m , P ,m and P ,m could be recovered exactly. Thus we wereable to construct three models: an all order model; a modelthat used the exact matter bispectrum and trispectrum; anda tree-level model.We applied the χ analysis using these three models andfor bispectra with k = 0 . h Mpc − . As expected, the ex-act model recovered the correct bias parameters. The modelwith the exact B mmm and P ,m , was in fact also exact for B bmm . For B bbm and B bbb the recovered parameters wereclose to the true values, but showed evolution with smooth-ing scale. Finally, for the tree-level model we showed thatthere was a significant evolution in the estimated bias pa-rameters with smoothing scale and with the type of statisticused.We conclude that estimates of nonlinear bias from thebispectrum that do not attempt to account for higher-ordercorrections, will most likely provide biased estimates for thebias parameters b and b . Robust modelling of nonlinearbias from bispectra, will, at the very least, require almostexact models for the matter bispectrum and trispectrum.Real space estimates of bias appear to be inconsistentwith Fourier space based ones. We believe that this owesprimarily to the mixing of large- and small-scale wavemodesin real space. We therefore recommend that perturbativemethods should strictly be applied in Fourier space. We alsorecommend that measurements focus on the bispectrum andthe associated cross statistics, rather than the reduced bis-pectra, since this appears very sensitive to nonlinearities inthe modelling and also shot-noise corrections.Finally, we emphasize the importance of smoothing inthe local model. Owing to the fact that the smoothing scaleassociated with the halo/galaxy distribution in question isnot known a priori , it must be treated as a nuisance param-eter and so marginalized over.An alternative strategy for recovering information from c (cid:13) , 000–000 J. E. Pollack, R. E. Smith & C. Porciani higher order statistics, which may be of interest for futureconsideration, is the use of ‘Gaussianizing transformations’or ‘density clipping’ (Neyrinck et al. 2009; Seo et al. 2011;Simpson et al. 2011). However, the theoretical connectionbetween what is measured and what is interpreted from suchapproaches still remains to be fully calculated.
ACKNOWLEDGEMENTS
We thank the anonymous referee for helpful suggestions. Wealso thank Tobias Baldauf, Martin Crocce, Roman Scocci-marro, Emiliano Sefussati, Ravi Sheth and Masahiro Takadafor useful discussions. We thank V. Springel for makingpublic
GADGET-2 and for providing his
B-FoF halo finder,and R. Scoccimarro for making public his code. JEPand CP were supported by funding provided through theSFB-Transregio 33 The Dark Universe by the DeutscheForschungsgemeinschaft. RES acknowledges support from aMarie Curie Reintegration Grant, the Alexander von Hum-boldt Foundation and partial support from the Swiss Na-tional Foundation under contract 200021-116696/1.
REFERENCES
Baldauf T., Seljak U., Senatore L., 2011, Journal of Cos-mology and Astro-Particle Physics, 4, 6Bardeen J. M., Bond J. R., Kaiser N., Szalay A. S., 1986,ApJ, 304, 15Bernardeau F., Colombi S., Gazta˜naga E., Scoccimarro R.,2002, Phys. Rep. , 367, 1Bouchet F. R., Colombi S., Hivon E., Juszkiewicz R., 1995,A&A, 296, 575Catelan P., Porciani C., Kamionkowski M., 2000, MNRAS,318, L39Coles P., 1993, MNRAS, 262, 1065Crocce M., Pueblas S., Scoccimarro R., 2006, MNRAS, 373,369Dalal N., Dor´e O., Huterer D., Shirokov A., 2008, PRD,77, 123514Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985,ApJ, 292, 371Dekel A., Lahav O., 1999, ApJ, 520, 24Dekel A., Rees M. J., 1987, Nature, 326, 455Desjacques V., Seljak U., Iliev I. T., 2009, MNRAS, 396,85Feldman H. A., Frieman J. A., Fry J. N., Scoccimarro R.,2001, Physical Review Letters, 86, 1434Fry J. N., 1994, ApJ, 421, 21Fry J. N., Gaztanaga E., 1993, ApJ, 413, 447Fry J. N., Scherrer R. J., 1994, ApJ, 429, 36Gazta˜naga E., Norberg P., Baugh C. M., Croton D. J.,2005, MNRAS, 364, 620Gazta˜naga E., Scoccimarro R., 2005, MNRAS, 361, 824Goroff M. H., Grinstein B., Rey S.-J., Wise M. B., 1986,ApJ, 311, 6Guo H., Jing Y. P., 2009, ApJ, 702, 425Guzik J., Bernstein G., Smith R. E., 2007, MNRAS, 375,1329Hartlap J., Simon P., Schneider P., 2007, A&A, 464, 399 Heavens A. F., Matarrese S., Verde L., 1998, MNRAS, 301,797Hivon E., Bouchet F. R., Colombi S., Juszkiewicz R., 1995,A&A, 298, 643Hockney R. W., Eastwood J. W., 1988, Computer simula-tion using particles. Bristol: Hilger, 1988Jain B., Bertschinger E., 1994, ApJ, 431, 495Jing Y. P., 2005, ApJ, 620, 559Joachimi B., Shi X., Schneider P., 2009, A&A, 508, 1193Juszkiewicz R., 1981, MNRAS, 197, 931Kaiser N., 1984, ApJL, 284, L9Komatsu E., The WMAP Team 2010, ArXiv e-printsKomatsu E., Dunkley J., Nolta M. R., Bennett C. L., GoldB., Hinshaw G., Jarosik N., Larson D., Limon M., Page L.,Spergel D. N., Halpern M., Hill R. S., Kogut A., MeyerS. S., Tucker G. S., Weiland J. L., Wollack E., WrightE. L., 2009, ApJS, 180, 330Makino N., Sasaki M., Suto Y., 1992, PRD, 46, 585Manera M., Gazta˜naga E., 2011, MNRAS, 415, 383Manera M., Sheth R. K., Scoccimarro R., 2010, MNRAS,402, 589Marin F., 2010, ArXiv e-printsMatarrese S., Verde L., Heavens A. F., 1997, MNRAS, 290,651Matsubara T., 2011, PRD, 83, 083518McBride C. K., Connolly A. J., Gardner J. P., ScrantonR., Newman J. A., Scoccimarro R., Zehavi I., SchneiderD. P., 2011, ApJ, 726, 13Mo H. J., Jing Y. P., White S. D. M., 1997, MNRAS, 284,189Mo H. J., White S. D. M., 1996, MNRAS, 282, 347Neyrinck M. C., Szapudi I., Szalay A. S., 2009, ApJL, 698,L90Nishimichi T., Ohmuro H., Nakamichi M., Taruya A., Ya-hata K., Shirata A., Saito S., Nomura H., Yamamoto K.,Suto Y., 2007, Publications of Astronomical Society ofJapan, 59, 1049Nishimichi T., Taruya A., Koyama K., Sabiu C., 2010,Journal of Cosmology and Astro-Particle Physics, 7, 2Norberg P., Baugh C. M., Gazta˜naga E., Croton D. J.,2009, MNRAS, 396, 19Padmanabhan N., White M., Norberg P., Porciani C., 2009,MNRAS, 397, 1862Peebles P. J. E., 1980, The large-scale structure of theuniverse. Research supported by the National ScienceFoundation. Princeton, N.J., Princeton University Press,1980. 435 p.Peebles P. J. E., Groth E. J., 1975, ApJ, 196, 1Pillepich A., Porciani C., Hahn O., 2010, MNRAS, 402, 191Porciani C., Giavalisco M., 2002, ApJ, 565, 24Roth N., Porciani C., 2011, MNRAS, 415, 829Scoccimarro R., 1998, MNRAS, 299, 1097Scoccimarro R., 2000, ApJ, 544, 597Scoccimarro R., Colombi S., Fry J. N., Frieman J. A., HivonE., Melott A., 1998, ApJ, 496, 586Scoccimarro R., Couchman H. M. P., Frieman J. A., 1999,ApJ, 517, 531Scoccimarro R., Feldman H. A., Fry J. N., Frieman J. A.,2001, ApJ, 546, 652Scoccimarro R., Sheth R. K., Hui L., Jain B., 2001, ApJ,546, 20Sefusatti E., 2009, PRD, 80, 123002 c (cid:13) , 000–000 odelling halo bias using the bispectrum Sefusatti E., Crocce M., Pueblas S., Scoccimarro R., 2006,PRD, 74, 023522Sefusatti E., Komatsu E., 2007, PRD, 76, 083004Seljak U., Sugiyama N., White M., Zaldarriaga M., 2003,PRD, 68, 083507Seljak U., Zaldarriaga M., 1996, ApJ, 469, 437Seo H.-J., Sato M., Dodelson S., Jain B., Takada M., 2011,ApJL, 729, L11+Sheth R. K., Lemson G., 1999, MNRAS, 304, 767Sheth R. K., Tormen G., 1999, MNRAS, 308, 119Simpson F., Berian James J., Heavens A. F., Heymans C.,2011, ArXiv e-printsSmith R. E., 2009, MNRAS, pp 1337–+Smith R. E., Hern´andez-Monteagudo C., Seljak U., 2009,PRD, 80, 063528Smith R. E., Marian L., 2011, MNRAS, p. 1484Smith R. E., Scoccimarro R., Sheth R. K., 2007, PRD, 75,063512Smith R. E., Sheth R. K., Scoccimarro R., 2008, PRD, 78,023523Springel V., 2005, MNRAS, 364, 1105Verde L., Heavens A. F., Percival W. J., Matarrese S., The2dFGRS Team 2002, MNRAS, 335, 432Vishniac E. T., 1983, MNRAS, 203, 345Warren M. S., Abazajian K., Holz D. E., Teodoro L., 2006,ApJ, 646, 881White S. D. M., Rees M. J., 1978, MNRAS, 183, 341 c (cid:13) , 000–000 J. E. Pollack, R. E. Smith & C. Porciani
APPENDIX A: BISPECTRUM ESTIMATION: ALGORITHM
Briefly, the algorithm that we employ is as follows: firstly the dark matter density field is computed by assigning the darkmatter particles to a cubical grid using the ‘Cloud in Cell’ (CIC) technique (Hockney & Eastwood 1988). Next the fastFourier transform (FFT) of the gridded density field is computed. Each Fourier mode is then corrected for convolution withthe Fourier mesh. We do this by dividing out from each mode, the Fourier transform of the window assignment function ofthe CIC scheme (Hockney & Eastwood 1988; Jing 2005): δ ( k ) = δ g ( k ) W CIC ( k ) ; W CIC ( k ) ≡ Y i =1 , (cid:20) sin( πk i / k Ny ) πk i / k Ny (cid:21) , (A1)where subscript g denotes gridded quantities, k Ny = πN g /L is the Nyquist frequency of the mesh and N g is the number ofFourier grid cells.The estimator for the bispectrum can be written (Scoccimarro et al. 1998): b B ( k , k , θ ) = V µ V B ( k , k , θ ) Z d q (2 π ) d q (2 π ) d q (2 π ) (2 π ) δ D ( q ) δ ( q ) δ ( q ) δ ( q ) , (A2)where V µ is the sample volume (in our case the simulation volume), the normalization factor, V B , can be written as(Sefusatti et al. 2006; Joachimi et al. 2009) V B ( k , k , µ ) ≡ Z d q (2 π ) d q (2 π ) d q (2 π ) (2 π ) δ D ( q ) ≈ π k k k (2 π ) (∆ k ) , (A3)and we write in shorthand δ D ( q ...n ) ≡ δ D ( q + . . . + q n ). A practical implementation of the above estimator may be achievedthrough (Smith et al. 2008): b B d ( k , k , θ ) = V µ N tri ( k , k , θ ) N tri ( k ,k ,θ ) X ( n , n ) R e [ δ ( k n ) δ ( k n ) δ ( k − n − n )] , (A4)where superscript d denotes discretized quantities; n i denotes an integer vector from the origin of the k -space to each meshpoint; ( n , n ) represents a pair of integer vectors, which lie in thin shells centred on k and k and whose angular separationlies in a narrow angular bin centred on θ , and for which k = − k − k . The upper limit of the sum N tri( k , k , θ )represents the total number of triangles that have such a configuration.The estimator for the bin-averaged reduced bispectrum, b Q , is written as: b Q d = b B d / b Q denom , d , (A5)where b Q denom , d is the estimator for the bin-averaged cyclical terms of the power spectrum generated from first computing thebin-averaged power-spectra, b P d i . Note that we estimate the power spectra that enter into this product in a slightly differentway than normal: we use only those modes that go into estimating the particular B triangle configuration to estimate b Q denom , d ( k , k , θ ). Hence, b P d i = V µ N tri ( k , k , θ ) N tri ( k ,k ,θ ) X ( n , n ) | δ ( k n i ) | , (A6)where i ∈ { , , } and where b P d3 is dependent on the angular bin, since, with | n | and | n | fixed, we still havecos θ = n · n / | n || n | and the closure criterion implies n = − n − n varies as function of θ . Therefore, b Q denom , d = b P d1 b P d2 + b P d2 b P d3 + b P d1 b P d3 . (A7)The estimates of B d and Q d are then corrected for discreteness, i.e. shot-noise. For the estimators of interest, thecorrections are (Peebles 1980; Smith et al. 2008): b P shot ≡ /n ; (A8) b B shot ≡ [ b P d1 + b P d2 + b P d3 ] /n − /n ; (A9) b Q denomshot ≡ b P d1 + b P d2 + b P d3 ] /n − /n . (A10)Shot-noise corrected estimates of the statistics are obtained: χ = χ d − χ shot , (A11)where χ ∈ { b P , b B, b Q denom } and where b Q = b B/ b Q denom . Note that the above recipe corrects some typos that are present inSmith et al. (2008). c (cid:13) , 000–000 odelling halo bias using the bispectrum APPENDIX B: HALO CROSS-BISPECTRA IN THE LOCAL MODEL
As was shown in Eq. (15), at quadratic order, the local model of nonlinear biasing can be written as: δ h ( k | R ) = b ( M ) δ ( k | R ) + b ( M )2 Z d q (2 π ) d q (2 π ) δ ( q | R ) δ ( q | R )(2 π ) δ D ( k − q − q ) . (B1)where δ ( q i | R ) ≡ δ ( q i ) W ( q i R ), is the filtered density. Using this model we may now proceed to calculate the halo auto- andhalo-mass cross-bispectra. B1 Halo-mass-mass bispectrum in the local model
Let us start with the simplest three-point cross-statistic, the halo-matter-matter bispectrum, this can be written: h δ h ( k | M, R ) δ ( k | R ) δ ( k | R ) i = b ( M ) h δ ( k | R ) δ ( k | R ) δ ( k | R ) i + b ( M )2 Z d q (2 π ) d q (2 π ) (2 π ) δ D ( k − q − q ) h δ ( q | R ) δ ( q | R ) δ ( k | R ) δ ( k | R ) i . (B2)Let us define the smoothed n -point spectra as: h δ ( k | R ) . . . δ ( k n | R ) i ≡ (2 π ) δ ( k + . . . + k n ) e P n ( k , . . . , k n | R ) , = (2 π ) δ ( k + . . . + k m ) W ( k R ) . . . W ( k n R ) P n ( k , . . . , k n ) (B3)where e P ≡ e P = W ( kR ) P , e P ≡ e B = W ( k R ) W ( k R ) W ( k R ) B , and where e P ≡ e T = W ( k R ) W ( k R ) W ( k R ) W ( k R ) T .We may now integrate over q to obtain h δ h ( k | M, R ) δ ( k | R ) δ ( k | R ) i = (2 π ) δ D ( k + k + k ) (cid:20) b ( M ) e B ( k , k , k ) + b ( M )2 Z d q (2 π ) e T ( q , k − q , k , k ) (cid:21) . (B4)On dividing the above expression by W ( k R ) W ( k R ) W ( k R ), then we find the halo-mass-mass bispectrum can be written as B hmm ( k , k , k ) = b ( M ) B mmm ( k , k , k ) + b ( M )2 Z d q (2 π ) f W q , k − q T mmmm ( q , k − q , k , k ) . (B5)We may symmetrize the above result by constructing the sum [ B hmm + B mhm + B mmh ] /
3, and this gives us: B hmm ( k , k , k ) = b ( M ) B mmm ( k , k , k ) + b ( M )6 Z d q (2 π ) hf W q , k − q T mmmm ( q , k − q , k , k ) + 2 cyc i (B6)On expanding B and T to fourth order in δ , the above expression can be approximated as, B (0)hmm ( k , k , k ) ≈ b ( M ) B (0)mmm ( k , k , k ) + b ( M )3 hf W k , k P (0) ( k ) P (0)mm ( k ) + 2 cyc i . (B7)Finally, in the large-scale limit k i →
0, or for arbitrarily small smoothing scales, k i R ≪
1, the above expression becomes, B (0)hmm ( k , k , k ) ≈ b ( M ) B (0)mmm ( k , k , k ) + b ( M )3 h P (0)mm ( k ) P (0)mm ( k ) + 2 cyc i (B8)The reduced bispectrum Q hmm is given by Q hmm ( k , k , k ) ≡ B hmm ( k , k , k ) P P hmm , (B9)where P P hmm = 23 [ P hm ( k ) P mm ( k ) + 2 cyc ] + 13 [ P hm ( k ) P hm ( k ) + 2 cyc ] . (B10)In order to calculate the reduced halo-mass cross-bispectrum, then we need to evaluate the halo-matter power spectrum.In the local model and up to quadratic order in the bias we have, P hm ( k ) = b ( M ) P mm ( k ) + b ( M )2 Z d q (2 π ) f W q , k − q B mmm ( q , k − q , − k ) . (B11)Using the above expression in Eq. (B10) we find P P hmm ( k , k , k ) = (cid:26) (cid:20) b ( M ) P mm ( k ) + b ( M )2 Z d q (2 π ) f W q , k − q B mmm ( q , k − q , − k ) (cid:21) P mm ( k ) + 2 cyc (cid:27) + 13 (cid:26)(cid:20) b ( M ) P mm ( k ) + b ( M )2 Z d q (2 π ) f W q , k − q B mmm ( q , k − q , − k ) (cid:21) × (cid:20) b ( M ) P mm ( k ) + b Z d q (2 π ) f W q , k − q B mmm ( q , k − q , − k ) (cid:21) + 2 cyc (cid:27) (B12)If we expand P P hmm to fourth order in the density then the above expression simplifies to: c (cid:13)000
1, the above expression becomes, B (0)hmm ( k , k , k ) ≈ b ( M ) B (0)mmm ( k , k , k ) + b ( M )3 h P (0)mm ( k ) P (0)mm ( k ) + 2 cyc i (B8)The reduced bispectrum Q hmm is given by Q hmm ( k , k , k ) ≡ B hmm ( k , k , k ) P P hmm , (B9)where P P hmm = 23 [ P hm ( k ) P mm ( k ) + 2 cyc ] + 13 [ P hm ( k ) P hm ( k ) + 2 cyc ] . (B10)In order to calculate the reduced halo-mass cross-bispectrum, then we need to evaluate the halo-matter power spectrum.In the local model and up to quadratic order in the bias we have, P hm ( k ) = b ( M ) P mm ( k ) + b ( M )2 Z d q (2 π ) f W q , k − q B mmm ( q , k − q , − k ) . (B11)Using the above expression in Eq. (B10) we find P P hmm ( k , k , k ) = (cid:26) (cid:20) b ( M ) P mm ( k ) + b ( M )2 Z d q (2 π ) f W q , k − q B mmm ( q , k − q , − k ) (cid:21) P mm ( k ) + 2 cyc (cid:27) + 13 (cid:26)(cid:20) b ( M ) P mm ( k ) + b ( M )2 Z d q (2 π ) f W q , k − q B mmm ( q , k − q , − k ) (cid:21) × (cid:20) b ( M ) P mm ( k ) + b Z d q (2 π ) f W q , k − q B mmm ( q , k − q , − k ) (cid:21) + 2 cyc (cid:27) (B12)If we expand P P hmm to fourth order in the density then the above expression simplifies to: c (cid:13)000 , 000–000 J. E. Pollack, R. E. Smith & C. Porciani
P P (0)hmm ( k , k , k ) ≈ b ( M )3 [2 + b ( M )] h P (0)mm ( k ) P (0)mm ( k ) + 2 cyc i . (B13)Hence we have, Q (0)hmm ( k , k , k ) ≈
32 + b ( M ) Q (0)mmm ( k , k , k ) + b ( M )2 b ( M ) + b ( M ) " f W k , k P (0)mm ( k ) P (0)mm ( k ) + 2 cyc P (0)mm ( k ) P (0)mm ( k ) + 2 cyc (B14)Finally, in the limit that k i R ≪ f W → Q (0)hmm ( k , k , k ) ≈
32 + b ( M ) Q (0)mmm ( k , k , k ) + b ( M )2 b ( M ) + b ( M ) . (B15) B2 Halo-halo-mass bispectrum in the local model
Again, using Eq. (B1), the halo-halo-mass bispectrum, symmetrized in the k i arguments, can be written: B hhm ( k , k , k ) = b ( M ) B mmm ( k , k , k ) + b ( M ) b ( M )3 Z d q (2 π ) hf W q , k − q T mmmm ( q , k − q , k , k ) + 2 cyc i + b ( M )12 Z d q (2 π ) d q (2 π ) hf W q , k − q f W q , k − q P , m ( q , k − q , q , k − q , k ) + 2 cyc i . (B16)If we use perturbation theory to expand P , B , T and P , and only keep terms that are 4 th order in the density field, then theabove expression can be approximated by: B (0)hhm ( k , k , k ) ≈ b ( M ) B (0)mmm ( k , k , k ) + 13 b ( M ) b ( M ) hf W k , k P (0)mm ( k ) P (0)mm ( k ) + 2 cyc i . (B17)In the large-scale limit k i R →
0, the above expression can be approximated as: B (0)hhm ( k , k , k ) ≈ b ( M ) B (0)mmm ( k , k , k ) + 13 b ( M ) b ( M ) h P (0)mm ( k ) P (0)mm ( k ) + 2 cyc i . (B18)The reduced bispectrum is given by, Q hhm ( k , k , k ) ≡ B hhm ( k , k , k ) P P hhm , (B19)where P P hhm = 23 [ P hh ( k ) P hm ( k ) + 2 cyc ] + 13 [ P hm ( k ) P hm ( k ) + 2 cyc ] . (B20)The halo-mass power spectrum is given by Eq. (B11) and the halo auto-power spectrum is given by: P hh ( k ) = b ( M ) P mm ( k ) + b ( M ) b ( M ) Z d q (2 π ) f W q , k − q B mmm ( q , k − q , − k )+ b ( M )4 Z d q (2 π ) d q (2 π ) f W q , k − q f W q , − k − q T mmmm ( q , k − q , q , − k − q ) . (B21)Expanding the above expression to lowest order in perturbation theory gives, P P (0)hhm ≈ b ( M )3 [2 b ( M ) + 1] h P (0)mm ( k ) P (0)mm ( k ) + 2 cyc i . (B22)Using the above expression, we find that the tree-level expression for the reduced bispectrum can be written: Q (0)hhm ( k , k , k ) ≈ b ( M ) + 1 Q (0)mmm ( k , k , k ) + 2 b ( M )2 b ( M ) + b ( M ) " f W k , k P (0)mm ( k ) P (0)mm ( k ) + 2 cyc P (0)mm ( k ) P (0)mm ( k ) + 2 cyc . (B23)In the large-scale limit k i R →
0, we again have f W → Q (0)hhm ( k , k , k ) ≈ b ( M ) + 1 Q (0)mmm ( k , k , k ) + 2 b ( M )2 b ( M ) + b ( M ) . (B24) B3 Halo-halo-halo bispectrum in the local model
Again using Eq. (B1), the halo-halo-halo bispectrum, symmetrized in the k i arguments, can be written: B hhh ( k , k , k ) = b ( M ) B mmm ( k , k , k ) + 12 b ( M ) b ( M ) Z d q (2 π ) hf W q , k − q T mmmm ( q , k − q , k , k ) + 2 cyc i + 14 b b Z d q (2 π ) d q (2 π ) hf W q , k − q f W q , k − q P , m ( q , k − q , q , k − q , k ) + 2 cyc i c (cid:13) , 000–000 odelling halo bias using the bispectrum + b Z d q (2 π ) d q (2 π ) d q (2 π ) f W q , k − q f W q , k − q f W q , k − q P , m ( q , k − q , q , k − q , q , k − q ) . (B25)If we use perturbation theory to expand P , B , T , P , and P , and keep only terms that are fourth order in the density field,then the above expression can be approximated by: B (0)hhh ( k , k , k ) ≈ b ( M ) B (0)mmm ( k , k , k ) + b ( M ) b ( M ) hf W k , k P (0)mm ( k ) P (0)mm ( k ) + 2 cyc i . (B26)In the large-scale limit k i R →
0, the above expression can be approximated as: B (0)hhh ( k , k , k ) ≈ b ( M ) B (0)mm ( k , k , k ) + b ( M ) b ( M ) h P (0)mm ( k ) P (0)mm ( k ) + 2 cyc i . (B27)The reduced halo-halo-halo bispectrum is given by, Q hhh ( k , k , k ) ≡ B hhh ( k , k , k ) P P hhh , (B28)where P P hhh = [ P hh ( k ) P hh ( k ) + 2 cyc ] , (B29)where the halo auto-power spectrum is given by Eq. (B21). On expanding the above expression to fourth order in the densitywe find, P P (0)hhh ≈ b ( M ) h P (0)mm ( k ) P (0)mm ( k ) + 2 cyc i . (B30)Using the above expression, we find that the tree-level expression for the reduced bispectrum can be written: Q (0)hhh ( k , k , k ) ≈ b ( M ) Q (0)mmm ( k , k , k ) + b ( M ) b ( M ) " f W k , k P (0)mm ( k ) P (0)mm ( k ) + 2 cyc P (0)mm ( k ) P (0)mm ( k ) + 2 cyc . (B31)In the large-scale limit k i R →
0, we again have f W → Q (0)hhh ( k , k , k ) ≈ b ( M ) Q (0)mmm ( k , k , k ) + b ( M ) b ( M ) . (B32) c (cid:13)000