Signature of primordial non-Gaussianity of phi^3-type in the mass function and bias of dark matter haloes
aa r X i v : . [ a s t r o - ph . C O ] N ov Signature of primordial non-Gaussianity of φ -type in the mass function and bias ofdark matter haloes Vincent Desjacques and Uroˇs Seljak , , Institute for Theoretical Physics, University of Z¨urich, Winterthurerstrasse 190, CH-8057 Z¨urich, Switzerland Physics and Astronomy Department, University of California, and Lawrence Berkeley National Laboratory,Berkeley, California 94720, USA IEU, Ewha University, Seoul, S. Korea ∗ We explore the effect of a cubic correction g NL φ on the mass function and bias of dark matterhaloes extracted from a series of large N-body simulations and compare it to theoretical predictions.Such cubic terms can be motivated in scenarios like the curvaton model, in which a large cubiccorrection can be produced while simultaneously keeping the quadratic f NL φ correction small. Thedeviation from the Gaussian halo mass function is in reasonable agreement with the theoreticalpredictions. The scale-dependent bias correction ∆ b κ ( k, g NL ) measured from the auto- and cross-power spectrum of haloes, is similar to the correction in f NL models, but the amplitude is lower thantheoretical expectations. Using the compilation of LSS data in Slosar et al. [JCAP, , 031 (2008)],we obtain for the first time a limit on g NL of − . × < g NL < +8 . × (at 95% CL). Thislimit will improve with the future LSS data by 1-2 orders of magnitude, which should test many ofthe scenarios of this type. PACS numbers: 98.65.-r, 98.80.Cq, 95.36.+x, 98.70.Vc
I. INTRODUCTION
In standard single field inflation, primordial curvatureperturbations are produced by the inflaton field as itslowly rolls down its potential ([1, 2, 3, 4]). Most of thesemodels predict a nearly scale-invariant spectrum of adi-abatic curvature fluctuations in agreement with cosmo-logical observations. In addition, very small deviationsfrom Gaussianity are expected [5, 6, 7]. Therefore, anyevidence for or against the detection of primordial non-Gaussianity would strongly constrain inflationary scenar-ios.Non-Gaussianity can be produced by nonlinearities inthe relation between the primordial curvature perturba-tion Φ (Here and henceforth, the usual Bardeen potentialin matter-dominated era) and the inflaton field, interac-tions of scalar fields, a modified dispersion relation or adeparture from the natural adiabatic vacuum state (see[8] for a review). Any non-Gaussianity that is generatedoutside the horizon induces a three-point function (or bis-pectrum) B Φ ( k , k , k ) that is peaked on squeezed tri-angles (i.e. k ≪ k ∼ k ) for realistic values of the scalarspectral index. The resulting non-Gaussianity dependsonly on the local value Φ( x ) of the Bardeen’s curvaturepotential and can thus be conveniently parametrised upto third order byΦ( x ) = φ ( x ) + f NL (cid:2) φ ( x ) − h φ i (cid:3) + g NL φ ( x ) , (1)where φ ( x ) is an isotropic Gaussian random field and f NL , g NL are dimensionless, phenomenological param-eters. While the quadratic term generates the irre- ∗ [email protected]; [email protected] ducible three-point function or bispectrum at leadingorder, the cubic term does so for the irreducible four-point function or trispectrum. These statistics can becomputed straightforwardly from a perturbative expan-sion of the homogeneous Robertson-Walker background[9, 10]. Convolved with the appropriate transfer func-tion (e.g. the radiation transfer function for the CMBtemperature anisotropy), they can be used to constrainthe value of the coupling parameters f NL and g NL . Nosignificant detection of primordial non-Gaussianity hasbeen reported from measurements of the three-point cor-relation function of the cosmic microwave background(CMB) anisotropies [11, 12, 13, 14, 15]. The tightestlimits are − < f NL <
80 at 95% confidence level [14].If O ( f NL ) ∼ O ( g NL ) then the cubic correction shouldalways be negligibly small compared to the quadraticone since curvature perturbations are typically O (10 − ).However, this condition is not satisfied by some mul-tifield inflationary models such as the curvaton sce-nario, in which a large g NL and a small f NL can be si-multaneously produced. In this model, curvature per-turbations are generated by an additional scalar field,the curvaton, whose energy density is negligible dur-ing inflation [16, 17, 18, 19]. Non-Gaussianity is gener-ated by curvaton self-interactions which effectively con-tribute a non-quadratic term to the curvaton potential[20, 21, 22, 23, 24]. While the value and the sign of g NL depend upon the exact form of the self-interaction term(which can dominate the mass term if the curvaton massis small enough and the curvaton vacuum expectationvalue during inflation is large enough [25]), it is generi-cally of magnitude | g NL | ∼ − for realistic modelsin which the ratio of the energy density of the curva-ton to the total energy density at time of decay is small.There are other realizations where one can have large g NL and small f NL [26, 27]. In ekpyrotic and cyclic models, f NL typically is of the order of a few tens while g NL isof the order of a few thousand [28]. If f NL were small,then the imprint of non-Gaussianity would be detectedonly in four-point statistics such as the CMB trispectrum[29, 30, 31]. Thus far, no observational limits have beenset on g NL by measuring the CMB trispectrum [32, 33].Nevertheless, since the current bound | f NL | .
100 impliesa relative contribution for the quadratic term of ∼ . | g NL | ∼ shouldalso be consistent with the data.Large-scale structures offer another venue to test forthe presence of primordial non-Gaussianity. Deviationfrom Gaussianity can significantly affect the high massend of the mass function [34, 35], the large-scale two-point correlation [36, 37], the bispectrum [38, 39, 40, 41]of dark matter haloes hosting the observed galaxies aswell as void abundances [42, 43] and topological mea-sures of the cosmic web [44, 45]. Recently, references[46, 47, 48] showed that the local quadratic coupling f NL φ induces a scale-dependent bias ∆ b κ ( k, f NL ) in thelarge-scale power spectrum of biased tracers,∆ b κ ( k, f NL ) = 3 f NL (cid:2) b ( M ) − (cid:3) δ c Ω m H k T ( k ) D ( z ) , (2)where b ( M ) is the linear bias parameter, H is the Hub-ble parameter, T ( k ) is the matter transfer function nor-malised to unity as k → D ( z ) is the growth factor nor-malised to (1+ z ) − in the matter era and δ c ∼ .
68 is thepresent-day (linear) critical density threshold. Reference[48] applied Eq. (2) to constrain the value of f NL us-ing a compilation of large-scale structure data and found − < f NL < +69 at 95% confidence. These limits arecomparable with those from the CMB, demonstrating thecompetitiveness of the method. Forthcoming all sky sur-veys should achieve constraints of the order of f NL ∼ f NL [52]. On the numerical side however,while simulations of structure formation have confirmedthe scaling ∆ b κ ( k, f NL ) with k [46, 53, 54, 55], the exactamplitude of the non-Gaussian bias correction remainssomewhat debatable.All numerical studies to date have only implementedthe quadratic term f NL φ . The purpose of this work is toquantify the impact of the cubic term g NL φ on the massfunction and bias of dark matter haloes extracted fromcosmological simulations and assess the ability of forth-coming measurements of the large-scale bias of galax-ies/quasars to constrain the size of a local cubic correc-tion. This paper is organized as follows. We begin witha brief description of the N-body simulations and illus-trate the extent to which the coupling g NL φ affects thematter power spectrum and the halo mass function (Sec.II). We pursue with the non-Gaussian halo bias (Sec.III), to which we derive analytically the scale-dependentand scale-independent contribution, ∆ b κ and ∆ b I , anddemonstrates the large suppression of the simulated ∆ b κ relative to theory. We then place limits on the coupling parameter g NL and forecast constraints from future large-scale surveys and CMB experiments (Sec. IV). We alsoshow that our findings consistently apply to more gen-eral models with non-zero f NL and g NL (Sec. V). Weconclude with a discussion of the results in Sec. VI. II. THE NON-GAUSSIAN SIMULATIONSA. Characteristics of the N-body runs
We utilize a series of large N-body simulations ofthe ΛCDM cosmology seeded with Gaussian and non-Gaussian initial conditions. The (dimensionless) powerspectrum of the Gaussian part φ ( x ) of the Bardeen po-tential is the usual power-law ∆ φ ( k ) ≡ k P φ ( k ) / (2 π ) = A φ ( k/k ) n s − . The non-Gaussianity is of the “local”form Φ = φ + g NL φ . We adopt the standard (CMB)convention in which Φ( x ) is primordial, and not extrap-olated to present epoch. It is important to note that thelocal transformation is performed before multiplicationby the matter transfer function. T ( k ) is computed with CMBFAST [56] for the WMAP5 best-fitting parameters[13] : h = 0 .
7, Ω m = 0 . b = 0 . n s = 0 .
96 anda normalisation of the Gaussian curvature perturbations A φ = 7 . × − at the pivot point k = 0 . − .This yields a density fluctuations amplitude σ ≈ . simulations, each of which has g NL = 0 , ± ,were run with the N-body code GADGET2 [57]. We usedthe same Gaussian random seed field φ in each set ofruns so as to minimise the sampling variance. We alsoexplored scenarios with non-zero f NL and g NL and ran 2realisations for each of the non-Gaussian models charac-terized by ( f NL , g NL ) = ( ± , − × ). In all cases, thebox size is 1600 h − Mpc with a force resolution of 0.04times the mean interparticle distance. The particle massof these simulations thus is 3 . × M ⊙ / h , enough toresolve haloes down to 10 M ⊙ / h .In the curvaton scenario, generic polynomial interac-tion terms of the form λm σ ( σ/m σ ) n (where λ is somecoupling strength and m σ is the curvaton mass) to thequadratic potential of the curvaton field σ yield | g NL | ≫ f NL is very small[21, 22]. One typically finds | g NL | ∼ O (10 ) − O (10 )when f NL varies in the range − < f NL < g NL adopted inour simulations are about an order of magnitude largerso as to produce an effect strong enough to be unam-biguously measured despite the small simulated volume.Furthermore, we have also considered positive and neg-ative values of g NL so as to assess the sensitivity of thenon-Gaussian bias to the sign of the coupling parameter.The simulations with ( f NL , g NL ) = ( − , − × ) maybe seen as a particular realisation of the curvaton modelin which the coupling constant λ is positive, and the non-quadratic term is very steep ( n ∼ −
10) but contributeslittle to the total curvaton potential.
B. Properties of the initial density field
In order to ensure that the initial conditions are suc-cessfully generated, we measure at the redshift of ourinitial conditions, z = 99, the (normalized) skewness S ( R, z ) = h δ R,z i /σ and kurtosis S ( R, z ) = ( h δ R,z i / − σ ) /σ of the density field δ R,z smoothed with a (spheri-cally symmetric) window function of characteristic radius R . We adopt a tophat filter throughout this paper. Notealso that σ ( R, z ) is the variance of smoothed density fluc-tuations at redshift z .In the weakly nonlinear regime, the skewness and kur-tosis of the density field may be written as the sum ofa part generated by gravitational clustering and a partinduced by primordial non-Gaussianity. For Zel’dovichinitial conditions [58] and Ω m ( z ) ≈
1, the contribu-tion generated by gravitational instabilities reads as[59, 60, 61, 62] S Zel3 ( R, z ) = 4 − ( n eff + 3) (3) S Zel4 ( R, z ) = 2729 −
503 ( n eff + 3) + 73 ( n eff + 3) , where n eff ( R ) is the effective spectral index at thesmoothing scale R , n eff ( R ) ≡ − d ln σ ( R, z ) d ln R − . (4)Note that the initial skewness and kurtosis given by theZel’dovich approximation differs from the exact valuespredicted by perturbation theory, to which they asymp-tote in the limit D ( z ) → ∞ [63]. In addition, the cubiccoupling g NL φ induces a nonzero kurtosis S Pri4 ( R, z ) ≡ g NL S (1)4 ( R, z ) at leading order which can be computedanalytically from the relation σ S (1)4 ( R, z ) = 4! Y i =1 Z d k i (2 π ) α R ( k i , z ) P φ ( k i ) ! × α R (cid:0) | k + · · · + k n − | , z (cid:1) , (5)where α R ( k, z ) = 23Ω m H D ( z ) k T ( k ) W R ( k ) (6)is evaluated at redshift z and W R ( k ) is the Fourier trans-form of the tophat function. When f NL = 0, primordialskewness is not generated at the first order and, there-fore, may be neglected.Fig. 1 displays the initial skewness (top panel) and kur-tosis (bottom panel) obtained upon distributing the darkmatter particles onto a regular 512 mesh (i.e. of cell size ≈ h − Mpc) using the cloud-in-cell (CIC) interpolationscheme. Symbols represent the numerical results aver-aged over the realisations. Because S Zel4 (which is the solecontribution to the kurtosis in the Gaussian case) variesconsiderably between the realisations, we only show the
FIG. 1: Skewness and kurtosis of the initial ( z = 99) densityfield as a function of smoothing radius. While the top panelshows the sum of the contributions arising from the Zel’dovichdynamics and from primordial non-Gaussianity, S = S Zel3 + S Pri3 , the bottom panel only shows the absolute value of theprimordial kurtosis, | S Pri4 | . Symbols represent the numericalresults averaged over the realisations. They have been slightlyshifted horizontally for clarity. Error bars show the scatteramong the realisations for the models with g NL = 10 . Solidlines indicate the theoretical expectations. absolute difference | S Pri4 | between the kurtosis in the non-Gaussian ( g NL = ± ) and the Gaussian ( g NL = 0)runs. As expected, the skewness is very similar amongthe Gaussian and non-Gaussian simulations. While S in the simulations agrees well with the skewness inducedby the Zel’dovich dynamics Eq. (3) (solid curve), it isgradually suppressed as the filtering radius R approachesthe cell size, presumably because of the finite resolutionwhich smoothes also the fluctuations. Note that, in thesimulations with nonzero g NL and f NL , there is a largeprimordial skewness which is of magnitude | S Pri3 | ∼ a fewon scale R . h − Mpc, consistent with theory. As canalso be seen, the absolute value of the primordial kurtosisincreases sharply with R , in very good agreement withthe theoretical prediction. C. The matter power spectrum
Non-Gaussian corrections to the primordial curvatureperturbation can renormalise the input power spectrumof fluctuations used to seed the simulations. Sinceour simulations implement its unrenormalised version∆ φ ( k ) = A φ ( k/k ) n s − , it is desirable to ascertain theeffect of the local coupling term on the simulated densitypower spectrum before discussing the halo mass functionand bias. For f NL models with | f NL | . g NL models with similar levelof non-Gaussianity.The cubic order term g NL φ renormalises the ampli-tude A φ of the power spectrum of initial curvature per-turbations to A φ → A φ + 6 g NL h φ i , where h φ i = Z d k (2 π ) P φ ( k ) . (7)For scale invariant initial conditions, h φ i has a logarith-mic divergence at large and small scales (see [64] for amore detailed discussion of perturbative corrections innon-Gaussian cosmologies). In practice, a low- and high- k cutoff are naturally provided by the finite box size andthe resolution of the simulations. Therefore, the effec-tive amplitude of density fluctuations in non-Gaussiansimulations with cubic coupling is σ + δσ with δσ = 3 g NL h φ i (8)= 3 g NL (cid:18) k k min (cid:19) − n s " − (cid:18) k min k max (cid:19) − n s A φ − n s . Recall that k = 0 . h Mpc − is our choice of normalisa-tion point, and k min and k max are the integration limitsset by the fundamental mode and the Nyquist frequencyof the periodic cubical box over which the initial condi-tions are generated. Equivalently, δσ = 3 g NL (cid:18) Lk π (cid:19) − n s (cid:2) − N n s − (cid:3) A φ − n s , (9)where N = 1024 is the number of mesh points along onedimension. This result becomes δσ = 3 g NL ln( N ) A φ inthe scale-invariant limit n s →
1. For the cosmologicalsetup considered here, the absolute deviation is δσ ≈ . (cid:16) g NL (cid:17) . (10)This correction is fairly large for the values of g NL adopted here and, therefore, must be taken into accountin the comparison between the theory and the simula-tions. As we will see below, this is especially importantwhen studying the high mass tail of the halo mass func-tion which is exponentially sensitive to the amplitude ofdensity fluctuations.The cubic coupling term g NL φ can also induce ascale-dependent correction to the matter power spec-trum which can be quantified by the fractional change β m ( k, g NL ) = P mm ( k, g NL ) /P mm ( k, g NL = 0) −
1. InFig. 2, symbols show the result of measuring β m ( k, g NL )from the snapshots at z = 0 and 2 after correction ofthe normalisation shift | δσ /σ | = 0 . k & . h Mpc − but the resulting deviation is broadlyconsistent with zero. We will thus neglect β m ( k, g NL )henceforth. D. The halo multiplicity function
Haloes were identified using the
MPI parallelised ver-sion of the
AHF halo finder [65] which is based on thespherical overdensity (SO) finder developed by [66]. Themain reason for using a SO finder is that it is moreclosely connected to the predictions of the spherical col-lapse model, on which most of the analytic formulae pre-sented in this paper are based. Namely, the virial mass M of a halo is defined by the radius at which the inneroverdensity exceeds ∆ vir ( z ) times the background den-sity ¯ ρ ( z ) [67, 68]. The value of the overdensity thresh-old ∆ vir ( z ) is obtained from the collapse of a sphericaltophat perturbation and has a dependence on redshiftthrough the matter density Ω m ( z ) [69, 70]. We discardpoorly resolved haloes and only study those containingat least 34 particles or, equivalently, with a mass largerthan M = 10 M ⊙ / h .Analytic arguments based on the Press-Schechter the-ory [71, 72] predict that the halo mass function n ( M, z )is entirely specified by the distribution νf ( ν ) of first-crossings, or multiplicity function νf ( ν ) = M n ( M, z )¯ ρ d ln Md ln ν . (11)The peak height ν ( M, z ) = δ c ( z ) /σ ( M ) is the typicalamplitude of fluctuations that produce haloes of mass M by redshift z . Here and henceforth, σ ( M ) denotesthe variance of the density field δ M smoothed on massscale M ∝ R and linearly extrapolated to present epoch,whereas δ c ( z ) ≈ . D (0) /D ( z ) is the critical linear over-density for (spherical) collapse at redshift z .Despite the lack of a compelling theoretical descrip-tion of the multiplicity function for Gaussian initial con-ditions, the fractional deviation from Gaussianity can bemodelled accurately using the Press-Schechter formalism.In this approach, the halo mass function n ( M, z ) is re-lated to the probability P ( > δ c , M ) that a region of mass M exceeds the critical density for collapse δ c ( z ) throughthe relation n ( M, z ) = − ρ/M ) dP/dM . The non-Gaussian fractional correction to the multiplicity func-tion then is R ( ν, g NL ) ≡ f ( ν, g NL ) /f ( ν,
0) = ( dP/dM )( >δ c , M, g NL ) / ( dP/dM )( > δ c , M, P ( > δ c , M, g NL ) can be computed once theprobability distribution function (PDF) of the smootheddensity field δ M , P ( δ M ), is known. Here, we will considerthe simple extensions proposed by [73] and [74], in which P ( δ M ) is generically expressed as the inverse transformof a cumulant generating function. Both extensions havebeen shown to give reasonable agreement with numericalsimulations of non-Gaussian cosmologies [53, 55, 75]. FIG. 2: Non-Gaussian fractional correction β m ( k, g NL ) = P mm ( k, g NL ) /P mm ( k, − g NL h φ i induced by the cubic coupling g NL φ . In [73], the saddle-point technique is applied directly to P ( δ M ). The resulting Edgeworth expansion is then usedto obtain P ( > δ c , M, g NL ). For f NL non-Gaussianity, ref-erence [53] found that the resulting non-Gaussian massfunction agrees well with the simulations. For g NL non-Gaussianity, neglecting cumulants other than the kurto-sis S Pri4 ( M ) (Hereafter, we drop the superscript for con-ciseness) and truncating the series expansion at S ( M ),the non-Gaussian fractional correction reads R LV ( ν, g NL ) ≈ (cid:26) σ S (cid:0) ν − ν − (cid:1) (12) − σ dS d ln ν (cid:0) ν − (cid:1)(cid:27) exp (cid:2) ν δσ (cid:3) = (cid:26) σ S (cid:0) ν − ν + 3 (cid:1) (13) − d ( σ S ) d ln ν (cid:0) ν − (cid:1)(cid:27) exp (cid:2) ν δσ (cid:3) after integration over regions above δ c ( z ). Note that wehave omitted writing the redshift dependence explicitly.Strictly speaking however, R ( ν, g NL ) depends distinctlyupon the variables M (or ν ) and z due to the presence of σ S ( M ). Our notation is motivated by the fact that themeasured non-Gaussian correction, as plotted in Fig.4,appears to depend mostly on the peak height. The ex-ponential factor in the right-hand side is the correctioninduced by the renormalisation of the amplitude of lin-ear density fluctuations, Eq. (10). For consistency, wehave also used the Press-Schechter multiplicity functionto derive this last term although a Sheth-Tormen massfunction [76] may be more appropriate.In [74], it is the level excursion probability P ( > δ c , M )that is calculated within the saddle-point approximation.Including only a cubic coupling g NL φ and truncating the resulting expression at the kurtosis, we find P ( > δ c , M, g NL ) (14) ≈ √ π σδ c (cid:18) g NL σ h φ i − S σ δ (cid:19) × exp (cid:26) − δ σ (cid:20) − g NL h φ i − S δ (cid:21)(cid:27) at first order in g NL . Note that we have already includedthe renormalisation of the fluctuation amplitude. Forrare events, σ ≪ ν ⋆ f ( ν ⋆ ) = M n ( M, z, g NL )¯ ρ d ln Md ln ν ⋆ . (15)for the non-Gaussian mass function, where ν ⋆ = δ ⋆ /σ , δ ⋆ = δ c p − δσ − S δ / and f is the same mul-tiplicity function as in the Gaussian case. Taking thederivative of the level excursion probability then gives( dP/dM )( > δ c , M, g NL )( dP/dM )( > δ c , M, ≈ exp (cid:20) S δ σ + ν δσ (cid:21) (16) × (cid:18) δ ⋆ δ c + 14! δ δ ⋆ dS d ln σ (cid:19) . The fractional change in the multiplicity function even-tually reads as R MVJ ( ν, g NL ) ≈ exp (cid:20) ν σ S + ν δσ (cid:21) (17) × (cid:26) − ν σ S − ν d ( σ S ) d ln ν (cid:27) after expanding δ ⋆ at the first order and ignoring the shiftin the normalisation amplitude, i.e. δ ± ⋆ ≈ ∓ S δ / σ S ≪ ν ≫
1, the two theoretical ex-pectations reduce to 1 + ν σ S /
24. However, they differin the coefficient of the ν σ S term, which is − / − / In local f NL models, this formula only involves the skewness S . As it is incorrectly quoted in some of the literature on non-Gaussian halo mass functions, let us write down its explicit ex-pression: R MVJ ( ν, f NL ) = exp » S δ σ – » δ δ ⋆ dS d ln σ + δ ⋆ δ c – , or, in terms of the peak height ν = δ c /σ , R MVJ ( ν, f NL ) ≈ exp » ν σS – » − ν σS − ν d ( σS ) d ln ν – , after expanding δ ⋆ = p − S δ c / FIG. 3: Variance σ (dotted), skewness σS (1)3 (dashed) andkurtosis σ S (1)4 (solid) of the smoothed linear density field δ M as a function of mass scale M . fore, we shall also consider the approximation R ( ν, g NL ) = exp (cid:20) ν σ S + ν δσ (cid:21) (18) × (cid:26) − ν σ S − ν d ( σ S ) d ln ν (cid:27) , which is designed to match better the Edgeworth expan-sion of [73] when the peak height is ν ∼ S ( M ) ≡ g NL S (1)4 ( M )of the smoothed density field δ M , which we compute an-alytically using the general formula (valid for n ≥ σ n − S (1) n ( M ) = n ! n − Y i =1 Z d k i (2 π ) α M ( k i ) P φ ( k i ) ! × α M (cid:0) | k + · · · + k n − | (cid:1) , (19)where α M ( k ) ≡ α R ( k, z = 0) in what follows. Over themass range probed by our simulations, 10 . M . × M ⊙ / h , the normalised kurtosis σ S (1)4 ( M ) isa monotonic decreasing function of M that varies in thenarrow range ∼ − × − for the top-hat filter assumedhere (see Fig. 3). Note also that the σ S term dominatesthe total contribution to the non-Gaussian correction eqs.(13), (17) and (18) when the peak height is ν & | d ( σ S ) /d ln ν | . . σ S ).The fractional correction is plotted in Fig. 4 for thehaloes extracted from the simulations at redshift z = 0 . z = 0. As we can see, the level of non- FIG. 4:
Top panel : Fractional correction to the Gaussian mul-tiplicity function of dark matter haloes as a function of thepeak height ν ( M, z ) for a coupling parameter g NL = ± .The dotted, dashed and solid curves show the theoretical pre-dictions eqs. (13), (17) and (18) at z = 0, respectively. Errorbars denote Poisson errors. For illustration, M = 10 M ⊙ / h corresponds to ν = 3 .
2, 5.2, 7.7 at redshift z = 0, 1 and2, respectively. Similarly, M = 10 M ⊙ / h and 10 M ⊙ / h correspond to ν = 1 .
9, 3, 4.5 and 1.2, 1.9, 2.9 respectively.
Bottom panel : a comparison with equation (18) evaluated at z = 0 and 2. Gaussianity in the halo multiplicity function is consis-tent with the theory. Our approximation (18) performsequally well regardless of the sign of g NL . It agrees betterwith the measurements than the formulae of [74] whichsignificantly overestimates the data for g NL = 10 , andthan that of [73] which is not always positive definite for g NL = − . The bottom panel shows that the discrep-ancy somewhat worsens at higher redshift, especially inthe case g NL = 10 . However, it is possible the agree-ment may be improved by adding higher order powers of σ S and higher order cumulants.To conclude this section, one should keep in mind thatall these extensions are based on Press-Schechter and,therefore, provide a bad fit to the Gaussian mass functionof haloes. In this respect, excursion set approaches maybe more promising since they seem to reproduce both theGaussian halo counts and the dependence on f NL [78, 79]. III. THE NON-GAUSSIAN BIAS SHIFTA. Theoretical considerations
In order to calculate the scale-dependent bias correc-tion induced by the g NL coupling term to the correlationof haloes of mass M collapsing at redshift z , we follow[47] and consider the two-point correlation ξ hh ( r ) of re-gions of the smoothed density field δ M above a threshold δ c ( z ) = ν ( z ) σ . The two-point correlation function of thislevel excursion set, which was originally derived by [37],can be expressed in the high threshold approximation as ξ hh ( r ) = − ∞ X n =2 n − X j =1 ν n σ − n j !( n − j )! (20) × ξ ( n ) (cid:18) x , · · · , x , x , · · · , x j times ( n − j ) times (cid:19)(cid:27) , where r = x − x . For the non-Gaussian model consid-ered here, the leading-order correction induced by non-zero three-point and four-point correlations of the densityfield reads∆ ξ hh = ν σ ξ (3) ( x , x , x ) + ν σ (cid:20) ξ (4) ( x , x , x , x )+ 14 ξ (4) ( x , x , x , x ) (cid:21) . (21)The non-Gaussian correction ∆ P hh to the power spec-trum of biased tracers is obtained by Fourier transform-ing this expression.In the case f NL = 0 and g NL = 0, only the four-point functions contribute at first order. It should alsobe noted that, at linear order, ξ (2) ( x , x ) amounts toa renormalisation of the linear bias and, therefore, doesnot contribute to the scale-dependent correction. Detailsof the calculation can be found in Appendix A. In short,the non-Gaussian correction ∆ P hh in the limit of longwavelength k ≪ ν ξ (4) ( x , x , x , x ) / σ ,∆ P hh ( k ) = g NL ν S (1)3 ( M ) α M ( k ) P φ ( k ) (22)= g NL b ( z ) δ c ( z ) S (1)3 ( M ) α M ( k ) P φ ( k ) , where we have used b L ( z ) = ν /δ c ( z ) as is appropriatefor high density peaks. The smoothing window that ap-pears in α M ( k ) effectively makes little difference becausewe are considering the limit where k − is much largerthan the smoothing radius, so we will omit it in the fol-lowing. For small non-Gaussianity, we can also write∆ P hh ≈ b L ∆ b κ P δ ( k ) where P δ ( k ) = α M P φ ( k ) is thepower spectrum of the smoothed density field. The scale-dependent bias correction ∆ b κ ( k, g NL ) can eventually be FIG. 5:
Top panel : Non-Gaussian bias correction com-puted from the halo-matter power spectrum of haloes ofmass
M > × M ⊙ / h extracted from the snapshotat z = 0 . P mh ( k, g NL ) /P mh ( k, − b ( k, g NL ) given by Eq. (26). The dashed, dotted and dotted-dashed curves show three separate contributions that arise atfirst order in g NL . Bottom panel : ∆ b ( k, g NL ) is replaced bythe theoretical model Eq. (27). The shaded region indicatesthe data points used to fit the parameters ǫ κ and ǫ I . Errorbars indicate the scatter among 5 realisations. recast into the form∆ b κ ( k, g NL ) = 12 g NL b L ( z ) δ c ( z ) S (1)3 ( M ) α − M ( k ) (23)= 34 g NL b L ( z ) δ c (0) D (0) D ( z ) S (1)3 ( M ) Ω m H k T ( k )= 14 g NL δ c ( z ) S (1)3 ( M ) ∆ b κ ( k, f NL = 1) , where ∆ b κ ( k, f NL ) is the scale-independent bias inducedby the quadratic coupling f NL φ , Eq. (2). We havealso assumed the Eulerian bias prescription b ( M ) =1 + b L ( M ).The change in the mean number density of haloes alsocreates a scale-independent shift which we denote by∆ b I ( g NL ). As shown in [53] for f NL models, the inclusionof this correction noticeably improves the agreement withthe simulations at wavenumber k . . h Mpc − . Using apeak-background split and considering the limit of smallnon-Gaussianity, this contribution reads∆ b I ( g NL ) = − σ ∂∂ν ln (cid:2) R ( ν, g NL ) (cid:3) ≈ − σ (cid:20) (cid:0) ν − ν (cid:1) σ S + 2 νδσ (24)+ 14! (cid:0) ν − ν (cid:1) d ( σ S ) d ln ν − ν d ( σ S ) d ln ν (cid:21) after truncating ∆ b I at first order in g NL This approxima-tion should perform reasonably well for moderate valuesof the peak height, ν .
4, for which the fractional changein the mass function, equation (18), matches well the nu-merical data. It is worth noticing that ∆ b I ( g NL ) has asign opposite to that of g NL because the bias decreaseswhen the mass function goes up. ∆ b I ( g NL ) also includesa correction induced by the renormalisation of σ . Inpractice, we estimate ∆ b I ( g NL ) for a given halo sampleby evaluating σ S and ν at the scale corresponding tothe average halo mass ¯ M of the sample, as it is unclear towhich extent R ( ν, g NL ) agrees with the data in the limit M ≫ M ⋆ . B. Comparison with the simulations
To assess the effect of primordial non-Gaussianity onthe halo bias, we will consider the ratios P mh ( k, g NL ) P mh ( k, − b ( k, g NL ) b ( M ) + 2 δσ σ (25) P hh ( k, g NL ) P hh ( k, − (cid:18) b ( k, g NL ) b ( M ) (cid:19) + 2 δσ σ − , where ∆ b ( k, g NL ) is generally the sum of a scale-dependent and a scale-independent term. One shouldbear in mind that the scale-independent shift 2 δσ /σ arises from the matter power spectrum and, therefore, isdistinct from the term − νδσ /σ appearing in ∆ b I . Fol-lowing [53], we shall also quantify the departure from thetheoretical scaling as a function of wavemode amplitudewith the ratio ∆ b s / ∆ b t , where ∆ b s is the non-Gaussianbias correction measured from the simulation and ∆ b t isEq. (27).We interpolate the dark matter particles and halo cen-tres onto a regular cubical mesh. The resulting darkmatter and halo fluctuation fields are then Fourier trans-formed to yield the matter-matter, halo-matter and halo-halo power spectra P mm ( k ), P mh ( k ) and P hh ( k ), respec-tively. These power spectra are measured for a range ofhalo masses and redshifts, covering the relevant range ofstatistical properties corresponding to the available datasets of galaxy or quasar populations with different lumi-nosities and bias. Note that these quantities are com-puted on a 512 grid, whose Nyquist wavenumber is suf-ficiently large ( ≈ h Mpc − ) to allow for an accuratemeasurement of the power in wavemodes of amplitude k . . h Mpc − . The halo power spectrum is corrected FIG. 6: Best-fitting ǫ κ and ∆ b I + ǫ I as a function of halo biasand g NL = ± , for two different mass cuts as indicated inthe Figure. In the top panel, the dotted curve is our best fitto ǫ κ ( b, g NL ). In the bottom panel, the solid and dashed linesshow the scale-independent shift ∆ b I predicted by a peak-background split, Eq. (24), for the low and high mass samples,respectively. for the shot-noise due to the discrete nature of dark mat-ter haloes, which we assume to be the standard Pois-son term 1 / ¯ n h . This discreteness correction is negligi-ble for P mm ( k ) due to the large number of dark matterparticles. Yet another important quantity is the linearhalo bias b ( M ) which must be measured accurately fromthe Gaussian simulations as it controls the magnitude ofthe scale-dependent shift. Here, we shall use the ratio P mh ( k ) /P mm ( k ) as a proxy for the halo bias since it isless sensitive to shot-noise.
1. An effective non-Gaussian bias correction
Summarizing the analytical considerations of Sec.III A, non-Gaussianity of the g NL type add a correction∆ b ( k, g NL ) to the bias b ( k ) of dark matter haloes whichis at leading order∆ b ( k, g NL ) = ∆ b κ ( k, g NL ) + ∆ b I ( g NL ) , (26)We found, however, that this theoretical expectationsignificantly overestimates the magnitude of the non-Gaussian bias shift measured from the simulations.This is exemplified in the top panel of Fig. 5, where P mh ( k, g NL ) /P mh ( k, − M > × M ⊙ / h identified at z = 0 .
5. Clearly, the
FIG. 7: Non-Gaussian bias correction measured in the simulation outputs at redshift 0 < z <
M > × M ⊙ / h . In each panel, the upper plot shows the ratio P hh ( k, g NL ) /P hh ( k, − P mh ( k, g NL ) /P mh ( k, − b s / ∆ b t . The shaded area indicates the domain where the deviation is less than 20 per cent. Theparameters ǫ κ and ǫ I are fitted individually to each sample. For illustration, ǫ κ takes the best-fit values 0.06 and 0.60 for thelowest and highest biased samples, respectively. FIG. 8: Same as Fig. 7 but for haloes in the mass range 10 < M < × M ⊙ / h . predicted scale-dependent correction ∆ b κ is much steeperthan measured from the halo samples. In order to im-prove the agreement with the numerical data, we modifythe above relation as follows :∆ b ( k, g NL ) = ǫ κ ∆ b κ ( k, g NL ) + (cid:2) ∆ b I ( g NL ) + ǫ I (cid:3) , (27)and treat ǫ κ and ǫ I as free parameters that we fit to ourmeasurements of the cross-power spectrum (weighted bythe scatter among 5 realisations) in the range 0 . ≤ k ≤ . h Mpc − where the scale-dependent effect islargest. The bottom panel of Fig. 5 shows the result-ing best-fit contributions ǫ κ ∆ b κ and ∆ b I + ǫ I for thehalo sample mentioned above. As seen in Fig. 6, ǫ κ and ǫ I appear to depend mainly upon the linear halo bias b ( M ) and the coupling parameter g NL , although depen-dencies on redshift or other halo observables are not ex-cluded (The data is too noisy for a reliable estimate ofthese). The most striking feature of Fig. 6 is the func-tional dependence of ǫ κ on b ( M ) and g NL . Firstly, ǫ κ is a monotonically increasing function of the bias andnever reaches unity, even for the most biased samplesfor which the high peak approximation should be valid. Secondly, ǫ κ is noticeably larger for g NL = − , suggest-ing thereby that second (and higher) order contributionsto the scale-dependent bias may be important. Further-more, for b . . g NL , in agreement with theoretical expectations from thepeak-background split (see Sec.III A). However, whereasfor b . b &
3, reaching up to 5-10 per cent of the linear halobias.Assuming ǫ κ is a function of b ( M ) and g NL only andasymptotes to a constant in the highly biased limit, wefind that the following parametrised form ǫ κ ( b, g NL ) = c − c g NL − c c b ) (28)captures reasonably well the increase of ǫ κ with halo biasfor 1 . < b ( M ) <
7. The best-fit values of the parameters1are c = 0 . ± . , c = (6 . ± . × − (29) c = 2 . ± . , c = 0 . ± . . We do not provide a fitting formula for ǫ I (or ∆ b I + ǫ I )since it is not directly measurable in real data.Since both the kurtosis of the initial density field andthe mass function of dark matter haloes are in fairly goodagreement with theoretical expectations, the discrepancyin the scale-dependent bias ∆ b κ ( k, g NL ) indicates that itis the high peak approximation considered here which isinaccurate, even in the limit b ( M ) ≫ b κ ( k, g NL ) under theassumption that the non-Gaussian correction to the halopower spectrum is small whereas, for most of the biasrange covered by our halo catalogues, the effect is in factalready very large at k . . h Mpc − . Note, however,that for the quadratic coupling f NL φ , this perturbativetreatment predicts an effect of the right magnitude [53].Explaining these findings clearly requires a better the-oretical understanding, which we leave for a future in-vestigation. Notwithstanding this, the phenomenologicalprescription (27) with parameters ǫ κ and ǫ I fitted to thedata provides, as we will see below, a good descriptionof the large-scale halo power spectrum in simulations of g NL models.
2. Non-Gaussian bias from auto- and cross-power spectra
We have measured auto- and cross-power spectra fora range of halo masses and redshifts spanning the range0 < z <
2. The ratios defined in Eq.(25) are shownin Figs 7 and 8 as a function of wavenumber for themass threshold
M > × M ⊙ / h and the mass bin10 < M < × M ⊙ / h , respectively. The fractionaldeviation ∆ b s / ∆ b t is also shown at the bottom of eachpanel. The shaded region indicates a departure less than20 per cent. Error bars denote the scatter around themean and, therefore, may underestimate the true errorsas they are computed from a small number of realisations.Note that, in order to reduce the impact of sampling vari-ance, we first compute the ratios P mh ( k, g NL ) /P mh ( k, P hh ( k, g NL ) /P hh ( k,
0) for each realisation before cal-culating the average.As we can see, once ǫ κ and ǫ I are fitted to the ratio ofcross-power spectra, the theoretical prediction Eq. (27)provides a reasonable description of the non-Gaussianbias in the halo power spectrum P hh , indicating thatnon-Gaussianity does not generate much stochasticityand the predicted scaling ∆ b κ ( k, g NL ) ∝ k − T ( k ) − ap-plies equally well for the auto- and cross-power spectrum.This was also found to be true in f NL models [53]. Theinclusion of a scale-independent correction ∆ b I + ǫ I sig-nificantly improve the agreement at k . . h Mpc − . For the highly biased samples b > k & . h Mpc − . Suchan effect is not seen in the simulations, but we expectlarge deviations from the relation (27) in that rangeof wavenumber, where second- and higher-order correc-tions induced by the cubic coupling g NL φ together withthe nonlinear bias created by the gravitational evolu-tion of matter density fluctuations may become impor-tant. Even though the data is noisier due to the lownumber density of haloes, it is worth noting that, forthe highly biased samples at k . . h Mpc − , thecross-power spectrum P mh ( k, g NL = − ) goes neg-ative while P hh ( k, g NL = − ) remains positive andincreases sharply, in agreement with the analytic pre-diction. Still, there is some evidence that the ratio P hh ( k, g NL = − ) /P hh ( k, − δσ /σ − ≈ − .
96 at the minimum.Fig. 8 further explore the effect in the low mass sam-ples, for which the z = 0 haloes with b ( M ) ≈ .
15 consti-tute an almost unbiased sample of the density field. Inthis case, the sign of the scale-dependent contribution isreversed, namely, the large-scale halo power spectrum insimulations of g NL = − is enhanced relative to thatof the Gaussian ones. This is in rough agreement withthe theory, which predicts a similar effect for b ( M ) < z = 1 .
39 shown in Fig. 8 corresponds closely to the quasarsample used by [48], for which z = 1 . b = 2 . IV. CONSTRAINTS ON THE COUPLINGPARAMETER g NL A. Constraints on g NL from current large-scalestructure data Reference [48] took advantage of the scale-dependenceof the bias to constrain f NL from a sample of highly bi-ased luminous red galaxies (LRGs) and quasars (QSOs).It is straightforward to translate their 2- σ limit − 8, which covers a large comovingvolume and is highly biased, b = 2 . 7. In light of ourresults (see Fig. 6), we expect the parameter ǫ κ ( b, g NL )to vary with g NL . However, in order to simplify theanalysis, we will assume that, at fixed b , ǫ κ ( b, g NL ) isgiven by the mean of ǫ κ ( b, g NL = ± ). For a sam-ple with bias b ∼ . 7, this implies ǫ κ ≃ . 4. Further-more, assuming the typical mass of QSO-hosting haloesis ∼ M ⊙ / h yields S (1)3 ( M ) ≃ . × − . Hence,the multiplicative factor (1 / δ c ( z ) ǫ κ S (1)3 ( M ) is approx-2imately ≃ . × − . Our limits on g NL thus are − . × < g NL < +8 . × (30)at 95% confidence level. The scale-independent correc-tion ∆ b I + ǫ I is not directly measured as it adds up tothe bias b which is fitted to the data. For the limits ob-tained here, | ∆ b I + ǫ I | should be much smaller than b andcan thus be ignored. Note also that, whereas the non-Gaussian bias scales as D ( z ) − in f NL models, we have∆ b ( k, g NL ) ∝ D ( z ) − for g NL non-Gaussianity, so one canachieve relatively larger gains from measurements of highredshift tracers. In fact, the extent to which one can im-prove the observational bounds will strongly depend onour ability to minimize the impact of sampling variancecaused by the random nature of the wavemodes, and theshot-noise caused by the discrete nature of the tracers.By comparing differently biased tracers of the same sur-veyed volume [51] and suitably weighting galaxies (e.g.by the mass of their host halo) [80, 81], it should bepossible to circumvent these problems and considerablyimprove the detection level. B. Predictions for future LSS surveys References [50, 51, 64, 82] applied the Fisher matrixformalism to forecast constraints on f NL from forthcom-ing galaxy redshift surveys. Here, we will simply try toestimate the detection limit for g NL . Following [50, 64],we consider a (nearly spherical) survey of volume V .Assuming the Fourier modes are still uncorrelated andGaussian distributed, the total signal-to-noise squaredreads (cid:18) SN (cid:19) ≈ V π Z k max k min dk k "(cid:18) b κ b (cid:19) − (31)in the limit where sampling variance dominates the er-rors. Here, k min ∼ π/V / is the smallest wavemodeaccessible and k max is not necessarily finite since the in-tegral does converge as one takes k max to infinity. Sub-stituting the expression Eq. (23) for the scale-dependentbias ∆ b κ ( k, g NL ) and setting T ( k ) ≡ (cid:18) SN (cid:19) ≈ Vπ ( k ⋆ ) (cid:18) k min − k max (cid:19) , (32)where k ⋆ ≃ . × − g NL ǫ κ (1 − /b ) D ( z ) S (1)3 − ! h Mpc − . (33)We have also assumed | g NL | . , such that | k ⋆ | is atmost of the order of k and the term linear in ∆ b κ /b dominates the signal. When k min ≪ k max , we can furthersimplify ( S/N ) to (cid:18) SN (cid:19) ≈ . × − g ǫ κ (cid:18) − b (cid:19) D ( z ) − × S (1)3 − ! (cid:18) V h − Gpc (cid:19) / . (34)Note the strong sensitivity of the signal-to-noise squaredto the growth factor D ( z ) (For f NL non-Gaussianity, thisdependence is only D ( z ) − ).To highlight the improvement one could achieve withfuture galaxy surveys, it is useful to first calculate the de-tection limit for the SDSS LRG sample centred at z ∼ . v ≈ h − Gpc . Assuming a lin-ear bias b = 2 and a skewness parameter S (1)3 ∼ × − appropriate for haloes of mass M ∼ − M ⊙ / h ,the minimum g NL detectable at the 1- σ level is ≃ for a correction factor ǫ κ = 0 . , with central redshift z = 0 . V = 6 h − Gpc , the minimum g NL wouldbe ∼ × for galaxies tracing haloes of similar massand bias. Finally, for a configuration like EUCLID with a V = 100 h − Gpc survey centred at z = 1 . ∼ . × . Clearly, theselimits are only indicative: they may be significantly im-proved by selecting highly biased, high redshift (single-or multi-)tracers. Nevertheless, this shows that futuregalaxy surveys should furnish interesting constraints onthe size of the cubic coupling g NL φ . C. Predictions for CMB temperature anisotropies The CMB trispectrum provides an alternative probe oflocal, non-quadratic correction to the Gaussian curvatureperturbations, so it is interesting to assess the sensitivityof this statistics to the nonlinear parameter g NL .The temperature anisotropy field is convenientlydecomposed into spherical harmonics, ∆ T ( ˆn ) /T = P lm a ml Y ml ( ˆn ). As shown in [29, 30], statistical isotropyand invariance under parity transformation ˆn → − ˆn im-plies that the 4-point correlation of the spherical har-monic coefficients a ml takes the form h a m l a m l a m l a m l i = X LM ( − M (cid:18) l l Lm m − M (cid:19) (35) × (cid:18) l l Lm m M (cid:19) Q l l l l ( L ) . Q l l l l ( L ) is the angular average trispectrum andbrackets are Wigner-3j symbols. Statistical homogene-ity also implies that Q l l l l ( L ) is independent of position.The connected part of the trispectrum, T l l l l ( L ), encodesinformation about non-Gaussianity and is obtained bysubtracting a Gaussian piece constructed from the powerspectra C l . Eq.(35) can be inverted with the aid of theorthogonality of the Wigner-3j symbols to form an esti-mator for the CMB trispectrum.The signal-to-noise for the CMB trispectrum T l l l l ( L )summed up to a certain l max is [30] (cid:18) SN (cid:19) ( < l max ) ≈ l max X l >l >l >l X L | T l l l l ( L ) | (2 L + 1) C l C l C l C l (36)when cosmic variance dominates the errors. Otherwise,one shall include a contribution from the power spec-trum of the detector noise to the C l . Galactic foregroundsubtraction on a fraction 1 − f sky would further reduce( S/N ) by a factor of f sky .Neglecting the ISW effect, the Sachs-Wolfe provides auseful order-of-magnitude estimate of the signal-to-noiseas long as l max does not exceed . 100 [7, 30, 31, 83]. Thecalculation is performed in Appendix B. We find that thesignal-to-noise can be recast into the compact form (cid:18) SN (cid:19) ( < l max ) = 92 g A φ (cid:26) Z +1 − dx s l max ( x ) t l max ( x )+ 12 Z +1 − dx r l max ( x ) s l max ( x ) (cid:27) (37)where the auxiliary functions r l ( x ), s l ( x ) and t l ( x ) aredefined as r l ( x ) = l X k =2 (2 k + 1) P k ( x ) (38) s l ( x ) = l X k =2 (2 k + 1) k ( k + 1) P k ( x ) (39) t l ( x ) = l X k =2 (2 k + 1) k ( k + 1) P k ( x ) . (40)Here, P l ( x ) are Legendre polynomials. Note that we haveexcluded the monopole and dipole from the summationsince these modes are unobservable. We have also as-sumed a nearly scale-invariant spectrum n s ≈ g NL = 1. Although this approximation breaksdown for l max & l max = 200 so as to extrapolate more robustly the l max -dependence to small angular resolution. A power-law fit to ( S/N ) in the range 50 ≤ l max ≤ 200 gives (cid:18) SN (cid:19) ( < l max ) ≃ . × − g (cid:18) A φ − (cid:19) l . . (41) FIG. 9: Signal-to-noise ratio squared for the CMB trispec-trum as a function of the maximum multipole l max . We haveassumed g NL = 1 and f sky = 1. Our results appear consistent with the findings of [30]shown in their Fig.2. However, our constant of pro-portionality is about 20-30 times larger, presumably be-cause they adopted a lower fluctuation amplitude (com-pare also their prediction for the f NL model with that of[31]). Adding the information encoded in temperature-polarization trispectra may enhance ( S/N ) by a factorof a few [29].Assuming the scaling Eq.(41) persists well beyond therange over which the Sachs-Wolfe effect dominates, theminimum g NL detectable at 1- σ level is g NL ≃ 20, 7.9,3.2, 1.9 and 1 . × for l max = 250, 500, 1000, 1500and 2000. A more realistic calculation should includethe full radiation transfer function, detector noise etc.In this respect, detailed calculations have shown that,for the quadratic coupling f NL φ , ( S/N ) of the CMBbispectrum and trispectrum closely follows the behaviourobtained in the Sachs-Wolfe approximation [30, 31]. Itseems reasonable, then, to expect that this is also truefor g NL φ .While our predictions are qualitative they show that,for the WMAP CMB temperature measurement (whichwe approximate as a noise-free experiment with l max ∼ | g NL | ≤ × at the 1- σ level. This is of the same orderas the limit we derived from the QSO sample analyzed by[48]. For a PLANCK-like experiment ( l max ∼ | g NL | ≤ . × at the 1- σ level. This is comparable to the detection http://map.gsfc.nasa.gov/ V. EFFECT OF NON-GAUSSIANITY WITHNON-ZERO f NL AND g NL In this Section, we examine the halo multiplicity func-tion and large-scale bias in numerical simulations ofstructure formation with non-zero coupling parameters( f NL , g NL ) = ( ± , − × ). We show that the resultsare consistent with those obtained from the simulationswith non-vanishing g NL solely. A. Mass function It is straightforward to calculate the fractional devia-tion from the Gaussian mass function, Eq.(18), to non-zero f NL and g NL . Again, we start with the MVJ formulaand neglect second order corrections such as ( σS ) etc.Adjusting the coefficient of the terms νσS and ν σ S to that of the small ν expansion obtained by [73], wearrive at R ( ν, NL) = exp (cid:20) ν σS + ν σ S + ν δσ (cid:21) × (cid:26) − ν σS − ν d ( σS ) d ln ν (42) − ν σ S − ν d ( σ S ) d ln ν (cid:27) , where the shorthand notation NL designates the combi-nation ( f NL , g NL ). In Fig. 10, this theoretical predictionis compared R ( ν, NL) measured in non-Gaussian simu-lations of ( f NL , g NL ) = ( ± , − × ). We accountfor the fact that the amplitude of density fluctuations isrenormalised by δσ ≈ . z = 2 when f NL = − f NL = 100, the positive and negativecontributions from the quadratic and cubic coupling, re-spectively, almost cancel each other and flatten the frac-tional deviation over most of the mass range probed bythe simulations.The scale-independent bias shift which arises from thechange in the mean number density of haloes can againbe estimated using the peak-background split. We find∆ b I (NL) ≈ − σ (cid:20) (cid:0) ν − ν (cid:1) σ S + 12 (cid:0) ν − (cid:1) σS + 14! (cid:0) ν − ν (cid:1) d ( σ S ) d ln ν − ν d ( σ S ) d ln ν + 16 (cid:0) ν − (cid:1) d ( σS ) d ln ν − d ( σS ) d ln ν +2 νδσ (cid:21) (43) at the first order in the nonlinear parameters f NL and g NL . B. Bias Having checked that the amount of non-Gaussianityin the mass function is also consistent with our simpletheoretical expectation when both f NL and g NL are non-zero, we now turn to the clustering of dark matter haloes.As shown in Appendix A, the non-Gaussian correction tothe halo power spectrum can be written down as∆ P hh ( k ) = 4 f NL b δ c ( z ) α M ( k ) P φ ( k )+ 4 f b δ ( z ) P φ ( k ) + (cid:18) g NL + 43 f (cid:19) × b δ ( z ) S (1)3 ( M ) α M ( k ) P φ ( k ) . (44)If we set g NL = 0 and keep only the first two terms in theright-hand side, then the non-Gaussian (Eulerian) halopower spectrum can be cast into the form P hh ( k ) = (cid:2) b ( M ) + f NL b φ ( k ) (cid:3) P δ ( k ) (45)where the scale-dependent bias parameter b φ ( k ) is b φ ( k ) = 2 (cid:2) b ( M ) − (cid:3) δ c ( z ) α − M ( k ) . (46)Note that reference [64] obtained this relation by con-sidering the halo power spectrum implied by a bias re-lation that is a local mapping of the density field. Inpractice, the term proportional to P φ ( k ) is negligibleas it contributes only at very small wavenumber k . . h − Mpc. The third term in the right-hand sideof Eq.(44) is derived in this paper for the first time. Inthe case g NL = 0, its magnitude relative to the termlinear in f NL is (1 / f NL δ c ( z ) S (1)3 ( M ), which is approx-imately ∼ . 03 at redshift z = 1 . M = 10 M ⊙ / h . Although its contribution becomes in-creasingly important at higher redshift, it is fairly smallfor the values of f NL considered here. Consequently, weshall employ the approximation∆ b ( k, NL) = ǫ κ ∆ b κ ( k, g NL ) + ∆ b κ ( k, f NL )+ (cid:2) ∆ b I (NL) + ǫ I (cid:3) . (47)to describe the non-Gaussian bias of dark matter haloes.The quadratic coupling f NL φ also affect the matterpower spectrum at leading order [38, 84], positive valuesof f NL increasing the small scale power. However, therelative size of this k -dependent correction, β m ( k, f NL ),is at a per cent level only in the weakly nonlinear regime k . . h Mpc − [53, 85] and fades rapidly as one goes tolarger scales. We will thus neglect it in what follows.In Fig. 11, the result of measuring ratios of auto- andcross-power spectra in the simulations with ( f NL , g NL ) =( ± , − × ) is shown at 0 < z < . M > × M ⊙ / h . We do not quote error bars5 FIG. 10: Fractional correction to the Gaussian multiplicityfunction of dark matter haloes as a function of the peakheight ν ( M, z ) for f NL = ± 100 and a cubic coupling param-eter g NL = − × . The solid and dashed curves show thetheoretical prediction Eq.(42) at z = 0 and 2, respectively.Error bars denote Poisson errors. since the data points are obtained by averaging over tworealisations only. The solid and dashed curves show thetheoretical prediction Eq. (47). The value of the mul-tiplicative factor ǫ κ ( b, g NL ) was obtained from the four-parameter formula Eq.(28), while ∆ b I + ǫ I was individ-ually fitted for each halo sample over the wavenumberrange 0 . < k < . h Mpc − . As can be seen, the the-oretical expectation Eq. (47) agrees reasonably well withthe numerical data. This demonstrates that the rangeof validity of the non-Gaussian bias formula Eq. (27) ex-tends to smaller values of g NL as well as models withnon-vanishing f NL and g NL .The lowest order, k -dependent corrections to the Gaus-sian bias induced by the quadratic and the cubic couplingare fully degenerated in the halo power spectrum as theyboth scale as α − M ( k ) ∝ k − T ( k ) − . For the values of f NL and g NL and the halo mass range considered here, theratio ∆ b κ ( k, g NL ) / ∆ b κ ( k, f NL ) increases approximatelyfrom 0.25 to 0.5 when the redshift increases from z = 0to 2. It is unclear whether higher order corrections couldhelp breaking such a degeneracy. A more promising al-ternative may be to measure the bispectrum of dark mat-ter haloes, which carries much more information aboutthe shape of the primordial three-point function than thepower spectrum of bias tracers [40, 41]. However, this isbeyond the scope of this paper. VI. DISCUSSION In this paper we explored the effect of a local cubic cou-pling g NL φ on the mass function and bias of dark matterhaloes. We derived analytical expressions for these statis-tics and tested them against the outcome of numericalsimulations.We showed that current theoretical predictions of thenon-Gaussian correction to the mass function reasonablyagree with the simulations. The LV formula [73] appearsto provide a better fit to the data than the MVJ formula[74], in agreement with some of the literature on the sub-ject [46, 86]. The two approximations can be combined toprovide an accurate description if one adjusts the low- ν expansion of the latter so as to match that of the former.We found that the magnitude of the non-Gaussianscale-dependent bias ∆ b κ ( k, g NL ) is suppressed relative toa theoretical prediction based on the statistics of highlyoverdense regions, even on linear scales k . . h Mpc − .This suppression is stronger for the low biased sam-ples b . b , for positive valuesof g NL . We were able to fit the measured halo biasat the expense of introducing two free parameters, ǫ κ and ǫ I , that depend mostly on the halo bias b ( M ) andthe coupling parameter g NL . These parameters quantifythe departure from the theoretical scale-dependent andscale-independent non-Gaussian bias correction, respec-tively. We provide a simple fitting formula for ǫ κ ( b, g NL ),Eq.(28), which should be used when analyzing observa-tional data. In non-Gaussian simulations of the f NL type,the data also hint at a (possibly f NL -dependent) suppres-sion of the non-Gaussian scale-dependent bias relative totheory for wavemodes k . . h Mpc − [53, 54, 55], butthe effect is much weaker than seen in our simulationsof g NL models. Clearly, these results require a bettertheoretical modelling of the non-Gaussian halo bias.Reference [55] argued that both the MVJ and LV an-alytic formula can be reconciled with measurements ofthe non-Gaussian fractional correction to the mass func-tion once non-spherical collapse is included. In practice,the critical density for collapse is replaced by δ c → √ qδ c ,where the value q = 0 . 75 is obtained from a fit to themass function measured in simulations [87]. Reference[88] claimed that such a relation is a consequence ofthe diffusive nature of the critical threshold for collapse.Their model predicts q ≃ . 8, in good agreement withthe findings of [55]. However, we found that substitut-ing δ c → √ qδ c in Eq.(18) only modestly improve theagreement with the data. Regarding the non-Gaussianbias, it is not obvious how one could justify the replace-ment δ c → √ qδ c given that the linear bias of our (Gaus-sian) halo samples converges towards the spherical col-lapse prediction ν /δ c for large peak height.A important ingredient is the choice of the halo iden-tification algorithm. While we used a spherical overden-sity (SO) finder, reference [88] considered a Friends-of-Friends (FoF) finder with a linking length b = 0 . 2. Thequestion of how the spherical overdensity masses can be6 FIG. 11: Non-Gaussian bias correction measured in the simulations with ( f NL , g NL ) = ( ± , − × ) for haloes of mass M > × M ⊙ / h . Error bars are not shown as the data points are averaged over two realisations solely. mapped onto friends-of-friends masses remains a matterof debate [89]. Clealy however, since the peak height de-pends on the halo mass M through the variance σ ( M ),any systematic difference will be reflected in the value of ν associated to a specific halo sample. This will in turnaffects the size of the fractional deviation from the Gaus-sian mass function at some specified peak height. Thesensitivity of the non-Gaussian mass function and biasto the halo finder will be presented elsewhere.The observational bound on f NL inferred by [48] fromthe clustering of a high redshift sample of quasars canbe straightforwardly translated into a limit on g NL since∆ b κ ( k, g NL ) also scales as k − T ( k ) − at low wavenumber.We have obtained − . × < g NL < +8 . × (95%) . (48)These are the first limits derived on g NL . While theyare too weak to provide interesting constraints on infla-tionary scenarios such as the curvaton model, future all-sky redshift surveys should improve them by a factor of ∼ Acknowledgements We thank Paolo Creminelli and Leonardo Senatore foruseful discussions. We acknowledge support from theSwiss National Foundation under contract No. 200021-116696/1 and WCU grant R32-2008-000-10130-0. APPENDIX A: NON-GAUSSIAN BIAS IN THEHIGH PEAK LIMIT In this Appendix, we detail the derivation of the scale-dependent bias correction induced by the g NL couplingto the two-point correlation of dark matter haloes. Wefollow [47] and approximate the latter by the two-point7correlation ξ hh ( r ) of regions of the smoothed density field δ M with a peak height ν ≫ 1. case g NL = 0 only In this case, only the four-point correlations of the den-sity field contribute at the first order. The leading-ordercorrection to the correlation of tresholded regions thusreads as∆ ξ hh = ν σ (cid:20) ξ (4) ( x , x , x , x ) + 14 ξ (4) ( x , x , x , x ) (cid:21) . (A1)The four-point correlations ξ (4) are Fourier transform ofthe trispectrum of the density field, T δ ( k , k , k , k ),with conservation of momentum enforced. Linearly ex-trapolating the density field to present epoch, the lattercan be expressed as T δ ( k , k , k , k ) = Y i =1 α M ( k i ) ! T Φ ( k , k , k , k ) , (A2) where the expression for the trispectrum of primordialcurvature perturbation T Φ ( k , k , k , k ) = 6 g NL (cid:2) P φ ( k ) P φ ( k ) P φ ( k )+(cyclic) (cid:3) (A3)follows straightforwardly from the Fourier mode relationΦ( k ) = φ ( k )+ g NL Z d k (2 π ) Z d k (2 π ) φ ( k ) φ ( k ) φ ( k − k − k ) . (A4)Density and curvature perturbations are related throughthe Poisson equation, whose scale-dependence is reflectedin the function α M ( k ) = α R ( k, z = 0) defined in Eq. (6).Combining these relations gives ξ (4) ( x , x , x , x ) = 6 g NL Y i =1 Z d k i (2 π ) α M ( k i ) P φ ( k i ) ! (cid:20) P φ ( k ) P φ ( k ) (cid:21) α M ( k ) e i k · r (A5) ξ (4) ( x , x , x , x ) = 6 g NL Y i =1 Z d k i (2 π ) α M ( k i ) P φ ( k i ) ! (cid:20) P φ ( k ) P φ ( k ) + P φ ( k ) P φ ( k ) (cid:21) α M ( k ) e i k · r , (A6)where we have defined k ij ··· l = k i + k j + · · · + k l forshorthand convenience. Since we will examine the effectof non-Gaussianity on Fourier space statistics only, we take the Fourier transform of the four-point functions.After some simplification, we arrive at Z d r ξ (4) ( x , x , x , x ) e − i k · r = 6 g NL α M ( k ) P φ ( k ) Z d k (2 π ) α M ( k ) P φ ( k ) Z d k (2 π ) α M ( k ) P φ ( k ) × α M (cid:0) | k + k | (cid:1) " P φ (cid:0) | k + k | (cid:1) P φ ( k ) (A7) Z d r ξ (4) ( x , x , x , x ) e − i k · r = 6 g NL Z d k (2 π ) α M ( k ) α M ( | k + k | ) P φ ( k ) P φ (cid:0) | k + k | (cid:1) Z d k (2 π ) α M ( k ) × α M (cid:0) | k + k | (cid:1) P φ ( k ) " P φ (cid:0) | k + k | (cid:1) P φ (cid:0) | k + k | (cid:1) + P φ (cid:0) | k + k | (cid:1) P φ ( k ) . (A8)For realistic values of the spectral index ( n s ∼ α M ( | k + k i | ) P φ ( | k + k i | ) appearing in the right- hand side of the above equalities formally diverges when-ever k + k = 0 due to the ultraviolet divergence of8 P φ ( k ). To cure this problem, one can set P φ ( k ) = 0for sufficiently small wavenumbers or excise a thin shellcentred at wavenumber k i from the integral. In the large- scale limit k ≪ k i , the ratio P φ ( | k + k i | ) /P φ ( k ) vanisheswhereas P φ ( | k + k i | ) /P φ ( k i ) tends towards unity. In thiscase, the above expressions reduce to Z d r ξ (4) ( x , x , x , x ) e − i k · r ≈ g NL σ S (1)3 ( M ) α M ( k ) P φ ( k )+ 6 g NL α M ( k ) Z d k (2 π ) α M ( k ) P φ ( k ) Z d k (2 π ) α M ( k ) P φ ( k ) α M ( k ) P φ ( k ) , (A9)and Z d r ξ (4) ( x , x , x , x ) e − i k · r ≈ g NL σ Z d k (2 π ) α M ( k ) P φ ( k ) . (A10)Only the first term in the right-hand side of Eq.(A9) isnot well behaved in the limit k → k n s − . The second scales as k , whilethe Fourier transform of ξ (4) ( x , x , x , x ) asymptotesto a constant. A similar decomposition also arises in f NL models. For this type of non-Gaussianity, the firstorder correction is furnished by the three-point function ξ (3) ( x , x , x ), whose Fourier transform can be split intothe familiar term 2 f NL σ α M ( k ) P φ ( k ), and a second piecegiven by12 f NL α M ( k ) Z d k (2 π ) α M ( k ) P φ ( k ) . (A11)In both quadratic and cubic local non-Gaussianity, theterm proportional to α M ( k ) can be neglected since, atthe pivot point k = k , its magnitude relative to theterm involving α M ( k ) P φ ( k ) is only O (0 . 01) and O (10 − ),respectively. Moreover, it decreases as one goes to largerscales. By contrast, it is not so obvious how to handle theterm (A10). In the non-Gaussian halo power spectrum,this term would appear multiplied by ν / (4 σ ),6 g NL b δ ( z ) σ Z d k (2 π ) α M ( k ) P φ ( k ) (A12) ≃ × − g NL b D (0) D ( z ) (cid:18) M M ⊙ / h (cid:19) . . The approximation (second line) holds for 10 ≤ M ≤ M ⊙ / h . For g NL = 10 and b L & 3, this can be muchlarger than the typical shot-noise correction applied to the halo power spectra we measure in the simulations(see below). For g NL = − , this will certainly pro-duce a halo power spectrum which is negative at suffi-ciently low wavenumber. It is plausible that higher ordercounter-terms in the expansion Eq.(21) renormalises itsvalue. However, such a calculation is beyond the scopeof this paper, so the simplest choice is to ignore thisterm in the following of the analysis. Hence, we canapproximate the non-Gaussian correction ∆ P hh in thelimit of long wavelength k ≪ ν ξ (4) ( x , x , x , x ) / σ , Eq. (22). 2. case g NL and f NL non-zero In addition to the first order trispectrum induced by g NL φ , Eq. (A3), the quadratic coupling f NL φ generatesthe bispectrum at leading order B Φ ( k , k , k ) = 2 f NL (cid:2) P ( k ) P ( k ) + (cyclic) (cid:3) , (A13)and an additional, albeit second order, contribution tothe trispectrum [90, 91, 92], T Φ ( k , k , k , k ) = 4 f (cid:2) P φ ( k ) P φ ( k ) P φ ( k )+11 permutations (cid:3) . (A14)The bispectrum (A13) induces a three-point contribution( ν /σ ) ξ (3) ( x , x , x ) to the power spectrum of biasedtracers which is calculated in [47]. Upon Fourier trans-formation, it reads as∆ P hh ( k ) = 4 f NL b δ c ( z ) α M ( k ) P φ ( k ) . (A15)We follow the steps outlined above to calculate the con-tribution from the second order trispectrum (A14). Af-ter some algebra, the Fourier transform of the four-pointcorrelations of the density field can be written down as9 Z d r ξ (4) ( x , x , x , x ) e − i k · r = 8 f α M ( k ) P φ ( k ) Z d k (2 π ) α M ( k ) P φ ( k ) Z d k (2 π ) α M ( k ) P φ ( k ) α M (cid:0) | k + k | (cid:1) × ( P φ ( k ) P φ ( k ) " P φ (cid:0) | k + k | (cid:1) P φ ( k ) + P φ (cid:0) | k + k | (cid:1) P φ ( k ) + P φ (cid:0) | k + k | (cid:1) P φ ( k )+ P φ (cid:0) | k + k | (cid:1) P φ (cid:0) | k + k | (cid:1) P φ ( k ) P φ ( k ) + P φ (cid:0) | k + k | (cid:1) P φ (cid:0) | k + k | (cid:1) P φ ( k ) P φ ( k ) ) (A16)and Z d r ξ (4) ( x , x , x , x ) e − i k · r = 4 f Z d k (2 π ) α M ( k ) α M (cid:0) | k + k | (cid:1) P φ ( k ) Z d k (2 π ) α M ( k ) α M (cid:0) | k + k | (cid:1) P φ ( k ) × ( P φ (cid:0) | k + k | (cid:1) " P φ (cid:0) | k + k | (cid:1) P φ ( k ) + P φ (cid:0) | k + k | (cid:1) P φ (cid:0) | k + k | (cid:1) P φ ( k ) P φ ( k ) + P φ ( k ) " P φ (cid:0) | k + k | (cid:1) P φ ( k ) + P φ (cid:0) | k + k | (cid:1) P φ (cid:0) | k + k | (cid:1) P φ ( k ) P φ ( k ) +2 P φ (cid:0) | k + k | (cid:1) (cid:20) P φ ( k ) P φ ( k ) + P φ ( k ) P φ ( k ) (cid:21)(cid:27) . (A17)In the large-scale limit k → 0, these expressions asymptote to Z d r ξ (4) ( x , x , x , x ) e − i k · r ≈ f σ S (1)3 ( M ) α M ( k ) P φ ( k ) + 8 f α M ( k ) Z d k (2 π ) α M ( k ) P φ ( k ) × Z d k (2 π ) α M ( k ) α M ( k ) (cid:2) P φ ( k ) P φ ( k ) + P φ ( k ) P φ ( k ) + P φ ( k ) (cid:3) (A18) Z d r ξ (4) ( x , x , x , x ) e − i k · r ≈ f σ P φ ( k ) + 16 f Z d k (2 π ) α M ( k ) Z d k (2 π ) α M ( k ) P φ ( k ) P φ ( k ) × (cid:2) P φ ( k ) + P φ ( k ) (cid:3) . (A19)Ignoring the second piece in the right-hand side of Eqs(A18) and (A19), the non-Gaussian correction to the halopower spectrum is then easily recast as Eq. (44). APPENDIX B: SIGNAL-TO-NOISE FOR THECMB TRISPECTRUM The formalism for the CMB trispectrum has been es-tablished in [29, 30]. The invariance of the 4-point har-monic function of the CMB temperature anisotropy fieldunder the 4! permutations of the coefficients a m i l i imposesconstraints on the CMB trispectrum T l l l l ( L ) which can be enforced by defining T l l l l ( L ) = P l l l l ( L ) + (2 L + 1) X L ′ (cid:20) ( − l + l × (cid:26) l l Ll l L ′ (cid:27) P l l l l ( L ′ ) + ( − L + L ′ × (cid:26) l l Ll l L ′ (cid:27) P l l l l ( L ′ ) (cid:21) , (B1)where curly brackets are Wigner-6j symbols, P l l l l ( L ) = t l l l l ( L ) + ( − P i l i t l l l l ( L ) + ( − L + l + l × t l l l l ( L ) + ( − L + l + l t l l l l ( L ) , (B2)and the reduced trispectrum t l l l l ( L ) is symmetric underthe exchange of its upper and lower indices and fullycharacterises the model.The expansion coefficients a ml are related to the pri-mordial curvature perturbation Φ( k ) through a ml = 4 π ( − i ) l Z d k (2 π ) Φ( k ) g T l ( k ) Y m⋆l ( ˆk ) , (B3)0where g Tl ( k ) is the radiation transfer function. The re-duced trispectrum can be calculated from this relationonce the four-point function T Φ ( k , k , k , k ) is speci-fied. For a local cubic coupling g NL φ , t l l l l ( L ) = Z ∞ dr r β l ( r ) β l ( r ) h l Ll h l Ll × [ µ l ( r ) β l ( r ) + β l ( r ) µ l ( r )] (B4)with β l ( r ) = 2 π Z ∞ dk k P φ ( k ) g Tl ( k ) j l ( kr ) (B5) µ l ( r ) = 2 π Z ∞ dk k g NL g Tl ( k ) j l ( kr ) , (B6)and h l Ll = 1 √ π π l Ll (cid:18) l l L (cid:19) . (B7)We also use the notation π l ··· lj = q (2 l + 1) × · · · × (2 l j + 1) . (B8)Note that most of the contribution to t l l l l ( L ) comes froma small volume centred at the comoving distance r ⋆ to thesurface of last scattering.The Sachs-Wolfe approximation g Tl ( k ) ≈ − j l ( kr ⋆ ) / l ≪ 100 provides a useful order-of-magnitude estimate [7, 30, 31, 83]. In this limit, wecan approximate µ l ( r ) as − g NL δ D ( r − r ⋆ ) / (3 r ⋆ ) since weassume g NL independent of wavenumber. Hence, the re-duced trispectrum simplifies to t l l l l ( L ) ≈ g NL C SW l C SW l (cid:0) C SW l + C SW l (cid:1) h l Ll h l Ll , (B9)Inserting this expression successively into eqs (B2) and(B1), the CMB trispectrum eventually reads as T l l l l ( L ) = 272 π g NL (2 L + 1) π l l l l (B10) × (cid:18) l l L (cid:19) (cid:18) l l L (cid:19) × h C SW l C SW l C SW l + (cyclic) i . where C SW l = 29 π Z ∞ dk k P φ ( k ) j l ( kr ⋆ ) ≈ πA φ l ( l + 1) . (B11)The last equality assumes a nearly scale-invariant spec-trum n s ≈ 1. The following relation between the Wigner-3j and 6j symbols (e.g., Appendix A of [29]), X l ′ (2 l ′ + 1)( − Σ+ l ′ − l − m − m ′ (cid:26) l l l l ′ l ′ l ′ (cid:27) × (cid:18) l l ′ l ′ m m ′ − m ′ (cid:19) (cid:18) l l ′ l ′ m m − m ′ (cid:19) = (cid:18) l l l m m − m (cid:19) (cid:18) l l ′ l ′ m m ′ − m ′ (cid:19) (B12) where Σ = l + l + l ′ + l ′ and the value of m is set by thetriangle condition, can be useful to derive Eq.(B10). Thesignal-to-noise summed up to multipole l max , Eq. (36),then becomes (cid:18) SN (cid:19) ( < l max ) = (cid:18) π (cid:19) g l max X l >l >l >l π l l l l × (cid:2) C SW l C SW l C SW l + (cyclic) (cid:3) C SW l C SW l C SW l C SW l (B13) × l max X L =0 (2 L + 1) (cid:18) l l L (cid:19) (cid:18) l l L (cid:19) . We can recast the sum over the diagonal modes L intoa manifestly symmetric form with the aid of the Gauntintegral12 Z +1 − dx P l ( x ) P l ( x ) P l ( x ) = (cid:18) l l l (cid:19) (B14)and the orthogonality relation ∞ X k =0 (2 k + 1) P k ( x ) P k ( y ) = 2 δ D ( x − y ) , (B15)where δ D is the Dirac delta. We find l max X L =0 (2 L + 1) (cid:18) l l L (cid:19) (cid:18) l l L (cid:19) (B16)= 12 Z +1 − dx P l ( x ) P l ( x ) P l ( x ) P l ( x ) . There is a strict equality because the Wigner-3j symbolsvanish for L > l + l . Eq. (37) for the signal-to-noise thenfollows by replacing the above equality into Eq.(B13) andsumming over all the 4! permutations of the quadruplet( l , l , l , l ). Although Eq. (37) becomes computation-ally expensive when l max ≫ 100 (because we are sum-ming over redundant configurations), we found that it isquite efficient for l max . l max by con-sidering only the contribution of the L = 1 mode inEq.(B13). Consequently, the product of the Wigner-3j symbols squared reduces to ∼ l l δ l − ,l δ l − ,l andyields ( S/N ) ∝ l . However, including all L modes asin Eq. (37) gives a steeper dependence, ( S/N ) ∝ l . (see Fig.9), due to the fact that the Wigner-3j symbolsdecay slowly with increasing L . This is quite apparent inthe classical limit l , l , L ≫ 1, where (cid:18) l l L (cid:19) (B17) ≈ π (cid:2) ( l + l ) − L (cid:3) − / (cid:2) L − ( l − l ) (cid:3) − / . /L for L ≫ l , l . By contrast, the second ordercontribution to the CMB trispectrum induced by thequadratic coupling f NL φ adds an additional multiplica-tive factor of ( C SW L ) in the summation over the L modes which increases the relative contribution of the low- L modes (since these now decay as L − ). This is the reasonwhy considering only L ≤ 10 modes as done in [31] stillprovides a good approximation to the signal-to-noise ofthe CMB trispectrum for the f NL model. [1] V. F. Mukhanov and G. V. Chibisov, Soviet Journal ofExperimental and Theoretical Physics Letters , 532(1981).[2] A. A. Starobinsky, Physics Letters B , 175 (1982).[3] S. W. Hawking, Physics Letters B , 295 (1982).[4] A. H. Guth and S. Pi, Physical Review Letters , 1110(1982).[5] T. J. Allen, B. Grinstein, and M. B. Wise, Physics LettersB , 66 (1987).[6] T. Falk, R. Rangarajan, and M. Srednicki,Phys. Rev. D , 4232 (1992).[7] A. Gangui, F. Lucchin, S. Matarrese, and S. Mollerach,ApJ , 447 (1994).[8] N. Bartolo, E. Komatsu, S. Matarrese, and A. Riotto,Phys. Rep. , 103 (2004).[9] V. Acquaviva, N. Bartolo, S. Matarrese, and A. Riotto,Nuclear Physics B , 119 (2003).[10] J. Maldacena, Journal of High Energy Physics , 13(2003).[11] E. Komatsu, A. Kogut, M. R. Nolta, C. L. Bennett, M.Halpern, G. Hinshaw, N. Jarosik, M. Limon, S. S. Meyer,L. Page, D. N. Spergel, G. S. Tucker, L. Verde, E. Wol-lack, and E. L. Wright, ”Astrophys. J. Supp.” , 119(2003).[12] P. Creminelli, L. Senatore, M. Zaldarriaga, and M.Tegmark, Journal of Cosmology and Astro-ParticlePhysics , 5 (2007).[13] E. Komatsu, J. Dunkley, M. R. Nolta, C. L. Bennett, B.Gold, G. Hinshaw, N. Jarosik, D. Larson, M. Limon, L.Page, D. N. Spergel, M. Halpern, R. S. Hill, A. Kogut,S. S. Meyer, G. S. Tucker, J. L. Weiland, E. Wollack, andE. L. Wright, ApJS , 330 (2009).[14] K. M. Smith, L. Senatore, and M. Zaldarriaga, ArXive-prints (2009).[15] A. Curto, E. Mart´ınez-Gonz´alez, P. Mukherjee, R. B.Barreiro, F. K. Hansen, M. Liguori, and S. Matarrese,Mon. Not. R. Astron. Soc. , 615 (2009).[16] D. H. Lyth, C. Ungarelli, and D. Wands, Phys. Rev. D ,023503 (2003).[17] N. Bartolo, S. Matarrese, and A. Riotto, Phys. Rev. D ,043503 (2004).[18] K. Enqvist and S. Nurmi, Journal of Cosmology andAstro-Particle Physics , 13 (2005).[19] K. A. Malik and D. H. Lyth, Journal of Cosmology andAstro-Particle Physics , 8 (2006).[20] M. Sasaki, J. V¨aliviita, and D. Wands, Phys. Rev. D ,103003 (2006).[21] K. Enqvist and T. Takahashi, Journal of Cosmology andAstro-Particle Physics , 12 (2008).[22] Q.-G. Huang and Y. Wang, Journal of Cosmology andAstro-Particle Physics , 25 (2008).[23] K. Ichikawa, T. Suyama, T. Takahashi, and M. Yam-aguchi, Phys. Rev. D , 023513 (2008).[24] P. Chingangbam and Q.-G. Huang, Journal of Cosmology and Astro-Particle Physics , 31 (2009).[25] Q.-G. Huang, Journal of Cosmology and Astro-ParticlePhysics , 5 (2008).[26] Q.-G. Huang, Journal of Cosmology and Astro-ParticlePhysics , 35 (2009).[27] C. T. Byrnes and G. Tasinato, Journal of Cosmology andAstro-Particle Physics , 16 (2009).[28] J.-L. Lehners and S. Renaux-Petel, Phys. Rev. D ,063503 (2009).[29] W. Hu, Phys. Rev. D , 083005 (2001).[30] T. Okamoto and W. Hu, Phys. Rev. D , 063008 (2002).[31] N. Kogo and E. Komatsu, Phys. Rev. D , 083007(2006).[32] E. Komatsu, ArXiv Astrophysics e-prints (2002).[33] M. Kunz, A. J. Banday, P. G. Castro, P. G. Ferreira, andK. M. G´orski, Astrophys. J. Lett. , L99 (2001).[34] F. Lucchin and S. Matarrese, ApJ , 535 (1988).[35] S. Colafrancesco, F. Lucchin, and S. Matarrese, ApJ ,3 (1989).[36] B. Grinstein and M. B. Wise, ApJ , 19 (1986).[37] S. Matarrese, F. Lucchin, and S. A. Bonometto, Astro-phys. J. Lett. , L21 (1986).[38] R. Scoccimarro, E. Sefusatti, and M. Zaldarriaga,Phys. Rev. D , 103513 (2004).[39] E. Sefusatti and E. Komatsu, Phys. Rev. D , 083004(2007).[40] E. Sefusatti, ArXiv e-prints (2009).[41] D. Jeong and E. Komatsu, ArXiv e-prints (2009).[42] M. Kamionkowski, L. Verde, and R. Jimenez, Journal ofCosmology and Astro-Particle Physics , 10 (2009).[43] T. Y. Lam and R. K. Sheth, Mon. Not. R. Astron. Soc. , 1743 (2009).[44] C. Hikage, E. Komatsu, and T. Matsubara, ApJ , 11(2006).[45] C. Hikage, P. Coles, M. Grossi, L. Moscardini, K. Dolag,E. Branchini, and S. Matarrese, Mon. Not. R. Astron.Soc. , 1613 (2008).[46] N. Dalal, O. Dor´e, D. Huterer, and A. Shirokov,Phys. Rev. D , 123514 (2008).[47] S. Matarrese and L. Verde, Astrophys. J. Lett. , L77(2008).[48] A. Slosar, C. Hirata, U. Seljak, S. Ho, and N. Padman-abhan, Journal of Cosmology and Astro-Particle Physics , 31 (2008).[49] N. Afshordi and A. J. Tolley, Phys. Rev. D , 123507(2008).[50] C. Carbone, L. Verde, and S. Matarrese, Astrophys. J.Lett. , L1 (2008).[51] U. Seljak, Physical Review Letters , 021302 (2009).[52] E. Sefusatti, M. Liguori, A. P. S. Yadav, M. G. Jackson,and E. Pajer, ArXiv e-prints (2009).[53] V. Desjacques, U. Seljak, and I. T. Iliev, Mon. Not. R.Astron. Soc. , 85 (2009).[54] A. Pillepich, C. Porciani, and O. Hahn, ArXiv e-prints (2008).[55] M. Grossi, L. Verde, C. Carbone, K. Dolag, E. Branchini,F. Iannuzzi, S. Matarrese, and L. Moscardini, ArXiv e-prints (2009).[56] U. Seljak and M. Zaldarriaga, ApJ , 437 (1996).[57] V. Springel, Mon. Not. R. Astron. Soc. , 1105 (2005).[58] Y. B. Zel’Dovich, A&A , 84 (1970).[59] J. N. Fry, ApJ , 499 (1984).[60] M. H. Goroff, B. Grinstein, S.-J. Rey, and M. B. Wise,ApJ , 6 (1986).[61] F. R. Bouchet, R. Juszkiewicz, S. Colombi, and R. Pellat,Astrophys. J. Lett. , L5 (1992).[62] F. Bernardeau, ApJ , 1 (1994).[63] R. Scoccimarro, Mon. Not. R. Astron. Soc. , 1097(1998).[64] P. McDonald, Phys. Rev. D , 123519 (2008).[65] S. R. Knollmann and A. Knebe, ArXiv e-prints (2009).[66] S. P. D. Gill, A. Knebe, and B. K. Gibson, Mon. Not. R.Astron. Soc. , 399 (2004).[67] M. S. Warren, P. J. Quinn, J. K. Salmon, and W. H.Zurek, ApJ , 405 (1992).[68] C. Lacey and S. Cole, Mon. Not. R. Astron. Soc. ,676 (1994).[69] V. R. Eke, S. Cole, and C. S. Frenk, Mon. Not. R. Astron.Soc. , 263 (1996).[70] G. L. Bryan and M. L. Norman, ApJ , 80 (1998).[71] W. H. Press and P. Schechter, ApJ , 425 (1974).[72] J. R. Bond, S. Cole, G. Efstathiou, and N. Kaiser,ApJ , 440 (1991).[73] M. Lo Verde, A. Miller, S. Shandera, and L. Verde, Jour-nal of Cosmology and Astro-Particle Physics , 14 (2008).[74] S. Matarrese, L. Verde, and R. Jimenez, ApJ , 10(2000). [75] M. Grossi, K. Dolag, E. Branchini, S. Matarrese, and L.Moscardini, Mon. Not. R. Astron. Soc. , 1261 (2007).[76] R. K. Sheth and G. Tormen, Mon. Not. R. Astron. Soc. , 119 (1999).[77] P. Valageas, ArXiv e-prints (2009).[78] T. Y. Lam and R. K. Sheth, ArXiv e-prints (2009).[79] T. Y. Lam, R. K. Sheth, and V. Desjacques, Mon. Not.R. Astron. Soc. , 1482 (2009).[80] A. Slosar, Journal of Cosmology and Astro-ParticlePhysics , 4 (2009).[81] U. Seljak, N. Hamaus, and V. Desjacques, Physical Re-view Letters , 091303 (2009).[82] P. McDonald and U. Seljak, ArXiv e-prints (2008).[83] E. Komatsu and D. N. Spergel, Phys. Rev. D , 063002(2001).[84] M. Grossi, E. Branchini, K. Dolag, S. Matarrese, and L.Moscardini, Mon. Not. R. Astron. Soc. , 438 (2008).[85] A. Taruya, K. Koyama, and T. Matsubara,Phys. Rev. D , 123534 (2008).[86] X. Kang, P. Norberg, and J. Silk, Mon. Not. R. Astron.Soc. , 343 (2007).[87] R. K. Sheth and G. Tormen, Mon. Not. R. Astron. Soc. , 61 (2002).[88] M. Maggiore and A. Riotto, ArXiv e-prints (2009).[89] Z. Luki´c, D. Reed, S. Habib, and K. Heitmann, ApJ ,217 (2009).[90] L. Boubekeur and D. H. Lyth, Phys. Rev. D , 021301(2006).[91] C. T. Byrnes, M. Sasaki, and D. Wands, Phys. Rev. D ,123519 (2006).[92] M.-X. Huang and G. Shiu, Phys. Rev. D74