PRIMORDIAL NON-GAUSSIANITY AND THE NRAO VLA SKY SURVEY

The investigation of primordial non-Gaussianity (NG) offers a powerful way of testing the generation mechanism of cosmological perturbations in the early universe. Although the standard single-field, slow-roll, canonical kinetic energy and adiabatic vacuum state inflation generate very small non-Gaussianity, any inflationary model that deviates from this may entail a larger level of it (Bartolo et al. 2004; Komatsu et al. 2009, and references therein).

Deviations from Gaussian initial conditions are commonly taken to be of the so-called local type and parameterized by the dimensionless parameter f_NL:

$$\begin{equation} \Phi =\phi +f_{\rm NL}(\phi ^2-\langle \phi ^2\rangle), \end{equation} \tag{ 1 }$$

where Φ denotes Bardeen's gauge-invariant potential and ϕ is a Gaussian random field. In this Letter, we use the cosmic microwave background (CMB) convention for the quoted f_NL values.

A method (Dalal et al. 2008; Matarrese & Verde 2008) for constraining non-Gaussianity from large scale structure (LSS) surveys exploits the fact that the clustering of extrema (i.e., dark matter halos where galaxies form) on large scales increases (decreases) for positive (negative) f_NL. In particular, a non-Gaussianity described by Equation (1) introduces a scale-dependent boost of the halo power spectrum proportional to 1/k² on large scales (k < 0.03 h  Mpc⁻¹), which evolves roughly as (1 + z).

Extragalactic radio sources are uniquely well suited to probe clustering on the largest scales: (1) radio surveys are unaffected by dust extinction which may introduce spurious features reflecting the inhomogeneous extinction due to Galactic dust; (2) due to their strong cosmological evolution, radio sources are very rare locally, so that radio samples are free from the profusion of local objects that dominate optically selected galaxy samples and tend to swamp very large-scale structures at cosmological distances; (3) thanks to the strong cosmological evolution, even shallow radio surveys reach out to substantial redshifts. The NRAO VLA Sky Survey (NVSS; Condon et al. 1998) offers the most extensive sky coverage (82% of the sky to a completeness limit of about 3 mJy at 1.4 GHz) with sufficient statistics to allow an accurate determination of the auto-correlation function (ACF), w(θ), on scales of up to several degrees (Blake & Wall 2002; Overzier et al. 2003).

When a realistic redshift distribution of the NVSS sources is adopted, the interpretation of the measured w(θ) in the framework of the standard hierarchical clustering scenario with Gaussian primordial perturbations requires an evolution of the bias factor radically different from that of optically selected QSOs (Negrello et al. 2006; Massardi et al. 2010), in stark contrast with the similar evolution of the luminosity function. In fact, the observed w(θ) is positive up to large (∼10°) angular scales, which, for the median source redshift (z_m ≃ 1), correspond to linear scales where the correlation function should be negative (see also Hernandez-Monteagudo 2009). Here, we explore whether the (scale-dependent) large-scale non-Gaussian halo bias could reproduce the observed shape of the NVSS source ACF, preserving the kinship with optically selected active galactic nuclei (AGNs).

We include in our analysis only NVSS sources brighter than 10 mJy, excluding the strip at |b| < 5°, where the catalog may be substantially affected by Galactic emissions. This ensures a uniform sky coverage (Blake & Wall 2002); the effect of possible residual large-scale gradients producing and offset of the ACF is discussed in Section 4. The NVSS source surface density at this threshold is 16.9 deg⁻². The redshift distribution has recently been determined by Brookes et al. (2008). Their sample, complete to a flux density of 7.2 mJy, comprises 110 sources with $S_{1.4\,\rm GHz}\ge 10\,$mJy, of which 78 (71%) have spectroscopic redshifts, 23 have redshift estimates via the K–z relation for radio sources, and nine were not detected in the K band and therefore have only a lower limit in z. Here, we have adopted the description given by de Zotti et al. (2010):

$$\begin{equation} dN/dz=1.29+32.37z-32.89z^2+11.13z^3-1.25z^4. \end{equation} \tag{ 2 }$$

The NVSS maps are pixelized using the HEALPix software package (Górski et al. 2005) with N_side = 64, corresponding to N_pix = 49,152 pixels with dimensions 092 × 092.

The ACF estimator $\hat{w}(\theta)$ reads

$$\begin{equation} \hat{w}(\theta)=\frac{1}{N_{\theta }}\sum _{i,j}\frac{(n_i-f_i\bar{n})(n_j-f_j\bar{n})}{\bar{n}^2} \;, \end{equation} \tag{ 3 }$$

where f_i and n_i are the coverage fraction and the number of radio sources in each pixel, respectively; $\bar{n}$ is the expectation value for the number of objects in the pixel (see Xia et al. 2009). The sum runs over all the pixels with a given angular separation θ. The equal weighting used here is nearly optimal because of the uniform NVSS sky coverage and because on large scales the noise is dominated by sample variance. For each angular bin centered around θ, N_θ is the number of pixel pairs separated by an angle within the bin, weighted by the coverage fractions. We used N_b = 9 angular bins in the range 1° ⩽ θ ⩽ 8° with a linear binning plus another estimate at θ = 40', since below 40' the correlation function is affected by multiple source components and above 8° the signal may be affected or even dominated by spurious density gradients (Blake & Wall 2002).

We estimated the covariance matrix of the data points using the jackknife re-sampling method (Scranton et al. 2002). We divide the data into M = 30 patches, then create M subsamples by neglecting each patch in turn, and in each subsample we measure the ACF. From the M estimates of the ACF functions, we estimate the diagonal (variance) and off-diagonal (covariance) elements of the covariance matrix for $\hat{w}(\theta)$. In Figure 1, we plot the observed NVSS ACF which is consistent with previous estimates using different approaches (e.g., Blake & Wall 2002). For comparison, we also show the best-fit theoretical ACF curve for the Gaussian case (f_NL = 0) assuming the redshift-independent minimal halo mass M_min. Clearly, the curve does not match the observed ACF data on large scales.

Zoom In Zoom Out Reset image size

Note that the integral over the full survey solid angle covered by the observationally determined ACF vanishes by construction. The estimated values go to zero for θ ≃ 30°, and the non-Gaussian model ACF shown in Figure 1 becomes negative approximately at the same θ. However, no special meaning should be attached to this coincidence since, as noted above, the observational estimate of w(θ) are unreliable for θ ≳ 8°.

3.1. Modeling the Effects of 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.

In the presence of 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 the Sheth–Tormen formula (Sheth & Tormen 1999), multiplied by a non-Gaussian correction factor (Matarrese et al. 2000; LoVerde et al. 2008):

$$\begin{eqnarray} &&\!\!\!\!\!\!R_{\rm NG}(M,z,f_{\rm NL})=1+\frac{\sigma ^2_{\rm M}}{6\delta _{\rm ec}(z)} \nonumber \\ &&\!\!\!\!\!\!\times \left[\!S_{\rm 3,M}\!\left(\!\frac{\delta ^4_{\rm ec}(z)}{\sigma ^4_{\rm M}}\!\,-\,\!2\frac{\delta ^2_{\rm ec}(z)}{\sigma ^2_{\rm M}}\!\,-\,\!1\right)\! +\!\frac{dS_{\rm 3,M}}{d\ln {\sigma _{\rm M}}}\!\left(\!\frac{\delta ^2_{\rm ec}(z)}{\sigma ^2_{\rm M}}\!\,-1\!\right)\!\right]\!, \end{eqnarray} \tag{ 4 }$$

where the normalized skewness of the density field S_3,M ∝ f_NL, and σ_M denotes the rms of the dark matter density field linearly extrapolated to z = 0 and smoothed on the scale R corresponding to a Lagrangian radius of a halo of mass M. Here, δ_ec(z) denotes the critical density for ellipsoidal collapse, which for high peaks is $\delta _{\rm ec}(z)\sim \delta _{\rm c}(z)\sqrt{q}$ (q = 0.75) and has been calibrated on N-body simulations (Grossi et al. 2009) and δ_c(z) = Δ_c(z)D(0)/D(z) where D(z) denotes the linear growth factor; Δ_c(z) ∼ 1.68 and evolves very weakly with redshift.

More importantly, the large-scale halo bias is also modified by the presence of non-Gaussianity (Dalal et al. 2008; Matarrese & Verde 2008; Grossi et al. 2009):

$$\begin{equation} b_{\rm NG}(z)-b_{\rm G}(z)\simeq 2(b_{\rm G}(z)-1)f_{\rm NL}\delta _{\rm ec}(z)\alpha _{\rm M}(k), \end{equation} \tag{ 5 }$$

where the factor α_M(k) encloses the scale and halo mass dependence. In practice, we find that, on large scales, α_M(k) ∝ 1/k² and is independent of the halo mass.

We start by assuming that the large-scale, linear halo bias for the Gaussian case is (Sheth & Tormen 1999)

$$\begin{eqnarray} b_{\rm G}&=&1+\frac{1}{D(z_{\rm o})}\left[\frac{q\delta _{\rm c}(z_{\rm f})}{\sigma ^2_{\rm M}}-\frac{1}{\delta _{\rm c}(z_{\rm f})}\right]\nonumber \\ &&+\frac{2p}{\delta _{\rm c}(z_{\rm f})D(z_{\rm o})}\left\lbrace 1+\left[\frac{q\delta ^2_{\rm c}(z_{\rm f})}{\sigma ^2_{\rm M}}\right]^p\right\rbrace ^{-1}, \end{eqnarray} \tag{ 6 }$$

where z_f is the halo formation redshift and z_o is the halo observation redshift. As we are interested in massive halos, we expect that z_f ≃ z_o. Here, q = 0.75 and p = 0.3 account for non-spherical collapse and are a fit to numerical simulations (see also Mo & White 1996; Mo et al. 1997; Scoccimarro et al. 2001). We will later relax this assumption.

Finally, the weighted effective halo bias is given by

$$\begin{equation} b_{\rm NG}^{\rm eff}(M_{\rm min},z,k,f_{\rm NL})=\frac{\int ^\infty _{M_{\rm min}}b_{\rm NG}n_{\rm NG}\,dM}{\int ^\infty _{M_{\rm min}}n_{\rm NG}\,dM}. \end{equation} \tag{ 7 }$$

Two things should be clear from Equation (5): (1) there is a degeneracy between b_G and f_NL (the same amount of non-Gaussian bias can be given by different pairs of b_G, f_NL values; strictly speaking, b_G is not a free parameter here, and the degeneracy is between M_min, which is a free parameter, and f_NL; however b_G is strongly dependent on M_min); (2) the 1/k² scale dependence means that large scales are mostly affected by f_NL and small scales are primarily affected by b_G. A positive f_NL enhances the amplitude of auto-correlation power spectra especially at large angular scales (ℓ < 200, θ>4°). This is the effect we shall use to constrain f_NL, and its impact on the ACF is clearly visible in Figure 1. In fact, in the Gaussian case, for the adopted redshift distribution and a redshift-independent M_min, the ACF is expected to become negative for θ>4°. This is also shown in Figure 2 where the ACF is plotted for several values of M_min.

Zoom In Zoom Out Reset image size

In general, the expected degeneracy between f_NL and b_G may be lifted in two ways: (1) on small scales the effect of f_NL is completely negligible but not that of b_G, (2) the redshift dependences of the two contributions are different. However, the angular correlation function encompasses the signal from different redshifts and different physical scales, complicating the separation of the f_NL and b_G signals.

To explore the effect of relaxing the assumption z_f ≃ z_o, we have also considered a model for the Gaussian bias given by (Matarrese et al. 1997; Moscardini et al. 1998, and references therein)

$$\begin{equation} b_{\rm G}(z)=b_1+\frac{b_2}{D^{\gamma }(z)} \,\,\,\,\,\,(0\le \gamma \le 2), \end{equation} \tag{ 8 }$$

with b₁ and b₂ being free parameters. Indeed, an "object-conserving" bias model corresponds to γ ≈ 1, while the bias of high-density peaks for objects that have just formed yields γ ≈ 2. In Figure 3, we show the effect of a Gaussian bias model on the auto-correlation power spectra and ACF when varying the power-law index γ; the larger the value of γ the larger is the large-scale non-Gaussian boost.

Zoom In Zoom Out Reset image size

Lin (2001) and Yoo et al. (2009) discussed the gauge dependence of matter power spectrum on very large scales (k < 0.003 h Mpc⁻¹). We find that this gauge-dependent effect on the matter power spectrum can be mimicked by that of a non-Gaussian halo bias model with f_NL ∼ 5. Here, we calculate the matter power spectrum in the conformal Newtonian gauge.

3.2. Implementation and Data Sets

We start by considering the case where the Gaussian bias is given by Equation (6). In Figure 4, we show the one-dimensional posterior probability distributions for M_min and f_NL after marginalizing over the other parameters. The right panel shows the degeneracy between M_min and f_NL. Note that the constraints on f_NL only come from the ACF data; external data sets are only used to set the underlying cosmology.

Zoom In Zoom Out Reset image size

We find that the current ACF implies f_NL>0 at the ∼3σ confidence level. The reason for that can be clearly seen in Figure 2 where the Gaussian model with the redshift-independent M_min yields an ACF dropping to zero for θ ≃ 4° and becoming negative on larger scales where the observed ACF is still positive. Non-Gaussianity of the local type adds power on large angular scales yielding a good fit to the observed data points.

The marginalized constraints on the non-Gaussianity parameter f_NL,

$$\begin{equation} f_{\rm NL}=62\pm 27 (1\sigma\, {\rm \mbox{CL}}), \end{equation} \tag{ 9 }$$

$$\begin{equation} (6)25<f_{\rm NL}<117(142)\, [95\% (99.7\%)\, {\rm \mbox{CL}}], \end{equation} \tag{ 10 }$$

are compatible with other previous estimates (Yadav & Wandelt 2008; Slosar et al. 2008; Curto et al. 2009; Smidt et al. 2009; Jimenez & Verde 2009; Smith et al. 2009; Rudjord et al. 2009) and in very good agreement with the recent WMAP7 estimate (Komatsu et al. 2010).

The minimal effective halo mass, M_min = 10^12.47±0.26 h⁻¹ M_☉ (1σ), turns out to be remarkably close to that found for optically selected QSOs: M_QSO = (3.0 ±  1.6) × 10¹² h⁻¹ M_☉ (Croom et al. 2005).

To explore whether a more general bias model (with γ allowed to conservatively vary even in an unphysical range γ < 1) than that of Equation (6) may reconcile a Gaussian model with the data, we repeated the analysis using the bias of Equation (8). We keep b₁ and b₂ fixed to 1.1 and 0.6. respectively, and we vary 0 < γ < 2 (the values for b₁ and b₂ have been chosen in order to provide a good fit to the data). In this case, the recovered central value of f_NL becomes a little higher and the 1σ error bar increases by about a factor of 2 (a smaller f_NL accommodates larger γ and larger bias).

We also perform a cross-check by using the published version of the ACF of Blake & Wall (2002; in this case zero covariance between data points is assumed, since no covariance matrix has been computed) and find f_NL = 58 ±  12(1σ  CL) (we rely on 12 data points in the range 05623 ⩽ θ ⩽ 770795). If we instead set the off-diagonal terms of our ACF estimate to zero, we get f_NL = 70 ± 15(1σ  CL). From these results, we can conclude that covariance between data points increases the error bar on f_NL by almost a factor of 2, and the Blake & Wall (2002) and our ACF measurements are in very good agreement with each other in terms of derived f_NL values. Another instructive cross-check is to see to what extent our conclusions are affected by applying an overall subtraction to all the ACF values of 10⁻⁴, in order to correct for a possible systematic offset that can contaminate the signal (Blake & Wall 2002). In this case, our constraints are weaker but consistent with the previous analysis and we get f_NL = 42 ± 30(1σ  CL). We have also checked that the correction proposed by Wands & Slosar (2009) to account for the infrared divergence of the non-Gaussian halo correlation function is negligibly small for our best-fit f_NL value.

To allow for the "integral constraint" (measurements probe the survey mean and not the ensemble mean), we have added to w(θ) a constant c and have marginalized over c, allowing this quantity to vary in the range [10⁻⁸, 10⁻⁴] (the upper limit cannot be larger since this is the theoretical variance expected on the scale of the survey scales). We find f_NL = 58 ±  28 (1σ), showing that, as expected given the large sky fraction covered by the NVSS survey, the best-fit value of f_NL is only marginally affected.

We use a Fisher matrix approach to estimate which scales contribute most to the signal for f_NL and log₁₀M_min from the ACF as a function of θ. Depending on where the signal is localized, this may give some insights into what systematic effects are important.

Figure 5 shows the Fisher-predicted error (for the ACF of a survey with NVSS characteristics) as a function of the maximum angle θ (the minimum is always set to 1°). The error-bar normalization is arbitrary and error bars of f_NL and log₁₀M_min at θ = 8° are set to be the same. The error on M_min stabilizes at θ ∼ 3°, indicating that the bias signal is mainly localized at separations smaller than 3°. The f_NL error decreases rapidly and stabilizes at larger θ, indicating that the non-Gaussian signal is localized in the ACF at θ = 2°–5°. Since the mean redshift of NVSS sources brighter than 10 mJy is about 1.2, 1° corresponds to a comoving size r ∼ 60 Mpc. The maximum non-Gaussianity signal thus comes from comoving scales in the range 120 Mpc < r < 300 Mpc, while the constraints on M_min come primarily from scales r < 180 Mpc.

Zoom In Zoom Out Reset image size

While previous analyses exploiting the NVSS to constrain primordial non-Gaussianity of local type have focused on the cross-correlation with the WMAP Internal Linear Combination map, we have shown that the angular correlation function alone is a very sensitive non-Gaussianity probe. The key point is that, given the redshift distribution of NVSS, which has been recently observationally determined, and the typical AGN halo mass estimated for optically selected QSOs, the standard ΛCDM cosmology with scale-independent bias would imply, on scales >4°, a negative ACF, but, on the contrary, it is observed to be positive. Careful analyses of the NVSS sample (e.g., Blake & Wall 2002) indicate that systematic offsets that may induce a spurious positive signal should be negligible for the sources with $S_{1.4\,\rm GHz}>10\,$mJy, selected for the present analysis. If so, the NVSS ACF may point at the presence of a small primordial non-Gaussianity that adds power to the largest scales. After marginalizing over all the other parameters, we find 25 < f_NL < 117 at the 95% confidence level, compatible with bounds derived by other studies. The minimum halo mass turns out to be M_min = 10^12.47±0.26 h⁻¹ M_☉ (1σ), remarkably close to the value found by Croom et al. (2005) for optically selected QSOs.

We have addressed the significance and the robustness of our findings by considering different bias models and by investigating the impact of gauge effects on large scales. Error bars were estimated by a jackknife re-sampling procedure, widely used in the literature (Scranton et al. 2003; Xia et al. 2009). It is known to be robust and accurate for the diagonal elements of the covariance matrix but not as widely tested and calibrated for the off-diagonal ones, which cannot be neglected because neighboring ACF data points are highly correlated. From Scranton et al. (2003), we infer that jackknife can underestimate parameter errors by up to 30%. If the error bars were to be increased by this (maximal) amount, f_NL would become compatible with zero at the ∼2σ confidence level. We conclude that our work should be seen as a "proof of principle," indicating that future surveys probing scales ∼100 Mpc at substantial redshifts can put stringent constraints on primordial non-Gaussianity (e.g., Carbone et al. 2008; Viel et al. 2009).

We thank the referee for constructive comments. We acknowledge the use of LAMBDA and the HEALPix package. Numerical analysis has been performed at the University of Cambridge High Performance Computing Service. Research is supported by the ASI Contract No. I/016/07/0 COFIS, the ASI/INAF Agreement I/072/09/0 for the Planck LFI, PRIN MIUR, MICCIN grant AYA2008-0353, FP7-IDEAS-Phys.LSS 240117, FP7-PEOPLE-2007-4-3-IRGn202182.