Constraints on Primordial Non-Gaussianity from Large Scale Structure Probes

In this paper we measure the angular power spectra $C_\ell$ of three high-redshift large-scale structure probes: the radio sources from the NRAO VLA Sky Survey (NVSS), the quasar catalogue of Sloan Digital Sky Survey Release Six (SDSS DR6 QSOs) and the MegaZ-LRG (DR7), the final SDSS II Luminous Red Galaxy (LRG) photometric redshift survey. We perform a global analysis of the constraints on the amplitude of primordial non-Gaussianity from these angular power spectra, as well as from their cross-correlation power spectra with the cosmic microwave background (CMB) temperature map. In particular, we include non-Gaussianity of the type arising from single-field slow roll, multifields, curvaton (local type), and those which effects on the halo clustering can be described by the equilateral template (related to higher-order derivative type non-Gaussianity) and by the enfolded template (related to modified initial state or higher-derivative interactions). When combining all data sets, we obtain limits of $f_{\rm NL}=48\pm20$, $f_{\rm NL}=50\pm265$ and $f_{\rm NL}=183\pm95$ at 68% confidence level for local, equilateral and enfolded templates, respectively. Furthermore, we explore the constraint on the cubic correction $g_{\rm NL}\phi^3$ on the bias of dark matter haloes and obtain a limit of $-1.2\times10^5


I. INTRODUCTION
The physical mechanisms responsible for the generation of primordial perturbations seeding present-day large-scale structure, may leave their imprint in the form of small deviations from a Gaussian distribution of the primordial perturbations. Searches for primordial non-Gaussianity can thereby provide key information on the origin and evolution of cosmological structures (e.g., ref. [1] and references therein). Although the standard single-field, slow-roll, canonical kinetic energy and adiabatic vacuum state inflation generates very small non-Gaussianity, any inflationary model that deviates from this may entail a larger level of it (refs. [2,3] and references therein).
Deviations from Gaussian initial conditions are often taken to be of the so-called local type and are parameterized by the constant dimensionless parameter f NL [4][5][6][7][8]: where Φ denotes Bardeen's gauge-invariant potential (evaluated deep in the matter era) and φ is a Gaussian random field 1 .
Recently, a method for constraining non-Gaussianity from large-scale structure surveys has been proposed [9,10], which exploits the fact that the clustering of dark matter halos -where galaxies form-on large scales is modified in a scale-dependent way by the presence of even small amount of non-Gaussianity. In particular, a non-Gaussianity described by eq. (1), introduces a scale-dependent boost (for f NL > 0 and a suppression for f NL < 0) of the halo power spectrum proportional to 1/k 2 on large scales (k < 0.03 h/Mpc), which evolves roughly as (1 + z). Large-Scale Structure (LSS) surveys covering large volumes are needed to access the scales where the signal arises (e.g., ref. [12] and references therein). The possibility of using high redshift data to constrain non-Gaussianity has also been addressed by means of hydrodynamic simulations of the Intergalactic Medium [13]. Among the many currently-available tracers of the LSS, the radio sources from the NRAO VLA Sky Survey (NVSS) [14] and the quasar catalogue of the Sloan Digital Sky Survey Release Six (SDSS DR6 QSOs) [15] are particularly interesting since they span large volumes extending out to substantial redshifts [16]. Indeed these source samples were shown to provide tight constraints on primordial non-Gaussianity by the cross-correlation measurements between CMB temperature fluctuations and the LSS number density which could be used to detect the late Integrated Sachs-Wolfe (ISW) [17] effect [18][19][20][21].
In recent papers [20,21] it was shown that the observed NVSS and SDSS DR6 QSOs Auto-Correlation Function (ACF) and Cross-Correlation Function (CCF) hint at a positive value of f NL at more than 2 σ. This is because the ACF is found to be still positive on angular scales θ > 4 • , which, for the median source redshift (z m ≃ 1), correspond to linear scales where the correlation function should be negative if the density fluctuation field is Gaussian. A positive f NL adds power on large angular scales, accounting for the observed ACF.
In this paper, we will revisit the constraints on the primordial non-Gaussianity from these two LSS tracers using the power spectrum C ℓ , instead of the correlation function w(θ). While it is true that the angular correlation function and the angular power spectrum are spherical harmonic pairs and thus in principle they should enclose the same physical information, each of the two measurements is affected by systematic errors in different ways. For example the correlation function is less affected by the survey mask than the power spectrum but error estimation for the power spectrum is much easier then for the correlation function. A direct comparison of parameter constraints obtained via the correlation function and the power spectrum analyses offer therefore a way to quantify possible systematic effects. We will also use the recently updated SDSS II LRG photometric redshift survey, MegaZ-LRG (DR7) [22] to constrain the primordial non-Gaussianity. Ref. [18] used the LRG angular power spectrum for a sample close to the SDSS DR3 to constrain f NL . Here we revisit the analysis using an improved bias model, improved LRG catalog, updated knowledge of the sources redshift distribution, up-to-date CMB maps and complementary analysis methods.
Furthermore, we explore the constraints on physically motivated primordial non-Gaussianity shapes that are different from the local case. We consider non-Gaussianities which effects on the halo clustering is well described by the equilateral and enfolded templates. We consider the angular power spectra of these three LSS tracers, as well as their cross-correlation power spectra with the CMB temperature fluctuation. Finally, we consider the cubic correction on the halo bias, which can be motivated in scenarios like the curvaton model, in which a large cubic correction can be produced while simultaneously keeping the f NL correction small [23,24].
The structure of the paper is as follows. In section II we review the effects of primordial non-Gaussianity on the power spectrum C ℓ . Section III contains the analysis of power spectrum C ℓ for NVSS radio sources, SDSS DR6 QSOs and MegaZ-LRG. In section IV we present the method used to derive constraints on f NL . Section V contains our main results. We conclude with a discussion in section VI.

II. EFFECTS OF PRIMORDIAL NON-GAUSSIANITY
The effects of non-Gaussianity on the source clustering properties arise because a non-zero f NL affects the halo mass function and enhances the halo clustering on large scales. The second effect is the dominant one on large scales.
In the presence of primordial non-Gaussianity, the mass function n NG (M, z, f NL ) can be written in terms of the Gaussian one n sim G (M, z), for which a good fit to the results of simulations is provided by e.g., the Sheth-Tormen formula [25], multiplied by a non-Gaussian correction factor [24,26,27] where the normalized skewness of the density field, S 3,M , is ∝ f NL , σ M denotes the rms of the dark matter density field linearly extrapolated to z = 0 and smoothed on scale R corresponding to the Lagrangian radius of a halo of mass M , and δ ec is the critical density for ellipsoidal collapse, calibrated on N-body simulations [30][31][32]. For high peaks (δ ec /σ M ≫ 1) and small f NL , δ ec is slightly smaller than the critical density for spherical collapse, δ c (z) = ∆ c (z)D(0)/D(z) where D(z) is the linear growth factor, and ∆ c (z) ∼ 1.68 evolves very weakly with redshift. The large-scale halo bias is also modified by the presence of non-Gaussianity [9-11, 30, 33, 34]: where the factor α M (k) encloses the scale and halo mass dependence. Here, we consider three types of non-Gaussianity. In practice, we find that, on large scales, α M (k) ∝ 1/k 2 , is independent of the halo mass for the local type. The factor α M (k) for the enfolded template is proportional to 1/k, while the equilateral template gives an almost scaleindependent α M (k). It is important to note here that the equilateral and enfolded templates were created to have factorizable expressions that gave a good description of the actual bispectra, on average, over all configurations. For the the halo bias effect, nearly squeezed configurations dominate. As discussed e.g., in Ref. [34] the templates may not necessarily offer a good approximation of the actual bispectra generated by higher-order derivatives or vacuum-state modifications. For all inflationary single-field models the bispectrum scales either as 1/k or as 1/k 3 in the squeezed limit (i.e. either like the equilateral template or like the local case) [34]. The orthogonal template instead scales like 1/k 2 in the squeezed limit. Nevertheless it is interesting to consider such an intermediate case: inflationary models that go beyond single field can show a range of scaling in the squeezed limit (e.g., ref. [35]).
We assume that the large-scale linear halo bias for the Gaussian case is [25]: where z f is the halo formation redshift, and z o is the halo observation redshift. As we are interested in massive haloes, we expect that z f ≃ z o . Here, q = 0.75 and p = 0.3 account for the non-spherical collapse and are a fit to numerical simulations (see also refs. [36][37][38]). Finally, the weighted effective halo bias is given by M min being the minimum halo mass hosting a source of the kind we are considering.

III. OBSERVED POWER SPECTRUM
A. Discrete Sources Maps In this subsection we describe three LSS probes, NVSS radio sources, high redshift SDSS DR6 QSOs and MegaZ DR7 LRGs sample. In figure 1 we show their redshift distribution, dN/dz(z). All distributions are normalized to unit integral.

NVSS Radio Sources
The NVSS [14] offers the most extensive sky coverage (82% of the sky to a completeness limit of about 3 mJy at 1.4 GHz). We start by confining our analysis to NVSS sources brighter than 10 mJy, since the surface density distribution of fainter sources suffers from declination-dependent fluctuations [39]. Density gradients in the NVSS catalog become increasingly unimportant as the source brightness threshold is increased. Figure 6 of ref. [39], figure 1 of ref. [40] and figure 3 of ref. [41] showed variations of NVSS source surface density as a function of declination for sources with different flux thresholds. Dim sources are strongly affected by density gradients, however this does not appreciably happen for the brighter sources. Furthermore, in our previous work we used the NVSS sources with different flux thresholds S ≥ 10 mJy and S ≥ 20 mJy to compute the ACF of NVSS source and its CCF with CMB temperature fluctuations and found that the obtained ACF and CCF results are stable [21]. Also we mask the stripe |b| < 5 • , where the catalog may be substantially affected by Galactic emissions. Here, we have explicitly checked that the choices of masked stripes with |b| < 10 • and |b| < 20 • have negligible impact on our final results. In order not to excessively reduce the sample, in our analysis we consider the NVSS sources brighter than 10 mJy and the masked strip with |b| < 5 • to reduce the contamination from Galactic emissions. The NVSS source surface density at this threshold is 16.9 deg −2 and the redshift distribution at this flux limit has been recently determined by ref. [42]. Their sample, complete to a flux density of 7.2 mJy, comprises 110 sources with S 1.4GHz ≥ 10 mJy, of which 78 (71%) have spectroscopic redshifts, 23 have redshift estimates via the K-z relation for radio sources, and 9 were not detected in the K band and therefore have a lower limit to z. We adopt here the smooth description of this redshift distribution given by ref. [43], which is shown in the red dashed line of figure 1: The mean redshift of this sample isz = 1.23 and the redshift range is 0 < z < 3.5.

SDSS DR6 Quasars
The SDSS DR6 QSO catalog released by Ref. [15] contains about 10 6 objects with photometric redshifts ranging from 0.065 to 6.075 over a total area of 8417 deg 2 (∼ 20% of the whole sky). We refer the reader to ref. [15] for a detailed description of the object selection with the non-parametric Bayesian classification kernel density estimator (NBC-KDE) algorithm. We use the electronically-published table that contains only objects with the "good" flag with values in the range [0, 6]. The higher the value, the more probable for the object to be a real QSO (see § 4.2 of ref. [15] for details). Furthermore we restrict ourselves to the "uvxts=1", i.e. to QSOs clearly showing a UV excess which should be a signature of a QSO spectrum (in this case we have N qso ≈ 6 × 10 5 QSOs). As we know, several systematics effects, including the galactic extinction by dust, sky brightness, number of point sources and poor seeing, could potentially affect both the observed auto-/cross-correlation power spectra. In our previous work we carefully checked for their contribution in calculations and in particular we considered extinction and point sources contamination. We found that the dominant systematic effect is the extinction [44,46]. Therefore, in order to minimize the effect of Galactic extinction on the observed QSO distribution, we use the extinction mask A g < 0.18 only in our following analysis [44][45][46]. These masks will remove about ∼20% of the considered area.
We fit the redshift distribution dN/dz of the DR6 QSO sample with a function of the form [44] (black solid line of figure 1): The best-fit values of the parameters are m = 2.00, β = 2.20, z 0 = 1.62; the mean redshift of the sample isz ∼ 1.49.

MegaZ DR7 LRG
We use the updated MegaZ LRG DR7 sample [22], which contains ∼ 1.5 × 10 6 galaxies from the SDSS DR7 in the redshift range 0.4 < z < 0.7 with limiting magnitude i < 20. To reduce stellar contamination, there is a variable of the MegaZ neural network estimator δ sg , defined such that δ sg = 1 if the object is a galaxy, and δ sg = 0 if it is a star. For a conservative analysis, we choose a cut δ sg > 0.2, which is reported to reduce stellar contamination below 2% while keeping 99.9% of the galaxies. In addition to the SDSS DR7 geometry mask, we also add two foreground masks to account for seeing (removing pixels with median seeing in the red band larger than 1.4 arcsec) and reddening (removing pixels with median extinction in the red band A r > 0.18) [47]. The redshift distribution function (blue dash-dotted line of figure 1) in this case is found directly from the photometric redshifts that are given in the catalogue [48], which is similar to a Gaussian function: where the mean redshift z 0 = 0.55 and the deviation σ = 0.055.

The pseudo-C ℓ estimator
A two dimensional density field σ(n) defined over the full sky can be decomposed in a series of spherical harmonics Y ℓm and their corresponding coefficients a ℓm : with On the full sky spherical harmonics are a complete and orthogonal basis set. However, in practice, a masked (incomplete) sky is observed. Thus the observedã ℓm coefficients on the partially cut sky become [49]: where W (n) is the position dependent, weight function imposed by the mask and, optionally, the adopted weighting scheme. From these coefficients we can obtain the observed pseudo power spectrumC ℓ , the pseudo-C ℓ estimator [49]: The pseudo power spectrumC ℓ , given by the weighted spherical harmonic transform of a masked map, is clearly different from the full sky true power spectrum C ℓ . The expectation value ofC ℓ is related to the true C ℓ by a convolution: where the coupling matrix G ℓℓ ′ , describing the mode coupling resulting from the weight function W (n) [50], can be expressed in terms of 3j symbols as: W ℓ is the angular power spectrum of the weight function: where the weight coefficientsw ℓm are:w For the LSS catalogues presented above, in practice, we construct the pixelized maps using the HEALPix software package [51] with high resolution N side = 512, yielding pixel areas of 6.87 ′ × 6.87 ′ and estimate their pseudo power spectrumĈ ℓ using the HEALPix function "anafast". We also subtract the shot noise contribution ∆Ω/N from the observed pseudo power spectrum, where ∆Ω is the surveyed area and N is the observed number of sources. The pseudo power spectrum is measured up to ℓ max = 400, at these multipoles the magnitude of the estimated error bars is dominated by shot-noise. In order to avoid the effects of gauge corrections on the power spectrum on very large scales [52,53], we set the minimal multipole ℓ min = 10. The effects of non-Gaussianity are large on very large scales. By ignoring ℓ < 10 we are neglecting the scales where the signal is maximal. However the theoretical predictions for the measured correlations is uncertain on very large scales as discussed in ref. [53]. In ref. [21] we estimated that these uncertainties could introduce a systematic error on f NL of 5. Here we take a more conservative approach and simply ignore the largest scales, even if this implies a reduction of the signal-to-noise.

Covariance matrix of the pseudo-C ℓ estimator
In the signal dominated limit, the inverse of the power spectrum covariance matrix F ℓℓ ′ is [54]: where C ℓ is the theoretical angular power spectrum, N ℓ = ∆Ω/N denotes the shot noise term. C ℓ and N ℓ are the necessary error contributions from both cosmic variance and shot noise, respectively. For the full sky survey, G ℓ1ℓ2 = δ ℓ1ℓ2 f sky and the covariance matrix becomes diagonal:

Cross-correlation power spectrum
Besides the angular power spectrum of each LSS tracer, we also consider the cross-correlation power spectrum among them. In figure 1, we show that the redshift distributions of these three tracers overlap. We can thus expect some cross-correlation signal among them. We also use the HEALPix software to estimate the cross-correlation power spectrumC XY ℓ between tracer X and Y :C whereã i ℓm are the observed coefficients of tracer i and they are computed on the common sky area only. Consequently, the inverse of the covariance matrix becomes: where the shot noise term N i ℓ = ∆Ω/N i is scale independent, ∆Ω is the common area of two surveys, N i is the observed number of sources in this common area.
Similarly, we also estimate the cross-correlation power spectrumC XT ℓ between the LSS number density map and the CMB map. While the Planck satellite is operating [55], here, we use the the 7-years Internal Linear Combination (ILC) map of the CMB provided by the Wilkinson Microwave Anisotropy Probe (WMAP) team [56]. Since an estimate of the uncertainty associated to this map is not provided, we have checked that using other CMB maps (e.g. the 5 years ILC map by the WMAP team) does not change the results in any significant way. For the cross-check, we also use the V and W band all-sky 7-years CMB maps provided by WMAP [57] and smooth these maps to the resolution of ILC map (60 arcmins). Using these maps, we perform all the calculations again and find that our results are undistinguishable. We adopt the WMAP KQ75 mask, excluding about 30% of the sky at low Galactic latitude, to avoid most of the residual Galactic contamination.
We still use the estimator of cross-correlation power spectrumC XT ℓ : whereã T ℓm are the observed coefficients of the CMB map ∆T (n), also computed on the common sky area only. The inverse of the covariance matrix is: where C TT ℓ is the theoretical angular power spectrum of CMB map, N TT ℓ is measurement error of the CMB temperature power spectrum, in which the main source is the pixel noise which accounts for the imperfections in measured temperature on the sky. Pixel noise for WMAP is reported as σ pix = σ 0 / √ N obs , where σ 0 is the noise per observation and N obs is the number of observations per given pixel [58]. After getting the theoretical C T ℓ and C XT ℓ , we include the correction of Gaussian beam B(ℓ) (using the HEALpix gaussbeam routine): . We neglect the pixel window function since its contribution at N side = 512 is negligible at the resolution of the WMAP ILC map which is 60 arcmins FWHM.

IV. METHOD & DATA ANALYSIS
We use the publicly available package CAMB − sources 3 [59] to calculate the theoretical power spectrum. In our analysis, we use the Halofit [60] built-in routine for non-linear correction to obtain the fully-evolved, nonlinear matter power spectrum P (k) at any epoch. We perform a global fitting of cosmological parameters, including f NL , for the data of section III including also the datasets described below, using the CosmoMC package [61], a Markov Chain Monte Carlo (MCMC) code, modified to calculate the theoretical ACF and CCF. We assume purely adiabatic initial conditions and a flat Universe, with no tensor contribution to primordial fluctuations. The following six cosmological parameters are allowed to vary with top-hat priors: the dark matter energy density parameter Ω c h 2 ∈ [0.01, 0.99], the baryon energy density parameter Ω b h 2 ∈ [0.005, 0.1], the primordial spectral index n s ∈ [0.5, 1.5], the primordial amplitude log[10 10 A s ] ∈ [2.7, 4.0], the ratio (multiplied by 100) of the sound horizon at decoupling to the angular diameter distance to the last scattering surface Θ s ∈ [0.5, 10], and the optical depth to reionization τ ∈ [0.01, 0.8]. The pivot scale is set at k s0 = 0.05 Mpc −1 and do not consider massive neutrinos and dynamical dark energy. Besides these six basic cosmological parameters, we have two more parameters related to the power spectrum C ℓ data: the non-Gaussianity parameter, the minimal halo mass M min . When we combine all the different probes together we have four parameters: three minimal halo masses for three LSS surveys and the non-Gaussianity parameter.
The model power spectrum C th ℓ is compared with the observed valuesC ℓ , respectively, through the Gaussian likelihood function [62]: where the Fisher matrix F ℓℓ ′ is the inverse of the power spectrum covariance matrix. While it is true that at very low ℓ a Gaussian likelihood is not a good approximation for CMB data, at high ℓ the Gaussian approximation works well even for CMB, for LSS it is customary to adopt a Gaussian likelihood.
The following cosmological data are also included in the fit: i) power spectra of CMB temperature and polarization anisotropies; ii) baryonic acoustic oscillations (BAOs) in the galaxy power spectra; iii) SNIa distance moduli.
To deal with the 7-years WMAP (WMAP7) CMB temperature and polarization power spectra we use the routines for computing the likelihood supplied by the WMAP team [63]. The WMAP7 data are used only to improve the constraints on the six basic cosmological parameters, not to constrain f NL .
The BAOs [64] can, in principle, measure not only the angular diameter distance, D A (z), but also the expansion rate of the Universe, H(z). However, the limited accuracy of current data only allows us to determine the ratio between the distance scale defined by ref. [65]: and the comoving sound horizon at the baryon-drag epoch, r s (z d ) (see ref. [66]). Accurate determinations of the distance ratio r s (z d )/D v (z), r s (z d ) have been obtained by Ref.
We adopt these values as Gaussian priors. The SNIa distance moduli provide the luminosity distance as a function of redshift D L (z) which provides strong constraints on the dark energy evolution. In this paper we will use the latest SNIa data sets from the Supernova Cosmology Project, "Union2 Compilation" which consists of 557 samples and spans the redshift range 0 z 1.55 [67]. In this data set, they improved the data analysis method by using and refining the approach of their previous work [68]. When comparing with the previous "Union Compilation", they extended the sample with the supernovae from refs. [67,[69][70][71]. The authors also provide the covariance matrix of data with and without systematic errors and, in order to be conservative, we include systematic errors in our calculations. In the calculation of the likelihood from SNIa we marginalize over the nuisance parameter as done in refs. [72,73].
Furthermore, we add a prior on the Hubble constant, H 0 = 74.2 ± 3.6 km/s/Mpc given by ref. [74]. Finally, in the analyses we set the minimal halo mass at M min > 10 12 h −1 M ⊙ consistent with the results of ref. [75].

V. NUMERICAL RESULTS
In this section, we present constraints on the primordial non-Gaussianity parameter f NL from the observed LSS power spectra, in combination with the external data sets introduced above: the WMAP7, BAO and SNIa. These extra datasets are only used to constrain the underlying cosmology. In table I and figure 4 we show the constraints on the primordial non-Gaussianity parameter from different data combinations, after marginalizing over all the other parameters. We also plot the observed power spectra data ℓC ℓ in figure 2 and figure 3, together with the corresponding theoretical best fit models for comparison. For illustration purposes, here we bin the observed power spectra with the bin size ∆ℓ = 10.
A. Local Type

Angular Power Spectrum
We begin by considering the constraint on the local type of non-Gaussianity f NL from the angular power spectrum only. In the first panel of figure 2, we plot the observed NVSS angular power spectrum data, which is consistent with the previous work [76]. When using the NVSSC ℓ data only, we obtain the marginalized constraint on the non-Gaussianity parameter: f NL = 78 ± 52 (1 σ) and −34 < f NL < 187 (2 σ), which is consistent with zero. The minimal halo mass, M min = 10 12.44±0.26 h −1 M ⊙ (1 σ), turns out to be remarkably close to our previous work [20]. We show the best fit model in this panel (black dashed line). The amplitude of the power spectrum on large scales (small ℓ) is enhanced to fit the data points, by the positive value of f NL . In our previous works [20,21], we used the autocorrelation function (ACF) of NVSS radio sources to constrain f NL and found that the NVSS ACF yields a positive f NL at more than 2 σ. The reason for this difference is that we only use the power spectrum data at ℓ ≥ 10, and the power spectrum approach does not include the information on the total number of NVSS radio sources N tot (which was folded in the ACF analysis). If we add a prior on the total number of sources N tot : 0.5 < N tot /N NVSS tot < 1.5, we obtain a limit of f NL = 62 ± 30 at 1 σ. This result is perfectly consistent with previous works [20,21].
Next, we use the estimated angular power spectrum of the SDSS DR6 QSO sample, which is shown in the second panel of figure 2. There is an obvious excess power at large scales, which favors a positive f NL : f NL = 62 ± 26 (1 σ) and 5 < f NL < 115 (2 σ). Similarly to our previous work [21], we also find evidence for positive f NL at more than 2 σ. However, ref. [18], using SDSS QSO data, obtained non-Gaussianity constraints consistent with zero at 95% confidence level. As we discussed before, ref. [18] used an extension of the DR3 QSO sample that include sources that subsequently were released with the DR6 sample: we use the more complete and better calibrated final official SDSS DR6 QSO catalog release [15]. Ref. [18] found that their QSO sample at z < 1.45 (QSO0) seem to suffer from contamination which they attribute to systematic calibration errors, discarded this sample from their analysis and only used the high redshift sample (QSO1). We have checked for this effect and find no evidence in the sample we use that the z < 1.45 QSOs have different contamination than the z > 1.45 [21]. In addition, we remove three narrow stripes in the southern hemisphere of the galaxy and set the minimal multipole ℓ min = 10, which could reduce the effect of star contamination further.
Finally, we consider the MegaZ DR7 LRG sample. Ref. [22] split this sample in four redshift bins according to each source photometric redshift estimate and found a large excess of power over the lowest multipoles (ℓ < 20) in the angular power spectrum, which is also found in their previous analysis using MegaZ DR4 LRG sample [77]. In our analysis, although we do not split the LRG sample in redshift bins, we still find this excess power on large scales, which is shown in the third panel of figure 2. There are several explanations for this excess power that range from systematic errors as due to contaminants to new physics [78]. Ref. [78] checked some possible systematic errors of this MegaZ LRG sample in details and found that systematics errors do not seem to be responsible for the excess power at large scales. We also perform a similar check and find no significant changes in the power spectrum at large scales. One of the possibilities is that this excess power could be induced by the primordial non-Gaussianity. Here, we use the observedC ℓ data of MegaZ DR7 LRG sample to constrain the primordial non-Gaussianity for the first time. We calculate the angular power spectrum on largest scales and account for the redshift distortion power as described in Ref. [79], although we set the minimal multipole ℓ min = 10. After marginalizing over other free parameters, we obtain the constraint: which is consistent with zero at 2 σ. We also obtain a limit of M min = 10 12.45±0.07 h −1 M ⊙ at 68% confidence level. Using eq. (3), we could obtain the limit of effective halo bias at the mean redshift of the survey: b(z = 0.55) = 1.93 ± 0.06 (1 σ), which is consistent with other works [80], although in our analysis we directly use the full sample. When we combine these three observed angular power spectra data together to constrain non-Gaussianity, we obtain a limit of f NL = 68 ± 22 at 68% confidence level and 22 < f NL < 108 at 95% confidence level, which implies the current angular power spectra of LSS tracers favor f NL > 0 at the ∼ 3 σ level, since non-Gaussianity adds power on large angular scales yielding a good fit to the observed data points.
As can be seen from figure 2 and table I, there is not a specific measurement that is driving the signal: all the different probes yield very consistent results and when combined the error bar on f NL is reduced to yield a > 2σ result.

Adding Cross-Correlation Power Spectrum
Comparing with the angular power spectrum, the constraining power of cross-correlation power spectrum on f NL is much weaker [12]. In the panels of figure 3, we plot the observed cross-correlation power spectraC XY ℓ among three LSS tracers andC XT ℓ between the CMB map and three LSS tracers. These cross-correlation data points are consistent with zero within 1 σ error bar but the best fit value f NL = 48 is an excellent fit to the data points, as shown in figure  3. However, cross-correlation data give a larger error bar for f NL than that from the angular power spectra data and are thus much less constraining.
First, we consider the single tracer case in combination with its cross-correlation with the CMB map. When we combine the angular power spectrum and its cross-correlation power spectrumC XT ℓ with CMB map to constrain the primordial non-Gaussianity, we obtain f NL = 74 ± 40, f NL = 59 ± 21 and f NL = 153 ± 95 at 68% confidence level for NVSS radio sources, SDSS DR6 QSOs and MegaZ DR7 LRGs, respectively. The results are consistent with those obtained from the angular power spectrum only, but those cross correlations does not seem to add much information.
If we consider all the angular and cross-correlation power spectra data sets of three LSS tracers together 4 , we obtain the constraint on the local type primordial non-Gaussianity: This result is compatible with previous estimates [18,33,[81][82][83][84][85][86], and with the WMAP7 limits [63]. The LSS data still implies a positive f NL at more than 2 σ. In figure 2 and figure 3, we plot the theoretical best fit model f NL = 48 of our non-Gaussianity calculations (red solid lines). These curves match the observed power spectra very well, especially on large scales.

B. Other Shapes
Generally, the correction to the standard halo bias due to the presence of primordial non-Gaussianity is ∆b h /b h = δ c β(k), where the function β(k) is given by [87]: where M R = 4πρ m R 3 /3 is the halo mass which is related to the scale R, B Φ (k) denotes the expression for the primordial bispectrum of the Bardeen potential Φ, P Φ its power spectrum and α ≡ k 2 1 + k 2 + 2k 1 kµ. In the local non-Gaussian case, the expression for the correction to the standard gaussian halo bias simplify [eq.
(3)]. In this subsection, we explore other two different types of primordial non-Gaussianity given by the Equilateral and Enfolded templates.
The equilateral type of non-Gaussianity [88,89], which can be generated from the inflationary models with higher-derivative operators of the inflaton, can be well described by the following template [90]: Then we find that the correction β(k) in the Equilateral non-Gaussian case is almost scale-independent (see figure 1 of ref. [87]). In principle there is a mass dependence (see discussion in ref. [34]), but the template and the actual physical bispectrum yield different effective f NL (denoted byf NL ) for the halo bias. Here, therefore, we consider a phenomenological scale-independent scaling and neglect the mass dependence. The interpretation of the resulting constraint should be interpreted as a constraint on an effective f NL , its connection to the physical bispectrum depends on the detailed model of inflation under consideration. Using all datasets, we obtain the constraint on the Equilateral type of primordial non-Gaussianity: which is consistent with the WMAP7 results [63]. The template of another type of non-Gaussianity "Enfolded", which arises in models with non-Bunch-Davies initial state [89,91] or in effective field theories of inflation with higher-derivative interactions [92,93] (see also ref. [94]), is given by ref. [95]: In figure 1 of ref. [87], it can be seen that the correction in this type is proportion to k −1 on large scales.
No single field inflationary model has this scaling in the squeezed limit: all inflationary single field models considered so far have only the local or the template equilateral scaling in this limit [34]. This is nevertheless an interesting intermediate case to consider. In fact non-single field models can yield all the intermediate scalings in the squeezed limit, see e.g., ref. [35]. We therefore consider a phenomenological scaling like k −1 which could cover this Enfolded template case and neglect the mass dependence. In this case, the whole data sets give the constraint on the effective non-Gaussianity parameter for the enfolded template: In the WMAP7 paper [63], a constraint on the Orthogonal type of non-Gaussianity [89] is given. The correction to the halo bias of this template is also proportional to k −1 on large scales, but have the opposite sign and a normalization factor off orth NL = −2f enf NL with a weak mass dependence which we neglect here. Therefore, based on the result of Enfolded template, in our simple analysis, we obtain the constraint on the non-Gaussianity parameter for the orthogonal template:f orth NL = −92 ± 47 (1 σ) and −179 <f orth NL < 6 (2 σ). The orthogonal type of non-Gaussianity in the squeezed limit however scales like the equilateral one and, as shown in ref. [34] its effect on the halo bias is almost degenerate with the equilateral non-Gaussianity. The enfolded type of non-Gaussianity arising for example from modified initial state, in the squeezed limit scales like the local template and as shown in ref. [34] its effect on the halo bias is almost degenerate with the local non-Gaussianity but with a different normalization f mod.in.state NL ≃ 1/8f local NL . This re-scaling can be used to re-interpret the local constraints into constraints for the modified initial state non-Gaussianity.

C. Cubic Correction gNL
Finally, we consider the effect of the cubic correction on the halo bias model. In this case, the Bardeen potential Φ can be conveniently parameterized up to third order by: where g NL is dimensionless, phenomenological parameter. Ref. [96] explored the effect of a cubic correction g NL φ 3 on the mass function and halo bias model in detail and gave the expression of the scale-dependent bias correction when f NL = 0 and g NL = 0: where the normalized skewness of the density field S 3,M = S 3,M /f NL and the parameter ǫ k ≃ 0.6 obtained by the fitting in ref. [96]. Then we obtain the constraint on g NL from all the data sets: The result is compatible with that obtained in refs. [96,97].

VI. SUMMARY & CONCLUSIONS
Searches for primordial non-Gaussianity can provide key information on the physical mechanism that generate the primordial perturbations and can constrain specific inflationary models. A relatively new and promising method to detect primordial non-Gaussianity exploits the fact that the clustering of dark matter halos (which host massive galaxies and QSOs) is modified on large-scales in a scale-dependent way. This effect is called non-Gaussian halo bias. To access these large scales, large-scale structure surveys covering large volumes, high redshift and tracers of massive halos are most suited. Among the currently available surveys, the NVSS and the SDSS QSO (in particular the DR6  catalog) are the largest. Previous work had detected excess large-scale power in their auto-correlation function yielding a detection of a positive value of the non-Gaussianity parameter f NL at > 2σ level. Here we have revisited the issue using the angular power spectrum rather than the correlation function. While the two statistics should enclose the same physical information they are affected by systematics effects in different ways: this complementary analysis offers thus a consistency check and a test for possible systematics. We have also used the recently updated MegaZ-LRG survey: a previous release of the same survey had also been found to have excess large-scale power, but mostly at ℓ < 10. In an effort to be as conservative as possible in our analysis, we excuded all multipoles at ℓ < 10: these multipoles are affected by theoretical uncertainties and systematics such as calibration, sample contamination etc. These also probe the largest scales, thus these are the multipoles where most of the non-Gaussian signal is likely to be concentrated. We combine the survey auto power spectra with the cross-power spectra among the surveys and of each of the survey with the CMB temperature map. In our analysis we marginalize over the LCDM cosmological parameters with a prior imposed by WMAP7 angular power spectra, BAO, supernovae and Hubble constant measurement. We have investigated the constraints on the parameter f NL , characterizing primordial non-Gaussianity not only for the popular local shape but also for the equilateral, enfolded and orthogonal templates. We have also considered constraints on a cubic non-Gaussianity correction (parameterized by the g NL parameter). We find no evidence for non-Gaussianity for any "shape" which halo bias effect on large scale scales less steeply than 1/k 2 and no evidence for non-zero g NL . However for the local type of non-Gaussianity we find a ∼ 2σ signal for positive f NL when all data-sets are combined.
The results are summarized in figure 4: no single data set drives the signal and the central f NL values recovered are fully in agreement with previous analyses. Due to the more conservative approach taken here (and thus larger error-bars), the only single survey that yield more than 2σ detection is the SDSS DR7 QSO survey. The NVSS survey central recovered f NL value is also in agreement with both the QSO measurement and previous analyses but the error-bar is larger. The combination of the data sets yields: f NL = 48 ± 20 and 5 < f NL < 84 at 1 and 2 σ. Thus the tantalizing hint for a positive local f NL found from the auto-correlation function remains despite the different analysis. We have investigated several possible sources of systematic errors such as the sources number density and the contamination by Galactic emissions for NVSS sources; different choices for the CMB temperature fluctuations templates; possible contamination of stars for the SDSS and LRG samples. Our main results are stable when these systematic errors are considered. Such a result would have profound implications for inflationary mechanisms. The improved statistical power of forthcoming Large-Scale Structure surveys is needed to fully resolve the issue. Moreover, it will be also important to use the CMB data which will be provided by Planck.