# Primordial Non-Gaussianity in the Large-Scale Structure of  the Universe

Authors: Vincent Desjacques; Uroš Seljak
Journal: Advances in Astronomy
 (2010)

---

## Abstract

Primordial non-Gaussianity is a potentially powerful discriminant of the physical mechanisms that
generated the cosmological fluctuations observed today. Any detection of significant non-Gaussianity would thus have profound implications for our understanding of cosmic structure formation. The large-scale mass distribution in the Universe is a sensitive probe of the nature of initial conditions. Recent theoretical progress together with rapid developments in observational techniques will enable us to critically confront predictions of inflationary scenarios and set constraints as competitive as those from the Cosmic Microwave Background. In this paper, we review past and current efforts in the search for primordial non-Gaussianity in the large-scale structure of the Universe.


---

## Body

## 1. Introduction

In generic inflationary models based on the slow roll of a scalar field, primordial curvature perturbations are produced by the inflaton field as it slowly rolls down its potential [1–4]. Most of these scenarios predict a nearly scale-invariant spectrum of adiabatic curvature fluctuations, a relatively small amount of gravity waves, and tiny deviations from Gaussianity in the primeval distribution of curvature perturbations [5–7]. Although the latest measurements of the cosmic microwave background (CMB) anisotropies favor a slightly red power spectrum [8], no significant detection of a B-mode or some level of primordial non-Gaussianity (NG) has thus far been reported from CMB observations.While the presence of aB-mode can only be tested with CMB measurements, primordial deviations from Gaussianity can leave a detectable signature in the distribution of CMB anisotropies and in the large-scale structure (LSS) of the Universe. Until recently, it was widely accepted that measurement of the CMB furnished the best probe of primordial non-Gaussianity [9]. However, these conclusions did not take into account the anomalous scale dependence of the galaxy power spectrum and bispectrum arising from primordial NG of the local fNLloc type [10, 11]. These theoretical results, together with rapid developments in observational techniques, will provide large amount of LSS data to critically confront predictions of non-Gaussian models. In particular, galaxy clustering should provide independent constraints on the magnitude of primordial non-Gaussianity as competitive as those from the CMB and, in the long run, may even give the best constraints.The purpose of this paper is to provide an overview of the search for a primordial non-Gaussian signal in the large-scale structure. We will begin by briefly summarizing how non-Gaussianity arises in inflationary models (Section2). Next, we will discuss the impact of primordial non-Gaussianity on the mass distribution in the low-redshift Universe (Section 3). The main body of this paper is Section 4, where we describe in detail a number of methods exploiting the abundance and clustering properties of observed tracers of the LSS to constrain the amount of initial non-Gaussianity. We conclude with a discussion of present and forecasted constraints achievable with LSS surveys (Section 5).


## 2. Models and Observables

Because they assume (i) a single dynamical field (the inflation), (ii) canonical kinetic energy terms (i.e., perturbations propagate at the speed of light), (iii) slow roll (i.e., the timescale over which the inflaton field changes are much larger than the Hubble rate), and (iv) an initial vacuum state, single-field slow-roll models lead to a small level of primordial non-Gaussianity [6, 7, 12]. The lowest-order statistics sensitive to non-Gaussian features in the initial distribution of scalar perturbations Φ(x) (we will adopt the standard CMB convention in which Φ(x) is Bardeen's curvature perturbation in the matter era) is the 3-point function or bispectrum BΦ(k1,k2,k3), which is a function of any triangle k1+k2+k3=0 (as follows from statistical homogeneity which we assume throughout this paper). It has been recently shown that, in the squeezed limit k3≪k1≈k2, the bispectrum of any single-field slow-roll inflationary model asymptotes to the local shape (defined in (3)) [13–15]. The corresponding nonlinear parameter predicted by these models is(1)fNLloc=512(1-ns)≈0.017(singlefield),
where ns is the tilt or spectral index of the power spectrum PΦ(k), which is accurately measured to be ns≈0.960±0.013 [8]. Therefore, any robust measurement of fNLloc well above this level would thus rule out single-field slow-roll inflation.


### 2.1. The Shape of the Primordial Bispectrum

Large, potentially detectable amount of Gaussianity can be produced when at least one of the assumptions listed above is violated, that is, by multiple scalar fields [16, 17], nonlinearities in the relation between the primordial scalar perturbations and the inflaton field [7, 12], interactions of scalar fields [18], a modified dispersion relation, or a departure from the adiabatic Bunch-Davies ground state [19]. Generation of a large non-Gaussian signal is also expected at reheating [20] and in the ekpyrotic scenario [21]. Each of these physical mechanisms leaves a distinct signature in the primordial 3-point function BΦ(k1,k2,k3), a measurement of which would thus provide a wealth of information about the physics of primordial fluctuations. Although the configuration shape of the primordial bispectrum can be extremely complex in some models, there are broadly three classes of shape characterizing the local, equilateral, and folded type of primordial non-Gaussianity [22, 23]. The magnitude of each template “X” is controlled by a dimensionless nonlinear parameter fNLX which we seek to constrain using CMB or LSS observations (instead of attempting a model-independent measurement of BΦ).Any non-Gaussianity generated outside the horizon induces a 3-point function that is peaked on squeezed or collapsed triangles for realistic values of the scalar spectral index. The resulting non-Gaussianity depends only on the local value of Bardeen's curvature potential and can thus be conveniently parameterized up to third order by [7, 9, 12, 24](2)Φ(x)=ϕ(x)+fNLlocϕ2(x)+gNLlocϕ3(x),
where ϕ(x) is an isotropic Gaussian random field and fNLloc, gNLloc are dimensionless, phenomenological parameters. Since curvature perturbations are of magnitude 𝒪(10-5), the cubic-order correction should always be negligibly small compared to the quadratic one when 𝒪(fNLloc)~𝒪(gNLloc). However, this condition is not satisfied by some multifield inflationary models such as the curvaton scenario, in which a large gNLloc and a small fNLloc can be simultaneously produced [17]. The quadratic term generates the 3-point function at leading order(3)BΦloc(k1,k2,k3)=2fNLloc[Pϕ(k1)Pϕ(k2)+(cyc.)],
where (cyc.) denotes all cyclic permutations of the indices and Pϕ(k) is the power spectrum of the Gaussian part ϕ(x) of the Bardeen potential. The cubic-order terms generate a trispectrum TΦ(k1,k2,k3,k4) at leading order.Equilateral type of non-Gaussianity, which arises in inflationary models with higher-derivative operators such as the DBI model, is well described by the factorizable form [25](4)BΦeq(k1,k2,k3)=6fNLeq[-(Pϕ(k1)Pϕ(k2)+(cyc.))-2(Pϕ(k1)Pϕ(k2)Pϕ(k3))2/3+(Pϕ1/3(k1)Pϕ2/3(k2)Pϕ(k3)+(perm.))].
It can be easily checked that the signal is largest in the equilateral configurations k1≈k2≈k3 and suppressed in the squeezed limit k3≪k1≈k2. Note that, in single-field slow-roll inflation, the 3-point function is a linear combination of the local and equilateral shapes [13].As a third template, we consider the folded or flattened shape which is maximized fork2≈k3≈k1/2 [26]:(5)BΦfol(k1,k2,k3)=6fNLfol[(Pϕ(k1)Pϕ(k2)+(cyc.))+3(Pϕ(k1)Pϕ(k2)Pϕ(k3))2/3-(Pϕ1/3(k1)Pϕ2/3(k2)Pϕ(k3)+(perm.))],
and approximate the non-Gaussianity due to modification of the initial Bunch-Davies vacuum in canonical single-field action (the actual 3-point function is not factorizable). As in the previous example, BΦfol is suppressed in the squeezed configurations. Unlike BΦeq, however, the folded shape induces a scale-dependent bias at large scales (see Section 4.3).


### 2.2. Statistics of the Linear Mass Density Field

Bardeen's curvature potentialΦ(x) is related to the linear density perturbation δ0(k,z) at redshift z through(6)δ0(k,z)=ℳ(k,z)Φ(k),
where(7)ℳ(k,z)=23k2T(k)D(z)ΩmH02.
Here, T(k) is the matter transfer function normalized to unity as k→0, Ωm is the present-day matter density, and D(z) is the linear growth rate normalized to 1+z. n-point correlators of the linear mass density field can thus be written as(8)〈δ0(k)⋯δ0(kn)〉=(∏i=1nℳ(ki))〈Φ(k1)⋯Φ(kn)〉.
Smoothing inevitably arises when comparing observations of the large-scale structure with theoretical predictions from, for example, perturbation theory (PT), which are valid only in the weakly nonlinear regime [27], or from the spherical collapse model which ignores the strongly nonlinear internal dynamics of the collapsing regions [28, 29]. For this reason, we introduce the smoothed linear density field(9)δR(k,z)=ℳ(k,z)WR(k)Φ(k)≡ℳR(k,z)Φ(k),
where WR(k) is a (spherically symmetric) window function of characteristic radius R or mass scale M that smoothes out the small-scale nonlinear fluctuations. We will assume a top-hat filter throughout. Furthermore, since M and R are equivalent variables, we shall indistinctly use the notations δR and δM in what follows.


### 2.3. Topological Defects Models

In addition to inflationary scenarios, there is a whole class of models, known as topological defect models, in which cosmological fluctuations are sourced by an inhomogeneously distributed component which contributes a small fraction of the total energy momentum tensor [30, 31]. The density field is obtained as the convolution of a discrete set of points with a specific density profile. Defects are deeply rooted in particle physics as they are expected to form at a phase transition. Since the early Universe may have plausibly undergone several phase transitions, it is rather unlikely that no defects at all were formed. Furthermore, high-redshift tracers of the LSS may be superior to CMB at finding non-Gaussianity sourced by topological defects [32]. However, CMB observations already provide stringent limits on the energy density of a defect component [8], so we shall only minimally discuss the imprint of these scenarios in the large-scale structure. Phenomenological defect models are, for instance,(10)δ(x)=ϕ(x)+αNL(ϕ2(x)-〈ϕ2〉),
in which the initial matter density δ(x) (rather than the curvature perturbation Φ(x)) contains a term proportional to the square of a Gaussian scalar field ϕ(x) [9], or the χ2 model (also known as isocurvature CDM model) in which δ(x)∝ϕ2(x) [33].


## 2.1. The Shape of the Primordial Bispectrum

Large, potentially detectable amount of Gaussianity can be produced when at least one of the assumptions listed above is violated, that is, by multiple scalar fields [16, 17], nonlinearities in the relation between the primordial scalar perturbations and the inflaton field [7, 12], interactions of scalar fields [18], a modified dispersion relation, or a departure from the adiabatic Bunch-Davies ground state [19]. Generation of a large non-Gaussian signal is also expected at reheating [20] and in the ekpyrotic scenario [21]. Each of these physical mechanisms leaves a distinct signature in the primordial 3-point function BΦ(k1,k2,k3), a measurement of which would thus provide a wealth of information about the physics of primordial fluctuations. Although the configuration shape of the primordial bispectrum can be extremely complex in some models, there are broadly three classes of shape characterizing the local, equilateral, and folded type of primordial non-Gaussianity [22, 23]. The magnitude of each template “X” is controlled by a dimensionless nonlinear parameter fNLX which we seek to constrain using CMB or LSS observations (instead of attempting a model-independent measurement of BΦ).Any non-Gaussianity generated outside the horizon induces a 3-point function that is peaked on squeezed or collapsed triangles for realistic values of the scalar spectral index. The resulting non-Gaussianity depends only on the local value of Bardeen's curvature potential and can thus be conveniently parameterized up to third order by [7, 9, 12, 24](2)Φ(x)=ϕ(x)+fNLlocϕ2(x)+gNLlocϕ3(x),
where ϕ(x) is an isotropic Gaussian random field and fNLloc, gNLloc are dimensionless, phenomenological parameters. Since curvature perturbations are of magnitude 𝒪(10-5), the cubic-order correction should always be negligibly small compared to the quadratic one when 𝒪(fNLloc)~𝒪(gNLloc). However, this condition is not satisfied by some multifield inflationary models such as the curvaton scenario, in which a large gNLloc and a small fNLloc can be simultaneously produced [17]. The quadratic term generates the 3-point function at leading order(3)BΦloc(k1,k2,k3)=2fNLloc[Pϕ(k1)Pϕ(k2)+(cyc.)],
where (cyc.) denotes all cyclic permutations of the indices and Pϕ(k) is the power spectrum of the Gaussian part ϕ(x) of the Bardeen potential. The cubic-order terms generate a trispectrum TΦ(k1,k2,k3,k4) at leading order.Equilateral type of non-Gaussianity, which arises in inflationary models with higher-derivative operators such as the DBI model, is well described by the factorizable form [25](4)BΦeq(k1,k2,k3)=6fNLeq[-(Pϕ(k1)Pϕ(k2)+(cyc.))-2(Pϕ(k1)Pϕ(k2)Pϕ(k3))2/3+(Pϕ1/3(k1)Pϕ2/3(k2)Pϕ(k3)+(perm.))].
It can be easily checked that the signal is largest in the equilateral configurations k1≈k2≈k3 and suppressed in the squeezed limit k3≪k1≈k2. Note that, in single-field slow-roll inflation, the 3-point function is a linear combination of the local and equilateral shapes [13].As a third template, we consider the folded or flattened shape which is maximized fork2≈k3≈k1/2 [26]:(5)BΦfol(k1,k2,k3)=6fNLfol[(Pϕ(k1)Pϕ(k2)+(cyc.))+3(Pϕ(k1)Pϕ(k2)Pϕ(k3))2/3-(Pϕ1/3(k1)Pϕ2/3(k2)Pϕ(k3)+(perm.))],
and approximate the non-Gaussianity due to modification of the initial Bunch-Davies vacuum in canonical single-field action (the actual 3-point function is not factorizable). As in the previous example, BΦfol is suppressed in the squeezed configurations. Unlike BΦeq, however, the folded shape induces a scale-dependent bias at large scales (see Section 4.3).


## 2.2. Statistics of the Linear Mass Density Field

Bardeen's curvature potentialΦ(x) is related to the linear density perturbation δ0(k,z) at redshift z through(6)δ0(k,z)=ℳ(k,z)Φ(k),
where(7)ℳ(k,z)=23k2T(k)D(z)ΩmH02.
Here, T(k) is the matter transfer function normalized to unity as k→0, Ωm is the present-day matter density, and D(z) is the linear growth rate normalized to 1+z. n-point correlators of the linear mass density field can thus be written as(8)〈δ0(k)⋯δ0(kn)〉=(∏i=1nℳ(ki))〈Φ(k1)⋯Φ(kn)〉.
Smoothing inevitably arises when comparing observations of the large-scale structure with theoretical predictions from, for example, perturbation theory (PT), which are valid only in the weakly nonlinear regime [27], or from the spherical collapse model which ignores the strongly nonlinear internal dynamics of the collapsing regions [28, 29]. For this reason, we introduce the smoothed linear density field(9)δR(k,z)=ℳ(k,z)WR(k)Φ(k)≡ℳR(k,z)Φ(k),
where WR(k) is a (spherically symmetric) window function of characteristic radius R or mass scale M that smoothes out the small-scale nonlinear fluctuations. We will assume a top-hat filter throughout. Furthermore, since M and R are equivalent variables, we shall indistinctly use the notations δR and δM in what follows.


## 2.3. Topological Defects Models

In addition to inflationary scenarios, there is a whole class of models, known as topological defect models, in which cosmological fluctuations are sourced by an inhomogeneously distributed component which contributes a small fraction of the total energy momentum tensor [30, 31]. The density field is obtained as the convolution of a discrete set of points with a specific density profile. Defects are deeply rooted in particle physics as they are expected to form at a phase transition. Since the early Universe may have plausibly undergone several phase transitions, it is rather unlikely that no defects at all were formed. Furthermore, high-redshift tracers of the LSS may be superior to CMB at finding non-Gaussianity sourced by topological defects [32]. However, CMB observations already provide stringent limits on the energy density of a defect component [8], so we shall only minimally discuss the imprint of these scenarios in the large-scale structure. Phenomenological defect models are, for instance,(10)δ(x)=ϕ(x)+αNL(ϕ2(x)-〈ϕ2〉),
in which the initial matter density δ(x) (rather than the curvature perturbation Φ(x)) contains a term proportional to the square of a Gaussian scalar field ϕ(x) [9], or the χ2 model (also known as isocurvature CDM model) in which δ(x)∝ϕ2(x) [33].


## 3. Evolution of the Matter Density Field with Primordial NG

In this section, we summarize a number of results relative to the effect of primordial NG on the mass density field. These will be useful to understand the complexity that arises when considering biased tracers of the density field (see Section4).


### 3.1. Setting Up Non-Gaussian Initial Conditions

Investigating the impact of non-Gaussian initial conditions (ICs) on the large-scale structure traced by galaxies and so forth requires simulations large enough so that many long-wavelength modes are sampled. At the same time, the simulations should resolve dark-matter halos hosting luminous red galaxies (LRGs) or quasars (QSOs), so that one can construct halo samples whose statistical properties mimic as closely as possible those of the real data. This favors the utilization of pure N-body simulations for which a large dynamical range can be achieved.The evolution of the matter density field with primordial non-Gaussianity has been studied in series of large cosmological N-body simulations seeded with Gaussian and non-Gaussian initial conditions; see, for example, [11, 34–43]. For generic non-Gaussian (scalar) random fields, we face the problem of setting up numerical simulations with a prescribed correlation structure [44]. While an implementation of the equilateral and folded bispectrum shape requires the calculation of several computationally demanding convolutions, the operation is straightforward for primordial NG described by a local mapping such as the χ2 or the fNLloc model. In the latter case, the local transformation Φ=ϕ+fNLlocϕ2 is performed before multiplication by the matter transfer function T(k) (computed with publicly available Boltzmann codes [45, 46]). The (dimensionless) power spectrum of the Gaussian part ϕ(x) of the Bardeen potential is the usual power-law Δϕ2(k)≡k3Pϕ(k)/(2π2)=Aϕ(k/k0)ns-1. Unless otherwise stated, we shall assume a normalization Aϕ=7.96×10-10 at the pivot point k0=0.02Mpc-1. To date, essentially all numerical studies of structure formation with inflationary non-Gaussianity have implemented the local shape solely, so we will focus on this model in what follows.Non-Gaussian corrections to the primordial curvature perturbations can renormalize the input (unrenormalized) power spectrum of fluctuations used to seed the simulations [47]. For the local fNLloc model with |fNLloc|≲100, renormalization effects are unlikely to be noticeable due to the limited dynamical range of current cosmological simulations. However, they can be significant, for example, in simulations of a local cubic coupling gNLlocϕ3 with a large primordial trispectrum [48]. The cubic-order term gNLlocϕ3 renormalizes the amplitude Aϕ of the power spectrum of initial curvature perturbations to Aϕ→Aϕ+6gNLloc〈ϕ2〉, where(11)〈ϕ2〉=∫d3k(2π)3Pϕ(k).
For scale-invariant initial conditions, 〈ϕ2〉 has a logarithmic divergence at large and small scales. In practice, however, the finite box size and the resolution of the simulations naturally furnish a low- and high-k cutoff. The effective correction to the amplitude of density fluctuations δσ8 in the gNLlocϕ3 model thus is(12)δσ8=3gNLloc(Lk02π)1-ns[1-Nns-1]Aϕ1-ns,
where N is the number of mesh points along one dimension and L is the simulation box length. For gNLloc=106, L=1.5h-1Gpc, and N=1024, for instance, we obtain δσ8≈0.015.To generate the initial particle distribution, the Zeldovich approximation is commonly used instead of the exact gravitational dynamics. This effectively corresponds to starting from non-Gaussian initial conditions [49]. Since the transition to the true gravitational dynamics proceeds rather gradually [50], one should ensure that the initial expansion factor is much smaller than that of the outputs analyzed. Alternatively, it is possible to generate more accurate ICs based on second-order Lagrangian perturbation theory (2LPT) [51]. At fixed initial expansion factor, they reduce transients such that the true dynamics is recovered more rapidly [52].


### 3.2. Mass Density Probability Distribution

In the absence of primordial NG, the probability distribution function (PDF) of the initial smoothed density field (the probability that a randomly placed cell of volumeV has some specific density) is Gaussian. Namely, all normalized or reduced smoothed cumulants SJ of order J≥3 are zero. An obvious signature of primordial NG would thus be an initially nonvanishing skewness S3=〈δR3〉c/〈δR2〉2 or kurtosis S4=〈δR4〉c/〈δR2〉3-3/〈δR2〉 [36, 53, 54]. Here, the subscript c denotes the connected piece of the n-point moment that cannot be simplified into a sum over products of lower-order moments. At third order, for instance, the cumulant of the smoothed density field is an integral of the 3-point function:(13)〈δR3〉c=∫d3k1(2π)3∫d3k2(2π)3∫d3k3(2π)3BR(k1,k2,k3,z),
where(14)BR(k1,k2,k3,z)=ℳR(k1,z)ℳR(k2,z)ℳR(k3,z)×BΦ(k1,k2,k3)
is the bispectrum of the smoothed linear density fluctuations at redshift z. Note that, owing to S3(R,z)∝D(z)-1, the product σS3(R) does not depend on redshift. Over the range of scale accessible to LSS observations, σS3(R) is a monotonic decreasing function of R that is of magnitude ~10-4 for the local, equilateral, and folded templates discussed above (Figure 1). Strictly speaking, all reduced moments should be specified to fully characterize the density PDF, but a reasonable description of the density distribution can be achieved with moments up to the fourth order.Figure 1
SkewnessσS3(1)(R) of the smoothed density field in unit of fNLX for the local, equilateral, and folded bispectrum shapes. The skewness for the equilateral and folded templates is a factor of ~3 smaller than in the local model. In any case, this implies that |σRS3(R)|≪1 on the scales probed by the large-scale structure for realistic values of the nonlinear coupling parameter, |fNLX|≲100. The shaded regions approximately indicate the range of scales probed by various LSS tracers. For the galaxy power spectrum and bispectrum, the upper limit sensitively depends upon the surveyed volume.Numerical and analytic studies generally find that a density PDF initially skewed towards positive values produces more overdense regions while a negatively skewed distribution produces larger voids. Gravitational instabilities also generate a positive skewness in the density field, reflecting the fact that the evolved density distribution exhibits an extended tail towards large overdensities [49, 55–59]. This gravitationally induced signal eventually dominates the primordial contribution such that, at fixed normalization amplitude, non-Gaussian models deviate more strongly from the Gaussian paradigm at high redshift. The time evolution of the normalized cumulants SJ can be worked out for generic Gaussian and non-Gaussian ICs using, for example, PT or the spherical collapse approximation. For Gaussian ICs, PT predicts the normalized cumulants to be time independent to the lowest nonvanishing order, with a skewness S3≈34/7, whereas, for non-Gaussian ICs, the linear contribution to the cumulants decays as SJ(R,z)=SJ(R,∞)/DJ-2(z) [60, 61].The persistence of the primordial hierarchical amplitudeSJ(R,∞) sensitively depends upon the magnitude of SN, N≥J, relative to unity. For example, an initially large nonvanishing kurtosis could source skewness with a time-dependence and amplitude similar to those induced by nonlinear gravitational evolution [60]. Although there is an infinite class of non-Gaussian models, we can broadly divide them into weakly and strongly non-Gaussian. In weak NG models, the primeval signal in the normalized cumulants is rapidly obliterated by gravity-induced non-Gaussianity. This is the case of hierarchical scaling models where n-point correlation functions satisfy ξn∝ξ2n-1 with ξ2≪1 at large-scales. By contrast, strongly NG initial conditions dominate the evolution of the normalized cumulants. This occurs when the hierarchy of correlation functions obeys the dimensional scaling ξn∝ξ2n/2, which arises in the particular case of χ2 initial conditions [62] or in defect models such as texture [37, 63, 64]. These scaling laws have been successfully confronted with numerical investigations of the evolution of cumulants [37, 38]. We note that the scaling of the contribution induced by gravity is, however, different for the kurtosis [65], suggesting that the latter is a better probe of the nature of initial conditions.Although gravitational clustering tends to erase the memory of initial conditions, numerical simulations of non-Gaussian initial conditions show that the occurrence of highly underdense and overdense regions is very sensitive to the presence of primordial NG. In fact, the imprint of primordial NG is best preserved in the negative tail of the PDFP(ρR) of the evolved (and smoothed) density field ρR [40]. A satisfactory description of this effect can be obtained from an Edgeworth expansion of the initial smoothed overdensity field [66]. At high densities ρR≫1, the non-Gaussian modification approximately scales as ρR3/5 whereas, at low densities ρR≃0, the deviation is steeper, ρR6/5. Taking into account the weak-scale dependence of σS3(R) further enhances this asymmetry.


### 3.3. Power Spectrum and Bispectrum

Primordial-non-Gaussianity also imprints a signature in Fourier-space statistics of the matter density field. Positive values offNLloc tend to increase the small-scale matter power spectrum Pδ(k) [10, 40, 67] and the large-scale matter bispectrum Bδ(k1,k2,k3) [10, 68]. In the weakly nonlinear regime where 1-loop PT applies, the Fourier mode of the density field for growing-mode initial conditions reads [56, 69](15)δ(k,z)=δ0(k,z)+1(2π)3∫d3q1d3q2δD(k-q1-q2)×F2(q1,q2)δ0(q1,z)δ0(q2,z).
The kernel F2(k1,k2)=5/7+μ(k1/k2+k2/k1)/2+2μ2/7, where μ is the cosine of the angle between k1 and k2, describes the nonlinear 2nd-order evolution of the density field. It is nearly independent of Ωm and ΩΛ and vanishes in the (squeezed) limit k1=-k2 as a consequence of the causality of gravitational instability. At 1-loop PT, (15) generates the mass power spectrum(16)Pδ(k,z)=PδG(k,z)+ΔPδNG(k,z)=P0(k,z)+[P(22)(k,z)+P(13)(k,z)]+ΔPδNG(k,z).
Here, P0(k,z) is the linear matter power spectrum at redshift z, P(22) and P(13) are the standard one-loop contributions in the case of Gaussian ICs [69, 70], and(17)ΔPδNG(k,z)=2∫d3q(2π)3F2(q,k-q)B0(-k,q,k-q,z)
is the correction due to primordial NG [67]. This terms scales as ∝D3(z), so the effect of non-Gaussianity is largest at low redshift. Moreover, because F2(k1,k2) vanishes in the squeezed limit, (17) is strongly suppressed at small wavenumbers, even in the local fNLloc model for which B0(-k,q,k-q,z) is maximized in the same limit (i.e., |k|→0). For fNLloc~𝒪(102), the magnitude of this correction is at a percent level in the weakly nonlinear regime k≲0.1hMpc-1 [41, 42, 71], in good agreement with the measurements (see Figure 2). Extensions of the renormalization group description of dark-matter clustering [72] to non-Gaussian initial density, and velocity perturbations can improve the agreement up to wavenumbers k≲0.25hMpc-1 [73, 74].Figure 2
Non-Gaussian fractional correctionβδ(k,z)=ΔPδNG(k,z)/PδG(k,z) to the matter power spectrum that originates from primordial non-Gaussianity of the local type. Results are shown at redshift z=0, 0.5, 1, and 2 for fNLloc=+100 (filled symbols) and fNLloc=-100 (empty symbols). The solid curves indicate the prediction from a 1-loop perturbative expansion.It is also instructive to compare measurements of the matter bispectrumBδ(k1,k2,k3) with perturbative predictions. To second order in PT, the matter bispectrum is the sum of a primordial contribution and two terms induced by gravitational instability [56, 75] (we will henceforth omit the explicit z-dependence for brevity):(18)Bδ(k1,k2,k3)=B0(k1,k2,k3)+[2F2(k1,k2)P0(k1)P0(k2)+(cyc.)]+∫d3q(2π)3[F2(q,k3-q)×T0(q,k3-q,k1,k2)+(cyc.)],
where T0(k1,k2,k3,k4) is the primordial trispectrum of the density field. A similar expression can also be derived for the matter trispectrum, which turns out to be less sensitive to gravitationally induced nonlinearities [76]. The reduced bispectrum Q3 is conveniently defined as(19)Q3(k1,k2,k3)=Bδ(k1,k2,k3)[Pδ(k1)Pδ(k2)+cyclic].
For Gaussian initial conditions, Q3 is independent of time and, at tree-level PT, is constant and equal to Q3(k,k,k)=4/7 for equilateral configurations [56]. For general triangles, moreover, it approximately retains this simple behavior, with a dependence on triangle shape through F2(k1,k2) [10]. Figure 3 illustrates the effect of primordial NG of the local fNLloc type on the shape dependence of Q3 for a particular set of triangle configurations. As can be seen, the inclusion of 1-loop corrections greatly improves the agreement with the numerical data [77]. An important feature that is not apparent in Figure 3 is the fact that the primordial part to the reduced-matter bispectrum scales as Q3∝1/ℳR(k) for approximately equilateral triangles (and under the assumption that fNLloc is scale independent) [10]. This anomalous scaling considerably raises the ability of the matter bispectrum to constrain primordial NG of the local fNLloc type. Unfortunately, neither the matter bispectrum nor the power spectrum is directly observable with the large-scale structure of the Universe. Temperature anisotropies in the redshifted 21 cm background from the prereionization epoch could in principle furnish a direct measurement of these quantities [78–80], but foreground contamination may severely hamper any detection. Weak lensing is another direct probe of the dark matter, although we can only observe it in projection along the line of sight [81].Figure 3
Reduced matter bispectrumQ3 as a function of the angle θ between k1 and k2 for a fixed k1=0.094hMpc-1 and k2=1.5k1. The panels show ratios between the non-Gaussian and Gaussian Q3 for fNLloc=+100 (top) and -100 (bottom). Dashed lines correspond to tree-level PT while continuous line indicates the 1-loop PT prediction.As we will see shortly, however, this large-scale anomalous scaling is also present in the bispectrum and power spectrum of observable tracers of the large-scale structure such as galaxies. It is this unique signature that will make future all-sky LSS surveys competitive with CMB experiments.


### 3.4. Velocities

Primordial non-Gaussianity also leaves a signature in the large-scale coherent bulk motions which, in the linear regime, are directly related to the linear density field [55]. The various non-Gaussian models considered by [82] tend to have lower-velocity dispersion and bulk flow than fiducial Gaussian model, regardless of the sign of the skewness. However, while the probability distribution of velocity components is sensitive to primordial NG of the local type, in defect models it can be very close to Gaussian, even when the density field is strongly non-Gaussian, as a consequence of the central limit theorem [83, 84]. In this regard, correlation of velocity differences could provide a better test of non-Gaussian initial conditions [85].To measure peculiar velocities, one must subtract the Hubble flow from the observed redshift. This requires an estimate of the distance which is only available for nearby galaxies and clusters (although, e.g., the kinetic Sunyaev-Zel'dovich (kSZ) effect could be used to measure the bulk motions of distant galaxy clusters [86]). So far, measurements of the local galaxy density and velocity fields [87] as well as reconstruction of the initial density PDF from galaxy density and velocity data [88] are consistent with Gaussian initial conditions.


## 3.1. Setting Up Non-Gaussian Initial Conditions

Investigating the impact of non-Gaussian initial conditions (ICs) on the large-scale structure traced by galaxies and so forth requires simulations large enough so that many long-wavelength modes are sampled. At the same time, the simulations should resolve dark-matter halos hosting luminous red galaxies (LRGs) or quasars (QSOs), so that one can construct halo samples whose statistical properties mimic as closely as possible those of the real data. This favors the utilization of pure N-body simulations for which a large dynamical range can be achieved.The evolution of the matter density field with primordial non-Gaussianity has been studied in series of large cosmological N-body simulations seeded with Gaussian and non-Gaussian initial conditions; see, for example, [11, 34–43]. For generic non-Gaussian (scalar) random fields, we face the problem of setting up numerical simulations with a prescribed correlation structure [44]. While an implementation of the equilateral and folded bispectrum shape requires the calculation of several computationally demanding convolutions, the operation is straightforward for primordial NG described by a local mapping such as the χ2 or the fNLloc model. In the latter case, the local transformation Φ=ϕ+fNLlocϕ2 is performed before multiplication by the matter transfer function T(k) (computed with publicly available Boltzmann codes [45, 46]). The (dimensionless) power spectrum of the Gaussian part ϕ(x) of the Bardeen potential is the usual power-law Δϕ2(k)≡k3Pϕ(k)/(2π2)=Aϕ(k/k0)ns-1. Unless otherwise stated, we shall assume a normalization Aϕ=7.96×10-10 at the pivot point k0=0.02Mpc-1. To date, essentially all numerical studies of structure formation with inflationary non-Gaussianity have implemented the local shape solely, so we will focus on this model in what follows.Non-Gaussian corrections to the primordial curvature perturbations can renormalize the input (unrenormalized) power spectrum of fluctuations used to seed the simulations [47]. For the local fNLloc model with |fNLloc|≲100, renormalization effects are unlikely to be noticeable due to the limited dynamical range of current cosmological simulations. However, they can be significant, for example, in simulations of a local cubic coupling gNLlocϕ3 with a large primordial trispectrum [48]. The cubic-order term gNLlocϕ3 renormalizes the amplitude Aϕ of the power spectrum of initial curvature perturbations to Aϕ→Aϕ+6gNLloc〈ϕ2〉, where(11)〈ϕ2〉=∫d3k(2π)3Pϕ(k).
For scale-invariant initial conditions, 〈ϕ2〉 has a logarithmic divergence at large and small scales. In practice, however, the finite box size and the resolution of the simulations naturally furnish a low- and high-k cutoff. The effective correction to the amplitude of density fluctuations δσ8 in the gNLlocϕ3 model thus is(12)δσ8=3gNLloc(Lk02π)1-ns[1-Nns-1]Aϕ1-ns,
where N is the number of mesh points along one dimension and L is the simulation box length. For gNLloc=106, L=1.5h-1Gpc, and N=1024, for instance, we obtain δσ8≈0.015.To generate the initial particle distribution, the Zeldovich approximation is commonly used instead of the exact gravitational dynamics. This effectively corresponds to starting from non-Gaussian initial conditions [49]. Since the transition to the true gravitational dynamics proceeds rather gradually [50], one should ensure that the initial expansion factor is much smaller than that of the outputs analyzed. Alternatively, it is possible to generate more accurate ICs based on second-order Lagrangian perturbation theory (2LPT) [51]. At fixed initial expansion factor, they reduce transients such that the true dynamics is recovered more rapidly [52].


## 3.2. Mass Density Probability Distribution

In the absence of primordial NG, the probability distribution function (PDF) of the initial smoothed density field (the probability that a randomly placed cell of volumeV has some specific density) is Gaussian. Namely, all normalized or reduced smoothed cumulants SJ of order J≥3 are zero. An obvious signature of primordial NG would thus be an initially nonvanishing skewness S3=〈δR3〉c/〈δR2〉2 or kurtosis S4=〈δR4〉c/〈δR2〉3-3/〈δR2〉 [36, 53, 54]. Here, the subscript c denotes the connected piece of the n-point moment that cannot be simplified into a sum over products of lower-order moments. At third order, for instance, the cumulant of the smoothed density field is an integral of the 3-point function:(13)〈δR3〉c=∫d3k1(2π)3∫d3k2(2π)3∫d3k3(2π)3BR(k1,k2,k3,z),
where(14)BR(k1,k2,k3,z)=ℳR(k1,z)ℳR(k2,z)ℳR(k3,z)×BΦ(k1,k2,k3)
is the bispectrum of the smoothed linear density fluctuations at redshift z. Note that, owing to S3(R,z)∝D(z)-1, the product σS3(R) does not depend on redshift. Over the range of scale accessible to LSS observations, σS3(R) is a monotonic decreasing function of R that is of magnitude ~10-4 for the local, equilateral, and folded templates discussed above (Figure 1). Strictly speaking, all reduced moments should be specified to fully characterize the density PDF, but a reasonable description of the density distribution can be achieved with moments up to the fourth order.Figure 1
SkewnessσS3(1)(R) of the smoothed density field in unit of fNLX for the local, equilateral, and folded bispectrum shapes. The skewness for the equilateral and folded templates is a factor of ~3 smaller than in the local model. In any case, this implies that |σRS3(R)|≪1 on the scales probed by the large-scale structure for realistic values of the nonlinear coupling parameter, |fNLX|≲100. The shaded regions approximately indicate the range of scales probed by various LSS tracers. For the galaxy power spectrum and bispectrum, the upper limit sensitively depends upon the surveyed volume.Numerical and analytic studies generally find that a density PDF initially skewed towards positive values produces more overdense regions while a negatively skewed distribution produces larger voids. Gravitational instabilities also generate a positive skewness in the density field, reflecting the fact that the evolved density distribution exhibits an extended tail towards large overdensities [49, 55–59]. This gravitationally induced signal eventually dominates the primordial contribution such that, at fixed normalization amplitude, non-Gaussian models deviate more strongly from the Gaussian paradigm at high redshift. The time evolution of the normalized cumulants SJ can be worked out for generic Gaussian and non-Gaussian ICs using, for example, PT or the spherical collapse approximation. For Gaussian ICs, PT predicts the normalized cumulants to be time independent to the lowest nonvanishing order, with a skewness S3≈34/7, whereas, for non-Gaussian ICs, the linear contribution to the cumulants decays as SJ(R,z)=SJ(R,∞)/DJ-2(z) [60, 61].The persistence of the primordial hierarchical amplitudeSJ(R,∞) sensitively depends upon the magnitude of SN, N≥J, relative to unity. For example, an initially large nonvanishing kurtosis could source skewness with a time-dependence and amplitude similar to those induced by nonlinear gravitational evolution [60]. Although there is an infinite class of non-Gaussian models, we can broadly divide them into weakly and strongly non-Gaussian. In weak NG models, the primeval signal in the normalized cumulants is rapidly obliterated by gravity-induced non-Gaussianity. This is the case of hierarchical scaling models where n-point correlation functions satisfy ξn∝ξ2n-1 with ξ2≪1 at large-scales. By contrast, strongly NG initial conditions dominate the evolution of the normalized cumulants. This occurs when the hierarchy of correlation functions obeys the dimensional scaling ξn∝ξ2n/2, which arises in the particular case of χ2 initial conditions [62] or in defect models such as texture [37, 63, 64]. These scaling laws have been successfully confronted with numerical investigations of the evolution of cumulants [37, 38]. We note that the scaling of the contribution induced by gravity is, however, different for the kurtosis [65], suggesting that the latter is a better probe of the nature of initial conditions.Although gravitational clustering tends to erase the memory of initial conditions, numerical simulations of non-Gaussian initial conditions show that the occurrence of highly underdense and overdense regions is very sensitive to the presence of primordial NG. In fact, the imprint of primordial NG is best preserved in the negative tail of the PDFP(ρR) of the evolved (and smoothed) density field ρR [40]. A satisfactory description of this effect can be obtained from an Edgeworth expansion of the initial smoothed overdensity field [66]. At high densities ρR≫1, the non-Gaussian modification approximately scales as ρR3/5 whereas, at low densities ρR≃0, the deviation is steeper, ρR6/5. Taking into account the weak-scale dependence of σS3(R) further enhances this asymmetry.


## 3.3. Power Spectrum and Bispectrum

Primordial-non-Gaussianity also imprints a signature in Fourier-space statistics of the matter density field. Positive values offNLloc tend to increase the small-scale matter power spectrum Pδ(k) [10, 40, 67] and the large-scale matter bispectrum Bδ(k1,k2,k3) [10, 68]. In the weakly nonlinear regime where 1-loop PT applies, the Fourier mode of the density field for growing-mode initial conditions reads [56, 69](15)δ(k,z)=δ0(k,z)+1(2π)3∫d3q1d3q2δD(k-q1-q2)×F2(q1,q2)δ0(q1,z)δ0(q2,z).
The kernel F2(k1,k2)=5/7+μ(k1/k2+k2/k1)/2+2μ2/7, where μ is the cosine of the angle between k1 and k2, describes the nonlinear 2nd-order evolution of the density field. It is nearly independent of Ωm and ΩΛ and vanishes in the (squeezed) limit k1=-k2 as a consequence of the causality of gravitational instability. At 1-loop PT, (15) generates the mass power spectrum(16)Pδ(k,z)=PδG(k,z)+ΔPδNG(k,z)=P0(k,z)+[P(22)(k,z)+P(13)(k,z)]+ΔPδNG(k,z).
Here, P0(k,z) is the linear matter power spectrum at redshift z, P(22) and P(13) are the standard one-loop contributions in the case of Gaussian ICs [69, 70], and(17)ΔPδNG(k,z)=2∫d3q(2π)3F2(q,k-q)B0(-k,q,k-q,z)
is the correction due to primordial NG [67]. This terms scales as ∝D3(z), so the effect of non-Gaussianity is largest at low redshift. Moreover, because F2(k1,k2) vanishes in the squeezed limit, (17) is strongly suppressed at small wavenumbers, even in the local fNLloc model for which B0(-k,q,k-q,z) is maximized in the same limit (i.e., |k|→0). For fNLloc~𝒪(102), the magnitude of this correction is at a percent level in the weakly nonlinear regime k≲0.1hMpc-1 [41, 42, 71], in good agreement with the measurements (see Figure 2). Extensions of the renormalization group description of dark-matter clustering [72] to non-Gaussian initial density, and velocity perturbations can improve the agreement up to wavenumbers k≲0.25hMpc-1 [73, 74].Figure 2
Non-Gaussian fractional correctionβδ(k,z)=ΔPδNG(k,z)/PδG(k,z) to the matter power spectrum that originates from primordial non-Gaussianity of the local type. Results are shown at redshift z=0, 0.5, 1, and 2 for fNLloc=+100 (filled symbols) and fNLloc=-100 (empty symbols). The solid curves indicate the prediction from a 1-loop perturbative expansion.It is also instructive to compare measurements of the matter bispectrumBδ(k1,k2,k3) with perturbative predictions. To second order in PT, the matter bispectrum is the sum of a primordial contribution and two terms induced by gravitational instability [56, 75] (we will henceforth omit the explicit z-dependence for brevity):(18)Bδ(k1,k2,k3)=B0(k1,k2,k3)+[2F2(k1,k2)P0(k1)P0(k2)+(cyc.)]+∫d3q(2π)3[F2(q,k3-q)×T0(q,k3-q,k1,k2)+(cyc.)],
where T0(k1,k2,k3,k4) is the primordial trispectrum of the density field. A similar expression can also be derived for the matter trispectrum, which turns out to be less sensitive to gravitationally induced nonlinearities [76]. The reduced bispectrum Q3 is conveniently defined as(19)Q3(k1,k2,k3)=Bδ(k1,k2,k3)[Pδ(k1)Pδ(k2)+cyclic].
For Gaussian initial conditions, Q3 is independent of time and, at tree-level PT, is constant and equal to Q3(k,k,k)=4/7 for equilateral configurations [56]. For general triangles, moreover, it approximately retains this simple behavior, with a dependence on triangle shape through F2(k1,k2) [10]. Figure 3 illustrates the effect of primordial NG of the local fNLloc type on the shape dependence of Q3 for a particular set of triangle configurations. As can be seen, the inclusion of 1-loop corrections greatly improves the agreement with the numerical data [77]. An important feature that is not apparent in Figure 3 is the fact that the primordial part to the reduced-matter bispectrum scales as Q3∝1/ℳR(k) for approximately equilateral triangles (and under the assumption that fNLloc is scale independent) [10]. This anomalous scaling considerably raises the ability of the matter bispectrum to constrain primordial NG of the local fNLloc type. Unfortunately, neither the matter bispectrum nor the power spectrum is directly observable with the large-scale structure of the Universe. Temperature anisotropies in the redshifted 21 cm background from the prereionization epoch could in principle furnish a direct measurement of these quantities [78–80], but foreground contamination may severely hamper any detection. Weak lensing is another direct probe of the dark matter, although we can only observe it in projection along the line of sight [81].Figure 3
Reduced matter bispectrumQ3 as a function of the angle θ between k1 and k2 for a fixed k1=0.094hMpc-1 and k2=1.5k1. The panels show ratios between the non-Gaussian and Gaussian Q3 for fNLloc=+100 (top) and -100 (bottom). Dashed lines correspond to tree-level PT while continuous line indicates the 1-loop PT prediction.As we will see shortly, however, this large-scale anomalous scaling is also present in the bispectrum and power spectrum of observable tracers of the large-scale structure such as galaxies. It is this unique signature that will make future all-sky LSS surveys competitive with CMB experiments.


## 3.4. Velocities

Primordial non-Gaussianity also leaves a signature in the large-scale coherent bulk motions which, in the linear regime, are directly related to the linear density field [55]. The various non-Gaussian models considered by [82] tend to have lower-velocity dispersion and bulk flow than fiducial Gaussian model, regardless of the sign of the skewness. However, while the probability distribution of velocity components is sensitive to primordial NG of the local type, in defect models it can be very close to Gaussian, even when the density field is strongly non-Gaussian, as a consequence of the central limit theorem [83, 84]. In this regard, correlation of velocity differences could provide a better test of non-Gaussian initial conditions [85].To measure peculiar velocities, one must subtract the Hubble flow from the observed redshift. This requires an estimate of the distance which is only available for nearby galaxies and clusters (although, e.g., the kinetic Sunyaev-Zel'dovich (kSZ) effect could be used to measure the bulk motions of distant galaxy clusters [86]). So far, measurements of the local galaxy density and velocity fields [87] as well as reconstruction of the initial density PDF from galaxy density and velocity data [88] are consistent with Gaussian initial conditions.


## 4. LSS Probe of Primordial Non-Gaussianity

Discrete and continuous tracers of the large-scale structure, such as galaxies, the Lyα forest, the 21 cm hydrogen line, and so forth, provide a distorted image of the matter density field. In CDM cosmologies, galaxies form inside overdense regions [89] and this introduces a bias between the mass and the galaxy distribution [90]. As a result, distinct samples of galaxies trace the matter distribution differently, the most luminous galaxies preferentially residing in the most massive DM halos. This biasing effect, which concerns most tracers of the large-scale structure, remains to be fully understood. Models of galaxy clustering assume, for instance, that the galaxy biasing relation only depends on the local mass density, but the actual mapping could be more complex [91, 92]. Because of biasing, tracers of the large-scale structure will be affected by primordial non-Gaussianity in a different way than the mass density field. In this section, we describe a number of methods exploiting the abundance and clustering properties of biased tracers to constrain the level of primordial NG. We focus on galaxy clustering as it provides the tightest limits on primordial NG (see Section 5).


### 4.1. Halo Finding Algorithm

Locating groups of bound particles, or DM halos, in simulations is central to the methods described below. In practice, we aim at extracting halo catalogs with statistical properties similar to those of observed galaxies or quasars. This, however, proves to be quite difficult because the relation between observed galaxies and DM halos is somewhat uncertain. Furthermore, there is freedom at defining a halo mass.An important ingredient is the choice of the halo identification algorithm. Halo finders can be broadly divided into two categories: friends-of-friends (FOF) finders [93] and spherical overdensity (SO) finders [94]. While the mass of an SO halo is defined by the radius at which the inner overdensity exceeds Δvir(z) (typically ~ a few hundred times the background density ρ¯(z)), the mass of an FOF halo is given by the number of particles within a linking length b from each other (b~0.15-0.2 in unit of mean interparticle distance). These definitions are somewhat arbitrary and may suit specific purposes only. In what follows, we shall present mainly results for SO halos as their mass estimate is more closely connected to the predictions of the spherical collapse model, on which most of the analytic formulae presented in this section are based. The question of how the spherical overdensity masses can be mapped onto friends-of-friends masses remains a matter of debate (e.g., [95]). Clearly, however, since the peak height ν(M,z) depends on the halo mass M through the variance σM (see below), any systematic difference will be reflected in the value of ν associated to a specific halo sample. As we will see shortly, this affects the size of the fractional deviation from the Gaussian mass function.Catalogs of mock galaxies with luminosities comparable to those of the targeted survey provide an additional layer of complication that can be used, among others, to assess the impact of observational errors on the non-Gaussian signal. However, most numerical studies of cosmic structure formation with primordial NG have not yet addressed this level of sophistication, so we will discuss results based on statistics of dark matter halos only.


### 4.2. Abundances of Voids and Bound Objects

It has long been recognized that departure from Gaussianity can significantly affect the abundance of highly biased tracers of the mass density field, as their frequency sensitively depends upon the tails of the initial density PDF [96–98]. The (extended) Press-Schechter approach has been extensively applied to ascertain the effect of primordial NG on the high mass tail of the mass function.


#### 4.2.1. Press-Schechter Approach

The Press-Schechter theory [99] and its extentions based on excursion sets [100–102] predict that the number density n(M,z) of halos of mass M at redshift z is entirely specified by a multiplicity function f(ν),(20)n(M,z)=ρ¯M2νf(ν)dlnνdlnM,
where the peak height or significance ν(M,z)=δc(z)/σM is the typical amplitude of fluctuations that produce those halos. Here and henceforth, σM denotes the variance of the initial density field δM smoothed on mass scale M∝R3 and linearly extrapolated to present epoch, whereas δc(z)≈1.68D(0)/D(z) is the critical linear overdensity for (spherical) collapse at redshift z. In the standard Press-Schechter approach, n(M,z) is related to the level excursion probability P(>δc,M) that the linear density contrast of a region of mass M exceeds δc(z),(21)νf(ν)=-2ρ¯MdPdM=2πνe-ν2/2,
where the last equality assumes Gaussian initial conditions. The factor of 2 is introduced to account for the contribution of low density regions embedded in overdensities at scale >M. In the extended Press-Schechter theory, δM evolves with M and νf(ν) is the probability that a trajectory is absorbed by the constant barrier δ=δc (as is appropriate in the spherical collapse approximation) on mass scale M. In general, the exact form of f(ν) depends on the barrier shape [103] and the filter shape [104]. Note also that ∫dνf(ν)=1, which ensures that all the mass is contained in halos.Despite the fact that the Press-Schechter mass function overpredicts (underpredicts) the abundance of low (high) mass objects, it can be used to estimate the fractional deviation from Gaussianity. In this formalism, the non-Gaussian fractional correction to the multiplicity function isR(ν,fNLX)≡f(ν,fNLX)/f(ν,0)=(dP/dM)(>δc,M,fNLX)/(dP/dM)(>δc,M,0), which is readily computed once the non-Gaussian density PDF P(δM) is known. In the simple extensions proposed by [105, 106], P(δM) is expressed as the inverse transform of a cumulant-generating function. In [106], the saddle-point technique is applied directly to P(δM). The resulting Edgeworth expansion is then used to obtain P(>δc,M). Neglecting cumulants beyond the skewness, one obtain (we momentarily drop the subscript M for convenience)(22)RLV(ν,fNLX)≈1+16σS3(ν3-3ν)-16d(σS3)dlnν(ν-1ν),
after integration over regions above δc(z). In [105], it is the level excursion probability P(>δc,M) that is calculated within the saddle-point approximation. This approximation better asymptotes to the exact large mass tail, which exponentially deviates from the Gaussian tail. To enforce the normalization of the resulting mass function, one may define ν⋆=δ⋆/σ with δ⋆=δc1-S3δc/3, and use [105, 107](23)ν⋆f(ν⋆)=M2nNG(M,z)ρ¯dlnMdlnν⋆.
The fractional deviation from the Gaussian mass function then becomes(24)RMVJ(ν,fNLX)≈exp(ν36σS3)[-σν26ν⋆dS3dlnν+ν⋆ν].
Both formulae have been shown to give reasonable agreement with numerical simulations of non-Gaussian cosmologies [41, 108, 109] (but note that [11, 110] have reached somewhat different conclusions). Expanding δ⋆ at the first order in fNLX shows that these two theoretical expectations differ in the coefficient of the νσS3 term. Therefore, it is interesting to consider also the approximation(25)R(ν,fNLX)≈exp(ν36σS3)[1-ν2σS3-ν6d(σS3)dlnν],
which is designed to match better the Edgeworth expansion of [106] when the peak height is ν~1 [48]. When the primordial trispectrum is large (i.e., when gNLloc~106), terms involving the kurtosis must be included [48, 105, 106, 111]. In this case, it is also important to take into account a possible renormalization of the fluctuation amplitude, σ8→σ8+δσ8 (12), to which the high mass tail of the mass function is exponentially sensitive [48].Figure4 shows the effect of primordial NG of the local fNLloc type on the mass function of SO halos identified with a redshift-dependent overdensity threshold Δvir(z) (motivated by spherical collapse in a ΛCDM cosmology [112]). Overall, the approximation (25) agrees better with the measurements than (24), which somewhat overestimates the data for fNLloc=100, and than (22), which is not always positive definite for fNLloc=-100. However, as can be seen in Figure 5, while the agreement with the data is reasonable for the SO halos, the theory strongly overestimates the effect in the mass function of FOF halos. Reference [109], who use a FOF algorithm with b=0.2, makes the replacement δc→δcq with q≃0.75 to match the Gaussian and non-Gaussian mass functions. A physical motivation of this replacement is provided by [113, 114], who demonstrate that the diffusive nature of the collapse barrier introduces a similar factor q=(1+DB)-1, regardless of the initial conditions. However, the value of the diffusion constant DB was actually measured from simulations that use a SO finder with constant Δvir=200 [115]. In the excursion set approach of [116], the value of q is obtained by normalizing the Gaussian mass function to simulation (i.e., it has nothing to do with the collapse dynamics) and, therefore, depends on the halo finder. Figure 5 demonstrates that this is also the case for the non-Gaussian correction R(ν,fNLloc): choosing q≃0.75 as advocated in [109] gives good agreement for FOF halos, but strongly underestimates the effect for SO halos, for which the best-fit q is close to unity. As we will see below, the strength of the non-Gaussian bias may also be sensitive to the choice of halo finder.Figure 4
Fractional deviation from the Gaussian mass function as a function of the peak heightν=δc/σ. Different symbols refer to different redshifts as indicated. The various curves are theoretical prediction at z=0 (see text). Halos were identified using a spherical overdensity (SO) finder with a redshift-dependent overdensity threshold Δvir(z) (with Δvir(z) increasing from ~200 at high redshift to attain ~350 at z=0). Error bars denote Poisson errors. For illustration, M=1015M⊙/h corresponds to a significance ν=3.2, 5.2, 7.7 at redshift z=0, 1 and 2, respectively. Similarly, M=1014M⊙/h and 1013M⊙/h correspond to ν=1.9, 3, 4.5 and 1.2, 1.9, 2.9, respectively.Figure 5
Fractional deviation from the Gaussian mass function as a function of the peak heightν=δc/σ. Different symbols refer to different redshifts as in Figure 4. The curves are the theoretical prediction (25) at z=0 with q=1 (solid) and q=0.75 (dotted). In the top panel, halos were identified using a spherical overdensity (SO) finder with a redshift-dependent overdensity threshold Δvir(z) whereas, in the bottom panel, a Friends-of-Friends (FOF) finding algorithm with linking length b=0.2 was used.More sophisticated formulations based on extended Press-Schechter (EPS) theory and/or modifications of the collapse criterion look promising since they can reasonably reproduce both the Gaussian halo counts and the dependence onfNLX [114, 117, 118]. The probability of first upcrossing can, in principle, be derived for any non-Gaussian density field and any choice of smoothing filter [119, 120]. For a general filter, the non-Markovian dynamics generates additional terms in the non-Gaussian correction to the mass function that arise from 3-point correlators of the smoothed density δM at different mass scales [114]. However, large error bars still make difficult to test for the presence of such subleading terms. For generic moving barriers B(σ) such as those appearing in models of triaxial collapse [121, 122], the leading contribution to the non-Gaussian correction approximately is [117](26)R(ν,fNLX)≈1+16σS3H3(B(σ)σ)≈1+16σS3q(qν3-3ν),
where H3(ν)≡ν3-3ν and the last equality assumes ν≫1. For SO halos, (26) with q~0.7 does not fit to the measured correction R(ν,fNLloc) better than (25). However, the ellipsoidal collapse barrier with q~0.7 matches better the Gaussian mass function for moderate peak height ν≲2 [118].Parameterizations of the fractional correction based on N-body simulations have also been considered. While [43] considers a fourth-order polynomial fit to account for values of fNLloc as large as 750, [11] exploits the fact that, for sufficiently small fNLloc, there is a one-to-one mapping between halos in Gaussian and non-Gaussian cosmologies. In both cases, the fitting functions are consistent with the simulations at the few percent level.


#### 4.2.2. Clusters Abundance

Rich clusters of galaxies trace the rare, high-density peaks in the initial conditions and thus offer the best probe of the high-mass tail of the multiplicity function. To infer the cluster mass function, the X-ray and millimeter windows are better suited than the optical-wave range because selection effects can be understood better (see, however, [123]).Following early theoretical [96, 97, 124–126] and numerical [35, 127–129] work on the effect of non-Gaussian initial conditions on the multiplicity function of cosmic structures, the abundance of clusters and X-ray counts in non-Gaussian cosmologies has received much attention in the literature. At fixed normalization of the observed abundance of local clusters, the protoclusters associated with high redshift (2<z<4) Lyα emitters are much more likely to develop in strongly non-Gaussian models than in the Gaussian paradigm [39, 110, 130]. Considering the redshift evolution of cluster abundances thus can break the degeneracy between the initial density PDF and the background cosmology. Simple extensions of the Press-Schechter formalism similar to those considered above have been shown to capture reasonably well the cluster mass function over a wide range of redshift for various non-Gaussian scenarios [131]. Scaling relations between the cluster mass, X-ray temperature and Compton y-parameter calibrated using theory, observations and detailed simulations of cluster formation [132, 133], can be exploited to predict the observed distribution functions of X-ray and SZ signals and assess the capability of cluster surveys to test the nature of the initial conditions [134–140].An important limitation of this method is that the impact of realistic models of primordial non-Gaussianity on cluster abundances is small compared to systematics in current and upcoming surveys [141–143]. Given the current uncertainties in the redshift evolution of clusters (one commonly assumes that clusters are observed at the epoch they collapse [142]), the selection effects in the calibration of X-ray and SZ fluxes with halo mass, the freedom in the definition of the halo mass, the degeneracy with the normalization amplitude σ8 (for positive fNLX, the mass function is more enhanced at the high mass end, and this is similar to an increase in the amplitude of fluctuations σ8 [144]) and the low number statistics, the prospects of using the cluster mass function only to place competitive limits on fNLX with the current data are small. A two-fold improvement in cluster mass calibration is required to provide constraints comparable to CMB measurements [143].


#### 4.2.3. Voids Abundance

The frequency of cosmic voids offers a probe of the low density tail of the initial PDF [145]. The Press-Schechter formalism can also be applied to ascertain the sensitivity of the voids abundance to non-Gaussian initial conditions. Voids are defined as regions of mass M whose density is less than some critical value δv≤0 or, alternatively, whose three eigenvalues of the tidal tensor [146] lie below some critical value λv≤0 [66, 118, 145, 147]. An important aspect in the calculation of the mass function of voids is the overcounting of voids located inside collapsing regions. This voids-in-clouds problem, as identified by [148]), can be solved within the excursion set theory by studying a two barriers problem: δc for halos and δv for voids. Including this effect reduces the frequency of the smallest voids [118]. Neglecting this complication notwithstanding, the differential number density of voids of radius R is [145, 147](27)dndR≈92π2π2|νv|R4e-νv2/2dln|νv|dlnM[1-16σS3H3(|νv|)],
where νv=δv/σM. While a positive fNLX produces more massive halos, it generates fewer large voids [118, 145]. Hence, the effect is qualitatively different from a simple rescaling of the normalization amplitude σ8. A joint analysis of both abundances of clusters and cosmic voids might thus provide interesting constraints on the shape of the primordial 3-point function. There are, however, several caveats to this method, including the fact that there is no unique way to define voids [145]. Clearly, voids identification algorithms will have to be tested on numerical simulations [149] before a robust method can be applied to real data.


### 4.3. Galaxy 2-Point Correlation

Before [90] showed that, in Gaussian cosmologies, correlations of galaxies and clusters can be amplified relative to the mass distribution, it was argued that primeval fluctuations have a non-Gaussian spectrum [150, 151] to explain the observed strong correlation of Abell clusters [152, 153]. Along these lines, [154] pointed out that primordial non-Gaussianity could significantly increase the amplitude of the two-point correlation of galaxies and clusters on large-scales. However, except from [155] who showed that correlations of high density peaks in non-Gaussian models are significantly stronger than in the Gaussian model with identical mass power spectrum, subsequent work focused mostly on abundances (Section 4.2) or higher order statistics such as the bispectrum (Section 4.4). It is only recently that [11] have demonstrated the strong scale-dependent bias arising in non-Gaussian models of the local fNLloc type.


#### 4.3.1. The Non-Gaussian Bias

In the original derivation of [11], the Laplacian is applied to the local mapping Φ=ϕ+fNLlocϕ2 in order to show that, upon substitution of the Poisson equation, the overdensity in the neighborhood of density peaks is spatially modulated by a factor proportional to the local value of ϕ. Taking into account the coherent motions induced by gravitational instabilities, the scale-dependent bias correction reads(28)Δbκ(k,fNLloc)=3fNLloc[b1(M)-1]δc(0)ΩmH02k2T(k)D(z),
where b1(M) is the linear, Eulerian bias of halos of mass M. The original result missed out a multiplicative factor of T(k)-1 which was introduced subsequently by [157] upon a derivation of (28) in the limit of high density peaks. The peak-background split approach [103, 158, 159] promoted by [156] shows that the scale-dependent bias applies to any tracer of the matter density field whose (Gaussian) multiplicity function depends on the local mass density only. In this approach, the Gaussian piece of the potential is decomposed into short- and long-wavelength modes, ϕ=ϕl+ϕs. The short-wavelength piece of the density field is then given by the convolution(29)δs=ℳ⋆Φs=ℳ⋆ϕs(1+2fNLlocϕl)+fNLlocℳ⋆ϕs2,
where ℳ is the transfer function (7). Ignoring the white-noise term, this provides an intuitive explanation of the effect in terms of a local rescaling of the small-scale amplitude of matter fluctuations or, equivalently, a local rescaling of the critical density threshold,(30)σs→σs(1+2fNLlocϕl(x)).
Assuming that the mass function depends only on the peak height ν=δc/σs, the long-wavelength part of the halo overdensity becomes [156] (see also [11, 71, 160])(31)δlh(x)=1n(ν)n(δc-δl(x)σs(1+2fNLlocϕl(x)))-1≈-1σs(δl(x)+2fNLlocϕl(x))dlnndν.
Upon a Fourier transformation and using the fact that, in the Gaussian case, δlh(k)=bLδl(k) with the Lagrangian bias bL=-σs-1dlnn/dν, we recover the non-Gaussian bias correction (28) provided that the tracers move coherently with the dark matter, that is, bL=b1(M)-1 [161]. As emphasized in [11], the scale-dependence arises from the fact that the non-Gaussian curvature perturbation Φ(x) is related to the density through the Poisson equation (6) (so that δl(k)=ℳ(k)ϕl(k)). There is no such effect in the (local) χ2 model, (10), nor in texture-seeded cosmologies [162] for instance.The derivation of [157], based on the clustering of regions of the smoothed density field δM above threshold δc(z), is formally valid for high density peaks only. However, it is general enough to apply to any shape of primordial bispectrum. The 2-point correlation function of that level excursion set, which was first derived by [124], can be expressed in the high threshold limit (ν≫1) as(32)ξ>ν(r)=-1+exp{∑n=2∞∑j=1n-1νnσ-nj!(n-j)!×ξR(n)(x1,…,x1,x2,…,x2jtimes(n-j)times,z=0)∑n=2∞},
where r=x1-x2. For most non-Gaussian models in which the primordial 3-point function is the dominant correction, this expansion can be truncated at the third order and Fourier transformed to yield the non-Gaussian correction ΔP>ν(k) to the power spectrum. Assuming a small level of primordial NG, we can also write ΔP>ν(k)≈2bLΔbκPR(k) where bL≈ν2/δc, and eventually obtain(33)Δbκ(k,fNLX)≡bϕ(k)ℱ(k,fNLX)=(2bLδc(z)ℳR(k,0))ℱ(k,fNLX).
The dependence on the shape of the 3-point function is encoded in the function ℱ(k,fNLX) [157, 163](34)ℱ(k,fNLX)=116π2σ2∫0∞dk1k12ℳR(k1,0)×∫-1+1dμℳR(α,0)BΦ(k1,α,k)PΦ(k),
where α2=k2+k12+2μkk1. Note that, for fNLloc<0, this first order approximation always breaks down at sufficiently small k because ΔP>ν(k)<0.Figure6 shows the non-Gaussian halo bias (33) induced by the local, equilateral and folded bispectrum [163]. In the local and folded non-Gaussianity, the deviation is negligible at k=0.1hMpc-1, but increases rapidly with decreasing wavenumber. Still, the large-scale correction is much smaller for the folded template, and nearly absent for the equilateral type, which make them much more difficult to detect with galaxy surveys [163]. To get insights into the behavior of Δbκ(k,fNLX) at large-scales, let us identify the dominant contribution to ℱ(k,fNLX) in the limit k≪1. Setting ℳR(α,0)≈ℳR(k1,0) and expanding Pϕ(α) at second order in k/k1, we find after some algebra(35)ℱ(k,fNLloc)≈fNLloc,ℱ(k,fNLeq)≈fNLeq[3ΣR(2(ns-4)3)k2(4-ns)/3+12(ns-4)ΣR(-2)k2]σR-2,ℱ(k,fNLfol)≈32fNLfolΣR(ns-43)k(4-ns)/3σR-2,
assuming no running scalar index, that is, dns/dlnk=0. The auxiliary function ΣR(n) is defined as(36)ΣR(n)=12π2∫0∞dkk(2+n)ℳR(k,0)2Pϕ(k).
Hence, we have ΣR(0)≡σR2. As can be seen in Figure 6, these approximations capture relatively well the large scale non-Gaussian bias correction induced by the equilateral and folded type of non-Gaussianity. For a nearly scale-invariant spectrum ns≈1, the effect scales as Δbκ∝k and Δbκ∝ const., respectively. Another important feature of the equilateral and folded non-Gaussian bias is the dependence on the mass scale M through the multiplicative factor σR-2. Indeed, choosing R=1h-1Mpc instead of R=5h-1Mpc as done in Figure 6 would suppress the effect by a factor of ~3. In the high peak limit, σR-2≈bL/δc(z) which cancels out the dependence on redshift but enhances the sensitivity to the halo bias, that is, Δbκ∝bL2 for the equilateral and folded shapes whereas Δbκ∝bL in the local model.Figure 6
Non-Gaussian bias correction (33) for the local, equilateral and folded primordial bispectrum as a function of wavenumber. Results are shown at z=0 for a smoothing radius R=5h-1Mpc and a nonlinear parameter fNLX=1. The dotted line represent an analytic approximation, (35), which is valid at large-scales k≪1. Note that while the magnitude of Δbκ(k,fNLloc) does not change with R, Δbκ(k,fNLeq) and Δbκ(k,fNLfol) strongly depends on the smoothing radius (see text).At this point, it is appropriate to mention a few caveats to these calculations. Firstly, (28) assumes that the tracers form after a spherical collapse, which may be a good approximation for the massive halos only. If one instead considers the ellipsoidal collapse dynamics, in which the evolution of a perturbation depends upon the three eigenvalues of the initial tidal shear, δc(0) should be replaced by its ellipsoidal counterparts δec(0) which is always larger than the spherical value [121]. In this model, the scale-dependent bias Δbκ is thus enhanced by a factor δec(0)/δc(0) [11, 160]. Secondly, (28) assumes that the biasing of the surveyed objects is described by the peak height ν only or, equivalently, the hosting halo mass M. However, this may not be true for quasars whose activity may be triggered by merger of halos [164, 165]. Reference [156] used the EPS formalism to estimate the bias correction Δbmerger induced by mergers(37)Δbmerger=δc-1,
so the factor b1(M)-1 should be replaced by b1(M)-1-δc-1≈b1(M)-1.6. The validity of this result should be evaluated with cosmological simulations of quasars formation. In this respect, semianalytic models of galaxy formation suggest that merger-triggered objects such as quasars do not cluster much differently than other tracers of the same mass [166]. However, this does not mean that the same should hold for the non-Gaussian scale-dependent bias. Still, since the recent merger model is an extreme case it seems likely that the actual bias correction is 0<Δbmerger<δc-1. Thirdly, the scale-dependent bias has been derived using the Newtonian approximation to the Poisson equation, so one may wonder whether general relativistic (GR) corrections to ℳR(k)-1 may suppress the effect on scales comparable to the Hubble radius. Reference [167] showed how large-scale primordial NG induced by GR corrections propagates onto small scales once cosmological perturbations reenter the Hubble radius in the matter dominated era. This effect generates a scale-dependent bias comparable, albeit of opposite sign to that induced by local NG [163]. More recently, [168] argues that there are no GR corrections to the non-Gaussian bias and that the scaling Δbκ∝k-2 applies down to smallish wavenumbers.We can also ask ourselves whether higher-order terms in the series expansion (32) furnish corrections to the non-Gaussian bias similar to (28). The quadratic coupling fNLlocϕ2 induces a second order correction to the halo power spectrum which reads [48](38)ΔPh(k)=43(fNLloc)2[b1(M)-1]2δc2(z)S3(1)(M)×ℳR(k,0)Pϕ(k).
Its magnitude relative to the term linear in fNLloc, (28), is approximately 0.03 at redshift z=1.8 and for a halo mass M=1013M⊙/h. Although its contribution becomes increasingly important at higher redshift, it is fairly small for realistic values of fNLloc. In local NG model, the power spectrum of biased tracers of the density field can also be obtained from a local Taylor series in the evolved (Eulerian) density contrast δ and the Gaussian part ϕ of the initial (Lagrangian) curvature perturbation [47, 71]. Using this approach, it can be shown that the halo power spectrum arising from the first order terms of the local bias expansion can be cast into the form [47](39)Ph(k)=[b1(M)+fNLlocbϕ(k)]2PR(k).
Hence, we also obtain a second order term proportional to (fNLloc)2ℳR-2PR(k)=(fNLloc)2Pϕ(k) which, however, contributes only at very small wavenumber k≲0.001h-1Mpc. All this suggests that (28) is the dominant contribution to the non-Gaussian bias in the wavenumber range 0.001≲k≲0.1hMpc-1.Finally, a non-Gaussian, scale-dependent bias correction can also arise in the local, deterministic bias ansatzδh(x)=b1δ(x)+b2δ(x)2/2+⋯ [169] if the initial density field is non-Gaussian. Here, bN is the Nth-order bias parameters (here again, the first-order bias is b1≡1+bL). In this approach, the correction is induced by the correlation b1b2〈δ(x1)δ2(x2)〉 between the linear and quadratic term in the galaxy biasing relation (which is in fact a collapsed or squeezed 3-point function) and thus reads [67, 68](40)Δbκ(k,fNLloc)=2fNLlocb2σR2ℳR(k,0)-1.
Even though b2σR2≈bLδc in the high-threshold limit ν≫1, b2σR2 behaves very differently than bLδc for moderate peak height because b2 is proportional to the second derivative of the mass function n(ν). So far however, (28) appears to describe reasonably well the numerical results for a wide range of halo bias.


#### 4.3.2. Comparison with Simulations

In order to fully exploit the potential of forthcoming large-scale surveys, a number of studies have tested the theoretical prediction against the outcome of large numerical simulations [11, 41–43, 71, 109].At the lowest order, there are two additional albeit relatively smaller corrections to the Gaussian bias which arise from the dependence of both the halo number densityn(M,z) and the matter power spectrum Pδ(k,z) on primordial NG [41]. Firstly, assuming the peak-background split holds, the change in the mean number density of halos induces a scale-independent offset which we denote ΔbI(fNLloc). In terms of the non-Gaussian fractional correction R(ν,fNLloc) to the mass function, this contribution is(41)ΔbI(fNLloc)=-1σ∂∂νln[R(ν,fNLloc)].
It is worth noticing that ΔbI(fNLloc) has a sign opposite to that of fNLloc, because the bias decreases when the mass function goes up. In practice, the approximation (25), which matches well the SO data for ν≲4, can be used for moderate values of the peak height. For FOF halos with linking length b=0.2, one should make the replacement δc→δcq with q≈0.75 in the calculation of the scale-independent offset. It is sensible to evaluate ΔbI(fNLloc) at a mass scale 〈M〉 equal to the average halo mass of the sample. Secondly, we also need to account for the change in the matter power spectrum (see Figure 2 in Section 3). Summarizing, local non-Gaussianity adds a correction Δb(k,fNLloc) to the bias b(k) of dark matter halos that reads [41](42)Δb(k,fNLloc)=Δbκ(k,fNLloc)+ΔbI(fNLloc)+b1(M)βδ(k,fNLloc)
at first order in fNLloc. As can be seen in Figure 7, the inclusion of these extra terms substantially improves the comparison between the theory and the simulations. Considering only the scale-dependent shift Δbκ leads to an apparent suppression of the effect in simulations relative to the theory. Including the scale-independent offset ΔbI considerably improves the agreement at wavenumbers k≲0.05hMpc-1. Finally, adding the scale-dependent term b1(M)βδ further adjusts the match at small scale k≳0.05hMpc-1 by making the non-Gaussian bias shift less negative. Along these lines, [71] find that the inclusion of ΔbI to the bias also improves the agreement with measurements of Δb(k,fNLloc) obtained for FOF halos, and show that taking into account second- and higher-order corrections could extend the validity of the theory up to scales k~0.1-0.3hMpc-1.Figure 7
Non-Gaussian bias correction (filled symbols) for halos of massM>2×1013M⊙/h extracted at z=0.5 from simulations of the local fNLloc model. The solid curve represents the theoretical model (42). The dashed, dotted-dashed-dotted, and dotted curves show the scale-dependent bias Δbκ, the scale-independent offset ΔbI and the contribution from the matter power spectrum b1(M)βδ that arise at first order in fNLloc (see Figure 2). The error bars indicate the scatter among 8 realizations of 10243 simulations with box size L=1600h-1Mpc.The non-Gaussian bias correction can be measured in the cross- and autopower spectrum of dark matter halos,Phδ(k) and Ph(k). To compute these quantities, dark matter particles and halo centers are interpolated onto a regular cubical mesh. The resulting dark matter and halo fluctuation fields, δm(k) and δh(k), are then Fourier transformed to yield the matter-matter, halo-matter and halo-halo power spectra Pδ(k), Phδ(k) and Ph(k), respectively. Ph(k) is then corrected for the shot noise, which is assumed to be 1/n¯h if dark matter halos are a Poisson sampling of some continuous field. This discreteness correction is negligible for Pδ(k) due to the large number of dark matter particles. On linear scales (k≲0.01hMpc-1), the halo bias b(k)=δh(k)/δm(k) approaches the constant value b1(M) which needs to be measured accurately as it controls the strength of the scale-dependent bias correction Δbκ. In this respect, the ratio Phδ(k)/Pδ(k) is a better proxy for the halo bias since it is less sensitive to shot-noise.Auto- and crosspower analyses may not agree with each other if the halos and dark matter do not trace each other on scalek≲0.01hMpc-1 where the non-Gaussian bias is large, that is, if there is stochasticity. Figure 8 shows Phδ(k) and Ph(k) averaged over 8 realizations of the models with fNLloc=0,±100 [41]. The same Gaussian random seed field ϕ was used in each set of runs so as to minimize the sampling variance. Measurements of the non-Gaussian bias correction obtained with the halo-halo or the halo-matter power spectrum are in a good agreement with each other, indicating that non-Gaussianity does not induce stochasticity and the predicted scaling (28) applies equally well for the auto- and cross-power spectrum. However, while a number of numerical studies of the fNLloc model have confirmed the scaling Δbκ(k,fNLloc)∝ℳR(k)-1 and the redshift dependence ∝D(z)-1 [11, 41, 43, 109], the exact amplitude of the non-Gaussian bias correction remains somewhat debatable. Reference [41] who use SO halos and [71] who use FOF halos find satisfactory agreement with the theory once the scale-independent offset ΔbI is included. By contrast, see [43], who use the same FOF halos as [71], argue that the scale-dependent piece Δbκ requires, among others, a multiplicative correction of the form (1-β1fNLloc), with β1~4×10-4>0. Similarly, [109] who also use FOF halos find that the theory is a good fit to the simulations only upon replacing bL by qbL in (33), with q≃0.75. Part of the discrepancy may be probably due to the fact that the last two references do not include ΔbI, which leads to an apparent suppression of the effect (see Figure 7). Another possible source of discrepancy may be choice of the halo finder which, as seen in Figure 5, has an impact on the strength of the non-Gaussian correction to the mass function. This possibility is investigated in Figure 9, which shows the non-Gaussian bias correction obtained with FOF halos. For this low biased sample, the scale-independent correction is |ΔbI|≲0.003 and can thus be neglected. The best-fit values of fNLloc are significantly below the input values of ±100, in agreement with the findings of [43, 109] (note, however, that this suppression is more consistent with δc being rescaled by qδc≈0.86δc and bL being unchanged). This indicates that the choice of halo finder may also affect the magnitude of the scale-dependent non-Gaussian bias. Discrepancies have also been observed between the theoretical and measured non-Gaussian bias corrections in non-Gaussian models of the local cubic-order coupling gNLlocϕ3 [48]. Understanding all these results clearly requires a better modeling of halo clustering.Figure 8
Halo-halo (solid curve) and halo-matter (dashed curve) power spectraPh(k) and Phδ(k) measured in simulations of the Gaussian and fNLloc=±100 models for halos of mass M>2×1013M⊙/h at redshift z=1. The error bars represent the scatter among 8 realizations. For fNLloc=-100, the crosspower spectrum is negative on scales k≲0.005hMpc-1, in good agreement with the theoretical prediction.Figure 9
Fractional correction to the Gaussian halo bias in thefNLloc=±100 and Gaussian models. In contrast to Figure 8, halos were identified with a FOF finder of linking length b=0.2. Only the wavemodes to the left of the vertical line were used to fit Δbκ(k,fNLloc). For this low biased sample, the scale-independent correction is |ΔbI|≲0.003 and can thus be ignored. The best-fit value of fNLloc and the corresponding 1σ error is quoted for each model (Figure taken from [170]).


#### 4.3.3. Redshift-Space Distortions

Peculiar velocities generate systematic differences between the spatial distribution of data in real and redshift space. These redshift-space distortions must be properly taken into account in order to extractfNLX from redshift surveys. On the linear scales of interest, the redshift-space power spectrum of biased tracers reads as [171, 172](43)Phs(k,μ)=[b12Pδ(k)+2b1fμ2Pδθ(k)+f2μ4Pθ(k)],
where Pδθ and Pθ are the density-velocity and velocity-divergence power spectra, μ is the cosine of the angle between the wavemode k and the line of sight and f is the logarithmic derivative of the growth factor. For Pθ, the 1-loop correction due to primordial NG is identical to (17) provided F2(k1,k2) is replaced by the kernel G2(k1,k2)=3/7+μ(k1/k2+k2/k1)/2+4μ2/7 describing the 2nd order evolution of the velocity divergence [62]. For Pδθ, this correction is(44)ΔPδθNG(k)=∫d3q(2π)3[F2(q,k-q)+G2(q,k-q)]×B0(-k,q,k-q).
Again, causality implies that G2(k1,k2) vanishes in the limit k1=-k2. For unbiased tracers with b1=1, the linear Kaiser relation is thus recovered at large scales k≲0.01hMpc-1 (this is consistent with the analysis of [173]). For biased tracers, we still expect the Kaiser formula to be valid, but the distortion parameter β should now be equal to β=f/(b1+Δbκ), where Δbκ(k,fNLX) is the scale-dependent bias induced by the primordial non-Gaussianity.


#### 4.3.4. Mitigating Cosmic Variance and Shot-Noise

Because of the finite number of large scale wavemodes accessible to a survey, any large-scale measurement of the power spectrum is limited by the cosmic (or sampling) variance caused by the random nature of the wavemodes. For discrete tracers such as galaxies, the shot noise is another source of error. Restricting ourselves to weak primordial NG, the relative error on the power spectrumP is σP/P≈1/N(1+σn2/P), where N is the number of independent modes measured and σn2 is the shot-noise [174]. Under the standard assumption of Poisson sampling, σn2 equals the inverse of the number density 1/n¯ and causes a scale-independent enhancement of the power spectrum. The extent to which one can improve the observational limits on the nonlinear parameters fNLX will strongly depend on our ability to minimize the impact of these two sources of errors. By comparing differently biased tracers of the same surveyed volume [175, 176] and suitably weighting galaxies (by the mass of their host halo for instance) [177, 178], it should be possible to circumvent these problems and considerably improve the detection level.Figure9 illustrates how the impact of sampling variance on the measurement of fNLloc can be mitigated. Namely, the data points show the result of taking the ratio Ph(k,fNLloc)/Pδ(k,fNLloc) for each set of runs with same Gaussian random seed field ϕ before averaging over the realizations. This procedure is equivalent to the multitracers method advocated by [175]. Here, Pδ can be thought as mimicking the power spectrum of a nearly unbiased tracer of the mass density field with high number density. Although, in practical applications, using the dark matter field works better [170], in real data Pδ should be replaced by a tracer of the same surveyed volume different than the one used to compute Ph. Figure 9 also shows that upon taking out most of the cosmic variance, there is some residual noise caused by the discrete nature of the dark matter halos. As shown recently [178] however, weighting the halos according to their mass can dramatically reduce the shot noise relative to the Poisson expectation, at least when compared against the dark matter. Applying such a weighting may thus significantly improve the error on the nonlinear parameter fNLloc, but this should be explored in realistic simulations of galaxies, especially because the halo mass M may not be easily measurable from real data [170]. This approach undoubtedly deserves further attention as it has the potential to substantially improve the extraction of the primordial non-Gaussian signal from galaxy surveys.To conclude this section, it is worth noting that, while the PDF of power valuesP(k) has little discriminatory power (for large surveyed volume, it converges towards the Rayleigh distribution as a consequence of the central limit theorem) [179], the covariance of power spectrum measurements (which is sensitive to the selection function, but also to correlations among the phase of the Fourier modes) may provide quantitative limits on certain type of non-Gaussian models [174, 180].


### 4.4. Galaxy Bispectrum and Higher-Order Statistics

Higher statistics of biased tracers, such as the galaxy bispectrum, are of great interest as they are much more sensitive to the shape of the primordial 3-point function than the power spectrum [10, 42, 68, 181, 182]. Therefore, they could break some of the degeneracies affecting the non-Gaussian halo bias (e. g., the leading order scale-dependent correction to the Gaussian bias induced by the local quadratic and cubic coupling are fully degenerated [48]).


#### 4.4.1. Normalized Cumulants of the Galaxy Distribution

The skewness of the galaxy count probability distribution function could provide constraints on the amount of non-Gaussianity in the initial conditions. As discussed in Section3, however, it is difficult to disentangle the primordial and gravitational causes of skewness in low redshift data unless the initial density field is strongly non-Gaussian. The first analyzes of galaxy catalogs in terms of count-in-cells densities all reached the conclusion that the skewness (and higher-order moments) of the observed galaxy count PDF is consistent with the value predicted by gravitational instability of initially Gaussian fluctuations [50, 57, 60, 183–185]. Back then however, most of the galaxy samples available were not large enough to accurately determine the SJ at large-scales [186]. Despite the two orders of magnitude increase in surveyed volume, these measurements are still sensitive to cosmic variance, that is, to the presence of massive super-clusters or large voids. Nevertheless, the best estimates of the first normalized cumulants SJ of the galaxy PDF strongly suggest that high order galaxy correlation functions indeed follow the hierarchical scaling predicted by the gravitational clustering of Gaussian ICs [187]. There is no evidence for strong non-Gaussianity in the initial density field as might by seeded by cosmic strings or textures [188].The genus statistics of constant density surfaces through the galaxy distribution measures the relative abundance of low and high density regions as a function of the smoothing scaleR and, therefore, could also be used as a diagnostic tool for primordial non-Gaussianity. For a Gaussian random field, the genus curve (i.e., the genus number as a function of the density contrast) is symmetric about δR=0 regardless the value of R. Primordial NG and nonlinear gravitational evolution can disrupt this symmetry [189]. The effect of non-Gaussian ICs on the topology of the galaxy distribution has been explored in a number of papers [35, 190–193]. For large values of R and realistic amount of primordial NG, the genus statistics can also be expanded in a series whose coefficients are the normalized cumulants SJ of the smoothed galaxy density field. In other words, the genus statistics essentially provides another measure of the (large-scale) cumulants. So far, measurements from galaxy data are broadly consistent with Gaussian initial conditions [194, 195].


#### 4.4.2. Galaxy Bispectrum

Most of the scale-dependence of the primordialn-point functions is integrated out in the normalized cumulants, which makes them weakly sensitive to primordial NG. However, while the effect of non-Gaussian initial conditions, galaxy bias, gravitational instabilities and so forth, are strongly degenerated in the SJ, they imprint distinct signatures in the galaxy bispectrum Bg(k1,k2,k3), an accurate measurement of which could thus constrain the shape of the primordial 3-point function.In the original derivation of [181], the large-scale (unfiltered) galaxy bispectrum in the fNLloc model is given by(45)Bg(k1,k2,k3)=b13B0(k1,k2,k3)+b12b2[P0(k1)P0(k2)+(cyc.)]+2b13[F2(k1,k2)P0(k1)P0(k2)+(cyc.)].
Again, b1 and b2 are the first- and second-order bias parameters that describe the galaxy biasing relation assumed local and deterministic [169]. The first term in the right-hand side is the primordial contribution which, for equilateral configurations and in the fNLloc model, scales as ℳR(k,z)-1 like in the matter bispectrum, (18). The two last terms are the contribution from nonlinear bias and the tree-level correction from gravitational instabilities, respectively. They have the smallest signal in squeezed configurations.As recognized by [68, 182], (45) misses an important term that may significantly enhance the sensitivity of the galaxy bispectrum to non-Gaussian initial conditions. This contribution is sourced by the trispectrum TR(k1,k2,k3,k4) of the smoothed mass density field(46)12b12b2∫d3q(2π)3TR(k1,k2,q,k3-q)+(2perms.).
At large-scale, this simplifies to the sum of the linearly evolved primordial trispectrum T0(k1,k2,k3,k4) and a coupling between the primordial bispectrum B0(k1,k2,k3) (linear in fNLX) and the second order PT corrections (through the kernel F2(k1,k2)). In the case of local non-Gaussianity and for equilateral configurations, the first piece proportional to T0 scales as (fNLloc)2k-4 times the Gaussian tree-level prediction, with the same redshift dependence. Hence, it is similar to the second order correction (fNLloc)2ℳR-2PR(k) that appears in the halo power spectrum (see (39)). The second piece linear in fNLX generates a signal at large-scales for essentially all triangle shapes in the local model as well as in the case of equilateral NG. This second contribution is maximized in the squeezed limit (where it is one order of magnitude larger than the result obtained by [181]) which helps disentangling it from the Gaussian terms. Note that a strong dependence on triangle shape is also present in other NG scenarios such as the χ2 model [62].This newly derived contributions are claimed to lead to more than one order of magnitude improvement in certain limits [182], but it is not yet clear whether these gains can be fully realized with upcoming galaxy surveys. To accurately predict the constraints that could be achieved with future measurements of the galaxy bispectrum, a comparison of these predictions with the halo bispectrum extracted from numerical simulations is highly desirable. To date, the only numerical study [42] has measured the halo bispectrum for some isosceles triangles (k1=k2). While the shape dependence is in reasonable agreement with the theory, the observed k-dependence appears to depart from the predicted scaling.


### 4.5. Intergalactic Medium and the Lyα Forest

Primordial non-Gaussianity also affects the intergalactic medium (IGM) as a positivefNLX enhances the formation of high-mass halos at early times and, therefore, accelerates reionization [196–198]. At lower redshift, small box hydrodynamical simulations of the Lyα forest indicate that non-Gaussian initial conditions could leave a detectable signature in the Lyα flux PDF, power spectrum and bispectrum [199]. However, while differences appear quite pronounced in the high transmissivity tail of the flux PDF (i.e., in underdense regions), the Lyα 1D flux power spectrum seems little affected. Given the small box size of these hydrodynamical simulations, it is worth exploring the effect in large N-body cosmological simulations using a semianalytic modeling of the Lyα forest [200], even though such an approach only provides a very crude approximation to the temperature-density diagram of the IGM in hydrodynamical simulations. Figure 10 shows the imprint of local type NG on the Lyα 3D flux power spectrum (which is not affected by projection effects) extracted at z=2 from a series of large simulations. The Lyα transmitted flux is calculated in the Gunn-Peterson approximation [201]. A clear signature similar to the non-Gaussian halo bias can be seen. As expected, it is of opposite sign since the Lyα forest is anti-biased relative to the mass density field (overdensities are mapped onto relatively low flux transmission).Figure 10
Ratio between thez=2 Lyα flux power spectrum extracted from simulations of Gaussian and non-Gaussian initial conditions. The mean transmission is set to F¯=0.8 and the power-law exponent α=1.65 (see text).To estimate the strength of the signal (see [200] for the details), one can assume that the (real space) optical depth τ(x) to Lyα absorption at comoving position x is approximately [202](47)τ(x)=τ¯[1+δg(x)]α,
where δg is the gas density, τ¯~1 is the optical depth at mean gas density and α~1-2 is some parameter that depends on the exact thermal history of the low density IGM. The above relation holds for the moderate overdensities δg≲10 that are responsible for most of the Lyα absorption features. To relate the gas density to the smoothed linear density field, we could make the simple ansatz δg≡δR [203]. In this linear approximation, however, the large-scale bias bF of the Lyα flux density field is much larger than that measured in detailed numerical simulations (e.g., bF2≃0.017 at z=3 [204]). Therefore, one may want to consider the lognormal mapping [205, 206](48)1+δg=exp(δR-σR22),
to better capture nonlinearities in the gas density field. Expanding exp(δR) at second order and noticing that in the presence of weak non-Gaussianity, the joint PDF P(δR(x1),δR(x2)) can generically be expanded into an Edgeworth series where the primordial 3-point function is the dominant correction, it is straightforward to compute the Lyα 3D flux power spectrum for nonzero fNLX. Upon a Fourier transformation, we arrive at (49)PF(k,fNLX)PF(k,0)=1-4gFσR2ℳR(k,z)-1ℱ(k,fNLX).
where gF is some auxiliary function of (τ¯,α,σR). This result is valid for any model of primordial NG characterized by an initial bispectrum. In the fNLloc model, the large-scale non-Gaussian Lyα bias scales as ΔbF(k,fNLX)≈-2gFσR2ℳR(k,z)-1∝k-2T(k)-1 like the non-Gaussian halo bias. Assuming that τ¯=0.7, σR=1.8 and α=1.65 yields a mean flux F¯≈0.8 and a ratio PF(k=0.01,fNLX)/PF(k=0.01,0)≈1∓0.13 for fNLloc=±100 comparable in magnitude to that seen in Figure 10. A detection of this effect, although challenging in particular because of continuum uncertainties, could be feasible with future data sets. Summarizing, the Lyα should provide interesting information on the non-Gaussian signal over a range of scale and redshift not easily accessible to galaxy and CMB observations [199, 200].


## 4.1. Halo Finding Algorithm

Locating groups of bound particles, or DM halos, in simulations is central to the methods described below. In practice, we aim at extracting halo catalogs with statistical properties similar to those of observed galaxies or quasars. This, however, proves to be quite difficult because the relation between observed galaxies and DM halos is somewhat uncertain. Furthermore, there is freedom at defining a halo mass.An important ingredient is the choice of the halo identification algorithm. Halo finders can be broadly divided into two categories: friends-of-friends (FOF) finders [93] and spherical overdensity (SO) finders [94]. While the mass of an SO halo is defined by the radius at which the inner overdensity exceeds Δvir(z) (typically ~ a few hundred times the background density ρ¯(z)), the mass of an FOF halo is given by the number of particles within a linking length b from each other (b~0.15-0.2 in unit of mean interparticle distance). These definitions are somewhat arbitrary and may suit specific purposes only. In what follows, we shall present mainly results for SO halos as their mass estimate is more closely connected to the predictions of the spherical collapse model, on which most of the analytic formulae presented in this section are based. The question of how the spherical overdensity masses can be mapped onto friends-of-friends masses remains a matter of debate (e.g., [95]). Clearly, however, since the peak height ν(M,z) depends on the halo mass M through the variance σM (see below), any systematic difference will be reflected in the value of ν associated to a specific halo sample. As we will see shortly, this affects the size of the fractional deviation from the Gaussian mass function.Catalogs of mock galaxies with luminosities comparable to those of the targeted survey provide an additional layer of complication that can be used, among others, to assess the impact of observational errors on the non-Gaussian signal. However, most numerical studies of cosmic structure formation with primordial NG have not yet addressed this level of sophistication, so we will discuss results based on statistics of dark matter halos only.


## 4.2. Abundances of Voids and Bound Objects

It has long been recognized that departure from Gaussianity can significantly affect the abundance of highly biased tracers of the mass density field, as their frequency sensitively depends upon the tails of the initial density PDF [96–98]. The (extended) Press-Schechter approach has been extensively applied to ascertain the effect of primordial NG on the high mass tail of the mass function.


### 4.2.1. Press-Schechter Approach

The Press-Schechter theory [99] and its extentions based on excursion sets [100–102] predict that the number density n(M,z) of halos of mass M at redshift z is entirely specified by a multiplicity function f(ν),(20)n(M,z)=ρ¯M2νf(ν)dlnνdlnM,
where the peak height or significance ν(M,z)=δc(z)/σM is the typical amplitude of fluctuations that produce those halos. Here and henceforth, σM denotes the variance of the initial density field δM smoothed on mass scale M∝R3 and linearly extrapolated to present epoch, whereas δc(z)≈1.68D(0)/D(z) is the critical linear overdensity for (spherical) collapse at redshift z. In the standard Press-Schechter approach, n(M,z) is related to the level excursion probability P(>δc,M) that the linear density contrast of a region of mass M exceeds δc(z),(21)νf(ν)=-2ρ¯MdPdM=2πνe-ν2/2,
where the last equality assumes Gaussian initial conditions. The factor of 2 is introduced to account for the contribution of low density regions embedded in overdensities at scale >M. In the extended Press-Schechter theory, δM evolves with M and νf(ν) is the probability that a trajectory is absorbed by the constant barrier δ=δc (as is appropriate in the spherical collapse approximation) on mass scale M. In general, the exact form of f(ν) depends on the barrier shape [103] and the filter shape [104]. Note also that ∫dνf(ν)=1, which ensures that all the mass is contained in halos.Despite the fact that the Press-Schechter mass function overpredicts (underpredicts) the abundance of low (high) mass objects, it can be used to estimate the fractional deviation from Gaussianity. In this formalism, the non-Gaussian fractional correction to the multiplicity function isR(ν,fNLX)≡f(ν,fNLX)/f(ν,0)=(dP/dM)(>δc,M,fNLX)/(dP/dM)(>δc,M,0), which is readily computed once the non-Gaussian density PDF P(δM) is known. In the simple extensions proposed by [105, 106], P(δM) is expressed as the inverse transform of a cumulant-generating function. In [106], the saddle-point technique is applied directly to P(δM). The resulting Edgeworth expansion is then used to obtain P(>δc,M). Neglecting cumulants beyond the skewness, one obtain (we momentarily drop the subscript M for convenience)(22)RLV(ν,fNLX)≈1+16σS3(ν3-3ν)-16d(σS3)dlnν(ν-1ν),
after integration over regions above δc(z). In [105], it is the level excursion probability P(>δc,M) that is calculated within the saddle-point approximation. This approximation better asymptotes to the exact large mass tail, which exponentially deviates from the Gaussian tail. To enforce the normalization of the resulting mass function, one may define ν⋆=δ⋆/σ with δ⋆=δc1-S3δc/3, and use [105, 107](23)ν⋆f(ν⋆)=M2nNG(M,z)ρ¯dlnMdlnν⋆.
The fractional deviation from the Gaussian mass function then becomes(24)RMVJ(ν,fNLX)≈exp(ν36σS3)[-σν26ν⋆dS3dlnν+ν⋆ν].
Both formulae have been shown to give reasonable agreement with numerical simulations of non-Gaussian cosmologies [41, 108, 109] (but note that [11, 110] have reached somewhat different conclusions). Expanding δ⋆ at the first order in fNLX shows that these two theoretical expectations differ in the coefficient of the νσS3 term. Therefore, it is interesting to consider also the approximation(25)R(ν,fNLX)≈exp(ν36σS3)[1-ν2σS3-ν6d(σS3)dlnν],
which is designed to match better the Edgeworth expansion of [106] when the peak height is ν~1 [48]. When the primordial trispectrum is large (i.e., when gNLloc~106), terms involving the kurtosis must be included [48, 105, 106, 111]. In this case, it is also important to take into account a possible renormalization of the fluctuation amplitude, σ8→σ8+δσ8 (12), to which the high mass tail of the mass function is exponentially sensitive [48].Figure4 shows the effect of primordial NG of the local fNLloc type on the mass function of SO halos identified with a redshift-dependent overdensity threshold Δvir(z) (motivated by spherical collapse in a ΛCDM cosmology [112]). Overall, the approximation (25) agrees better with the measurements than (24), which somewhat overestimates the data for fNLloc=100, and than (22), which is not always positive definite for fNLloc=-100. However, as can be seen in Figure 5, while the agreement with the data is reasonable for the SO halos, the theory strongly overestimates the effect in the mass function of FOF halos. Reference [109], who use a FOF algorithm with b=0.2, makes the replacement δc→δcq with q≃0.75 to match the Gaussian and non-Gaussian mass functions. A physical motivation of this replacement is provided by [113, 114], who demonstrate that the diffusive nature of the collapse barrier introduces a similar factor q=(1+DB)-1, regardless of the initial conditions. However, the value of the diffusion constant DB was actually measured from simulations that use a SO finder with constant Δvir=200 [115]. In the excursion set approach of [116], the value of q is obtained by normalizing the Gaussian mass function to simulation (i.e., it has nothing to do with the collapse dynamics) and, therefore, depends on the halo finder. Figure 5 demonstrates that this is also the case for the non-Gaussian correction R(ν,fNLloc): choosing q≃0.75 as advocated in [109] gives good agreement for FOF halos, but strongly underestimates the effect for SO halos, for which the best-fit q is close to unity. As we will see below, the strength of the non-Gaussian bias may also be sensitive to the choice of halo finder.Figure 4
Fractional deviation from the Gaussian mass function as a function of the peak heightν=δc/σ. Different symbols refer to different redshifts as indicated. The various curves are theoretical prediction at z=0 (see text). Halos were identified using a spherical overdensity (SO) finder with a redshift-dependent overdensity threshold Δvir(z) (with Δvir(z) increasing from ~200 at high redshift to attain ~350 at z=0). Error bars denote Poisson errors. For illustration, M=1015M⊙/h corresponds to a significance ν=3.2, 5.2, 7.7 at redshift z=0, 1 and 2, respectively. Similarly, M=1014M⊙/h and 1013M⊙/h correspond to ν=1.9, 3, 4.5 and 1.2, 1.9, 2.9, respectively.Figure 5
Fractional deviation from the Gaussian mass function as a function of the peak heightν=δc/σ. Different symbols refer to different redshifts as in Figure 4. The curves are the theoretical prediction (25) at z=0 with q=1 (solid) and q=0.75 (dotted). In the top panel, halos were identified using a spherical overdensity (SO) finder with a redshift-dependent overdensity threshold Δvir(z) whereas, in the bottom panel, a Friends-of-Friends (FOF) finding algorithm with linking length b=0.2 was used.More sophisticated formulations based on extended Press-Schechter (EPS) theory and/or modifications of the collapse criterion look promising since they can reasonably reproduce both the Gaussian halo counts and the dependence onfNLX [114, 117, 118]. The probability of first upcrossing can, in principle, be derived for any non-Gaussian density field and any choice of smoothing filter [119, 120]. For a general filter, the non-Markovian dynamics generates additional terms in the non-Gaussian correction to the mass function that arise from 3-point correlators of the smoothed density δM at different mass scales [114]. However, large error bars still make difficult to test for the presence of such subleading terms. For generic moving barriers B(σ) such as those appearing in models of triaxial collapse [121, 122], the leading contribution to the non-Gaussian correction approximately is [117](26)R(ν,fNLX)≈1+16σS3H3(B(σ)σ)≈1+16σS3q(qν3-3ν),
where H3(ν)≡ν3-3ν and the last equality assumes ν≫1. For SO halos, (26) with q~0.7 does not fit to the measured correction R(ν,fNLloc) better than (25). However, the ellipsoidal collapse barrier with q~0.7 matches better the Gaussian mass function for moderate peak height ν≲2 [118].Parameterizations of the fractional correction based on N-body simulations have also been considered. While [43] considers a fourth-order polynomial fit to account for values of fNLloc as large as 750, [11] exploits the fact that, for sufficiently small fNLloc, there is a one-to-one mapping between halos in Gaussian and non-Gaussian cosmologies. In both cases, the fitting functions are consistent with the simulations at the few percent level.


### 4.2.2. Clusters Abundance

Rich clusters of galaxies trace the rare, high-density peaks in the initial conditions and thus offer the best probe of the high-mass tail of the multiplicity function. To infer the cluster mass function, the X-ray and millimeter windows are better suited than the optical-wave range because selection effects can be understood better (see, however, [123]).Following early theoretical [96, 97, 124–126] and numerical [35, 127–129] work on the effect of non-Gaussian initial conditions on the multiplicity function of cosmic structures, the abundance of clusters and X-ray counts in non-Gaussian cosmologies has received much attention in the literature. At fixed normalization of the observed abundance of local clusters, the protoclusters associated with high redshift (2<z<4) Lyα emitters are much more likely to develop in strongly non-Gaussian models than in the Gaussian paradigm [39, 110, 130]. Considering the redshift evolution of cluster abundances thus can break the degeneracy between the initial density PDF and the background cosmology. Simple extensions of the Press-Schechter formalism similar to those considered above have been shown to capture reasonably well the cluster mass function over a wide range of redshift for various non-Gaussian scenarios [131]. Scaling relations between the cluster mass, X-ray temperature and Compton y-parameter calibrated using theory, observations and detailed simulations of cluster formation [132, 133], can be exploited to predict the observed distribution functions of X-ray and SZ signals and assess the capability of cluster surveys to test the nature of the initial conditions [134–140].An important limitation of this method is that the impact of realistic models of primordial non-Gaussianity on cluster abundances is small compared to systematics in current and upcoming surveys [141–143]. Given the current uncertainties in the redshift evolution of clusters (one commonly assumes that clusters are observed at the epoch they collapse [142]), the selection effects in the calibration of X-ray and SZ fluxes with halo mass, the freedom in the definition of the halo mass, the degeneracy with the normalization amplitude σ8 (for positive fNLX, the mass function is more enhanced at the high mass end, and this is similar to an increase in the amplitude of fluctuations σ8 [144]) and the low number statistics, the prospects of using the cluster mass function only to place competitive limits on fNLX with the current data are small. A two-fold improvement in cluster mass calibration is required to provide constraints comparable to CMB measurements [143].


### 4.2.3. Voids Abundance

The frequency of cosmic voids offers a probe of the low density tail of the initial PDF [145]. The Press-Schechter formalism can also be applied to ascertain the sensitivity of the voids abundance to non-Gaussian initial conditions. Voids are defined as regions of mass M whose density is less than some critical value δv≤0 or, alternatively, whose three eigenvalues of the tidal tensor [146] lie below some critical value λv≤0 [66, 118, 145, 147]. An important aspect in the calculation of the mass function of voids is the overcounting of voids located inside collapsing regions. This voids-in-clouds problem, as identified by [148]), can be solved within the excursion set theory by studying a two barriers problem: δc for halos and δv for voids. Including this effect reduces the frequency of the smallest voids [118]. Neglecting this complication notwithstanding, the differential number density of voids of radius R is [145, 147](27)dndR≈92π2π2|νv|R4e-νv2/2dln|νv|dlnM[1-16σS3H3(|νv|)],
where νv=δv/σM. While a positive fNLX produces more massive halos, it generates fewer large voids [118, 145]. Hence, the effect is qualitatively different from a simple rescaling of the normalization amplitude σ8. A joint analysis of both abundances of clusters and cosmic voids might thus provide interesting constraints on the shape of the primordial 3-point function. There are, however, several caveats to this method, including the fact that there is no unique way to define voids [145]. Clearly, voids identification algorithms will have to be tested on numerical simulations [149] before a robust method can be applied to real data.


## 4.2.1. Press-Schechter Approach

The Press-Schechter theory [99] and its extentions based on excursion sets [100–102] predict that the number density n(M,z) of halos of mass M at redshift z is entirely specified by a multiplicity function f(ν),(20)n(M,z)=ρ¯M2νf(ν)dlnνdlnM,
where the peak height or significance ν(M,z)=δc(z)/σM is the typical amplitude of fluctuations that produce those halos. Here and henceforth, σM denotes the variance of the initial density field δM smoothed on mass scale M∝R3 and linearly extrapolated to present epoch, whereas δc(z)≈1.68D(0)/D(z) is the critical linear overdensity for (spherical) collapse at redshift z. In the standard Press-Schechter approach, n(M,z) is related to the level excursion probability P(>δc,M) that the linear density contrast of a region of mass M exceeds δc(z),(21)νf(ν)=-2ρ¯MdPdM=2πνe-ν2/2,
where the last equality assumes Gaussian initial conditions. The factor of 2 is introduced to account for the contribution of low density regions embedded in overdensities at scale >M. In the extended Press-Schechter theory, δM evolves with M and νf(ν) is the probability that a trajectory is absorbed by the constant barrier δ=δc (as is appropriate in the spherical collapse approximation) on mass scale M. In general, the exact form of f(ν) depends on the barrier shape [103] and the filter shape [104]. Note also that ∫dνf(ν)=1, which ensures that all the mass is contained in halos.Despite the fact that the Press-Schechter mass function overpredicts (underpredicts) the abundance of low (high) mass objects, it can be used to estimate the fractional deviation from Gaussianity. In this formalism, the non-Gaussian fractional correction to the multiplicity function isR(ν,fNLX)≡f(ν,fNLX)/f(ν,0)=(dP/dM)(>δc,M,fNLX)/(dP/dM)(>δc,M,0), which is readily computed once the non-Gaussian density PDF P(δM) is known. In the simple extensions proposed by [105, 106], P(δM) is expressed as the inverse transform of a cumulant-generating function. In [106], the saddle-point technique is applied directly to P(δM). The resulting Edgeworth expansion is then used to obtain P(>δc,M). Neglecting cumulants beyond the skewness, one obtain (we momentarily drop the subscript M for convenience)(22)RLV(ν,fNLX)≈1+16σS3(ν3-3ν)-16d(σS3)dlnν(ν-1ν),
after integration over regions above δc(z). In [105], it is the level excursion probability P(>δc,M) that is calculated within the saddle-point approximation. This approximation better asymptotes to the exact large mass tail, which exponentially deviates from the Gaussian tail. To enforce the normalization of the resulting mass function, one may define ν⋆=δ⋆/σ with δ⋆=δc1-S3δc/3, and use [105, 107](23)ν⋆f(ν⋆)=M2nNG(M,z)ρ¯dlnMdlnν⋆.
The fractional deviation from the Gaussian mass function then becomes(24)RMVJ(ν,fNLX)≈exp(ν36σS3)[-σν26ν⋆dS3dlnν+ν⋆ν].
Both formulae have been shown to give reasonable agreement with numerical simulations of non-Gaussian cosmologies [41, 108, 109] (but note that [11, 110] have reached somewhat different conclusions). Expanding δ⋆ at the first order in fNLX shows that these two theoretical expectations differ in the coefficient of the νσS3 term. Therefore, it is interesting to consider also the approximation(25)R(ν,fNLX)≈exp(ν36σS3)[1-ν2σS3-ν6d(σS3)dlnν],
which is designed to match better the Edgeworth expansion of [106] when the peak height is ν~1 [48]. When the primordial trispectrum is large (i.e., when gNLloc~106), terms involving the kurtosis must be included [48, 105, 106, 111]. In this case, it is also important to take into account a possible renormalization of the fluctuation amplitude, σ8→σ8+δσ8 (12), to which the high mass tail of the mass function is exponentially sensitive [48].Figure4 shows the effect of primordial NG of the local fNLloc type on the mass function of SO halos identified with a redshift-dependent overdensity threshold Δvir(z) (motivated by spherical collapse in a ΛCDM cosmology [112]). Overall, the approximation (25) agrees better with the measurements than (24), which somewhat overestimates the data for fNLloc=100, and than (22), which is not always positive definite for fNLloc=-100. However, as can be seen in Figure 5, while the agreement with the data is reasonable for the SO halos, the theory strongly overestimates the effect in the mass function of FOF halos. Reference [109], who use a FOF algorithm with b=0.2, makes the replacement δc→δcq with q≃0.75 to match the Gaussian and non-Gaussian mass functions. A physical motivation of this replacement is provided by [113, 114], who demonstrate that the diffusive nature of the collapse barrier introduces a similar factor q=(1+DB)-1, regardless of the initial conditions. However, the value of the diffusion constant DB was actually measured from simulations that use a SO finder with constant Δvir=200 [115]. In the excursion set approach of [116], the value of q is obtained by normalizing the Gaussian mass function to simulation (i.e., it has nothing to do with the collapse dynamics) and, therefore, depends on the halo finder. Figure 5 demonstrates that this is also the case for the non-Gaussian correction R(ν,fNLloc): choosing q≃0.75 as advocated in [109] gives good agreement for FOF halos, but strongly underestimates the effect for SO halos, for which the best-fit q is close to unity. As we will see below, the strength of the non-Gaussian bias may also be sensitive to the choice of halo finder.Figure 4
Fractional deviation from the Gaussian mass function as a function of the peak heightν=δc/σ. Different symbols refer to different redshifts as indicated. The various curves are theoretical prediction at z=0 (see text). Halos were identified using a spherical overdensity (SO) finder with a redshift-dependent overdensity threshold Δvir(z) (with Δvir(z) increasing from ~200 at high redshift to attain ~350 at z=0). Error bars denote Poisson errors. For illustration, M=1015M⊙/h corresponds to a significance ν=3.2, 5.2, 7.7 at redshift z=0, 1 and 2, respectively. Similarly, M=1014M⊙/h and 1013M⊙/h correspond to ν=1.9, 3, 4.5 and 1.2, 1.9, 2.9, respectively.Figure 5
Fractional deviation from the Gaussian mass function as a function of the peak heightν=δc/σ. Different symbols refer to different redshifts as in Figure 4. The curves are the theoretical prediction (25) at z=0 with q=1 (solid) and q=0.75 (dotted). In the top panel, halos were identified using a spherical overdensity (SO) finder with a redshift-dependent overdensity threshold Δvir(z) whereas, in the bottom panel, a Friends-of-Friends (FOF) finding algorithm with linking length b=0.2 was used.More sophisticated formulations based on extended Press-Schechter (EPS) theory and/or modifications of the collapse criterion look promising since they can reasonably reproduce both the Gaussian halo counts and the dependence onfNLX [114, 117, 118]. The probability of first upcrossing can, in principle, be derived for any non-Gaussian density field and any choice of smoothing filter [119, 120]. For a general filter, the non-Markovian dynamics generates additional terms in the non-Gaussian correction to the mass function that arise from 3-point correlators of the smoothed density δM at different mass scales [114]. However, large error bars still make difficult to test for the presence of such subleading terms. For generic moving barriers B(σ) such as those appearing in models of triaxial collapse [121, 122], the leading contribution to the non-Gaussian correction approximately is [117](26)R(ν,fNLX)≈1+16σS3H3(B(σ)σ)≈1+16σS3q(qν3-3ν),
where H3(ν)≡ν3-3ν and the last equality assumes ν≫1. For SO halos, (26) with q~0.7 does not fit to the measured correction R(ν,fNLloc) better than (25). However, the ellipsoidal collapse barrier with q~0.7 matches better the Gaussian mass function for moderate peak height ν≲2 [118].Parameterizations of the fractional correction based on N-body simulations have also been considered. While [43] considers a fourth-order polynomial fit to account for values of fNLloc as large as 750, [11] exploits the fact that, for sufficiently small fNLloc, there is a one-to-one mapping between halos in Gaussian and non-Gaussian cosmologies. In both cases, the fitting functions are consistent with the simulations at the few percent level.


## 4.2.2. Clusters Abundance

Rich clusters of galaxies trace the rare, high-density peaks in the initial conditions and thus offer the best probe of the high-mass tail of the multiplicity function. To infer the cluster mass function, the X-ray and millimeter windows are better suited than the optical-wave range because selection effects can be understood better (see, however, [123]).Following early theoretical [96, 97, 124–126] and numerical [35, 127–129] work on the effect of non-Gaussian initial conditions on the multiplicity function of cosmic structures, the abundance of clusters and X-ray counts in non-Gaussian cosmologies has received much attention in the literature. At fixed normalization of the observed abundance of local clusters, the protoclusters associated with high redshift (2<z<4) Lyα emitters are much more likely to develop in strongly non-Gaussian models than in the Gaussian paradigm [39, 110, 130]. Considering the redshift evolution of cluster abundances thus can break the degeneracy between the initial density PDF and the background cosmology. Simple extensions of the Press-Schechter formalism similar to those considered above have been shown to capture reasonably well the cluster mass function over a wide range of redshift for various non-Gaussian scenarios [131]. Scaling relations between the cluster mass, X-ray temperature and Compton y-parameter calibrated using theory, observations and detailed simulations of cluster formation [132, 133], can be exploited to predict the observed distribution functions of X-ray and SZ signals and assess the capability of cluster surveys to test the nature of the initial conditions [134–140].An important limitation of this method is that the impact of realistic models of primordial non-Gaussianity on cluster abundances is small compared to systematics in current and upcoming surveys [141–143]. Given the current uncertainties in the redshift evolution of clusters (one commonly assumes that clusters are observed at the epoch they collapse [142]), the selection effects in the calibration of X-ray and SZ fluxes with halo mass, the freedom in the definition of the halo mass, the degeneracy with the normalization amplitude σ8 (for positive fNLX, the mass function is more enhanced at the high mass end, and this is similar to an increase in the amplitude of fluctuations σ8 [144]) and the low number statistics, the prospects of using the cluster mass function only to place competitive limits on fNLX with the current data are small. A two-fold improvement in cluster mass calibration is required to provide constraints comparable to CMB measurements [143].


## 4.2.3. Voids Abundance

The frequency of cosmic voids offers a probe of the low density tail of the initial PDF [145]. The Press-Schechter formalism can also be applied to ascertain the sensitivity of the voids abundance to non-Gaussian initial conditions. Voids are defined as regions of mass M whose density is less than some critical value δv≤0 or, alternatively, whose three eigenvalues of the tidal tensor [146] lie below some critical value λv≤0 [66, 118, 145, 147]. An important aspect in the calculation of the mass function of voids is the overcounting of voids located inside collapsing regions. This voids-in-clouds problem, as identified by [148]), can be solved within the excursion set theory by studying a two barriers problem: δc for halos and δv for voids. Including this effect reduces the frequency of the smallest voids [118]. Neglecting this complication notwithstanding, the differential number density of voids of radius R is [145, 147](27)dndR≈92π2π2|νv|R4e-νv2/2dln|νv|dlnM[1-16σS3H3(|νv|)],
where νv=δv/σM. While a positive fNLX produces more massive halos, it generates fewer large voids [118, 145]. Hence, the effect is qualitatively different from a simple rescaling of the normalization amplitude σ8. A joint analysis of both abundances of clusters and cosmic voids might thus provide interesting constraints on the shape of the primordial 3-point function. There are, however, several caveats to this method, including the fact that there is no unique way to define voids [145]. Clearly, voids identification algorithms will have to be tested on numerical simulations [149] before a robust method can be applied to real data.


## 4.3. Galaxy 2-Point Correlation

Before [90] showed that, in Gaussian cosmologies, correlations of galaxies and clusters can be amplified relative to the mass distribution, it was argued that primeval fluctuations have a non-Gaussian spectrum [150, 151] to explain the observed strong correlation of Abell clusters [152, 153]. Along these lines, [154] pointed out that primordial non-Gaussianity could significantly increase the amplitude of the two-point correlation of galaxies and clusters on large-scales. However, except from [155] who showed that correlations of high density peaks in non-Gaussian models are significantly stronger than in the Gaussian model with identical mass power spectrum, subsequent work focused mostly on abundances (Section 4.2) or higher order statistics such as the bispectrum (Section 4.4). It is only recently that [11] have demonstrated the strong scale-dependent bias arising in non-Gaussian models of the local fNLloc type.


### 4.3.1. The Non-Gaussian Bias

In the original derivation of [11], the Laplacian is applied to the local mapping Φ=ϕ+fNLlocϕ2 in order to show that, upon substitution of the Poisson equation, the overdensity in the neighborhood of density peaks is spatially modulated by a factor proportional to the local value of ϕ. Taking into account the coherent motions induced by gravitational instabilities, the scale-dependent bias correction reads(28)Δbκ(k,fNLloc)=3fNLloc[b1(M)-1]δc(0)ΩmH02k2T(k)D(z),
where b1(M) is the linear, Eulerian bias of halos of mass M. The original result missed out a multiplicative factor of T(k)-1 which was introduced subsequently by [157] upon a derivation of (28) in the limit of high density peaks. The peak-background split approach [103, 158, 159] promoted by [156] shows that the scale-dependent bias applies to any tracer of the matter density field whose (Gaussian) multiplicity function depends on the local mass density only. In this approach, the Gaussian piece of the potential is decomposed into short- and long-wavelength modes, ϕ=ϕl+ϕs. The short-wavelength piece of the density field is then given by the convolution(29)δs=ℳ⋆Φs=ℳ⋆ϕs(1+2fNLlocϕl)+fNLlocℳ⋆ϕs2,
where ℳ is the transfer function (7). Ignoring the white-noise term, this provides an intuitive explanation of the effect in terms of a local rescaling of the small-scale amplitude of matter fluctuations or, equivalently, a local rescaling of the critical density threshold,(30)σs→σs(1+2fNLlocϕl(x)).
Assuming that the mass function depends only on the peak height ν=δc/σs, the long-wavelength part of the halo overdensity becomes [156] (see also [11, 71, 160])(31)δlh(x)=1n(ν)n(δc-δl(x)σs(1+2fNLlocϕl(x)))-1≈-1σs(δl(x)+2fNLlocϕl(x))dlnndν.
Upon a Fourier transformation and using the fact that, in the Gaussian case, δlh(k)=bLδl(k) with the Lagrangian bias bL=-σs-1dlnn/dν, we recover the non-Gaussian bias correction (28) provided that the tracers move coherently with the dark matter, that is, bL=b1(M)-1 [161]. As emphasized in [11], the scale-dependence arises from the fact that the non-Gaussian curvature perturbation Φ(x) is related to the density through the Poisson equation (6) (so that δl(k)=ℳ(k)ϕl(k)). There is no such effect in the (local) χ2 model, (10), nor in texture-seeded cosmologies [162] for instance.The derivation of [157], based on the clustering of regions of the smoothed density field δM above threshold δc(z), is formally valid for high density peaks only. However, it is general enough to apply to any shape of primordial bispectrum. The 2-point correlation function of that level excursion set, which was first derived by [124], can be expressed in the high threshold limit (ν≫1) as(32)ξ>ν(r)=-1+exp{∑n=2∞∑j=1n-1νnσ-nj!(n-j)!×ξR(n)(x1,…,x1,x2,…,x2jtimes(n-j)times,z=0)∑n=2∞},
where r=x1-x2. For most non-Gaussian models in which the primordial 3-point function is the dominant correction, this expansion can be truncated at the third order and Fourier transformed to yield the non-Gaussian correction ΔP>ν(k) to the power spectrum. Assuming a small level of primordial NG, we can also write ΔP>ν(k)≈2bLΔbκPR(k) where bL≈ν2/δc, and eventually obtain(33)Δbκ(k,fNLX)≡bϕ(k)ℱ(k,fNLX)=(2bLδc(z)ℳR(k,0))ℱ(k,fNLX).
The dependence on the shape of the 3-point function is encoded in the function ℱ(k,fNLX) [157, 163](34)ℱ(k,fNLX)=116π2σ2∫0∞dk1k12ℳR(k1,0)×∫-1+1dμℳR(α,0)BΦ(k1,α,k)PΦ(k),
where α2=k2+k12+2μkk1. Note that, for fNLloc<0, this first order approximation always breaks down at sufficiently small k because ΔP>ν(k)<0.Figure6 shows the non-Gaussian halo bias (33) induced by the local, equilateral and folded bispectrum [163]. In the local and folded non-Gaussianity, the deviation is negligible at k=0.1hMpc-1, but increases rapidly with decreasing wavenumber. Still, the large-scale correction is much smaller for the folded template, and nearly absent for the equilateral type, which make them much more difficult to detect with galaxy surveys [163]. To get insights into the behavior of Δbκ(k,fNLX) at large-scales, let us identify the dominant contribution to ℱ(k,fNLX) in the limit k≪1. Setting ℳR(α,0)≈ℳR(k1,0) and expanding Pϕ(α) at second order in k/k1, we find after some algebra(35)ℱ(k,fNLloc)≈fNLloc,ℱ(k,fNLeq)≈fNLeq[3ΣR(2(ns-4)3)k2(4-ns)/3+12(ns-4)ΣR(-2)k2]σR-2,ℱ(k,fNLfol)≈32fNLfolΣR(ns-43)k(4-ns)/3σR-2,
assuming no running scalar index, that is, dns/dlnk=0. The auxiliary function ΣR(n) is defined as(36)ΣR(n)=12π2∫0∞dkk(2+n)ℳR(k,0)2Pϕ(k).
Hence, we have ΣR(0)≡σR2. As can be seen in Figure 6, these approximations capture relatively well the large scale non-Gaussian bias correction induced by the equilateral and folded type of non-Gaussianity. For a nearly scale-invariant spectrum ns≈1, the effect scales as Δbκ∝k and Δbκ∝ const., respectively. Another important feature of the equilateral and folded non-Gaussian bias is the dependence on the mass scale M through the multiplicative factor σR-2. Indeed, choosing R=1h-1Mpc instead of R=5h-1Mpc as done in Figure 6 would suppress the effect by a factor of ~3. In the high peak limit, σR-2≈bL/δc(z) which cancels out the dependence on redshift but enhances the sensitivity to the halo bias, that is, Δbκ∝bL2 for the equilateral and folded shapes whereas Δbκ∝bL in the local model.Figure 6
Non-Gaussian bias correction (33) for the local, equilateral and folded primordial bispectrum as a function of wavenumber. Results are shown at z=0 for a smoothing radius R=5h-1Mpc and a nonlinear parameter fNLX=1. The dotted line represent an analytic approximation, (35), which is valid at large-scales k≪1. Note that while the magnitude of Δbκ(k,fNLloc) does not change with R, Δbκ(k,fNLeq) and Δbκ(k,fNLfol) strongly depends on the smoothing radius (see text).At this point, it is appropriate to mention a few caveats to these calculations. Firstly, (28) assumes that the tracers form after a spherical collapse, which may be a good approximation for the massive halos only. If one instead considers the ellipsoidal collapse dynamics, in which the evolution of a perturbation depends upon the three eigenvalues of the initial tidal shear, δc(0) should be replaced by its ellipsoidal counterparts δec(0) which is always larger than the spherical value [121]. In this model, the scale-dependent bias Δbκ is thus enhanced by a factor δec(0)/δc(0) [11, 160]. Secondly, (28) assumes that the biasing of the surveyed objects is described by the peak height ν only or, equivalently, the hosting halo mass M. However, this may not be true for quasars whose activity may be triggered by merger of halos [164, 165]. Reference [156] used the EPS formalism to estimate the bias correction Δbmerger induced by mergers(37)Δbmerger=δc-1,
so the factor b1(M)-1 should be replaced by b1(M)-1-δc-1≈b1(M)-1.6. The validity of this result should be evaluated with cosmological simulations of quasars formation. In this respect, semianalytic models of galaxy formation suggest that merger-triggered objects such as quasars do not cluster much differently than other tracers of the same mass [166]. However, this does not mean that the same should hold for the non-Gaussian scale-dependent bias. Still, since the recent merger model is an extreme case it seems likely that the actual bias correction is 0<Δbmerger<δc-1. Thirdly, the scale-dependent bias has been derived using the Newtonian approximation to the Poisson equation, so one may wonder whether general relativistic (GR) corrections to ℳR(k)-1 may suppress the effect on scales comparable to the Hubble radius. Reference [167] showed how large-scale primordial NG induced by GR corrections propagates onto small scales once cosmological perturbations reenter the Hubble radius in the matter dominated era. This effect generates a scale-dependent bias comparable, albeit of opposite sign to that induced by local NG [163]. More recently, [168] argues that there are no GR corrections to the non-Gaussian bias and that the scaling Δbκ∝k-2 applies down to smallish wavenumbers.We can also ask ourselves whether higher-order terms in the series expansion (32) furnish corrections to the non-Gaussian bias similar to (28). The quadratic coupling fNLlocϕ2 induces a second order correction to the halo power spectrum which reads [48](38)ΔPh(k)=43(fNLloc)2[b1(M)-1]2δc2(z)S3(1)(M)×ℳR(k,0)Pϕ(k).
Its magnitude relative to the term linear in fNLloc, (28), is approximately 0.03 at redshift z=1.8 and for a halo mass M=1013M⊙/h. Although its contribution becomes increasingly important at higher redshift, it is fairly small for realistic values of fNLloc. In local NG model, the power spectrum of biased tracers of the density field can also be obtained from a local Taylor series in the evolved (Eulerian) density contrast δ and the Gaussian part ϕ of the initial (Lagrangian) curvature perturbation [47, 71]. Using this approach, it can be shown that the halo power spectrum arising from the first order terms of the local bias expansion can be cast into the form [47](39)Ph(k)=[b1(M)+fNLlocbϕ(k)]2PR(k).
Hence, we also obtain a second order term proportional to (fNLloc)2ℳR-2PR(k)=(fNLloc)2Pϕ(k) which, however, contributes only at very small wavenumber k≲0.001h-1Mpc. All this suggests that (28) is the dominant contribution to the non-Gaussian bias in the wavenumber range 0.001≲k≲0.1hMpc-1.Finally, a non-Gaussian, scale-dependent bias correction can also arise in the local, deterministic bias ansatzδh(x)=b1δ(x)+b2δ(x)2/2+⋯ [169] if the initial density field is non-Gaussian. Here, bN is the Nth-order bias parameters (here again, the first-order bias is b1≡1+bL). In this approach, the correction is induced by the correlation b1b2〈δ(x1)δ2(x2)〉 between the linear and quadratic term in the galaxy biasing relation (which is in fact a collapsed or squeezed 3-point function) and thus reads [67, 68](40)Δbκ(k,fNLloc)=2fNLlocb2σR2ℳR(k,0)-1.
Even though b2σR2≈bLδc in the high-threshold limit ν≫1, b2σR2 behaves very differently than bLδc for moderate peak height because b2 is proportional to the second derivative of the mass function n(ν). So far however, (28) appears to describe reasonably well the numerical results for a wide range of halo bias.


### 4.3.2. Comparison with Simulations

In order to fully exploit the potential of forthcoming large-scale surveys, a number of studies have tested the theoretical prediction against the outcome of large numerical simulations [11, 41–43, 71, 109].At the lowest order, there are two additional albeit relatively smaller corrections to the Gaussian bias which arise from the dependence of both the halo number densityn(M,z) and the matter power spectrum Pδ(k,z) on primordial NG [41]. Firstly, assuming the peak-background split holds, the change in the mean number density of halos induces a scale-independent offset which we denote ΔbI(fNLloc). In terms of the non-Gaussian fractional correction R(ν,fNLloc) to the mass function, this contribution is(41)ΔbI(fNLloc)=-1σ∂∂νln[R(ν,fNLloc)].
It is worth noticing that ΔbI(fNLloc) has a sign opposite to that of fNLloc, because the bias decreases when the mass function goes up. In practice, the approximation (25), which matches well the SO data for ν≲4, can be used for moderate values of the peak height. For FOF halos with linking length b=0.2, one should make the replacement δc→δcq with q≈0.75 in the calculation of the scale-independent offset. It is sensible to evaluate ΔbI(fNLloc) at a mass scale 〈M〉 equal to the average halo mass of the sample. Secondly, we also need to account for the change in the matter power spectrum (see Figure 2 in Section 3). Summarizing, local non-Gaussianity adds a correction Δb(k,fNLloc) to the bias b(k) of dark matter halos that reads [41](42)Δb(k,fNLloc)=Δbκ(k,fNLloc)+ΔbI(fNLloc)+b1(M)βδ(k,fNLloc)
at first order in fNLloc. As can be seen in Figure 7, the inclusion of these extra terms substantially improves the comparison between the theory and the simulations. Considering only the scale-dependent shift Δbκ leads to an apparent suppression of the effect in simulations relative to the theory. Including the scale-independent offset ΔbI considerably improves the agreement at wavenumbers k≲0.05hMpc-1. Finally, adding the scale-dependent term b1(M)βδ further adjusts the match at small scale k≳0.05hMpc-1 by making the non-Gaussian bias shift less negative. Along these lines, [71] find that the inclusion of ΔbI to the bias also improves the agreement with measurements of Δb(k,fNLloc) obtained for FOF halos, and show that taking into account second- and higher-order corrections could extend the validity of the theory up to scales k~0.1-0.3hMpc-1.Figure 7
Non-Gaussian bias correction (filled symbols) for halos of massM>2×1013M⊙/h extracted at z=0.5 from simulations of the local fNLloc model. The solid curve represents the theoretical model (42). The dashed, dotted-dashed-dotted, and dotted curves show the scale-dependent bias Δbκ, the scale-independent offset ΔbI and the contribution from the matter power spectrum b1(M)βδ that arise at first order in fNLloc (see Figure 2). The error bars indicate the scatter among 8 realizations of 10243 simulations with box size L=1600h-1Mpc.The non-Gaussian bias correction can be measured in the cross- and autopower spectrum of dark matter halos,Phδ(k) and Ph(k). To compute these quantities, dark matter particles and halo centers are interpolated onto a regular cubical mesh. The resulting dark matter and halo fluctuation fields, δm(k) and δh(k), are then Fourier transformed to yield the matter-matter, halo-matter and halo-halo power spectra Pδ(k), Phδ(k) and Ph(k), respectively. Ph(k) is then corrected for the shot noise, which is assumed to be 1/n¯h if dark matter halos are a Poisson sampling of some continuous field. This discreteness correction is negligible for Pδ(k) due to the large number of dark matter particles. On linear scales (k≲0.01hMpc-1), the halo bias b(k)=δh(k)/δm(k) approaches the constant value b1(M) which needs to be measured accurately as it controls the strength of the scale-dependent bias correction Δbκ. In this respect, the ratio Phδ(k)/Pδ(k) is a better proxy for the halo bias since it is less sensitive to shot-noise.Auto- and crosspower analyses may not agree with each other if the halos and dark matter do not trace each other on scalek≲0.01hMpc-1 where the non-Gaussian bias is large, that is, if there is stochasticity. Figure 8 shows Phδ(k) and Ph(k) averaged over 8 realizations of the models with fNLloc=0,±100 [41]. The same Gaussian random seed field ϕ was used in each set of runs so as to minimize the sampling variance. Measurements of the non-Gaussian bias correction obtained with the halo-halo or the halo-matter power spectrum are in a good agreement with each other, indicating that non-Gaussianity does not induce stochasticity and the predicted scaling (28) applies equally well for the auto- and cross-power spectrum. However, while a number of numerical studies of the fNLloc model have confirmed the scaling Δbκ(k,fNLloc)∝ℳR(k)-1 and the redshift dependence ∝D(z)-1 [11, 41, 43, 109], the exact amplitude of the non-Gaussian bias correction remains somewhat debatable. Reference [41] who use SO halos and [71] who use FOF halos find satisfactory agreement with the theory once the scale-independent offset ΔbI is included. By contrast, see [43], who use the same FOF halos as [71], argue that the scale-dependent piece Δbκ requires, among others, a multiplicative correction of the form (1-β1fNLloc), with β1~4×10-4>0. Similarly, [109] who also use FOF halos find that the theory is a good fit to the simulations only upon replacing bL by qbL in (33), with q≃0.75. Part of the discrepancy may be probably due to the fact that the last two references do not include ΔbI, which leads to an apparent suppression of the effect (see Figure 7). Another possible source of discrepancy may be choice of the halo finder which, as seen in Figure 5, has an impact on the strength of the non-Gaussian correction to the mass function. This possibility is investigated in Figure 9, which shows the non-Gaussian bias correction obtained with FOF halos. For this low biased sample, the scale-independent correction is |ΔbI|≲0.003 and can thus be neglected. The best-fit values of fNLloc are significantly below the input values of ±100, in agreement with the findings of [43, 109] (note, however, that this suppression is more consistent with δc being rescaled by qδc≈0.86δc and bL being unchanged). This indicates that the choice of halo finder may also affect the magnitude of the scale-dependent non-Gaussian bias. Discrepancies have also been observed between the theoretical and measured non-Gaussian bias corrections in non-Gaussian models of the local cubic-order coupling gNLlocϕ3 [48]. Understanding all these results clearly requires a better modeling of halo clustering.Figure 8
Halo-halo (solid curve) and halo-matter (dashed curve) power spectraPh(k) and Phδ(k) measured in simulations of the Gaussian and fNLloc=±100 models for halos of mass M>2×1013M⊙/h at redshift z=1. The error bars represent the scatter among 8 realizations. For fNLloc=-100, the crosspower spectrum is negative on scales k≲0.005hMpc-1, in good agreement with the theoretical prediction.Figure 9
Fractional correction to the Gaussian halo bias in thefNLloc=±100 and Gaussian models. In contrast to Figure 8, halos were identified with a FOF finder of linking length b=0.2. Only the wavemodes to the left of the vertical line were used to fit Δbκ(k,fNLloc). For this low biased sample, the scale-independent correction is |ΔbI|≲0.003 and can thus be ignored. The best-fit value of fNLloc and the corresponding 1σ error is quoted for each model (Figure taken from [170]).


### 4.3.3. Redshift-Space Distortions

Peculiar velocities generate systematic differences between the spatial distribution of data in real and redshift space. These redshift-space distortions must be properly taken into account in order to extractfNLX from redshift surveys. On the linear scales of interest, the redshift-space power spectrum of biased tracers reads as [171, 172](43)Phs(k,μ)=[b12Pδ(k)+2b1fμ2Pδθ(k)+f2μ4Pθ(k)],
where Pδθ and Pθ are the density-velocity and velocity-divergence power spectra, μ is the cosine of the angle between the wavemode k and the line of sight and f is the logarithmic derivative of the growth factor. For Pθ, the 1-loop correction due to primordial NG is identical to (17) provided F2(k1,k2) is replaced by the kernel G2(k1,k2)=3/7+μ(k1/k2+k2/k1)/2+4μ2/7 describing the 2nd order evolution of the velocity divergence [62]. For Pδθ, this correction is(44)ΔPδθNG(k)=∫d3q(2π)3[F2(q,k-q)+G2(q,k-q)]×B0(-k,q,k-q).
Again, causality implies that G2(k1,k2) vanishes in the limit k1=-k2. For unbiased tracers with b1=1, the linear Kaiser relation is thus recovered at large scales k≲0.01hMpc-1 (this is consistent with the analysis of [173]). For biased tracers, we still expect the Kaiser formula to be valid, but the distortion parameter β should now be equal to β=f/(b1+Δbκ), where Δbκ(k,fNLX) is the scale-dependent bias induced by the primordial non-Gaussianity.


### 4.3.4. Mitigating Cosmic Variance and Shot-Noise

Because of the finite number of large scale wavemodes accessible to a survey, any large-scale measurement of the power spectrum is limited by the cosmic (or sampling) variance caused by the random nature of the wavemodes. For discrete tracers such as galaxies, the shot noise is another source of error. Restricting ourselves to weak primordial NG, the relative error on the power spectrumP is σP/P≈1/N(1+σn2/P), where N is the number of independent modes measured and σn2 is the shot-noise [174]. Under the standard assumption of Poisson sampling, σn2 equals the inverse of the number density 1/n¯ and causes a scale-independent enhancement of the power spectrum. The extent to which one can improve the observational limits on the nonlinear parameters fNLX will strongly depend on our ability to minimize the impact of these two sources of errors. By comparing differently biased tracers of the same surveyed volume [175, 176] and suitably weighting galaxies (by the mass of their host halo for instance) [177, 178], it should be possible to circumvent these problems and considerably improve the detection level.Figure9 illustrates how the impact of sampling variance on the measurement of fNLloc can be mitigated. Namely, the data points show the result of taking the ratio Ph(k,fNLloc)/Pδ(k,fNLloc) for each set of runs with same Gaussian random seed field ϕ before averaging over the realizations. This procedure is equivalent to the multitracers method advocated by [175]. Here, Pδ can be thought as mimicking the power spectrum of a nearly unbiased tracer of the mass density field with high number density. Although, in practical applications, using the dark matter field works better [170], in real data Pδ should be replaced by a tracer of the same surveyed volume different than the one used to compute Ph. Figure 9 also shows that upon taking out most of the cosmic variance, there is some residual noise caused by the discrete nature of the dark matter halos. As shown recently [178] however, weighting the halos according to their mass can dramatically reduce the shot noise relative to the Poisson expectation, at least when compared against the dark matter. Applying such a weighting may thus significantly improve the error on the nonlinear parameter fNLloc, but this should be explored in realistic simulations of galaxies, especially because the halo mass M may not be easily measurable from real data [170]. This approach undoubtedly deserves further attention as it has the potential to substantially improve the extraction of the primordial non-Gaussian signal from galaxy surveys.To conclude this section, it is worth noting that, while the PDF of power valuesP(k) has little discriminatory power (for large surveyed volume, it converges towards the Rayleigh distribution as a consequence of the central limit theorem) [179], the covariance of power spectrum measurements (which is sensitive to the selection function, but also to correlations among the phase of the Fourier modes) may provide quantitative limits on certain type of non-Gaussian models [174, 180].


## 4.3.1. The Non-Gaussian Bias

In the original derivation of [11], the Laplacian is applied to the local mapping Φ=ϕ+fNLlocϕ2 in order to show that, upon substitution of the Poisson equation, the overdensity in the neighborhood of density peaks is spatially modulated by a factor proportional to the local value of ϕ. Taking into account the coherent motions induced by gravitational instabilities, the scale-dependent bias correction reads(28)Δbκ(k,fNLloc)=3fNLloc[b1(M)-1]δc(0)ΩmH02k2T(k)D(z),
where b1(M) is the linear, Eulerian bias of halos of mass M. The original result missed out a multiplicative factor of T(k)-1 which was introduced subsequently by [157] upon a derivation of (28) in the limit of high density peaks. The peak-background split approach [103, 158, 159] promoted by [156] shows that the scale-dependent bias applies to any tracer of the matter density field whose (Gaussian) multiplicity function depends on the local mass density only. In this approach, the Gaussian piece of the potential is decomposed into short- and long-wavelength modes, ϕ=ϕl+ϕs. The short-wavelength piece of the density field is then given by the convolution(29)δs=ℳ⋆Φs=ℳ⋆ϕs(1+2fNLlocϕl)+fNLlocℳ⋆ϕs2,
where ℳ is the transfer function (7). Ignoring the white-noise term, this provides an intuitive explanation of the effect in terms of a local rescaling of the small-scale amplitude of matter fluctuations or, equivalently, a local rescaling of the critical density threshold,(30)σs→σs(1+2fNLlocϕl(x)).
Assuming that the mass function depends only on the peak height ν=δc/σs, the long-wavelength part of the halo overdensity becomes [156] (see also [11, 71, 160])(31)δlh(x)=1n(ν)n(δc-δl(x)σs(1+2fNLlocϕl(x)))-1≈-1σs(δl(x)+2fNLlocϕl(x))dlnndν.
Upon a Fourier transformation and using the fact that, in the Gaussian case, δlh(k)=bLδl(k) with the Lagrangian bias bL=-σs-1dlnn/dν, we recover the non-Gaussian bias correction (28) provided that the tracers move coherently with the dark matter, that is, bL=b1(M)-1 [161]. As emphasized in [11], the scale-dependence arises from the fact that the non-Gaussian curvature perturbation Φ(x) is related to the density through the Poisson equation (6) (so that δl(k)=ℳ(k)ϕl(k)). There is no such effect in the (local) χ2 model, (10), nor in texture-seeded cosmologies [162] for instance.The derivation of [157], based on the clustering of regions of the smoothed density field δM above threshold δc(z), is formally valid for high density peaks only. However, it is general enough to apply to any shape of primordial bispectrum. The 2-point correlation function of that level excursion set, which was first derived by [124], can be expressed in the high threshold limit (ν≫1) as(32)ξ>ν(r)=-1+exp{∑n=2∞∑j=1n-1νnσ-nj!(n-j)!×ξR(n)(x1,…,x1,x2,…,x2jtimes(n-j)times,z=0)∑n=2∞},
where r=x1-x2. For most non-Gaussian models in which the primordial 3-point function is the dominant correction, this expansion can be truncated at the third order and Fourier transformed to yield the non-Gaussian correction ΔP>ν(k) to the power spectrum. Assuming a small level of primordial NG, we can also write ΔP>ν(k)≈2bLΔbκPR(k) where bL≈ν2/δc, and eventually obtain(33)Δbκ(k,fNLX)≡bϕ(k)ℱ(k,fNLX)=(2bLδc(z)ℳR(k,0))ℱ(k,fNLX).
The dependence on the shape of the 3-point function is encoded in the function ℱ(k,fNLX) [157, 163](34)ℱ(k,fNLX)=116π2σ2∫0∞dk1k12ℳR(k1,0)×∫-1+1dμℳR(α,0)BΦ(k1,α,k)PΦ(k),
where α2=k2+k12+2μkk1. Note that, for fNLloc<0, this first order approximation always breaks down at sufficiently small k because ΔP>ν(k)<0.Figure6 shows the non-Gaussian halo bias (33) induced by the local, equilateral and folded bispectrum [163]. In the local and folded non-Gaussianity, the deviation is negligible at k=0.1hMpc-1, but increases rapidly with decreasing wavenumber. Still, the large-scale correction is much smaller for the folded template, and nearly absent for the equilateral type, which make them much more difficult to detect with galaxy surveys [163]. To get insights into the behavior of Δbκ(k,fNLX) at large-scales, let us identify the dominant contribution to ℱ(k,fNLX) in the limit k≪1. Setting ℳR(α,0)≈ℳR(k1,0) and expanding Pϕ(α) at second order in k/k1, we find after some algebra(35)ℱ(k,fNLloc)≈fNLloc,ℱ(k,fNLeq)≈fNLeq[3ΣR(2(ns-4)3)k2(4-ns)/3+12(ns-4)ΣR(-2)k2]σR-2,ℱ(k,fNLfol)≈32fNLfolΣR(ns-43)k(4-ns)/3σR-2,
assuming no running scalar index, that is, dns/dlnk=0. The auxiliary function ΣR(n) is defined as(36)ΣR(n)=12π2∫0∞dkk(2+n)ℳR(k,0)2Pϕ(k).
Hence, we have ΣR(0)≡σR2. As can be seen in Figure 6, these approximations capture relatively well the large scale non-Gaussian bias correction induced by the equilateral and folded type of non-Gaussianity. For a nearly scale-invariant spectrum ns≈1, the effect scales as Δbκ∝k and Δbκ∝ const., respectively. Another important feature of the equilateral and folded non-Gaussian bias is the dependence on the mass scale M through the multiplicative factor σR-2. Indeed, choosing R=1h-1Mpc instead of R=5h-1Mpc as done in Figure 6 would suppress the effect by a factor of ~3. In the high peak limit, σR-2≈bL/δc(z) which cancels out the dependence on redshift but enhances the sensitivity to the halo bias, that is, Δbκ∝bL2 for the equilateral and folded shapes whereas Δbκ∝bL in the local model.Figure 6
Non-Gaussian bias correction (33) for the local, equilateral and folded primordial bispectrum as a function of wavenumber. Results are shown at z=0 for a smoothing radius R=5h-1Mpc and a nonlinear parameter fNLX=1. The dotted line represent an analytic approximation, (35), which is valid at large-scales k≪1. Note that while the magnitude of Δbκ(k,fNLloc) does not change with R, Δbκ(k,fNLeq) and Δbκ(k,fNLfol) strongly depends on the smoothing radius (see text).At this point, it is appropriate to mention a few caveats to these calculations. Firstly, (28) assumes that the tracers form after a spherical collapse, which may be a good approximation for the massive halos only. If one instead considers the ellipsoidal collapse dynamics, in which the evolution of a perturbation depends upon the three eigenvalues of the initial tidal shear, δc(0) should be replaced by its ellipsoidal counterparts δec(0) which is always larger than the spherical value [121]. In this model, the scale-dependent bias Δbκ is thus enhanced by a factor δec(0)/δc(0) [11, 160]. Secondly, (28) assumes that the biasing of the surveyed objects is described by the peak height ν only or, equivalently, the hosting halo mass M. However, this may not be true for quasars whose activity may be triggered by merger of halos [164, 165]. Reference [156] used the EPS formalism to estimate the bias correction Δbmerger induced by mergers(37)Δbmerger=δc-1,
so the factor b1(M)-1 should be replaced by b1(M)-1-δc-1≈b1(M)-1.6. The validity of this result should be evaluated with cosmological simulations of quasars formation. In this respect, semianalytic models of galaxy formation suggest that merger-triggered objects such as quasars do not cluster much differently than other tracers of the same mass [166]. However, this does not mean that the same should hold for the non-Gaussian scale-dependent bias. Still, since the recent merger model is an extreme case it seems likely that the actual bias correction is 0<Δbmerger<δc-1. Thirdly, the scale-dependent bias has been derived using the Newtonian approximation to the Poisson equation, so one may wonder whether general relativistic (GR) corrections to ℳR(k)-1 may suppress the effect on scales comparable to the Hubble radius. Reference [167] showed how large-scale primordial NG induced by GR corrections propagates onto small scales once cosmological perturbations reenter the Hubble radius in the matter dominated era. This effect generates a scale-dependent bias comparable, albeit of opposite sign to that induced by local NG [163]. More recently, [168] argues that there are no GR corrections to the non-Gaussian bias and that the scaling Δbκ∝k-2 applies down to smallish wavenumbers.We can also ask ourselves whether higher-order terms in the series expansion (32) furnish corrections to the non-Gaussian bias similar to (28). The quadratic coupling fNLlocϕ2 induces a second order correction to the halo power spectrum which reads [48](38)ΔPh(k)=43(fNLloc)2[b1(M)-1]2δc2(z)S3(1)(M)×ℳR(k,0)Pϕ(k).
Its magnitude relative to the term linear in fNLloc, (28), is approximately 0.03 at redshift z=1.8 and for a halo mass M=1013M⊙/h. Although its contribution becomes increasingly important at higher redshift, it is fairly small for realistic values of fNLloc. In local NG model, the power spectrum of biased tracers of the density field can also be obtained from a local Taylor series in the evolved (Eulerian) density contrast δ and the Gaussian part ϕ of the initial (Lagrangian) curvature perturbation [47, 71]. Using this approach, it can be shown that the halo power spectrum arising from the first order terms of the local bias expansion can be cast into the form [47](39)Ph(k)=[b1(M)+fNLlocbϕ(k)]2PR(k).
Hence, we also obtain a second order term proportional to (fNLloc)2ℳR-2PR(k)=(fNLloc)2Pϕ(k) which, however, contributes only at very small wavenumber k≲0.001h-1Mpc. All this suggests that (28) is the dominant contribution to the non-Gaussian bias in the wavenumber range 0.001≲k≲0.1hMpc-1.Finally, a non-Gaussian, scale-dependent bias correction can also arise in the local, deterministic bias ansatzδh(x)=b1δ(x)+b2δ(x)2/2+⋯ [169] if the initial density field is non-Gaussian. Here, bN is the Nth-order bias parameters (here again, the first-order bias is b1≡1+bL). In this approach, the correction is induced by the correlation b1b2〈δ(x1)δ2(x2)〉 between the linear and quadratic term in the galaxy biasing relation (which is in fact a collapsed or squeezed 3-point function) and thus reads [67, 68](40)Δbκ(k,fNLloc)=2fNLlocb2σR2ℳR(k,0)-1.
Even though b2σR2≈bLδc in the high-threshold limit ν≫1, b2σR2 behaves very differently than bLδc for moderate peak height because b2 is proportional to the second derivative of the mass function n(ν). So far however, (28) appears to describe reasonably well the numerical results for a wide range of halo bias.


## 4.3.2. Comparison with Simulations

In order to fully exploit the potential of forthcoming large-scale surveys, a number of studies have tested the theoretical prediction against the outcome of large numerical simulations [11, 41–43, 71, 109].At the lowest order, there are two additional albeit relatively smaller corrections to the Gaussian bias which arise from the dependence of both the halo number densityn(M,z) and the matter power spectrum Pδ(k,z) on primordial NG [41]. Firstly, assuming the peak-background split holds, the change in the mean number density of halos induces a scale-independent offset which we denote ΔbI(fNLloc). In terms of the non-Gaussian fractional correction R(ν,fNLloc) to the mass function, this contribution is(41)ΔbI(fNLloc)=-1σ∂∂νln[R(ν,fNLloc)].
It is worth noticing that ΔbI(fNLloc) has a sign opposite to that of fNLloc, because the bias decreases when the mass function goes up. In practice, the approximation (25), which matches well the SO data for ν≲4, can be used for moderate values of the peak height. For FOF halos with linking length b=0.2, one should make the replacement δc→δcq with q≈0.75 in the calculation of the scale-independent offset. It is sensible to evaluate ΔbI(fNLloc) at a mass scale 〈M〉 equal to the average halo mass of the sample. Secondly, we also need to account for the change in the matter power spectrum (see Figure 2 in Section 3). Summarizing, local non-Gaussianity adds a correction Δb(k,fNLloc) to the bias b(k) of dark matter halos that reads [41](42)Δb(k,fNLloc)=Δbκ(k,fNLloc)+ΔbI(fNLloc)+b1(M)βδ(k,fNLloc)
at first order in fNLloc. As can be seen in Figure 7, the inclusion of these extra terms substantially improves the comparison between the theory and the simulations. Considering only the scale-dependent shift Δbκ leads to an apparent suppression of the effect in simulations relative to the theory. Including the scale-independent offset ΔbI considerably improves the agreement at wavenumbers k≲0.05hMpc-1. Finally, adding the scale-dependent term b1(M)βδ further adjusts the match at small scale k≳0.05hMpc-1 by making the non-Gaussian bias shift less negative. Along these lines, [71] find that the inclusion of ΔbI to the bias also improves the agreement with measurements of Δb(k,fNLloc) obtained for FOF halos, and show that taking into account second- and higher-order corrections could extend the validity of the theory up to scales k~0.1-0.3hMpc-1.Figure 7
Non-Gaussian bias correction (filled symbols) for halos of massM>2×1013M⊙/h extracted at z=0.5 from simulations of the local fNLloc model. The solid curve represents the theoretical model (42). The dashed, dotted-dashed-dotted, and dotted curves show the scale-dependent bias Δbκ, the scale-independent offset ΔbI and the contribution from the matter power spectrum b1(M)βδ that arise at first order in fNLloc (see Figure 2). The error bars indicate the scatter among 8 realizations of 10243 simulations with box size L=1600h-1Mpc.The non-Gaussian bias correction can be measured in the cross- and autopower spectrum of dark matter halos,Phδ(k) and Ph(k). To compute these quantities, dark matter particles and halo centers are interpolated onto a regular cubical mesh. The resulting dark matter and halo fluctuation fields, δm(k) and δh(k), are then Fourier transformed to yield the matter-matter, halo-matter and halo-halo power spectra Pδ(k), Phδ(k) and Ph(k), respectively. Ph(k) is then corrected for the shot noise, which is assumed to be 1/n¯h if dark matter halos are a Poisson sampling of some continuous field. This discreteness correction is negligible for Pδ(k) due to the large number of dark matter particles. On linear scales (k≲0.01hMpc-1), the halo bias b(k)=δh(k)/δm(k) approaches the constant value b1(M) which needs to be measured accurately as it controls the strength of the scale-dependent bias correction Δbκ. In this respect, the ratio Phδ(k)/Pδ(k) is a better proxy for the halo bias since it is less sensitive to shot-noise.Auto- and crosspower analyses may not agree with each other if the halos and dark matter do not trace each other on scalek≲0.01hMpc-1 where the non-Gaussian bias is large, that is, if there is stochasticity. Figure 8 shows Phδ(k) and Ph(k) averaged over 8 realizations of the models with fNLloc=0,±100 [41]. The same Gaussian random seed field ϕ was used in each set of runs so as to minimize the sampling variance. Measurements of the non-Gaussian bias correction obtained with the halo-halo or the halo-matter power spectrum are in a good agreement with each other, indicating that non-Gaussianity does not induce stochasticity and the predicted scaling (28) applies equally well for the auto- and cross-power spectrum. However, while a number of numerical studies of the fNLloc model have confirmed the scaling Δbκ(k,fNLloc)∝ℳR(k)-1 and the redshift dependence ∝D(z)-1 [11, 41, 43, 109], the exact amplitude of the non-Gaussian bias correction remains somewhat debatable. Reference [41] who use SO halos and [71] who use FOF halos find satisfactory agreement with the theory once the scale-independent offset ΔbI is included. By contrast, see [43], who use the same FOF halos as [71], argue that the scale-dependent piece Δbκ requires, among others, a multiplicative correction of the form (1-β1fNLloc), with β1~4×10-4>0. Similarly, [109] who also use FOF halos find that the theory is a good fit to the simulations only upon replacing bL by qbL in (33), with q≃0.75. Part of the discrepancy may be probably due to the fact that the last two references do not include ΔbI, which leads to an apparent suppression of the effect (see Figure 7). Another possible source of discrepancy may be choice of the halo finder which, as seen in Figure 5, has an impact on the strength of the non-Gaussian correction to the mass function. This possibility is investigated in Figure 9, which shows the non-Gaussian bias correction obtained with FOF halos. For this low biased sample, the scale-independent correction is |ΔbI|≲0.003 and can thus be neglected. The best-fit values of fNLloc are significantly below the input values of ±100, in agreement with the findings of [43, 109] (note, however, that this suppression is more consistent with δc being rescaled by qδc≈0.86δc and bL being unchanged). This indicates that the choice of halo finder may also affect the magnitude of the scale-dependent non-Gaussian bias. Discrepancies have also been observed between the theoretical and measured non-Gaussian bias corrections in non-Gaussian models of the local cubic-order coupling gNLlocϕ3 [48]. Understanding all these results clearly requires a better modeling of halo clustering.Figure 8
Halo-halo (solid curve) and halo-matter (dashed curve) power spectraPh(k) and Phδ(k) measured in simulations of the Gaussian and fNLloc=±100 models for halos of mass M>2×1013M⊙/h at redshift z=1. The error bars represent the scatter among 8 realizations. For fNLloc=-100, the crosspower spectrum is negative on scales k≲0.005hMpc-1, in good agreement with the theoretical prediction.Figure 9
Fractional correction to the Gaussian halo bias in thefNLloc=±100 and Gaussian models. In contrast to Figure 8, halos were identified with a FOF finder of linking length b=0.2. Only the wavemodes to the left of the vertical line were used to fit Δbκ(k,fNLloc). For this low biased sample, the scale-independent correction is |ΔbI|≲0.003 and can thus be ignored. The best-fit value of fNLloc and the corresponding 1σ error is quoted for each model (Figure taken from [170]).


## 4.3.3. Redshift-Space Distortions

Peculiar velocities generate systematic differences between the spatial distribution of data in real and redshift space. These redshift-space distortions must be properly taken into account in order to extractfNLX from redshift surveys. On the linear scales of interest, the redshift-space power spectrum of biased tracers reads as [171, 172](43)Phs(k,μ)=[b12Pδ(k)+2b1fμ2Pδθ(k)+f2μ4Pθ(k)],
where Pδθ and Pθ are the density-velocity and velocity-divergence power spectra, μ is the cosine of the angle between the wavemode k and the line of sight and f is the logarithmic derivative of the growth factor. For Pθ, the 1-loop correction due to primordial NG is identical to (17) provided F2(k1,k2) is replaced by the kernel G2(k1,k2)=3/7+μ(k1/k2+k2/k1)/2+4μ2/7 describing the 2nd order evolution of the velocity divergence [62]. For Pδθ, this correction is(44)ΔPδθNG(k)=∫d3q(2π)3[F2(q,k-q)+G2(q,k-q)]×B0(-k,q,k-q).
Again, causality implies that G2(k1,k2) vanishes in the limit k1=-k2. For unbiased tracers with b1=1, the linear Kaiser relation is thus recovered at large scales k≲0.01hMpc-1 (this is consistent with the analysis of [173]). For biased tracers, we still expect the Kaiser formula to be valid, but the distortion parameter β should now be equal to β=f/(b1+Δbκ), where Δbκ(k,fNLX) is the scale-dependent bias induced by the primordial non-Gaussianity.


## 4.3.4. Mitigating Cosmic Variance and Shot-Noise

Because of the finite number of large scale wavemodes accessible to a survey, any large-scale measurement of the power spectrum is limited by the cosmic (or sampling) variance caused by the random nature of the wavemodes. For discrete tracers such as galaxies, the shot noise is another source of error. Restricting ourselves to weak primordial NG, the relative error on the power spectrumP is σP/P≈1/N(1+σn2/P), where N is the number of independent modes measured and σn2 is the shot-noise [174]. Under the standard assumption of Poisson sampling, σn2 equals the inverse of the number density 1/n¯ and causes a scale-independent enhancement of the power spectrum. The extent to which one can improve the observational limits on the nonlinear parameters fNLX will strongly depend on our ability to minimize the impact of these two sources of errors. By comparing differently biased tracers of the same surveyed volume [175, 176] and suitably weighting galaxies (by the mass of their host halo for instance) [177, 178], it should be possible to circumvent these problems and considerably improve the detection level.Figure9 illustrates how the impact of sampling variance on the measurement of fNLloc can be mitigated. Namely, the data points show the result of taking the ratio Ph(k,fNLloc)/Pδ(k,fNLloc) for each set of runs with same Gaussian random seed field ϕ before averaging over the realizations. This procedure is equivalent to the multitracers method advocated by [175]. Here, Pδ can be thought as mimicking the power spectrum of a nearly unbiased tracer of the mass density field with high number density. Although, in practical applications, using the dark matter field works better [170], in real data Pδ should be replaced by a tracer of the same surveyed volume different than the one used to compute Ph. Figure 9 also shows that upon taking out most of the cosmic variance, there is some residual noise caused by the discrete nature of the dark matter halos. As shown recently [178] however, weighting the halos according to their mass can dramatically reduce the shot noise relative to the Poisson expectation, at least when compared against the dark matter. Applying such a weighting may thus significantly improve the error on the nonlinear parameter fNLloc, but this should be explored in realistic simulations of galaxies, especially because the halo mass M may not be easily measurable from real data [170]. This approach undoubtedly deserves further attention as it has the potential to substantially improve the extraction of the primordial non-Gaussian signal from galaxy surveys.To conclude this section, it is worth noting that, while the PDF of power valuesP(k) has little discriminatory power (for large surveyed volume, it converges towards the Rayleigh distribution as a consequence of the central limit theorem) [179], the covariance of power spectrum measurements (which is sensitive to the selection function, but also to correlations among the phase of the Fourier modes) may provide quantitative limits on certain type of non-Gaussian models [174, 180].


## 4.4. Galaxy Bispectrum and Higher-Order Statistics

Higher statistics of biased tracers, such as the galaxy bispectrum, are of great interest as they are much more sensitive to the shape of the primordial 3-point function than the power spectrum [10, 42, 68, 181, 182]. Therefore, they could break some of the degeneracies affecting the non-Gaussian halo bias (e. g., the leading order scale-dependent correction to the Gaussian bias induced by the local quadratic and cubic coupling are fully degenerated [48]).


### 4.4.1. Normalized Cumulants of the Galaxy Distribution

The skewness of the galaxy count probability distribution function could provide constraints on the amount of non-Gaussianity in the initial conditions. As discussed in Section3, however, it is difficult to disentangle the primordial and gravitational causes of skewness in low redshift data unless the initial density field is strongly non-Gaussian. The first analyzes of galaxy catalogs in terms of count-in-cells densities all reached the conclusion that the skewness (and higher-order moments) of the observed galaxy count PDF is consistent with the value predicted by gravitational instability of initially Gaussian fluctuations [50, 57, 60, 183–185]. Back then however, most of the galaxy samples available were not large enough to accurately determine the SJ at large-scales [186]. Despite the two orders of magnitude increase in surveyed volume, these measurements are still sensitive to cosmic variance, that is, to the presence of massive super-clusters or large voids. Nevertheless, the best estimates of the first normalized cumulants SJ of the galaxy PDF strongly suggest that high order galaxy correlation functions indeed follow the hierarchical scaling predicted by the gravitational clustering of Gaussian ICs [187]. There is no evidence for strong non-Gaussianity in the initial density field as might by seeded by cosmic strings or textures [188].The genus statistics of constant density surfaces through the galaxy distribution measures the relative abundance of low and high density regions as a function of the smoothing scaleR and, therefore, could also be used as a diagnostic tool for primordial non-Gaussianity. For a Gaussian random field, the genus curve (i.e., the genus number as a function of the density contrast) is symmetric about δR=0 regardless the value of R. Primordial NG and nonlinear gravitational evolution can disrupt this symmetry [189]. The effect of non-Gaussian ICs on the topology of the galaxy distribution has been explored in a number of papers [35, 190–193]. For large values of R and realistic amount of primordial NG, the genus statistics can also be expanded in a series whose coefficients are the normalized cumulants SJ of the smoothed galaxy density field. In other words, the genus statistics essentially provides another measure of the (large-scale) cumulants. So far, measurements from galaxy data are broadly consistent with Gaussian initial conditions [194, 195].


### 4.4.2. Galaxy Bispectrum

Most of the scale-dependence of the primordialn-point functions is integrated out in the normalized cumulants, which makes them weakly sensitive to primordial NG. However, while the effect of non-Gaussian initial conditions, galaxy bias, gravitational instabilities and so forth, are strongly degenerated in the SJ, they imprint distinct signatures in the galaxy bispectrum Bg(k1,k2,k3), an accurate measurement of which could thus constrain the shape of the primordial 3-point function.In the original derivation of [181], the large-scale (unfiltered) galaxy bispectrum in the fNLloc model is given by(45)Bg(k1,k2,k3)=b13B0(k1,k2,k3)+b12b2[P0(k1)P0(k2)+(cyc.)]+2b13[F2(k1,k2)P0(k1)P0(k2)+(cyc.)].
Again, b1 and b2 are the first- and second-order bias parameters that describe the galaxy biasing relation assumed local and deterministic [169]. The first term in the right-hand side is the primordial contribution which, for equilateral configurations and in the fNLloc model, scales as ℳR(k,z)-1 like in the matter bispectrum, (18). The two last terms are the contribution from nonlinear bias and the tree-level correction from gravitational instabilities, respectively. They have the smallest signal in squeezed configurations.As recognized by [68, 182], (45) misses an important term that may significantly enhance the sensitivity of the galaxy bispectrum to non-Gaussian initial conditions. This contribution is sourced by the trispectrum TR(k1,k2,k3,k4) of the smoothed mass density field(46)12b12b2∫d3q(2π)3TR(k1,k2,q,k3-q)+(2perms.).
At large-scale, this simplifies to the sum of the linearly evolved primordial trispectrum T0(k1,k2,k3,k4) and a coupling between the primordial bispectrum B0(k1,k2,k3) (linear in fNLX) and the second order PT corrections (through the kernel F2(k1,k2)). In the case of local non-Gaussianity and for equilateral configurations, the first piece proportional to T0 scales as (fNLloc)2k-4 times the Gaussian tree-level prediction, with the same redshift dependence. Hence, it is similar to the second order correction (fNLloc)2ℳR-2PR(k) that appears in the halo power spectrum (see (39)). The second piece linear in fNLX generates a signal at large-scales for essentially all triangle shapes in the local model as well as in the case of equilateral NG. This second contribution is maximized in the squeezed limit (where it is one order of magnitude larger than the result obtained by [181]) which helps disentangling it from the Gaussian terms. Note that a strong dependence on triangle shape is also present in other NG scenarios such as the χ2 model [62].This newly derived contributions are claimed to lead to more than one order of magnitude improvement in certain limits [182], but it is not yet clear whether these gains can be fully realized with upcoming galaxy surveys. To accurately predict the constraints that could be achieved with future measurements of the galaxy bispectrum, a comparison of these predictions with the halo bispectrum extracted from numerical simulations is highly desirable. To date, the only numerical study [42] has measured the halo bispectrum for some isosceles triangles (k1=k2). While the shape dependence is in reasonable agreement with the theory, the observed k-dependence appears to depart from the predicted scaling.


## 4.4.1. Normalized Cumulants of the Galaxy Distribution

The skewness of the galaxy count probability distribution function could provide constraints on the amount of non-Gaussianity in the initial conditions. As discussed in Section3, however, it is difficult to disentangle the primordial and gravitational causes of skewness in low redshift data unless the initial density field is strongly non-Gaussian. The first analyzes of galaxy catalogs in terms of count-in-cells densities all reached the conclusion that the skewness (and higher-order moments) of the observed galaxy count PDF is consistent with the value predicted by gravitational instability of initially Gaussian fluctuations [50, 57, 60, 183–185]. Back then however, most of the galaxy samples available were not large enough to accurately determine the SJ at large-scales [186]. Despite the two orders of magnitude increase in surveyed volume, these measurements are still sensitive to cosmic variance, that is, to the presence of massive super-clusters or large voids. Nevertheless, the best estimates of the first normalized cumulants SJ of the galaxy PDF strongly suggest that high order galaxy correlation functions indeed follow the hierarchical scaling predicted by the gravitational clustering of Gaussian ICs [187]. There is no evidence for strong non-Gaussianity in the initial density field as might by seeded by cosmic strings or textures [188].The genus statistics of constant density surfaces through the galaxy distribution measures the relative abundance of low and high density regions as a function of the smoothing scaleR and, therefore, could also be used as a diagnostic tool for primordial non-Gaussianity. For a Gaussian random field, the genus curve (i.e., the genus number as a function of the density contrast) is symmetric about δR=0 regardless the value of R. Primordial NG and nonlinear gravitational evolution can disrupt this symmetry [189]. The effect of non-Gaussian ICs on the topology of the galaxy distribution has been explored in a number of papers [35, 190–193]. For large values of R and realistic amount of primordial NG, the genus statistics can also be expanded in a series whose coefficients are the normalized cumulants SJ of the smoothed galaxy density field. In other words, the genus statistics essentially provides another measure of the (large-scale) cumulants. So far, measurements from galaxy data are broadly consistent with Gaussian initial conditions [194, 195].


## 4.4.2. Galaxy Bispectrum

Most of the scale-dependence of the primordialn-point functions is integrated out in the normalized cumulants, which makes them weakly sensitive to primordial NG. However, while the effect of non-Gaussian initial conditions, galaxy bias, gravitational instabilities and so forth, are strongly degenerated in the SJ, they imprint distinct signatures in the galaxy bispectrum Bg(k1,k2,k3), an accurate measurement of which could thus constrain the shape of the primordial 3-point function.In the original derivation of [181], the large-scale (unfiltered) galaxy bispectrum in the fNLloc model is given by(45)Bg(k1,k2,k3)=b13B0(k1,k2,k3)+b12b2[P0(k1)P0(k2)+(cyc.)]+2b13[F2(k1,k2)P0(k1)P0(k2)+(cyc.)].
Again, b1 and b2 are the first- and second-order bias parameters that describe the galaxy biasing relation assumed local and deterministic [169]. The first term in the right-hand side is the primordial contribution which, for equilateral configurations and in the fNLloc model, scales as ℳR(k,z)-1 like in the matter bispectrum, (18). The two last terms are the contribution from nonlinear bias and the tree-level correction from gravitational instabilities, respectively. They have the smallest signal in squeezed configurations.As recognized by [68, 182], (45) misses an important term that may significantly enhance the sensitivity of the galaxy bispectrum to non-Gaussian initial conditions. This contribution is sourced by the trispectrum TR(k1,k2,k3,k4) of the smoothed mass density field(46)12b12b2∫d3q(2π)3TR(k1,k2,q,k3-q)+(2perms.).
At large-scale, this simplifies to the sum of the linearly evolved primordial trispectrum T0(k1,k2,k3,k4) and a coupling between the primordial bispectrum B0(k1,k2,k3) (linear in fNLX) and the second order PT corrections (through the kernel F2(k1,k2)). In the case of local non-Gaussianity and for equilateral configurations, the first piece proportional to T0 scales as (fNLloc)2k-4 times the Gaussian tree-level prediction, with the same redshift dependence. Hence, it is similar to the second order correction (fNLloc)2ℳR-2PR(k) that appears in the halo power spectrum (see (39)). The second piece linear in fNLX generates a signal at large-scales for essentially all triangle shapes in the local model as well as in the case of equilateral NG. This second contribution is maximized in the squeezed limit (where it is one order of magnitude larger than the result obtained by [181]) which helps disentangling it from the Gaussian terms. Note that a strong dependence on triangle shape is also present in other NG scenarios such as the χ2 model [62].This newly derived contributions are claimed to lead to more than one order of magnitude improvement in certain limits [182], but it is not yet clear whether these gains can be fully realized with upcoming galaxy surveys. To accurately predict the constraints that could be achieved with future measurements of the galaxy bispectrum, a comparison of these predictions with the halo bispectrum extracted from numerical simulations is highly desirable. To date, the only numerical study [42] has measured the halo bispectrum for some isosceles triangles (k1=k2). While the shape dependence is in reasonable agreement with the theory, the observed k-dependence appears to depart from the predicted scaling.


## 4.5. Intergalactic Medium and the Lyα Forest

Primordial non-Gaussianity also affects the intergalactic medium (IGM) as a positivefNLX enhances the formation of high-mass halos at early times and, therefore, accelerates reionization [196–198]. At lower redshift, small box hydrodynamical simulations of the Lyα forest indicate that non-Gaussian initial conditions could leave a detectable signature in the Lyα flux PDF, power spectrum and bispectrum [199]. However, while differences appear quite pronounced in the high transmissivity tail of the flux PDF (i.e., in underdense regions), the Lyα 1D flux power spectrum seems little affected. Given the small box size of these hydrodynamical simulations, it is worth exploring the effect in large N-body cosmological simulations using a semianalytic modeling of the Lyα forest [200], even though such an approach only provides a very crude approximation to the temperature-density diagram of the IGM in hydrodynamical simulations. Figure 10 shows the imprint of local type NG on the Lyα 3D flux power spectrum (which is not affected by projection effects) extracted at z=2 from a series of large simulations. The Lyα transmitted flux is calculated in the Gunn-Peterson approximation [201]. A clear signature similar to the non-Gaussian halo bias can be seen. As expected, it is of opposite sign since the Lyα forest is anti-biased relative to the mass density field (overdensities are mapped onto relatively low flux transmission).Figure 10
Ratio between thez=2 Lyα flux power spectrum extracted from simulations of Gaussian and non-Gaussian initial conditions. The mean transmission is set to F¯=0.8 and the power-law exponent α=1.65 (see text).To estimate the strength of the signal (see [200] for the details), one can assume that the (real space) optical depth τ(x) to Lyα absorption at comoving position x is approximately [202](47)τ(x)=τ¯[1+δg(x)]α,
where δg is the gas density, τ¯~1 is the optical depth at mean gas density and α~1-2 is some parameter that depends on the exact thermal history of the low density IGM. The above relation holds for the moderate overdensities δg≲10 that are responsible for most of the Lyα absorption features. To relate the gas density to the smoothed linear density field, we could make the simple ansatz δg≡δR [203]. In this linear approximation, however, the large-scale bias bF of the Lyα flux density field is much larger than that measured in detailed numerical simulations (e.g., bF2≃0.017 at z=3 [204]). Therefore, one may want to consider the lognormal mapping [205, 206](48)1+δg=exp(δR-σR22),
to better capture nonlinearities in the gas density field. Expanding exp(δR) at second order and noticing that in the presence of weak non-Gaussianity, the joint PDF P(δR(x1),δR(x2)) can generically be expanded into an Edgeworth series where the primordial 3-point function is the dominant correction, it is straightforward to compute the Lyα 3D flux power spectrum for nonzero fNLX. Upon a Fourier transformation, we arrive at (49)PF(k,fNLX)PF(k,0)=1-4gFσR2ℳR(k,z)-1ℱ(k,fNLX).
where gF is some auxiliary function of (τ¯,α,σR). This result is valid for any model of primordial NG characterized by an initial bispectrum. In the fNLloc model, the large-scale non-Gaussian Lyα bias scales as ΔbF(k,fNLX)≈-2gFσR2ℳR(k,z)-1∝k-2T(k)-1 like the non-Gaussian halo bias. Assuming that τ¯=0.7, σR=1.8 and α=1.65 yields a mean flux F¯≈0.8 and a ratio PF(k=0.01,fNLX)/PF(k=0.01,0)≈1∓0.13 for fNLloc=±100 comparable in magnitude to that seen in Figure 10. A detection of this effect, although challenging in particular because of continuum uncertainties, could be feasible with future data sets. Summarizing, the Lyα should provide interesting information on the non-Gaussian signal over a range of scale and redshift not easily accessible to galaxy and CMB observations [199, 200].


## 5. Current Limits and Prospects

As the importance of primordial non-Gaussianity relative to the non-Gaussianity induced by gravitational clustering and galaxy bias increases towards high redshift, the optimal strategy to constrain the nonlinear coupling parameter(s) with LSS is to use large-scale, high-redshift observations [32].


### 5.1. Existing Constraints on Primordial NG

The non-Gaussian halo bias presently is the only LSS method that provides a robust limit on the magnitude of a primordial 3-point function of the local shape. It is a broadband effect that can be easily measured with photometric redshifts. The authors of [156] have applied (28) to constrain the value of fNLloc using a compilation of large-scale clustering data. Their constraint arise mostly from the QSO sample at median redshift z=1.8, which covers a large comoving volume and is highly biased, b1=2.7. They obtain(50)-29<fNLloc<+69,
at 95% confidence level. These limits are competitive with those from CMB measurements, -10<fNLloc<+74 [207]. It is straightforward to translate this 2-σ limit into a constraint on the cubic order coupling gNLloc since the non-Gaussian scale-dependent bias Δbκ(k,gNLloc) has the same functional form as Δbκ(k,fNLloc) [48]. Assuming fNLloc=0, one obtains(51)-3.5×105<gNLloc<+8.2×105.
These limits are comparable with those inferred from the analysis of CMB data.Measurements of the galaxy bispectrum in several redshift catalogs have shown evidence for a configuration shape dependence in agreement with that predicted from gravitational instability, ruling outχ2 initial conditions at the 95% C.L. [208, 209]. Recent analyses of the SDSS LRGs catalogue indicate that the shape dependence of the reduced 3-point correlation Q3~ξ3/(ξ2)2 is also consistent with Gaussian ICs [210], although a primordial (hierarchical) non-Gaussian contribution in the range Q3~0.5-3 cannot be ruled out [211]. Other LSS probes of primordial non-Gaussianity, such as the abundance of massive clusters, are still too affected by systematics to furnish tight constraints on the shape and magnitude of a primordial 3-point function. Still, the observation of a handful of unexpectedly massive high-redshift clusters has been interpreted as evidence of a substantial degree of primordial NG [212–214].


### 5.2. Future Prospects

Improving the current limits will further constrain the physical mechanisms for the generation of cosmological perturbations.The non-Gaussian halo bias also leaves a signature in cross-correlation statistics of weak cosmic shear (galaxy-galaxy and galaxy-CMB) [215, 216] and in the integrated Sachs-Wolfe (ISW) effect [156, 160, 217]. Measurements of the lensing bispectrum could also constrain a number of non-Gaussian models [218]. However, galaxy clustering will undoubtedly offer the most promising LSS diagnostic of primordial non-Gaussianity. The detectability of a local primordial bispectrum has been assessed in a series of papers. It is expected that future all-sky galaxy surveys will achieve constraints of the order of ΔfNLloc~1 assuming that all systematics are reasonably under control [47, 106, 156, 160, 217, 219–221]. Realistic models of cubic type non-Gaussianity [48], modifications of the initial vacuum state or horizon-scale GR corrections [163] should also be tested with future measurement of the galaxy power spectrum.Upcoming observations of high redshift clusters will provide increased leverage on measurement of primordial non-Gaussianity with abundances and possibly put limits on any nonlinear parameterfNLX at the level of a few tens [139]. Combining the information provided by the evolution of the mass function and power spectrum of galaxy clusters should yield constraints with a precision ΔfNLloc~10 for a wide field survey covering half of the sky [214]. Alternatively, using the full covariance of cluster counts (which is sensitive to the non-Gaussian halo bias) can furnish constraints of ΔfNLloc~1-5 for a Dark Energy Survey-type experiment [222, 223].As emphasized in “Section4” however, the exact magnitude of the non-Gaussian bias is still uncertain partly due to the freedom at the definition of the halo mass and the uncertainty in the correspondence between simulated quantities and observables. Understanding this type of systematics will be crucial to set reliable constraints on a primordial non-Gaussian component. To fully exploit the potential of future galaxy surveys, it will also be essential to extend the theoretical and numerical analyses to other bispectrum shapes than the local template used so far. Ultimately, the gain that can be achieved will critically depend on our ability to minimize the impact of sampling variance and shot-noise. In this regards, multi-tracers methods combined with optimal weighting schemes should deserve further attention as they hold the promise to become the most accurate method to extract the primordial non-Gaussian signal from galaxy surveys [175–178].


## 5.1. Existing Constraints on Primordial NG

The non-Gaussian halo bias presently is the only LSS method that provides a robust limit on the magnitude of a primordial 3-point function of the local shape. It is a broadband effect that can be easily measured with photometric redshifts. The authors of [156] have applied (28) to constrain the value of fNLloc using a compilation of large-scale clustering data. Their constraint arise mostly from the QSO sample at median redshift z=1.8, which covers a large comoving volume and is highly biased, b1=2.7. They obtain(50)-29<fNLloc<+69,
at 95% confidence level. These limits are competitive with those from CMB measurements, -10<fNLloc<+74 [207]. It is straightforward to translate this 2-σ limit into a constraint on the cubic order coupling gNLloc since the non-Gaussian scale-dependent bias Δbκ(k,gNLloc) has the same functional form as Δbκ(k,fNLloc) [48]. Assuming fNLloc=0, one obtains(51)-3.5×105<gNLloc<+8.2×105.
These limits are comparable with those inferred from the analysis of CMB data.Measurements of the galaxy bispectrum in several redshift catalogs have shown evidence for a configuration shape dependence in agreement with that predicted from gravitational instability, ruling outχ2 initial conditions at the 95% C.L. [208, 209]. Recent analyses of the SDSS LRGs catalogue indicate that the shape dependence of the reduced 3-point correlation Q3~ξ3/(ξ2)2 is also consistent with Gaussian ICs [210], although a primordial (hierarchical) non-Gaussian contribution in the range Q3~0.5-3 cannot be ruled out [211]. Other LSS probes of primordial non-Gaussianity, such as the abundance of massive clusters, are still too affected by systematics to furnish tight constraints on the shape and magnitude of a primordial 3-point function. Still, the observation of a handful of unexpectedly massive high-redshift clusters has been interpreted as evidence of a substantial degree of primordial NG [212–214].


## 5.2. Future Prospects

Improving the current limits will further constrain the physical mechanisms for the generation of cosmological perturbations.The non-Gaussian halo bias also leaves a signature in cross-correlation statistics of weak cosmic shear (galaxy-galaxy and galaxy-CMB) [215, 216] and in the integrated Sachs-Wolfe (ISW) effect [156, 160, 217]. Measurements of the lensing bispectrum could also constrain a number of non-Gaussian models [218]. However, galaxy clustering will undoubtedly offer the most promising LSS diagnostic of primordial non-Gaussianity. The detectability of a local primordial bispectrum has been assessed in a series of papers. It is expected that future all-sky galaxy surveys will achieve constraints of the order of ΔfNLloc~1 assuming that all systematics are reasonably under control [47, 106, 156, 160, 217, 219–221]. Realistic models of cubic type non-Gaussianity [48], modifications of the initial vacuum state or horizon-scale GR corrections [163] should also be tested with future measurement of the galaxy power spectrum.Upcoming observations of high redshift clusters will provide increased leverage on measurement of primordial non-Gaussianity with abundances and possibly put limits on any nonlinear parameterfNLX at the level of a few tens [139]. Combining the information provided by the evolution of the mass function and power spectrum of galaxy clusters should yield constraints with a precision ΔfNLloc~10 for a wide field survey covering half of the sky [214]. Alternatively, using the full covariance of cluster counts (which is sensitive to the non-Gaussian halo bias) can furnish constraints of ΔfNLloc~1-5 for a Dark Energy Survey-type experiment [222, 223].As emphasized in “Section4” however, the exact magnitude of the non-Gaussian bias is still uncertain partly due to the freedom at the definition of the halo mass and the uncertainty in the correspondence between simulated quantities and observables. Understanding this type of systematics will be crucial to set reliable constraints on a primordial non-Gaussian component. To fully exploit the potential of future galaxy surveys, it will also be essential to extend the theoretical and numerical analyses to other bispectrum shapes than the local template used so far. Ultimately, the gain that can be achieved will critically depend on our ability to minimize the impact of sampling variance and shot-noise. In this regards, multi-tracers methods combined with optimal weighting schemes should deserve further attention as they hold the promise to become the most accurate method to extract the primordial non-Gaussian signal from galaxy surveys [175–178].



---

*Source: 908640-2010-08-04.xml*