Inhomogeneous non-Gaussianity

We propose a method to probe higher-order correlators of the primordial density field through the inhomogeneity of local non-Gaussian parameters, such as f_NL, measured within smaller patches of the sky. Correlators between n-point functions measured in one patch of the sky and k-point functions measured in another patch depend upon the (n+k)-point functions over the entire sky. The inhomogeneity of non-Gaussian parameters may be a feasible way to detect or constrain higher-order correlators in local models of non-Gaussianity, as well as to distinguish between single and multiple-source scenarios for generating the primordial density perturbation, and more generally to probe the details of inflationary physics.


I. INTRODUCTION
Primordial non-Gaussianity is a window onto the early universe, which can be used to distinguish between different models for the origin of cosmological fluctuations [1][2][3][4][5]. For this reason, it is important to characterize as completely as possible the properties of observable quantities in a non-Gaussian distribution. In this paper, we discuss under what conditions local non-Gaussian observables can be inhomogeneous to a measurable level.
Consider subdividing the observable sky of size H −1 into a large number of small patches of size ℓ ≪ H −1 , and measuring cosmological observables on each of these patches. The effective values of the cosmological parameters can differ from one patch to the other, since each small part of the sky experiences a different local background. In particular primordial curvature fluctuations within a small patch are defined with respect to a local background which can be decomposed into a cosmological background solution of the homogeneous field equations, plus the cumulative effect of random modes with wavelengths larger than the size ℓ of the patch (but smaller than the size H −1 of the observable sky). Since different patches of our observable universe are affected by different contributions from the random long wavelength modes, each patch is defined on a different background configuration. As a consequence, quantities characterizing the properties of n-point functions of curvature fluctuations can be different in one patch or the other. Indeed the effect of long-wavelength modes on the small-scale gravitational collapse of dark matter halos is described by the peak-background split in astrophysics and is used to explain the biased clustering of dark matter halos [6]. In a non-Gaussian field the long-wavelength modes also modulate the amplitude of the short-wavelength density perturbations which can lead to a distinctive scale-dependent bias which is a powerful probe of primordial non-Gaussianity [7]. Of course there is also a cosmic variance between patches, even if there are no non-zero higher order correlators. However this effect will be small as long as there are a large number of independent measurements in every patch.
Correlators between n-point functions measured in one patch, and k-point functions measured in another patch, depend on (n + k)-point functions in the entire sky. As an example, we will show that the autocorrelation of f NL between different patches depends on a class of parameters controlling 5-and 6-point functions, that in the single-source case reduce to squares of trispectrum parameters. This implies that inhomogeneities in quantities characterizing non-Gaussian observables can be measurable, since at present we have only weak or even no constraints on parameters controlling higher-order point functions in the full sky. Moreover, although the optimal way to constrain higher n-point functions in the full sky would be through a direct analysis of that correlator, in practice this is exceedingly time consuming and computer intensive. This problem will become significantly more pressing when the Planck satellite data becomes avalailable in early 2013 [8], since it has a much higher resolution than the WMAP satellite [9]. Only during the last few years have the first constraints on the trispectrum been made [10][11][12][13][14], and there are not yet any constraints on higher n-point functions. Therefore, our method of relating n-point functions to inhomogeneities of lower-order functions should also provide a feasible and practical way of constraining higher-point functions in the full sky, even though the method is suboptimal. There is a clear theoretical interest in constraining higher-order correlators: they may be much larger than the correlators which have already been constrained, and their detection would provide important additional information for characterizing the origin of primordial fluctuations. Refs. [15][16][17] consider the amplitude of correlators of higher order than the trispectrum for local models of non-Gaussianity. See [18] for an analysis of non-Gaussianity not characterized by a local shape, that investigates 10-point functions and beyond.
The inhomogeneities of non-Gaussian quantities that we analyze in this paper are distinct from the anisotropies of inflationary observables, arising when background vector fields (perhaps of curvaton nature) are turned on [19][20][21][22][23][24][25][26][27]. We are going to consider primordial non-Gaussianity due to scalar fields that do not induce statistical anisotropy. Conversely, generic models of inflation with vector fields may produce statistical anisotropy but do not induce inhomogeneities of inflationary parameters. The inhomogeneities that we consider here are also different from the scale-dependence of the non-linearity parameters, that describe a variation with size rather than with position of patches [28].
The plan of this paper is as follows: We first focus on single-source scenarios, in which only one field generates the curvature perturbation. In Sec. II A we consider correlators of ζ and the power spectrum, in Sec. II B of f NL and in Sec. II C possible tests of single versus multi-source scenarios. In Sec. III we focus on the more general formula of multi-source models, before going on to analyze in Sec. IV when the inhomogeneities can be large enough to be observable. Finally we conclude in Sec. V.

II. STATISTICAL INHOMOGENEITIES IN SINGLE-SOURCE MODELS
Let us start by considering single-source models where the curvature perturbation ζ arises from fluctuations of a single scalar field. The single-source scenarios that we have in mind do not necessarily correspond to single-field models of inflation. We have in mind set-ups in which more than one scalar field may be present in the system, but only one of them, which we dub σ, is responsible for generating the primordial curvature perturbation. Examples are the curvaton model [29,30], or the modulated reheating scenario [31,32]. These models are particularly interesting since they are capable of generating large non-Gaussianity of the local type, potentially observable by the Planck satellite.
The statistics of the curvature perturbation in general depend on the size of the patch we consider. This is due to long-wavelength fluctuations which contribute to local background quantities in any patch smaller than the horizon. The curvature perturbation within a patch of size ℓ is given by [33][34][35][36][37] where σ ℓ denotes the background field value in the patch. It consists of the classical homogenous solutionσ(t), and of fluctuations δσ k with wavelengths greater than the patch size, k < a ℓ −1 (here a denotes the scale factor). N (σ) = H dt denotes the number of e-foldings from an initial time t i , soon after horizon exit during inflation of the modes of interest, to some final time (for example, during the primordial radiation-dominated era) at which ζ has frozen to its final, constant value. The fluctuations δ ℓ σ(x) include modes ranging from k = a ℓ −1 to k i = a i H i , the latter corresponding to modes that exit the horizon at the initial time t i 1 : We use a top-hat window function to select the modes within this momentum interval. Different choices of the window function should not significantly alter our results. δ ℓ σ(x) consists of superhorizon modes at the initial time, δσ k , which can subsequently be treated as a classical, Gaussian random field with the two-point function given by 2 Here the brackets denote ensemble averages. We are interested in computing ζ in patches smaller than the observable universe, ℓ ≪ H −1 , see Fig 1. The local background field value in a patch centered at some position x is given by where σ H −1 denotes the background field in the entire observable universe. ∆ ℓ σ x is comprised of fluctuations δσ k with wavelengths greater than the patch size, k < a ℓ −1 , which do not average out when computing spatial averages over the patch ℓ. For a patch ℓ located at a position x, ∆ ℓ σ x therefore acts as a constant, local background. We emphasize that we use the label x to indicate quantities evaluated in a patch at position x, and we consider them not as functions of the coordinate x. Denoting the spatial average over the patch by . . . ℓ , we in general obtain the non-vanishing result which depends on the location x of the patch. On the other hand, when ℓ ≪ H −1 , the spatial average of ∆ ℓ σ x , when computed in the entire observable universe H −1 , vanishes: ∆ ℓ σ x H −1 = 0. This is because ∆ ℓ σ x in our approximation depends linearly on the fluctuations δσ k , which have a vanishing spatial average over the full sky. By the ergodic theorem, and provided that the quantity ℓ · H is sufficiently small, the average of ∆ ℓ σ x computed over the full sky coincides with the ensemble average, ∆ ℓ σ x H −1 = ∆ ℓ σ x = 0. The consequences of long-wavelength fluctuations have been extensively studied in the context of the infrared growth of inflationary correlators. It has been shown that, for sufficiently local observations, a suitable shift of background quantities is able to remove infrared divergences associated with adiabatic field fluctuations [38][39][40][41], see also [42][43][44]. On the other hand, infrared-enhanced effects associated with long wavelength isocurvature fluctuations cannot be removed by shifts in background quantities, and have the opportunity to provide sizable contributions to inflationary observables [45]. In our analysis, we will consider correlators between cosmological observables measured on distinct subhorizon patches within the observable universe, in set-ups where longwavelength modes generate sizable statistical inhomogeneities. This effect was pointed out already in [46], but for the specific case of single-field slow-roll inflation where inhomogeneities are slow-roll suppressed. Our treatment applies to generic single and multiple source models, and we find that the inhomogeneities can become large in scenarios characterized by observable non-Gaussianity. We will only consider the inhomogeneities arising from scalar perturbations. Tensor modes, included in the analysis of [38,39], should only generate subleading, slow roll suppressed corrections to our results which we neglect.
After this general discussion aimed to define the set-up, we move on discussing in more detail small scale statistical inhomogeneities of quantities characterizing the properties of the curvature perturbation.

A. Correlators of one and two-point functions
The two-point function of the curvature perturbation (1) in a patch of size ℓ ≪ H −1 is given by If we re-express this with respect to the background field in the entire observable universe, σ H −1 , using Eq. (4) this becomes where x denotes the location of the patch and y 1 , y 2 are coordinates inside the patch. Recall that brackets ℓ denote spatial averages computed over the region ℓ. The long-wavelength contribution ∆ ℓ σ x breaks translational invariance: the two-point function depends not only on the separation of the points, |y 1 − y 2 |, but also on the location of the patch x. The long-wavelength modes therefore generate statistical inhomogeneity. For the case of 2-pt functions, this fact was also pointed out in [46] and, in more generality, in [48]. Taking a Fourier transformation of (7) within the region ℓ, i.e., integrating over the y-coordinates and treating x as a constant, we find a relation between the spectra of ζ ℓ and ζ H −1 where the full sky power spectrum is and f NL describes the amplitude of the primordial bispectrum relative to the square of the power spectrum, which is given in the δN -approach by [37] f NL = 5 6 ∆ ℓ P ζ (x) measures deviations of the local spectrum, measured in a patch of size ℓ, from the global power spectrum characterizing 2-pt functions of perturbations over the entire observable universe. Since ∆ ℓ P ζ (x) is proportional to ∆ ℓ σ x , it is a Gaussian field with zero mean over the full sky for ℓ ≪ H −1 where . . . H −1 denotes the spatial average computed over the entire observable universe, which we assume coincides with the ensemble average . . . . Hence we will drop the label H −1 from the angle brackets in what follows. We denote ∆x = |x 1 − x 2 |, k ℓ = a ℓ −1 and k hor = a H. The integral can be evaluated as with ci(x) the cosine integral function. Notice that, as expected, the correlators between two-point functions evaluated in different patches of the universe are proportional to 6 5 f NL 2 = τ NL , which in single-source models is a parameter characterizing the four-point function measured in the entire sky. 3 In the limit ∆x → 0, the result takes a particularly simple form corresponding to the well-known logarithm associated with IR-enhancements. For where γ E ≈ 0.58 is the Euler constant. A comparison between the exact result (12) and the estimate (14) is shown in Figure 2.
The magnitude of the statistical inhomogeneities is described by the variance of ∆ ℓ P ζ . Using COBE normalization P ζ = 2.4 × 10 −9 and the bound |f NL | 10 2 we find hor − k −1 ℓ , corresponding to the possible distance between non-overlapping spheres of coordinate radius k −1 ℓ /2 within the observable universe.
For patches ℓ corresponding to a few percent fraction of the observable universe, the logarithm contributes with a factor of order unity, |ln (ℓH)| ≃ O(1). 4 The statistical inhomogeneity seen on these scales could therefore be at a few percent level. This observation was already made in [46] (see also [47], and the more general discussion of [48]). Observational constraints were discussed in [49]. It is also interesting to consider the correlators of the quantity , measuring the difference in the amplitudes of the curvature perturbation evaluated in the small patch ℓ and in the full sky H −1 . To first order in perturbations we find The correlator between ∆ ℓ ζ(x)'s evaluated on different patches reads that is proportional to the two point function measured in the full sky. Also, we have that is proportional to the three point function.
In this paper we will focus on the two-point correlators of n-point functions measured in separate patches. It would also be possible to extend our analysis to higher point correlators of the same quantities between different patches. This would require extending our analysis to non-linear orders in ∆ ℓ σ x and, depending on the model, also extending the definition of ∆ ℓ σ x , given in Eq. (4), to non-linear orders in perturbations. This should be technically straightforward although it would lead to longer expressions. We will not pursue this here, also since the results will be suppressed by higher powers of the power spectrum. Measuring the higher-order correlators of observables in different patches would also become computationally more demanding, something which our analysis, constraining the higher-point correlators, alleviates.

B. Correlators of three-point functions
Because inhomogeneity of the power spectrum is already restricted by observational bounds on the primordial bispectrum, it is more interesting to study the correlators of non-Gaussian observables such as f NL itself. The non-linearity parameter f NL in a patch ℓ, located at x, can be expressed as The prime denotes a derivative with respect to σ H −1 and for single-source models we find where g NL is a parameter describing the primordial trispectrum relative to the power spectrum cubed, and is given in the δN -approach by [50] g NL = 25 54 For brevity we suppress the argument σ H −1 in what follows, and write f NL (σ H −1 ) ≡ f NL etc meaning that quantities without argument are evaluated in the full sky.
To leading order in ∆ ℓ σ x , the difference between f NL measured in the entire observable universe H −1 and in a patch of size ℓ is then given by and the two-point function of ∆ ℓ f NL (x) reads We then learn that ∆ ℓ f 2 NL is proportional to a combination of terms containing powers of g NL and f NL . As we will discuss in the next sections, those terms are related, by consistency relations valid in the single-source limit, to parameters controlling five and six point functions.
Using the amplitude of spectrum of curvature perturbations, P ζ = 2.4 × 10 −9 , we find For |g NL | f 2 NL the inhomogeneities are unobservably small, (∆ ℓ f NL ) 2 1, as f NL is constrained by |f NL | 10 2 [9]. On the other hand, if |g NL | ≫ f 2 NL , we discuss this possibility further in the next subsection, the inhomogeneities can become significant as the observational constraint for g NL is rather weak, |g NL | 10 6 [10][11][12][13][14]. In this limit, the last two terms in eq. (24), proportional to powers of f NL , give only subleading contributions and they can therefore be neglected. The function F (∆x, k hor , k ℓ ) contributes with a factor of order unity, and thus we arrive at the result For |g NL | ∼ 10 6 , the variation of f NL measured in different patches of size ℓ < H −1 becomes quite large, (∆ ℓ f NL ) 2 ∼ 10 4 , implying that f NL becomes inhomogeneous on small scales 5 . This is a generic feature of all single-source models with large g NL . 6 It is illuminating to compare the result with the observational accuracy for detecting f NL using a small fraction of the sky. For example, the BOOMERanG experiment found −670 < f NL < 30 at 65% CL observing a 3% region of the sky [51]. For patches of the same size, we find the variance (∆ ℓ f NL ) 2 ≈ 1.4 × 10 −8 g 2 NL using Eq. (24). For |g NL | ∼ 10 6 , the variation of f NL measured on different 3%-of-the-sky patches is therefore of the same order of magnitude as the BOOMERanG 1-σ accuracy for detecting f NL .
In general, when estimating n-point functions on a limited part of the sky, observational errors increase (at least) proportionally to 1/ f sky (f sky denotes the fraction of the sky which is observed), since there is less data available in a small patch. Therefore when dividing the sky into, say, one hundred equal patches, the error bar on every observation will grow by at least a factor of ten compared to the full sky measurement. On the other hand, we then correlate, over the full sky, the measurements of n and k-point functions evaluated in distinct small patches. This implies that, for the aim of constraining the values of higher-order point functions in the full sky, we are able to regain much of the accuracy we had previously lost.
There is however an important second effect to take into account: when making measurements of n-point functions in a patch smaller than the entire sky, a smaller number of correlators can be built. For example, for n = 3, when estimating f NL in a patch the range of scales available is reduced, thereby limiting our ability to analyze the squeezed limit of triangles 7 . This effect is expected to depend only logarithmically on the size of the patches. Indeed, the signal to noise ratio for measurements of local non-Gaussianity scales as [52,53] where ℓ max and ℓ min denote the largest and smallest multipoles accessible by the experiment. For example, for patches covering one per cents of the sky, ℓ min /ℓ max ≪ 1, and the patch size, which determines ℓ min , enters only in the logarithm, reducing the accuracy only by a factor of order three with respect to full sky. We will not develop this issue further in our work, but to account for these logarithmic effects we simply increase by one order of magnitude the values of our expected constraints on higher order n-point functions in the full sky. We will see that, even within this conservative estimate, we can find interesting constraints on these quantities. To conclude this section, let us point out that alternative techniques based on needlets analysis of CMB data are particularly well suited for testing non-Gaussianity in selected regions of the sky. This fact has already been used for investigating inhomogeneity of non-Gaussianity in different, large regions of the sky in [54,55], with the main aim of investigating foreground contaminations and directional variations of f NL . Although no significant hints of anisotropies have been reported in those studies, it would be very interesting to apply the same techniques to a collection of smaller regions of the sky, to instead test inhomogeneities of f NL along the lines we suggest here.
C. New perspectives for discriminating single from multiple source models We have learned from the previous analysis that a large g NL leads to sizable inhomogeneities of f NL . On the other hand, other inflationary observables can be affected by a particularly large g NL . We discuss in this section this possibility, showing that the conditions allowing to obtain large inhomogeneities for f NL offer new perspectives for distinguishing single from multiple-source scenarios.
Loop corrections [56] are known to influence observables by providing logarithmic contributions that in some scenarios can be large enough to dominate over tree-level quantities [57,58]. In our set-up, τ NL is the observable that is most sensitive to loop corrections when g NL is large. Expressing it as one finds that for |g NL | ≫ |f NL | 2 , the dominant part of the one-loop correction to τ NL is given by assuming that the fourth order derivative N ′′′′ is constrained according to footnote 5. We have used approximation methods similar to [59][60][61] in evaluating the loop integral. The cut-off scale l should be a bit larger than the box in which we are making our measurements, such that ln(kℓ) 1: because we are making a leading log approximation one cannot consider the limit of the logarithm going to zero. Comparing equations (23) and (27), we observe that the inhomogeneities of f NL are directly related to the fraction of τ NL generated by loops This consistency relation, which holds for |g NL | ≫ |f NL | 2 , is specific for the single-source case. As we will learn in Sec IV A, in multiple source models there are additional quantities contributing to the right hand side of Eq. (28), which in principle make it possible to discriminate between single and multiple source scenarios.
Although for some models one cannot realise |g NL | ≫ f 2 NL , especially for those with quadratic potentials, for a review see [62], it is certainly possible to realise |g NL | ≫ f 2 NL in general early universe models. This can happen for example in the interacting curvaton scenario [63], since the value of f NL oscillates depending on the initial field value of the curvaton, and the points where f NL = 0 do not coincide with points where g NL is suppressed.
Also an isocurvature field direction during inflation with a quartic potential and zero VEV gives rise to g NL but not f NL [64]. In general finite volume effects will also generate a non-zero f NL in horizon size patches [65,66], but in some patches f NL will still be small while g NL is not. Some models of multifield inflation can also give rise to |g NL | ≫ f 2 NL at the end of inflation, see Eq. (82) of [67] (see also [68]). Interestingly, large loop corrections to τ NL also break the well-known consistency relation τ NL = (6f NL /5) 2 , characteristic for single-source models. This relation is indeed not protected against loop corrections, as can be nicely demonstrated by considering models where the curvature perturbation takes the form ζ = N ′ δφ + N ′′′ δφ 3 /6, with no higher order terms. For this class of models f NL = 0 to all orders in perturbation theory, within the same patch that this ansatz for ζ is valid in. At tree-level τ NL also vanishes, τ tree NL = 0, but the unique one-loop corrections generate a non-vanishing result, τ NL ≈ (54/25) 2 g 2 NL P ζ ln(kℓ) + O(g 3 NL P 2 ζ ) 10 4 . There is no two-loop correction and the three-loop correction is constrained to be much less than unity due to the observational bound on g NL . The relation τ NL = (6f NL /5) 2 is therefore clearly violated by the loop corrections. To the best of our knowledge, this has not been demonstrated previously. However the Suyama-Yamaguchi (SY) inequality, τ NL (6f NL /5) 2 , which was originally proved at tree level in [69], and the possible effects of loop corrections were discussed in [61,62,70], has recently been shown to remain true against loop corrections at all orders in a completely model independent sense [71] (although there are subtelties about the definition of the local correlators beyond tree level [61]).
The previous example is by no means a special case: for generic single-source models we also find that the relation τ NL = (6f NL /5) 2 can easily be broken by loops. Assuming there are no large loop corrections with higher than three-point vertices, and concentrating on the limit |g NL | ≫ f 2 NL , the one-loop corrected f NL is given by f NL ≈ f tree NL (1 + 13 g NL P ζ ln(kℓ)). Using that |g NL | 10 6 , we conclude that the loop corrections to f NL are small, at most at one per cent level, since 13 g NL P ζ ln(kℓ) 10 −2 . On the other hand, the one-loop corrections to τ NL of Eq. (27) can be large, significantly altering the tree level relation τ tree NL = (6f tree NL /5) 2 . Combining the results so far we find the general one-loop corrected single-source relation between τ NL and f NL given by The second term on the right hand side arises from the loop corrections and can be as large as 10 4 , irrespectively of the value of f NL . It may therefore give the dominant contribution to τ NL . In conclusion, the condition that leads to sizeable inhomogeneities for f NL , that is having a large g NL , also leads to breaking the single-source consistency relation between f NL and τ NL . On the other hand, even though the loop corrections can cause large deviations from the tree-level single-source result τ tree NL = (6f tree NL /5) 2 , the relation between τ NL and f NL still serves as a useful discriminator between single and multiple source models. Indeed as we will show in section IV, the result (27) is converted into a lower bound for τ 1−loop NL in the multiple source case, which also implies that (29) becomes a lower bound for τ NL . It would be interesting to generalize this analysis allowing for large loop corrections with higher than three-point vertices and see if the same conclusion holds even in that case.

D. Correlators of four-point functions
Considering correlators involving four-point functions measured in small patches, we can constrain higherorder n-point functions in the full sky, with n ≥ 5. In the single-source case, the only new parameter which we can independently constrain (even if going beyond 2-point correlators) is the quantity This is associated with the fourth order term in the standard expansion of the curvature perturbation In multiple-source models, more possibilities arise as we are going to discuss in section III; in this subsection, we focus on correlators that specifically depend on h NL . These involve measurements of g NL in small patches, and read 8 where in the right hand side of each correlator we only include the terms proportional to h NL and that are weighted by the largest coefficient, under the assumption that the actual values of f NL , and g NL saturate their observational bounds. For simplicity, we also drop the factor F defined in Eq. (12) which is of order unity and does not modify our discussion. Approximate 1σ constraints with WMAP for ζ, log P ζ , f NL , g NL are respectively 10 −5 , 10 −1 , 10 , 10 6 . Considering around one hundred patches of the sky, the measurement errors in each patch will be about ten times larger than these. On the other hand, when making correlations over the full sky, we should be able to recover the accuracy we loose when focussing on a small patch. As discussed at the end of subsection II B, in order to take into account effects associated with the fact that only a limited class of momentum space configurations fits into the small patches, we increase the expected error bars on the correlations above by one order of magnitude. The resulting error bars are consequently of order 10 2 , 10 6 , 10 8 and 10 13 , respectively for eqs. (32)- (35). Considering the constraints these terms could give on h NL , we see the correlation between ζ and g NL , eq. (32), would lead to a competitive constraint, of order |h NL | few × 10 11 .
Moreover, this correlation is sensitive to the sign of h NL . A comparable constraint can be obtained from the autocorrelation of g NL , eq. (35), although this quantity cannot determine the sign of h NL 9 . Note that the assumption |h NL | ≪ 10 5 |g NL |, which justifies neglecting the second order terms in eq. (19), breaks down when |g NL | ∼ 10 6 and the bound |h NL | few×10 11 is saturated. To work out the precise numerical factors in this limit one should therefore include the second-order terms and repeat the analysis leading to eqs. (32) - (35). This however should not affect the order of magnitude of the bound |h NL | 10 11 .
How useful is this constraint on h NL ? As we did in the previous sections with respect to g NL , we should compare them with constraints obtained from loop contributions to inflationary observables that are sensitive to h NL . There is a "dressed-vertex" loop diagram of the bispectrum [72] which gives f 1−loop NL ∼ h NL P ζ , which gives effectively a constraint of the same order as (36), unless there is a cancellation between loop terms and the tree level term contributing to f NL . One would need to include numerical factors to see which constraint is the most competitive.
If one renormalises the δN coefficients to avoid the dressed vertices, as suggested in [72], one should shift the derivative of N introducingÑ ′′ = N ′′ + N ′′′′ δφ 2 /2, which absorbs this term. On the other hand, in this way, the fine-tuning needed to compensate the loop corrections does not go away, it just becomes harder to spot. If the loops associated with dressed vertices do cancel, then the largest loop term involving h NL becomes g 1−loop NL ∼ f NL h NL P ζ , which implies |h NL | 10 14 , which is a much weaker constraint than we forecast cound be found by considering the autocorrelation of ∆ l g NL .

III. MULTIPLE SOURCE MODELS
The extension of the previous analysis to the multiple-source case suggests that our method is able to probe many of the parameters controlling higher-order point functions.
Taking ensemble averages of correlators among n-point and k-point parameters calculated in small patches, we are probing parameters associated to (n + k)-point functions evaluated in the entire sky. More precisely, for n ≥ k the correlators of local n-and k-point functions probe those parameters of the full-sky (n + k)-point function that contain at most n:th order derivatives of N , the number of e-foldings. This is the generalization of the observation made in [71] for the power spectrum. In the single-source limit, some of the (n + k)-point function parameters coincide with the squares of parameters associated with lower point functions (analogously to the tree-level consistency relation τ NL = (6/5f NL ) 2 holding between four and three point functions in the single source case). This is the reason for which we found that the two-point correlators of ∆ ℓ f NL 's calculated in different small patches are proportional to g 2 NL . On the other hand, in the multiple source case, our method in principle allows us to probe individually the various parameters controlling the higher-point functions.
In this section, we discuss this topic in more detail. Recently, a diagrammatic approach, based on the δN -formalism, has been adopted in [16,73] to determine the parameters characterizing the five and six point functions. While the three-point function is characterized by one parameter f NL and the four-point function by the two parameters τ NL and g NL , the five and six point functions are characterized by three and six parameters, respectively. In this section, we borrow the notation of [16] for defining parameters associated to 5 and 6-point functions (see [16] for more details, but notice that we do not adopt their normalization condition N c N c = 1).
Two point function: Three point function: Four point function: Five point function: Six point function: Notice that the parameters associated with 3 and 4 point functions carry the conventional numerical coefficients. The aim of this section is to show that ensemble averages of correlators between quantities involving two and three-point functions depend on parameters associated with five and six point functions. This suggests a method for probing the quantities appearing in eqs. (40) and (41). Generalizing the notation of the previous section, we can write in the multiple-field case (summations over repeated indices are understood) where the lower-case upper indices a, b, ... denote the different scalar fields involved. Analogous expansions hold for ζ, P ζ , τ NL and g NL . We assume that the metric in field space is flat, and that (We will neglect the slow-roll suppressed σ-dependence of H throughout this paper). Then we have Denoting and the same for ζ, ln P ζ , τ NL and g NL , we obtain the following results (the function F ≡ F (∆x, k hor , k ℓ ) is given in Eq. (12)): These ensemble averages are able to probe combinations of parameters controlling 5 and 6 point functions, as explained above. The previous combinations allow one to individually test each of the parameters characterizing the five-point function. Indeed, a combination of Eqs. (51), (52) and (55) lead to the following expressions for these parameters evaluated in the full sky In order to individually test the quantities associated with 6-pt functions one can consider correlations among τ NL and g NL measured in small patches -although the expressions for these quantities depend also on 7 and 8-pt functions. However the parameter g 6 , containing fifth order derivatives, can not be probed by considering the inhomogeneities of trispectrum to leading order in ∆ ℓ σ. In order to probe g (4) 6 , one should consider the inhomogeneities of the 5-point function or include the next-to-leading-order contributions in ∆ ℓ σ. However these in general are suppressed by a further factor of P ζ .
It is straightforward to check that, in the single-source limit, one recovers the results of the previous sections. This is due to consistency relations associating parameters in eqs (40) and (41) with parameters in (38) and (39). In the multiple source case, the parameters of eqs. (40) and (41) are not directly constrained by the observational bounds on f NL , τ NL and g NL . The inhomogeneities we are calculating could therefore be very large in certain classes of models -even at the level that might be detectable by re-analysing the WMAP data.

IV. WHEN IS THE INHOMOGENEITY OF NON-GAUSSIAN PARAMETERS LARGE?
In this section, we investigate in more detail when the autocorrelations of f NL and g NL measured in distinct patches can be large, and discuss the differences between single and multiple-source scenarios. We focus our attention on a representative, two-field case, in which the curvature perturbation is expanded as An expression of this form can be obtained by doing a rotation in field space to the direction, φ, generating the first-order contribution to the power spectrum of ζ; in this way only one linear term in scalar perturbations appear in Eq. (61). One could instead choose a different field basis in which the second order derivatives are diagonalised (i.e. N 12 = 0), as was chosen in [70], but this would be at the expense of introducing a second linear term which leads to more cumbersome expressions for the non-linearity parameters. We note that, in general, either choice of eliminating one term from the δN expression (61) will lead to a non-trivial field-metric, because the required field space rotation will depend on the scale we are focussing in. However since we can only observe a limited range of scales (around 5 e-foldings) this should typically only introduce small corrections, unless there happen to be large numerical factors which override the expected slow-roll suppression. We will not develop this interesting issue further, but simply adopt the expression (61) as a specific Ansatz for the curvature perturbations. Within this set-up, the non-linearity parameters controlling quantities up to the 6-pt function are: A. Inhomogeneous autocorrelation of fNL We have learned that, in a single-source set-up, the autocorrelation of f NL can be large provided that we saturate the observational limit for the full sky value of g NL . Let us investigate what happens instead in the multiple source case. In order to get an inhomogeneity for the parameter f NL large enough to be potentially observable with Planck, we need the combination in brackets of (54) to be of order O(10 12 ).
Given the present observational constraints, |f NL | 10 2 , τ NL 10 4 , |g NL | 10 6 [9][10][11][12][13][14], one may check that the only term which could be so large is the one proportional to τ (2) 6 , that provides Notice that the first term in the right hand side with g NL is the same leading contribution we found in the singlesource case, see Eq. (23). The second term represents a pure multiple-source contribution and its magnitude is independent of g NL . While the form of τ in Eq. (66) is specific for the two-source model (61) that we are considering as an example, the conclusion that the inhomogeneities of f NL are not uniquely determined by the magnitude of g NL holds for generic multiple-source scenarios, as we will show below.
The result (66) shows that large inhomogeneities of f NL are generated by cubic derivatives of N , which also generate loop corrections to τ NL . In the limit f 2 NL ≪ |g NL |, the dominant part of one-loop corrections to τ NL is given by which implies that This has to be compared with what we found in the single-source case, Eq. It is straightforward to show that this conclusion holds for generic multiple-source models, thus opening interesting possibilities for discriminating between single and multiple-source scenarios. Defining a unit vector u (k) a = δ ak , the Cauchy-Schwarz inequality for the inner product of the vectors u (k) b N a N abc and N c leads to where only the repeated lower case indices are summed over. The one-loop corrections to τ NL containing two three-point vertices are given by where we have neglected the logarithm associated with loop integrals. This represents the dominant contribution to the loop corrections, provided that |g NL | ≫ f 2 NL , and that there are no large loops associated with higher than three-point vertices. Similarly we can write Comparing the above expressions to the inequality (69), and making use of the inequality (54/25) 2 g 2 NL τ (2) 6 derived in [16], we find 54 25 In the limit |g NL | ≫ f 2 NL , using Eq. (54) we therefore find the general result The inequality is saturated in the single-source case (28) where the magnitude of the inhomogeneities is set by τ 1−loop NL . In multiple source scenarios, τ 1−loop NL can be greater than (∆ ℓ f NL ) 2 . The inhomogeneities of f NL could therefore provide an interesting new tool for discriminating between single and multiple source models. While the result (73) only applies in the limit |g NL | ≫ f 2 NL and assuming there are no large loops with higher than three-point vertices, it could nevertheless offer an intriguing new window for probing inflationary physics.
Similar information can be obtained by considering the structure of the trispectrum. The inequality (72) implies 54 25 which is again saturated for the single-source case. Under the assumptions stated above, there are no significant loop corrections to g NL , and τ NL is dominated by the loop correction as we are considering the limit τ tree NL ∼ f 2 NL ≪ |g NL |. We therefore see that both the inhomogeneities of f NL and the ratio of the two parameters g NL and τ NL have the opportunity to provide interesting new tools for distinguishing between single and multiple source models.

B. Inhomogeneities of gNL
We can apply the same procedure to analyze the autocorrelator of g NL . This quantity is given by Based on a measurement accuracy of |g NL | ∼ 10 6 (which will be improved by up to two orders of magnitude with Planck [14]) we require the term in brackets to be O(10 22 ) in order to be able to probe the inhomogeneity with WMAP data. With Planck data values as small as O(10 18 ) could be relevant.
The new non-linearity parameters which enter at this order are: Notice that only two new quantities are important for studying the new correlator of Eq. (75). These are N 1111 and N 1112 , which are two out of the five fourth derivatives of N . Unless one of these two quantities is large, the correlator (75) does not give competitive information compared to the correlators we considered in the previous sections. If one of the fourth derivatives is assumed to be very large, then the dominant term is simply τ 8 . This implies that we can constrain N 1111 /N 4 1 and N 1112 /N 4 1 to the level of 10 15 with WMAP and 10 13 with Planck.
On the other hand, |f 1−loop NL | ≃ |N 1111 |/N 4 1 P ζ 10 2 provides the constraint |N 1111 |/N 4 1 10 11 , so barring an accurate cancellation between this term and the tree level f tree NL we do not find interesting constraints from this contribution. However for N 1112 there is no 1-loop constraint from f NL , and instead the tightest constraint comes from τ 1−loop NL = N 1112 N 12 P ζ /N 6 1 10 4 . This does not constitute a dominant constraint compared to the expected sensitivity with the Planck satellite unless |N 12 |/N 2 1 is larger than unity. In conclusion, we see some similarities between studying the autocorrelation of ∆ ℓ f NL and here the autocorrelation of ∆ ℓ g NL . In both cases one can constrain two higher order derivatives, the one which is present in the single-source case, i.e. N 111 and N 1111 which are proportional to g NL and h NL respectively, but this constraint is not very competitive compared to the natural constraint we get from consider the loop contribution to lower order correlators unless there is a chance cancellation between the loop and the tree level terms. For the 'crossderivatives', N 112 and N 1112 , the loop constraints are weaker (especially in the case of N 112 ) and we can achieve a tighter probe by considering the inhomogeneity of non-Gaussian correlators.

V. CONCLUSIONS
In this paper, we investigated under which conditions inflationary parameters can be inhomogeneous to an observable level, focussing in particular on observables controlling local non-Gaussianity. We have demonstrated that if we subdivide the entire sky into a large number of small patches, the value of cosmological observables associated with the properties of the curvature fluctuation can differ from patch to patch. In particular, we have shown that correlators between n-point functions of curvature fluctuations as measured in one patch, and k-point functions as measured in another patch, depend on (n + k)−point functions in the entire sky. This implies that the expected degree of inhomogeneity in observable quantities can be quantified in a rather model independent manner. In interesting cases it results large enough to be measurable, since at present we have only weak constraints on parameters controlling higher-order point functions in the full sky. Consequently, inhomogeneities of non-Gaussian parameters can also be seen as feasible method for probing or constraining higher-point functions.
We have analyzed in detail the degree of inhomogeneity of local non-Gaussian observables, first in the singlesource case, in which only one scalar field contributes to the generation of primordial curvature perturbation (as in curvaton models), then in multiple-source set-ups. In the case of single-source models, we have shown that autocorrelators of f NL evaluated in different patches depend on the value of g NL in the full sky. If g NL turns out to be large enough to saturate its present day bound, we should expect variations of f NL of order one hundred from patch to patch, large enough to be observable. Autocorrelators of g NL , on the other hand, can be used to set constraints on h NL , the parameter characterizing 5-point functions: present day data are accurate enough to be able to set an upper bound |h NL | few × 10 11 . In the case of multiple source models, we have shown that correlators between parameters controlling n-point functions in different patches, with n ≤ 4, can be useful for probing individually the several new parameters that characterize five and six point functions.
We have pointed out interesting connections between the degree of inhomogeneities and loop corrections to non-Gaussian observables. Typically, models that lead to sizable inhomogeneities are also characterized by large loop corrections to inflationary observables. We have used this fact to determine consistency relations between quantities denoting respectively inhomogeneities and loop contributions to non-Gaussian observables. These consistency relations take the form of inequalities in multiple source models, that are saturated in the single-source limit: consequently, they offer new observational perspectives for distinguishing between single and multiple source set-ups.
The conclusion of our theoretical analysis is that a sizable degree of inhomogeneity in non-Gaussian observables is allowed by present day bounds on g NL and τ NL . We have then described observational prospects for probing inhomogeneities of non-Gaussian observables, discussing the accuracy we should expect for determining correlators of n-point functions measured in different patches of the skies. We discussed geometrical effects that reduce the accuracy, and conservatively took them into account in our estimates. We have also pointed out that alternative techniques based on needlets analysis of CMB data are particularly well suited for testing non-Gaussianity in selected regions of the sky. It would be very interesting to apply those techniques to study inhomogeneities.