ournal of C osmology and A stroparticle P hysics J Constraining primordial non-Gaussianity with high-redshift probes

. We present an analysis of the constraints on the amplitude of primordial non-Gaussianity of local type described by the dimensionless parameter f NL . These constraints are set by the auto-correlation functions (ACFs) of two large scale structure probes, the radio sources from NRAO VLA Sky Survey (NVSS) and the QSO catalogue of Sloan Digital Sky Survey Release Six (SDSS DR6 QSOs), as well as by their cross-correlation functions (CCFs) with the cosmic microwave background (CMB) temperature map (Integrated Sachs-Wolfe eﬀect). Several systematic eﬀects that may aﬀect the observational estimates of the ACFs and of the CCFs are investigated and conservatively accounted for. Our approach exploits the large-scale scale-dependence of the non-Gaussian halo bias. The derived constraints on f NL coming from the NVSS CCF and from the QSO ACF and CCF are weaker than those previously obtained from the NVSS ACF, but still consistent with them. Finally, we obtain the constraints on f NL = 53 ± 25 (1 σ ) and f NL = 58 ± 24 (1 σ ) from NVSS data and SDSS DR6 QSO data, respectively.


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 and φ is a Gaussian random field. In the literature there are two conventions: in the large scale structure (LSS) convention Φ is linearly extrapolated to z = 0, while in the cosmic microwave background (CMB) convention it is evaluated deep in the matter era. Thus, f LSS NL = [g(z = ∞)/g(z = 0)]f CMB NL ∼ 1.3f CMB NL , where g(z) denotes the Λ-induced linear growth suppression factor. In this paper we will use the CMB convention.
A new method [9,10] for constraining non-Gaussianity from large-scale structure surveys exploits the fact that the clustering of dark matter halos -where galaxies form-on large scales increases (decreases) for positive (negative) f NL . In particular, a non-Gaussianity described by eq. (1.1), introduces a scale-dependent boost of the halo power spectrum proportional JCAP08(2010)013 to 1/k 2 on large scales (k < 0.03 h/Mpc), which evolves roughly as (1 + z). LSS surveys covering large volumes are needed to access the scales where the signal arises (e.g., ref. [11] and references therein). Among the many currently-available tracers of the LSS, the radio sources from NRAO VLA Sky Survey (NVSS) [12] and the QSO catalogue of Sloan Digital Sky Survey Release Six (SDSS DR6 QSOs) [13] are particularly interesting since they span large volumes extending out to substantial redshifts [14]. Indeed these source samples were shown to provide tight constraints on primordial non-Gaussianity [15][16][17].
Extragalactic radio sources are uniquely well suited to probe clustering on the largest scales. Radio surveys are in fact unaffected by dust extinction which may introduce in the observed large-scale distribution spurious features reflecting the inhomogeneous extinction due to Galactic dust. Moreover, 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; thanks to the strong cosmological evolution, even relatively shallow radio surveys reach out to substantial redshifts. The NVSS [12] 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 angular correlation function, w(θ) on scales of up to several degrees [18,19].
In a recent paper [17] it was shown that the observed NVSS Auto-Correlation Function (ACF) hints at a positive value of f NL at about the 3 σ confidence level. 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.
A cross-check for a positive f NL can be provided by an enhancement, compared to the Gaussian case, of the Cross-Correlation Function (CCF) of the CMB with LSS probes (late-time Integrated Sachs-Wolfe (ISW) [20] effect). Earlier analyses [15] did not find any evidence for f NL > 0 from the NVSS-CMB CCF, although the NVSS-CMB cross-spectrum showed an anomalous point, suggestive of some unaccounted systematic effect in the spatial distribution of NVSS sources and/or in the adopted CMB map. In this paper we will revisit this issue investigating possible systematics that may affect the results.
We will also revisit the ACF of SDSS DR6 QSOs [13] and their CCF. This was previously analyzed by ref. [21] to exploit the high redshift regime for early dark energy models and for the ISW effect. Here we use up-to-date CMB maps to derive constraints on primordial non-Gaussianity, following the approach of ref. [15]. Optically selected QSOs are also well suited to test primordial non-Gaussianity, as they probe large-volumes and high redshifts and are not seriously affected by dust; their major contaminating source being stars from our own galaxy. Ref. [15] used an extension of the SDSS DR3 QSOs sample to constrain f NL . Here we revisit the analysis using an improved bias model, improved QSO catalog, updated knowledge of the sources redshift distribution, up-to-date CMB maps and complementary analysis methods.
The structure of the paper is as follows. In section 2 we review the effects of primordial non-Gaussianity on the ACF and on the CCF, and the theory of the late-time ISW effect. Section 3 contains the analysis of ACF and CCF for NVSS radio sources and SDSS DR6 QSOs. In section 4 we present the method used to derive constraints on f NL . Section 5 contains our main results. We conclude with a discussion and comparison with related work in section 6.

JCAP08(2010)013 2 Theoretical framework
In this section we briefly present the equations describing the ISW effect, and review the impact of primordial non-Gaussianity on cosmic observables relevant to our analysis, namely the halo mass function and the halo bias.

Integrated Sachs-Wolfe effect
The temperature anisotropy due to the ISW effect can be expressed as an integral of the time derivative of the gravitational potential Φ over conformal time η A CMB photon falling into a gravitational potential well gains energy, while loses energy when it climbs out of it. These effects exactly cancel out if the potential Φ is time independent, such as in the matter dominated era, when the fractional density contrast δ m is proportional to the scale factor a, so thatΦ = 0 and no ISW is produced. However, when dark energy (or curvature) becomes important at later times, the potential evolves as the photon passes through it andΦ = 0, producing additional CMB anisotropies. The late-time ISW effect can thereby be a powerful tool for probing dark energy and its evolution. However, the main contribution of the ISW effect to CMB anisotropies occurs on large scales that are strongly affected by cosmic variance. This problem can be overcome by cross-correlating the CMB temperature fluctuations with the distribution of extragalactic sources.
The number density contrast in a given directionn 1 is: is the scale-independent bias factor relating the density contrast of visible objects, δ g , to the mass density contrast δ m (δ g = b g δ m ), and dN/dz is the normalized redshift distribution of the survey. The ISW temperature anisotropy in a directionn 2 is: The ACF of a complete sample of sources, C gg (θ), and its CCF with a CMB map, C gT (θ), as a function of the angular separation θ betweenn 1 andn 2 , can be written as: where the P ℓ [cos(θ)] are the Legendre polynomials, θ s is the smoothing scale of the CMB map.
Here C gg ℓ and C gT ℓ are the auto-correlation and cross-correlation power spectra, respectively, JCAP08(2010)013 given by: in terms of the present day matter power spectrum, P (k), and of the functions I g ℓ (k) and I ISW ℓ (k):

9)
j ℓ (x) being the spherical Bessel functions, and χ the comoving distance. We use the publicly available package CAMB − sources 1 [22] to calculate the theoretical ACF and CCF. JCAP08(2010)013

Effects of primordial non-Gaussianity
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 [23], multiplied by a non-Gaussian correction factor [24][25][26] 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 [29]. 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,10,29]: where the factor α M (k) encloses the scale and halo mass dependence. In practice, we find that, on large scales, α M (k) ∝ 1/k 2 , is independent of the halo mass. We assume that the large-scale linear halo bias for the Gaussian case is [23]: 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. [30][31][32]). Finally, the weighted effective halo bias is given by

13)
M min being the minimum halo mass hosting a source of the kind we are considering. Two things should be stressed when using eq. (2.11): a degeneracy between b G and f NL is expected at a given scale (the same amount of non-Gaussian bias can be given by different pairs of b G and 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 ); the 1/k 2 scale-dependence implies that large-scales are primarily affected by f NL , while small scales are mainly affected by b G (see a quantitative discussion for this in ref. [17]).
From the two left-hand panels of figure 1, one can see that a positive f NL enhances the amplitudes of the auto-correlation and cross-correlation power spectra especially on large angular scales (ℓ < 100). Consequently, the amplitudes of the ACF and CCF are also enhanced (right-hand panels of figure 1). The ACF remains positive up to angular scales JCAP08(2010)013 θ > 4 • where, for the adopted redshift distribution and a redshift-independent M min , the ACF is expected to become negative. The lower left-hand panel of figure 1 shows that the effects of non-Gaussianity on the cross-correlation power spectrum become very large for the lowest multipoles. Actually, most of the non-Gaussian contribution to the CCF comes from ℓ < 10, corresponding to angular scales 18 • . Its detection therefore requires uniform sky coverage up to very large scales, but this is very difficult to achieve observationally. For example, in the case of NVSS sources, very small modulations of the source surface density on very large scales, such as those of interest here, are easily swamped by systematic variations of the survey effective depth correlated with the varying observing conditions [18]. Also, the "integral constraint" [33] needs to be dealt carefully. In fact, the ACF and CCF estimators involve differences between local values of some quantity and their mean value over the considered area (see eqs. (3.2) and (3.3)). But if such area is not much larger than the correlation scale, the computed mean may differ from the true mean by a constant offset c, not known a priori, that must therefore be considered as a free parameter.

Observed ACF & CCF data
In this section we describe two LSS probes, NVSS radio sources and high redshift SDSS DR6 QSOs, and present observational determinations of their ACF and CCF with the CMB.

NVSS radio sources
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 [18]. Density gradients in the NVSS catalog become increasingly unimportant as the source brightness threshold is increased (see refs. [18,34]). Also we mask the strip |b| < 5 • , where the catalog may be substantially affected by 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. [35]. 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. [36]: The mean redshift of this sample is z = 1.23. From the NVSS catalog we construct a pixelized map using the HEALPix software package [37] with resolution N side = 64, yielding pixel areas of 0.92 Our ACF estimatorŵ(θ) reads: where n i is the number of sources in the i-th pixel, f i is the un-masked fraction of the pixel area,n is the expectation value for the number of objects in the pixel area [21]. The sum runs over all pixel pairs whose centers are separated by an angle θ. 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. We use N b = 9 angular bins, 1 • wide. The

JCAP08(2010)013
first one is centered at θ = 40 ′ (on smaller angular scales the dominant contribution to the angular correlation function comes from multiple components of individual radio galaxies); the other 8 are centered at 1 • · · · 8 • (on still larger scales the signal may be affected or even dominated by spurious density gradients [18,19]). The estimated NVSS ACF was previously found to be positive over the full range of angular scales we consider [17] (see also ref. [18]), although the integral of w(θ) over the full survey solid angle vanishes by construction. As pointed out by ref. [38], for a Gaussian distribution of primordial density fluctuations and a realistic redshift distribution for NVSS sources, the w(θ) must vanish at θ ≃ 4 • and become negative on larger scales unless the redshift dependence of the bias parameter for radio sources is drastically different from that of optical QSOs. Consistency with the clustering properties of the latter sources can be recovered allowing for a small amount of non-Gaussianity (f NL = 62 ± 27 (1 σ) [17]).
Since the NVSS covers most of the sky it can safely be assumed to provide a "fair sample" of the universe, so that the derived angular correlation function is not affected by the "integral constraint" significantly. Nevertheless, in ref. [17] we have explored whether the indications of an excess (compared to the Gaussian case) positive contribution to the ACF on large scales may be due to a unexpectedly large positive offset c (a negative c would obviously exacerbate the discrepancy with expectations from the Gaussian case). To this end we have added to w(θ) a constant c and have marginalized over it, allowing this quantity to vary in the range [10 −8 , 10 −4 ], where the lower limit means essentially zero (we work on a logarithmic scale) while the upper limit corresponds to the case where the offset accounts for the full clustering signal on the largest scales (see ref. [38]). The best-fit value of c is about ≃ 10 −5 , i.e. negligibly small, confirming that the data do not show indications that w(θ) is appreciably affected by the integral constraint. On the other hand, c is weakly constrained by the data, and values much larger than the adopted upper limit are formally allowed. Such large values of c would imply that w true (θ) = w estimated (θ)−c is negative on scales θ > 4 • , and the indications of non-Gaussianity would no longer be statistically significant. Such values however are not physically plausible because c is necessarily negligibly small when, as in the NVSS case, the survey area is so much larger than any other relevant scale (see section 9.4.2 of ref. [33]), even in the non-Gaussian case. Allowing also for the correction proposed by ref. [39] to account for the infrared divergence of the non-Gaussian halo correlation function the constraints on f NL , become f NL = 58 ± 28 (1 σ) [17].
A potentially trickier issue are large-scale surface density gradients due to instrumental effects that may spuriously enhance the ACF estimates. This issue has been extensively discussed by ref. [18] (see also ref. [19]) who concluded that this effect is negligible if we restrict ourselves to a flux limit of 10 mJy, as we did. An upper limit to the magnitude of the effect is an offset of 10 −4 , that, as noted above, would account for the estimated clustering signal on the largest scales. Having allowed c to take on this value, we have automatically allowed for the corresponding decrease of the statistical significance of f NL > 0. As a further test of the possibility that the large-scale clustering is spuriously enhanced by surface density gradients due to instrumental effects we have redone the analysis restricting ourselves to sources above 20 mJy, for which such effects are negligible on all scales, so that c = 0. In spite of the poorer statistics we find indications of f NL > 0 at approximately the same significance level as for the S ≥ 10 mJy sample: f NL = 73±32 (1 σ error) and 19 < f NL < 139 (95% confidence level).

JCAP08(2010)013
To measure the CCF between the NVSS number density map and the CMB map, we use the following estimator:ĉ where T i is the CMB temperature in the i-th pixel andT is the mean (monopole) value for the CMB temperature in the unmasked area. The results presented here rely on the 7 years Internal Linear Combination (ILC) map of the CMB provided by the Wilkinson Microwave Anisotropy Probe (WMAP) team [40]. We have checked that the using other CMB maps (the 5 years ILC map by the WMAP team and the improved (cleaner) map by ref. [41]) does not change the results in any significant way. We adopt the WMAP KQ75 mask, excluding about 30% of the sky at low Galactic latitude, to avoid most of the residual Galactic contamination. This mask is then combined with the NVSS sky coverage; regions of the NVSS catalogue that may introduce spurious features because of missing snapshot observations and over-dense regions associated with bright of extended sources (table 1 of ref. [19]) are also masked. As shown in section 2.2, the effect on the CCF of primordial non-Gaussianity shows up primarily on very large angular scales, at variance with the ACF for which the non-Gaussian signal is localized at θ = 2 • −5 • , as shown by ref. [17]. Thus the residual large-scale systematic fluctuations in the NVSS source density with declination due to the varying projection of the beam and the change in observing configuration, that ref. [18] found to be still present at the 10 mJy flux limit, may seriously bias the CCF estimate. While it is true that these fluctuations should be, on average, uncorrelated with the CMB signal and therefore, on average, should not contribute to the CCF signal, we deal here with a single realization of the Universe, and few independent modes; in practice, the jackknife procedure we employ to estimate the errors is not suited to deal properly with such an effect. We therefore restrict ourselves to S 1.4GHz ≥ 20 mJy, where such fluctuations are negligible. The normalized redshift distribution for the sample of ref. [35] does not change significantly if the flux limit is increased to 20 mJy, and therefore can still be described by eq. (3.1). The ACF for this sample does not differ in any significant way from that of the S 1.4GHz ≥ 10 mJy on the scales of interest.
The CCF is computed in N b = 16 angular bins, spaced by 1 • , with 0 • ≤ θ ≤ 15 • . The covariance matrix of the data points is estimated using the jackknife re-sampling method [42], dividing the unmasked area into M = 30 patches. From the 30 ACF and the CCF estimates obtained neglecting a patch in turn, we compute the diagonal (variance) and off-diagonal (covariance) elements of the covariance matrix. Our estimate of the NVSS-CMB CCF, shown in figure 2, is fully consistent with previous estimates using different approaches [14,34,[43][44][45].

SDSS DR6 QSOs
The SDSS DR6 QSO catalog released by ref. [13] 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. [13] 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 section4.2 of ref. [13] 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 JCAP08(2010)013 In the latter case, A amp is a free parameter and the fit is obtained for A amp = 1.14. It is indeed very close to the best-fit non-Gaussian curve, but there is no physical justification for A amp = 1.
spectrum (in this case we have N qso ≈ 6 × 10 5 QSOs). In order to minimize the effect of Galactic extinction on the observed QSO distribution, we mask regions with A g ≥ 0.18. We fit the redshift distribution dN/dz of the DR6 QSO sample with a function of the form: 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.
In spite of the high efficiency of the selection algorithm adopted to define the SDSS DR6 QSO catalog, some contamination from UV-excess stars is unavoidable. Following ref. [46], we model the observed ACF,ĉ tt , as the sum of contributions from QSOs,ĉ qq , and from stars, c ss , plus an offset, ǫ, arising from cross-terms: where a is the efficiency of the QSO catalog, i.e. the fraction of true QSOs. As shown by ref. [46], the KDE classification technique is efficient enough for the offset ǫ to be safely neglected. The contributions from stars and QSOs can be disentangled exploiting their different dependencies on θ. In fact refs. [46] and [21] showed thatĉ ss (θ) keeps almost flat up to large angular scales, and is therefore expected to dominate the signal on scales of several degrees. We have checked that the contribution from contaminating stars to the QSO-CMB CCF can safely be neglected.
As in the NVSS case we used the HEALPix software to pixelize the QSO map at the N side = 64 resolution and the jackknife re-sampling method to estimate the covariance matrix among data points. Our estimates of the ACF for the SDSS DR6 QSO catalogue and of the JCAP08(2010)013 QSO-CMB CCF, shown in figure 3, are consistent with previous analyses [21,43]. Note that we do not split the SDSS QSO in different redshift bins using the objects' photometric redshift estimates, but we consider the projected ACF (and CCF) of the full sample.

Method & data analysis
We perform a global fitting of cosmological parameters, including f NL , for the data of section 3 including also the datasets described below, using the Besides these six basic cosmological parameters, we have three more parameters related to the ACF and CCF data: the non-Gaussianity parameter f NL , the minimal halo mass M min , two offset parameters c, one for the ACF and one for the CCF, allowing for the effect of the integral constraint. The ranges for f NL and M min are the same as in ref. [17]. Those for the offsets were chosen to be large enough to account for the large-scale signals.

JCAP08(2010)013
In the QSO case we have, in addition, the efficiency of QSO classification a. Several authors also treated as an additional free parameter the ISW amplitude A amp defined bȳ c gT (θ) = A amp c gT (θ), wherec gT and c gT are the observed and theoretical CCF. We note that values of A amp = 1 have no physical meaning, and are indicative of some unrecognized problem in the interpretation of the data. Therefore we fix A amp = 1 in all our calculations, except to estimate the significance of the ISW signal, in order to compare with previous works in which this parameter was used.
The model ACF c gg (θ) and CCF c gT (θ) are compared with the observed valuesĉ gg (θ) and CCFĉ gT (θ), respectively, through the Gaussian likelihood function: where C ij and C ′ ij are the observed ACF and CCF covariance matrices. 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 [48]. The WMAP7 data are exploited only to improve the constraints on the six basic cosmological parameters, not to constrain f NL .
The BAOs [49] 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. [50]: We adopt these values as Gaussian priors. The SNIa data yield the luminosity distance as a function of redshift which provides strong constraints on the dark energy evolution. We use the Union compilation data (307 samples) from the Supernova Cosmology project [52], which include the samples of SNIa from the (Supernovae Legacy Survey) SNLS and from the ESSENCE survey, and span the redshift range 0 z 1.55. In the calculation of the likelihood from SNIa we marginalize over the nuisance parameter as done in refs. [53,54].
Furthermore, we add a prior on the Hubble constant, H 0 = 74.2 ± 3.6 km/s/Mpc given by ref. [55]. Finally, in the analyses of NVSS sources and SDSS QSOs we set the minimal halo mass at M min > 10 12 h −1 M ⊙ consistent with the results of ref. [56].

NVSS radio sources
To check the significance of the detection of the late-time ISW effect from the NVSS-CMB CCF, following ref. [21], we compute the best fit value of the amplitude A amp in the Gaussian case (f NL = 0) and the jackknife error. We obtain: A amp = 1.14 ± 0.50 (1 σ error) , (5.1) so that the significance is ≃ A amp /σ A , i.e. ≃ 2.3 σ, broadly consistent with earlier analyses [43][44][45] and very well compatible with its physical value (unlike ref. [14]).
The constraints on f NL , obtained setting A amp ≡ 1 and marginalizing over all the other parameters, are summarized in table 1. The NVSS CCF turns out to be only weakly sensitive to the minimum halo mass. Therefore we fix it to the best fit value from NVSS ACF data JCAP08(2010)013  Figure 5. Cross-correlation between the ILC CMB and a Galactic emission template comprising synchrotron, free-free and dust contributions (black squares). The Galactic template is our best guess of the residual Galactic contamination in the ILC obtained as described in ref. [58]. Error bars are jackknife-estimated.
alone, M min = 10 12.47 h −1 M ⊙ . While the NVSS ACF yields a positive f NL at more than 2 σ confidence level, f NL = 58 ± 28 (1 σ), as previously found by ref. [17], the NVSS-CMB CCF provides much weaker constraints: f NL = 29 ± 48 (1 σ). The best fit value for the offset is c ≃ −0.05 for S 1.4GHz ≥ 20mJy. To check the effect of residual systematic fluctuations in the NVSS source density, we have repeated the calculation for S 1.4GHz ≥ 10 mJy. We find that the amplitude of the CCF decreases, consistent with the cross-correlation being somewhat blurred (a similar effect can be gleaned also from figure 11 of ref. [34]), and the constraints on f NL become f NL = 6 ± 37 (1 σ).
Other systematic effects may affect the CCF estimates (for a parallel analysis on the ACF see ref. [17]). A first possibility is that the CMB map is contaminated by residual foreground emissions that increase the noise and thus swamp the signal we are looking for. The contamination from point sources below the detection limit is increasingly diluted for larger and larger pixel sizes. However we find that the CCF estimate does not change in any appreciable way when we use larger pixel sizes (N side = 32 and N side = 16). Note that increasing the pixel size we also decrease the relative contribution of CMB fluctuations generated at the recombination, which have much higher amplitude than those due to the ISW effect.
The impact of residual contamination from diffuse Galactic emissions is assessed by cross-correlating the CMB map with the standard dust, free-free and synchrotron templates, used also by the WMAP team (e.g., ref. [57]), smoothed to 1 deg resolution. No indications of significant cross-correlations are found. However, residual foreground contamination will be due a certain mixture of the different Galactic emissions, thus it is not expected to strongly correlate with single foreground templates. We thus build a typical foreground residual map to be correlated with the CMB map, as described in the following.
Ref. [58] shows how to predict the residual contamination after the application of a linear component separation technique. The method requires to have a model of the data and to know the weights of the linear mixture for component separation.
Our data model is composed as follows: for the frequency dependencies of the components is based on model M2 of ref. [59], obtained from an analysis of the WMAP 3 years data; for the spatial morphology is based on the previously mentioned foreground templates. The component-separation linear-mixture weights are the ILC weights adopted to build the CMB map, available from the LAMBDA site.
The foreground residual map we obtain is found to have a statistically significant negative CCF with the ILC map (figure 5), indicating that the separation of the Galactic emission from the CMB map is not perfect. The amplitude of this spurious cross-correlation is however too small to affect significantly the NVSS-CMB CCF.
If we combine with the NVSS ACF and CCF data we obtain the following constraint on f NL : f NL = 53 ± 25 (1 σ error), (5.2) or 10 < f NL < 106 at 95% confidence level. This result is compatible with previous estimates [15,[60][61][62][63][64][65][66], and with the WMAP7 limits [48]. In figure 6 we plot the two-dimensional constraints on (f NL , c ACF ) and (f NL , c CCF ) from the NVSS ACF and CCF data. As noted above, a positive offset means that the true ACF or CCF is lower than the estimated one, thus weakening or spoiling indications for non-Gaussianity. On the contrary, a negative offset points to larger values of f NL . As discussed in section 3.1, the integral constraint should have a negligible effect on the NVSS ACF but an offset of up to 10 −4 may be induced by large scale density gradients of instrumental origin. This value has been adopted as the upper boundary for the allowed range for c ACF .

SDSS DR6 QSOs
As in section 5.1, to assess the significance of the detection of a cross-correlation between the distribution of SDSS DR6 QSOs and the CMB map we compute the best fit value of A amp and its error keeping all the cosmological parameters fixed at the WMAP7 values. We find: the NVSS area and its effective angular size is not much larger than the scale over which positive correlations are induced by the f NL value indicated by the NVSS ACF data. As for the NVSS case we added constants c, as free parameters, to the ACF and to the CCF, and marginalized over them. In the ACF case, the effect of this constant is negligible. Its best fit value is c ACF = −2 · 10 −4 , and we find f NL = 42 ± 28 (1 σ); if we set c ACF = 0 we get f NL = 32 ± 19 (1 σ). As expected, the impact of the constant offset is much more relevant for the CCF. We find c CCF = −0.15 and f NL = 60 ± 42 (1 σ); setting c CCF = 0 we get f NL = 0 ± 10 (1 σ). Clearly marginalizing over the c CCF parameter has a big effect on the recovered f NL value. The larger the sky coverage of the survey the better this constant is determined and the less is its impact on f NL . Combining the ACF and CCF data we find f NL = 58 ± 24 (1 σ), c ACF = −3 · 10 −4 , c CCF = −0.13. The value of f NL is nicely consistent with that obtained from the analysis of NVSS data. The best fit value of the selection efficiency of the QSO catalogue, kept as a free parameter, is a = 97.1% ± 1%, consistent with the estimate by ref. [13] who claim an efficiency of over 97% for the UVX sub-sample considered here.
Finally, in figure 8 we plot the two-dimensional constraints on (f NL , 10 4 c ACF ) and (f NL , c CCF ) from the SDSS QSO ACF and CCF data. Both c ACF and c CCF are anti-correlated with f NL . A negative c ACF and c CCF enhances the estimate of f NL .

Discussion & conclusions
All the datasets we considered here yield broadly consistent results. The NVSS ACF drives the positive signal 16 < f NL < 114 at 95% confidence level. The SDSS QSO ACF and CCF individually yield results consistent with f NL = 0. This effect arises because of degeneracies in a high-dimensional parameter space. In particular, there is a degeneracy in the M min − f NL plane, with large values of M min corresponding to low f NL values, which, when one marginalizes over M min lowers the mean f NL value. The combination ACF+CCF however seems to break this degeneracy (ACF depends on b 2 while CCF depends on b) cutting out the large M min tail and thus increasing the recovered f NL value. In addition the relatively large and negative c CCF skews the maximum posterior f NL value for the CCF to larger positive values. Before drawing our conclusions we compare our findings with previous work in the literature where very similar data-sets were used.

Comparison with previous works
The work presented here is closely related to refs. [15,16] and to the NVSS-CMB cross power spectrum as obtained by ref. [14]. The most immediate difference in the analysis is that we work with the ACF and CCF while previous works has used the auto and cross power spectra. Correlation functions and power spectra are spherical harmonics transform pairs and should thus be fully equivalent. The two approaches however are effected by systematic effects in slightly different ways: whenever possible we try to directly compare the two approaches. In addition, compared with previous analyses, we marginalize over a suite of possible systematic effects, modeled by the constant offsets c ACF and c CCF and over a possible stellar contamination and QSO efficiency parameter. Ref. [16] used the NVSS-WMAP 5 year map cross power spectrum as provided by ref. [14] to obtain f NL = 236 ± 127 (1 σ). The main differences with our NVSS CCF analysis are: i) ref. [16] used a constant bias, here we use eq. (2.12), we include the non-Gaussianity correction in the mass function, thus have a free parameter given by the minimum host halo mass for the sources; ii) refs. [14,16] used WMAP 5 years map, while here we use up-todate cleaned CMB maps; iii) ref. [14] use NVSS sources much dimmer (≥ 2.5 mJy) than our sample (≥ 20 mJy). They have a higher source density but we have seen that density gradients are most important for dim sources. Finally, we do not attempt to model the redshift distribution of NVSS sources but use directly the one that has been measured; this results in our errors on f NL to be smaller, but the two results are fully consistent.
Ref. [15] obtained non-Gaussianity constraints using both the NVSS and SDSS QSO data. For their NVSS analysis they also used it only in cross correlation with the CMB (they do not use the NVSS sources power spectrum) and, like ref. [16], used the cross-correlation data and errors of ref. [14]. In their analysis the combinations: bias times source redshift distribution, b·dn/dz, and the bias as function of redshift, b(z), as well as f NL are determined JCAP08(2010)013 from the correlation function data themselves, here instead we use the actual redshift distribution and use the minimum halo host mass as a parameter, rather than b(z = 0). In their reported ISW results several data sets are combined: the SDSS QSOs, the NVSS, 2MASS etc. although they report their results to be dominated by NVSS. Our NVSS CCF constraints on f NL are consistent with theirs although with smaller error-bars due to the fact that we do not need to estimate the source redshift distribution from the correlations themselves.
For their SDSS QSO analysis ref. [15] considered both the auto and cross correlation signals. Ref. [15] used an extension of the DR3 QSO sample that include many of the sources that subsequently were released with the DR6 sample: we use the more complete and better calibrated final official SDSS DR6 QSO catalog release [13]. Ref. [15] found that their QSO sample at z < 1.45 seem to suffer from contamination which they described to systematic calibration errors, and discarded this sample from their analysis. 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 or that the two samples give different estimates of f NL . Our sample is more than double the size than theirs.
In our analysis we do not split the QSO in redshift bins according to each source photometric redshift estimate as ref. [15] do, but we consider the projected correlation of the full sample. Ref. [15] had therefore to estimate the source redshift distribution from the (cross)correlation signals themselves for each of their sub-samples; we avoid this by using the independently-estimated full sample source redshift distribution. We also leave the stellar contamination as a free parameter to marginalize over, while ref. [15] kept it fixed at the fiducial value.
The two analyses also use different bias models: ref. [15] used a functional form b(z) ∼ (1 + z) 5 or b(z) ∼ 1/D(z) with different constant offset depending on assumptions about "recent merging" activity. We use the extended Press-Schecter bias formulation of eq. (2.12).
We have however verified that if we use their bias model, fix the QSO efficiency (i.e. fraction of stellar contamination in the catalog) to the fiducial value and do not marginalize over the constant offset we recover even central f NL values fully consistent with theirs (f NL = 6 ± 18 (1 σ)). We therefore conclude that the largest effect driving the different central f NL values we recover is the treatment of systematic effects -stellar contamination and offset.
We should stress here that for the data in common between the different analyses, the results we find are fully consistent (i.e., at the 1 σ confidence level) with both refs. [16] and [15] results. Our positive f NL signal is driven by the NVSS auto-correlation function and by the combination of SDSS DR6 QSO ACF and CCF where the combination breaks the degeneracy with M min pushing he marginalized f NL values slightly up.
Better knowledge of the bias of the sample and of systematic effects such as stellar contamination, calibration and selection effects are clearly needed to make further progress on this front; this is especially relevant to future data-sets where the large volume surveyed will further reduce the statistical errors.

Conclusions
We have investigated the constraints on the parameter f NL , characterizing primordial non-Gaussianity of the local type, coming from the currently available data sets best suited for this purpose, namely the observationally determined auto-correlation functions (ACFs) of NVSS radio sources and the SDSS-DR6 QSOs, and the cross-correlation functions (CCFs) of the surface density distributions of these sources with the WMAP 7 years CMB map. The mean redshift of NVSS sources isz ≃ 1.23, and that of SDSS-DR6 QSOs isz ≃ 1.49.
A signature of f NL > 0 should also be present in the CCF between tracers of large scale structure and the map of CMB anisotropies. Most of the signal, however, comes from low multipoles (ℓ < 10) corresponding to very large angular scales where the signal may be blurred by systematic effects. In the case of NVSS sources we have indeed found that the amplitude of the CCF with the CMB increases somewhat as we increase the adopted flux limit from S 1.4GHz = 10 mJy, which is found to be adequate for the ACF analysis but is still affected by small spurious large scale surface density gradients, to S 1.4GHz = 20 mJy where such gradients become negligible. Other possible systematic effects, such as contamination of the CMB map by point sources below the detection limit or by residual Galactic emissions have been shown to be of minor importance. The NVSS-CMB CCF however turned out to provide rather loose constraints on f NL (f NL = 29 ± 48 (1 σ)).
The limited area covered by the SDSS-DR6 QSO catalog implies that the CCF is liable to a quite substantial effect of the integral constraint, at variance with the NVSS case, where the area is 4 times larger (the ACFs of both NVSS sources and SDSS QSOs are found to be insensitive to this effect). We have allowed for the integral constraint by adding to the estimated CCF a constant offset, treated as a free parameter, and marginalizing over it. In this way we obtained f NL = 60 ± 42 (1 σ).
Finally, we obtain the constraints on f NL = 53 ± 25 (1 σ) and f NL = 58 ± 24 (1 σ) from NVSS data and QSO data, respectively. The tantalizing hints of a positive f NL reported by ref. [17] have survived the present more extensive analysis, involving other data sets. In order to reconcile the data with f NL = 0, we need to show that either the NVSS catalog is affected by many-times larger systematic effects (spurious large-scale correlations introduced by density gradients) than known and quantified so far, or that the redshift evolution of the bias of these sources is radically different from that of their radio-quiet, optical counterpart.
Alternatively, one may question the accuracy of jackknife error-estimation especially when there are many non-zero off-diagonal elements in the covariance matrix. From ref. [67], we infer that jackknife can underestimate parameter errors by up to 30%. If the error bars were to be increased, f NL would become more compatible with zero. Much tighter and robust constraints on f NL should be provided by the forthcoming surveys with ASKAP [68] and MeerKAT [69], the SKA pathfinders, and by much larger area QSO surveys.