Non-Gaussian CMB and LSS statistics beyond polyspectra

TL;DR

This paper explores advanced methods to detect primordial non-Gaussianity in cosmic microwave background and large-scale structure data, focusing on a generalized local ansatz that captures non-Gaussian features beyond traditional polyspectra.

Contribution

It introduces a generalized local ansatz for primordial curvature perturbations, enabling analysis of non-Gaussianity beyond polyspectra and linking it to inflationary potential shapes.

Findings

Derived probability distribution function classes for data analysis

Proposed strategies for reconstructing non-Gaussianity from observations

Connected non-Gaussian features to inflationary landscape potential

Abstract

Cosmic inflation may have led to non-Gaussian initial conditions that cannot be fully parametrised by 3- and/or 4-point functions. In this work, we discuss various strategies to search for primordial non-Gaussianity beyond polyspectra with the help of cosmological data. Our starting point is a generalised local ansatz for the primordial curvature perturbation $\zeta$ of the form $\zeta = \zeta_{\rm G} + \mathcal{F}_{\rm NG} (\zeta_{\rm G})$, where $\zeta_{\rm G}$ is a Gaussian random field and $\mathcal{F}_{\rm NG}$ is an arbitrary function parametrising non-Gaussianity that, in principle, could be reconstructed from data. Noteworthily, in the case of multi-field inflation, the function $\mathcal{F}_{\rm NG}$ can be shown to be determined by the shape of tomographic sections of the landscape potential responsible for driving inflation. We discuss how this generalised local ansatz leads to a probability distribution functional that may be used to extract information about inflation from current and future observations. In particular, we derive various classes of probability distribution functions suitable for the statistical analysis of the cosmic microwave background and large-scale structure.

Equations

$$\zeta (\mathbf{x})={\zeta }_{\mathrm{G}}(\mathbf{x})+{\mathcal{F}}_{\mathrm{NG}}\left({\zeta }_{\mathrm{G}}(\mathbf{x}),\nabla \right),$$

$$\zeta_{\bm{k}}=\zeta_{\bm{k}}^{G}+(2\pi)^{3}\sum_{n=2}\frac{1}{n!}\int_{{\bm{k}}_{1}}...\int_{{\bm{k}}_{n}}\delta^{(3)}\Big{(}{\bm{k}}-\sum_{i=1}^{n}{\bm{k}}_{i}\Big{)}F_{n}\left({\bm{k}}_{1},...,{\bm{k}}_{n}\right)\zeta^{G}_{{\bm{k}}_{1}}...\zeta^{G}_{{\bm{k}}_{n}},$$

$$F_2({\mathbf{k}}_1,{\mathbf{k}}_2)=\frac{1}{2}(\eta -ϵ)+\frac{ϵ}{2}\frac{k_1+k_2}{|{\mathbf{k}}_1+{\mathbf{k}}_2|}+2ϵ\frac{k_1^2+k_2^2}{|{\mathbf{k}}_1+{\mathbf{k}}_2|(k_1+k_2+|{\mathbf{k}}_1+{\mathbf{k}}_2|)},$$

$$F_n\propto g^{n-1}.$$

$${\mathcal{F}}_{\mathrm{NG}}(\zeta )\propto \frac{\partial }{\partial \zeta }\Delta V\left(\psi (\zeta )\right),\mathrm{with}\psi (\zeta )\equiv \frac{H}{2\pi A_s^{1/2}}\zeta ,$$

$$\zeta (\mathbf{x})={\zeta }_{\mathrm{G}}(\mathbf{x})+{\mathcal{F}}_{\mathrm{NG}}\left[{\zeta }_{\mathrm{G}}\right](\mathbf{x}),\mathrm{with}{\mathcal{F}}_{\mathrm{NG}}\left[{\zeta }_{\mathrm{G}}\right](\mathbf{x})={\int }_{\mathbf{k}}{\int }_{\mathbf{y}}{e}^{i\mathbf{k}\cdot (\mathbf{x}-\mathbf{y})}F\left({\zeta }_{\mathrm{G}}(\mathbf{y})\right),$$

$$⟨\zeta ({\mathbf{x}}_1)\cdots \zeta ({\mathbf{x}}_n)⟩=\int D\zeta \mathrm{P}[\zeta ]\zeta ({\mathbf{x}}_1)\cdots \zeta ({\mathbf{x}}_n),$$

$${\mathrm{P}}_{\mathrm{G}}[{\zeta }_{\mathrm{G}}]=\mathcal{N}\exp \left[-\frac{1}{2}{\int }_{\mathbf{x}}{\int }_{\mathbf{y}}{\zeta }_{\mathrm{G}}(\mathbf{x}){\Sigma }^{-1}(\mathbf{x},\mathbf{y}){\zeta }_{\mathrm{G}}(\mathbf{y})\right]=\mathcal{N}\exp \left[-\frac{1}{2}{\int }_{\mathbf{k}}\frac{{\zeta }_{\mathrm{G}}(\mathbf{k}){\zeta }_{\mathrm{G}}(-\mathbf{k})}{P_{\zeta }(k)}\right],$$

$${\Sigma }^{-1}(\mathbf{x},\mathbf{y})\equiv {\int }_{\mathbf{k}}\frac{e^{i\mathbf{k}\cdot (\mathbf{x}-\mathbf{y})}}{P_{\zeta }(k)}.$$

$$\det \left(\frac{\delta \zeta (\mathbf{y})}{\delta {\zeta }_{\mathrm{G}}(\mathbf{x})}\right)=\exp \left[\mathrm{tr}\left(\ln \left(\frac{\delta \zeta (\mathbf{y})}{\delta {\zeta }_{\mathrm{G}}(\mathbf{x})}\right)\right)\right]$$

$$\mathrm{P}[\zeta ]=\mathcal{N}\exp \left[-\frac{1}{2}{\int }_{\mathbf{x}}{\int }_{\mathbf{y}}G[\zeta ](\mathbf{x}){\Sigma }^{-1}(\mathbf{x},\mathbf{y})G[\zeta ](\mathbf{y})-{\int }_{\mathbf{x}}{\int }_{\mathbf{k}}\ln \left\{1+\frac{\mathrm{d}F}{\mathrm{d}\zeta }({\zeta }_{\mathrm{G}}(\mathbf{x}))\right\}\right].$$

$$\mathrm{P}[\zeta ]={\mathrm{P}}_{\mathrm{G}}[\zeta ]\times \exp \left[{\int }_{\mathbf{x}}{\int }_{\mathbf{y}}\zeta (\mathbf{x}){\Sigma }^{-1}(\mathbf{x},\mathbf{y}){\mathcal{F}}_{\mathrm{NG}}(\zeta (\mathbf{y}))-{\int }_{\mathbf{x}}\frac{\delta }{\delta \zeta (\mathbf{x})}{\int }_{\mathbf{k}}{\int }_{\mathbf{y}}{e}^{i\mathbf{k}\cdot (\mathbf{x}-\mathbf{y})}{\mathcal{F}}_{\mathrm{NG}}(\zeta (\mathbf{y}))\right].$$

$$\mathrm{P}[\zeta ]={\mathrm{P}}_{\mathrm{G}}[\zeta ]\times \exp \left[{\int }_{\mathbf{x}}{\int }_{\mathbf{k}}{\int }_{\mathbf{y}}{e}^{i\mathbf{k}\cdot (\mathbf{x}-\mathbf{y})}\left(\frac{\zeta (\mathbf{x})}{P_{\zeta }(k)}-\frac{\delta }{\delta \zeta (\mathbf{x})}\right){\mathcal{F}}_{\mathrm{NG}}[\zeta ](\mathbf{y})\right].$$

$${\mathrm{P}}_F[\zeta ]$$

$${\zeta }_W(\mathbf{x})\equiv [W\star \zeta ](\mathbf{x})={\int }_{\mathbf{k}}{\int }_{\mathbf{y}}{e}^{i\mathbf{k}\cdot (\mathbf{x}-\mathbf{y})}W(k)\zeta (\mathbf{y}).$$

$$\begin{array}{rl}{\mathrm{P}}_W[{\zeta }_W]=& {\mathcal{N}}_W\exp \left[-\frac{1}{2}{\int }_{\mathbf{k}}\frac{{\zeta }_W(\mathbf{k}){\zeta }_W(-\mathbf{k})}{P_W(k)}\right]\\ & \times \left[1-{\int }_{\mathbf{y}}{\int }_{\mathbf{k}}{e}^{-i\mathbf{k}\cdot \mathbf{y}}W(k)\left(\frac{\delta }{\delta {\zeta }_W(-\mathbf{k})}-\frac{{\zeta }_W(\mathbf{k})}{P_W(k)}\right){\mathcal{F}}_{\mathrm{NG}}\left[{\int }_q\frac{e^{-i\mathbf{q}\cdot \mathbf{y}}{\zeta }_W(\mathbf{q})}{W(q)}\right](\mathbf{y})\right],\end{array}$$

$$W(k)=\{\begin{array}{ll}1& \text{if}k\in (k_L,k_S)\\ 0& \text{if}k\notin (k_L,k_S),\end{array}$$

$${\mathrm{P}}_W[\zeta ]=\int D{\zeta }_{|\mathbf{k}|\notin (k_L,k_S)}\mathrm{P}[\zeta ].$$

$$\begin{array}{rl}{\int }_{\mathbf{k}}{e}^{-i\mathbf{k}\cdot \mathbf{y}}\left(\frac{\delta }{\delta \zeta (-\mathbf{k})}-\frac{\zeta (\mathbf{k})}{P_{\zeta }(k)}\right){\mathcal{F}}_{\mathrm{NG}}[\zeta ](\mathbf{y})& ={\int }_{|\mathbf{k}|\in (k_L,k_S)}{e}^{-i\mathbf{k}\cdot \mathbf{y}}\left(\frac{\delta }{\delta \zeta (-\mathbf{k})}-\frac{\zeta (\mathbf{k})}{P_{\zeta }(k)}\right){\mathcal{F}}_{\mathrm{NG}}[\zeta ](\mathbf{y})\\ & +{\int }_{|\mathbf{k}|\notin (k_L,k_S)}{e}^{-i\mathbf{k}\cdot \mathbf{y}}\left(\frac{\delta }{\delta \zeta (-\mathbf{k})}-\frac{\zeta (\mathbf{k})}{P_{\zeta }(k)}\right){\mathcal{F}}_{\mathrm{NG}}[\zeta ](\mathbf{y}).\end{array}$$

$$\left(\frac{\delta }{\delta \zeta (-\mathbf{k})}-\frac{\zeta (\mathbf{k})}{P_{\zeta }(k)}\right)\int D{\zeta }_{|\mathbf{k}|\notin (k_L,k_S)}\exp \left[-\frac{1}{2}{\int }_{\mathbf{k}}\frac{\zeta (\mathbf{k})\zeta (-\mathbf{k})}{P_{\zeta }(k)}\right]{\int }_{\mathbf{y}}{\int }_{|\mathbf{k}|\in (k_L,k_S)}{e}^{-i\mathbf{k}\cdot \mathbf{y}}{\mathcal{F}}_{\mathrm{NG}}[\zeta ](\mathbf{y}),$$

$$\bar{F}({\zeta }_W(\mathbf{x}))\equiv \int D{\zeta }_{k\notin (k_L,k_S)}\exp \left[-\frac{1}{2}{\int }_{|\mathbf{k}|\notin (k_L,k_S)}\frac{\zeta (\mathbf{k})\zeta (-\mathbf{k})}{P_{\zeta }(k)}\right]F(\zeta (\mathbf{x}))={\int }_{-\infty }^{\infty }\mathrm{d}\bar{\zeta }\frac{e^{-\frac{({\zeta }_W(\mathbf{x})-\bar{\zeta }{)}^2}{2{\sigma }_{\mathrm{out}}^2}}}{\sqrt{2\pi {\sigma }_{\mathrm{out}}^2}}F(\bar{\zeta }),$$

$${\bar{\mathcal{F}}}_{\mathrm{NG}}[\zeta ](\mathbf{x})={\int }_{\mathbf{y}}{\int }_{|\mathbf{k}|\in (k_L,k_S)}{e}^{i\mathbf{k}\cdot (\mathbf{x}-\mathbf{y})}\bar{F}(\zeta (\mathbf{y})),$$

$${\sigma }_{\mathrm{out}}^2={\int }_{|\mathbf{k}|\notin (k_L,k_S)}{P}_{\zeta }(k),$$

$${\mathrm{P}}_{{\mathcal{F}}_{\mathrm{NG}},W}[\zeta ]={\mathrm{P}}_{{\bar{\mathcal{F}}}_{\mathrm{NG}}}[{\zeta }_{|\mathbf{k}|\in (k_L,k_S)}],$$

$$\bar{F}(\zeta )=\exp \left[\frac{{\sigma }_{\mathrm{out}}^2}{2}\frac{{\partial }^2}{\partial {\zeta }^2}\right]F(\zeta ).$$

$${\zeta }_W(\mathbf{x})={\zeta }_W^G(\mathbf{x})+{\bar{\mathcal{F}}}_{\mathrm{NG}}\left[{\zeta }_W^G\right](\mathbf{x}),$$

$$\langle\zeta({{\bm{k}}}_{1})\cdots\zeta({{\bm{k}}}_{n})\rangle=\frac{\delta^{n}Z[J]}{(i\delta J(-{{\bm{k}}}_{1}))\cdots(i\delta J(-{{\bm{k}}}_{n}))}\bigg{|}_{J=0},$$

$$Z[J]=\int D\zeta \mathrm{P}[\zeta ]e^{i{\int }_{\mathbf{k}}\zeta (\mathbf{k})J(-\mathbf{k})}.$$

$$\begin{array}{rl}Z[& J]=\exp \left[-\frac{1}{2}{\int }_{\mathbf{k}}J(\mathbf{k})J(-\mathbf{k})P_{\zeta }(k)\right]\times \\ & \left(1-{\int }_{\mathbf{x}}\frac{{\int }_{\mathbf{k}}{e}^{i\mathbf{k}\cdot \mathbf{x}}J(-\mathbf{k})}{{\int }_{\mathbf{k}}{e}^{i\mathbf{k}\cdot \mathbf{x}}J(-\mathbf{k})P_{\zeta }(k)}{\int }_{\bar{\zeta }}\frac{\exp \left[-\frac{{\left(\bar{\zeta }-i{\int }_{\mathbf{k}}{e}^{i\mathbf{k}\cdot \mathbf{x}}J(-\mathbf{k})P_{\zeta }(k)\right)}^2}{2{\sigma }_{\zeta }^2}\right]}{\sqrt{2\pi }{\sigma }_{\zeta }}\left({\sigma }_{\zeta }^2\frac{\partial }{\partial \bar{\zeta }}-\bar{\zeta }\right)F(\bar{\zeta })\right).\end{array}$$

$$⟨{\zeta }_W({\mathbf{x}}_1)\cdots {\zeta }_W({\mathbf{x}}_n){⟩}_c=f_{n-1}\sum _{i=1}^n{\int }_{\mathbf{x}}\frac{{\int }_{\mathbf{k}}W(k)e^{i\mathbf{k}\cdot ({\mathbf{x}}_i-\mathbf{x})}}{{\int }_{\mathbf{k}}{e}^{i\mathbf{k}\cdot ({\mathbf{x}}_i-\mathbf{x})}W(k)P_{\zeta }(k)}\left(\prod _{i=1}^n{\int }_{\mathbf{k}}W(k)P_{\zeta }(k)e^{i\mathbf{k}\cdot ({\mathbf{x}}_i-\mathbf{x})}\right),$$

Full text

1. 1 Introduction

2. 1.1 Main results

3. 1.2 Structure of the paper

4. 2 $\zeta$ Statistics

5. 2.1 Probability distribution functional: “bare” theory

6. 2.2 Probability distribution functional: running and renormalising

7. 2.3 Partition function and $n$-point correlators

8. 2.4 Fixed-point probability distribution functions

9. 2.4.1 1-Point probability distribution function

10. 2.4.2 2-Point probability distribution function

11. 2.5 Linear transfer functions and their limitations

12. 3 $\Theta$ Statistics

13. 3.1 1-Point temperature distribution function

14. 3.2 2-Point temperature distribution function

15. 4 $\delta$ Statistics

16. 4.1 Halo mass function

17. 4.2 Scale dependent halo bias

18. 5 Concluding remarks

19. A Generalised ansatz from quantum fluctuations during inflation

20. A.1 A concrete example: multi-field inflation

21. B The partition function

22. C Details of the fixed-point PDFs

23. C.1 2-Point PDF: general case

24. C.2 2-Point PDF: map to the sphere $S^{2}$

1 Introduction

Our observable universe is consistent with an extremely simple set of initial conditions. For all practical purposes, the observed cosmological inhomogeneities were seeded by a primordial curvature fluctuation $\zeta$ distributed according to a Gaussian profile parametrised by an almost scale invariant power spectrum [1, 2, 3, 4, 5]. The confirmation of this state of affairs by future surveys would reinforce our confidence in the single-field slow-roll inflation paradigm, that is, the idea that $\zeta$ was the consequence of quantum perturbations of a single scalar fluid (the inflaton) that evolved adiabatically during inflation [6, 7, 8, 9, 10]. Tiny deviations from Gaussianity, due to small nonlinear self-interactions affecting $\zeta$, are known to emerge in single-field inflation but these are predicted to be too small to be observed in the near future. On the other hand, large non-Gaussianity (within current bounds) may emerge in models of inflation beyond the canonical single-field paradigm, resulting from possible nongravitational self-interactions of $\zeta$ and/or interactions between $\zeta$ and other degrees of freedom. The observation of non-Gaussianity (NG) would therefore offer a unique opportunity to characterise the class of interactions that affected $\zeta$ during inflation, allowing us to pin down certain fundamental aspects about the period of inflation and, consequently, have a glimpse on the structure of the ultra-violet (UV) framework where it is realised.

While current cosmic microwave background (CMB) observations show no evidence of primordial non-Gaussianity, future large-scale structure (LSS) surveys, such as Lsst [11], Euclid [12], Spherex [13] and Ska [14], promise to revitalise its search. The proliferation of modes due to the three-dimensional volume probed by galaxy surveys is expected to yield constraints on primordial non-Gaussianity that might not only complement current bounds from CMB observations, but even surpass them. Among the most prominent effects of non-Gaussianity on the matter distribution is the celebrated observation that a nonzero skewness of the probability distribution function of $\zeta$ leads to an enhanced abundance of collapsed structures and a scale dependent correction in the halo bias [15], a result which has brought LSS surveys in the front line of the search for non-Gaussianity. Furthermore, there has been an intense effort to study how UV physics can show up in the matter power spectrum and bispectrum [16, 17, 18, 19].

Given that inflation might not be unique111Other scenarios, such as the ekpyrotic [20, 21] and bouncing universe [22, 23], have been proposed as models able to reproduce nearly Gaussian initial conditions. in explaining the observed inhomogeneities [24], in order to confront future data with theory, we may adopt an agnostic perspective about the details involved in the description of the pre-Big-Bang dynamics. In that case, we should agree that the main outcome from inflation, or any other framework intending to explain the initial conditions of our universe, consists of a relation giving us back the profile of $\zeta({\bm{x}})$ written in terms of a purely Gaussian random field $\zeta_{\rm G}({\bm{x}})$. Such a relation must be of the form

[TABLE]

where ${\cal F}_{\rm NG}$ represents a nonlinear function of the field $\zeta_{\rm G}({\bm{x}})$ and spatial gradients $\nabla$ acting on it. In Fourier space, the previous relation may be reexpressed as the following expansion in powers of $\zeta_{\rm G}$, starting at quadratic order:

[TABLE]

where $\int_{\bm{k}}\equiv(2\pi)^{-3}\int{\rm d}^{3}k$, and where $F_{n}\left({\bm{k}}_{1},...,{\bm{k}}_{n}\right)$ are functions of the momenta, symmetric under their permutations. This relation is sufficiently general to describe any form of primordial non-Gaussianity, and may be formally obtained in a generic manner from a quantum field theoretical framework (see App. A for details). The $F_{n}$ functions parametrise the non-Gaussian deviations generated by nonlinear interactions to which $\zeta$ were subject and, in the case of inflation, may be deduced from a particular model by studying the evolution of $\zeta_{\bm{k}}$ from sub-horizon scales up until the end of inflation (e.g., using the in-in formalism). In fact, any $n$-point correlation function for $\zeta$ may be computed out of (1.2). For instance, at tree-level, the bispectrum parametrising the amplitude of the 3-point function is found to be given by $B({\bm{k}}_{1},{\bm{k}}_{2},{\bm{k}}_{3})=\left[P_{\zeta}(k_{1})P_{\zeta}(k_{2})F_{2}({\bm{k}}_{1},{\bm{k}}_{2})+\textrm{perm}\right]$, where $P_{\zeta}$ is the power spectrum of the Gaussian field $\zeta^{G}_{{\bm{k}}}$. In the particular case of single-field slow-roll inflation, up to first order in the slow-roll parameters, the bispectrum is recovered as long as $F_{2}$ is given by (up to terms that vanish upon imposing momentum conservation)

[TABLE]

where $\epsilon$ and $\eta$ are the usual slow-roll parameters describing the steady evolution of the FLRW background during inflation. The effective field theory (EFT) of inflation approach [25] to study models beyond the canonical single-field paradigm will also yield a specific form of $F_{2}({\bm{k}}_{1},{\bm{k}}_{2})$, in which the sound speed of $\zeta$ plays an important role. Moreover, the well-known local ansatz [26, 27, 28, 29, 30, 31, 32, 33], related to the presence of multi-field dynamics, corresponds to another particular instance of this relation, where $F_{2}=\frac{6}{5}f_{\rm NL}^{\rm local}$.

One could thus take upon the challenge of directly reconstructing the form of ${\cal F}_{\rm NG}$ in Eq. (1.1) —or equivalently, the functions $F_{n}$ appearing in (1.2)— out of cosmological data. This would constitute a bottom-up approach to determine the properties of the model that gave origin to the initial conditions. To guide such a reconstruction, one could consider restricting the functions $F_{n}$ according to certain rules dictated by the symmetries of the alleged bulk model that led to (1.2) at the end of the pre-Big Bang period. For instance, scale invariance of the spectra is equivalent to the invariance of $F_{n}$ under the simultaneous rescaling of all the momenta: $F_{n}(\lambda{\bm{k}}_{1},...,\lambda{\bm{k}}_{n})=F_{n}({\bm{k}}_{1},...,{\bm{k}}_{n})$. Furthermore, the validity of soft theorems [29, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43] (under certain circumstances) would require some relations among $F_{n}$ of different order in the limit where one or more of the momenta vanish.

An objection to this program (the direct reconstruction of ${\cal F}_{\rm NG}$) is that perturbation theory applied to the study of the evolution of $\zeta$ implies that the expansion involved in the writing of Eq. (1.2) is hierarchical. That is, given a small dimensionless coupling constant $g$ parametrising the self-interactions experienced by $\zeta$ during inflation, the $F_{n}$ functions are naturally expected to satisfy

[TABLE]

For instance, in the case of single-field slow-roll inflation $g$ happens to be of order $\epsilon$ and $\eta$ and hence, non-Gaussianity is expected to be well-parametrised by the bispectrum. In noncanonical single-field models described by the EFT of inflation, where the $\zeta$ fluctuations propagate with a reduced sound speed $c_{s}<1$, the coupling $g$ is enhanced by a factor $c_{s}^{-2}$, but in order to trust perturbation theory, one still requires that it stay sufficiently suppressed. Based on this argument, we could say that most efforts to characterise non-Gaussianity so far have focused on a truncated version of (1.2), where only $F_{2}({\bm{k}}_{1},{\bm{k}}_{2})$ and $F_{3}({\bm{k}}_{1},{\bm{k}}_{2},{\bm{k}}_{3})$ (which at tree-level give the bispectrum and trisprectrum) are taken under consideration, with the Planck data implying weak constraints on the form of $F_{2}({\bm{k}}_{1},{\bm{k}}_{2})$, with the help of the so-called local, equilateral, orthogonal and folded templates [4].

In the present article, we wish to argue in favour of reconstructing the full function ${\cal F}_{\rm NG}$ of Eq. (1.1) from CMB and LSS observations without necessarily assuming a hierarchical dependence of the functions $F_{n}$ on a given coupling constant $g$. We posit that the search for non-Gaussianity focused on low $n$-point correlation functions may miss the existence of richer types of statistics [44, 45, 46, 47, 48, 49, 50] that we may be unable to predict by following standard perturbation theory techniques. For example, in App. A we show that in certain classes of multi-field models, the function ${\cal F}_{\rm NG}$ is found to be proportional to the gradient of the landscape potential $\Delta V(\psi)$ controlling the dynamics of the isocurvature field $\psi$. More to the point, we show that

[TABLE]

where $H$ is the Hubble expansion rate during inflation and $A_{s}$ is the amplitude of the power spectrum of $\zeta$ (constrained by Planck to be $10^{9}A_{s}=2.105\pm 0.030$ at the CMB pivot scale $k=0.05$ Mpc*-1* [2]). Hence, the shape of ${\cal F}_{\rm NG}$ provides a tomographic view of the landscape potential away from the inflationary trajectory. Certainly, one cannot discard the possibility of having potentials $\Delta V(\psi)$ characterised by a rich structure (i.e. with consecutive minima separated by field distances smaller than the Hubble expansion rate during inflation) in which case the bispectrum would not constitute an efficient tool to parametrise non-Gaussianity [48, 49]. A simple example of this, related to the presence of an axionic isocurvature direction, is ${\cal F}_{\rm NG}(\zeta_{\rm G})\propto\sin(\zeta_{\rm G}/f_{\zeta})$ with $f_{\zeta}^{2}<\sigma_{\zeta}^{2}$, where $\sigma_{\zeta}^{2}$ is the variance of the Gaussian field $\zeta_{\rm G}$. In this case, the function ${\cal F}_{\rm NG}(\zeta_{\rm G})$ remains bounded, but varies vigorously within the relevant range determined by $\sigma_{\zeta}$.

Despite of having a well motivated construction justifying the need of studying non-Gaussianity beyond polyspectra, we wish to keep the discussion as general as possible and, thus, in the main part of this work we will not refer to any specific origin of this parametrisation. Instead, we will focus our attention on the consequences —for various cosmological observables— implied by having initial conditions parametrised by a function ${\cal F}_{\rm NG}(\zeta)$ not necessarily respecting a hierarchical structure. In order to make the discussion more tractable, we will focus on the specific case of local NG, for which ${\cal F}_{\rm NG}$ in Eq. (1.1) does not involve gradients and hence the coefficients $F_{n}$ in the expansion (1.2) are independent of momenta222See Ref. [51] for a discussion on how the usual local ansatz involving $f_{\rm NL}$ (where the momenta are disregarded) can in principle be observationally distinguished from an ansatz where the function $F_{2}$ satisfies the consistency relation. The methods used in [51] should be valid to analyse our proposed generalised version of the local ansatz.. This implies a simplified version of (1.1) given by $\zeta({\bm{x}})=\zeta_{\rm G}({\bm{x}})+{\cal F}_{\rm NG}\left(\zeta_{\rm G}({\bm{x}})\right)$. However, since in practice observations are restricted to a finite range of wavelengths, we will consider a version of (1.1) recovered from (1.2) but with every momenta restricted to a finite range (including those appearing in the integrals). This relation is given by

[TABLE]

where $\int_{\bm{y}}\equiv\int{\rm d}^{3}y$. In the previous expression, $\int_{{\bm{k}}}\equiv(2\pi)^{-3}\int{\rm d}^{3}k$ is restricted to a given range of scales specified by convenience333For now, we do not need to specify the range of momenta and the reader may take the range as the entire momentum space, in which case ${\cal F}_{\rm NG}=F$. In Sec. 2.2, we will consider the problem of reducing the range of momenta (in a Wilsonian manner) and explore how the function $F$ appearing in Eq. (1.5) runs with the IR and UV cutoffs employed to restrict the momenta. (with $\zeta({\bm{x}})$ and $\zeta_{\rm G}({\bm{x}})$ also restricted in that manner). In this version of the ansatz, $F\left(\zeta_{\rm G}\right)$ is the function of the Gaussian random field $\zeta_{\rm G}$ that determines the values of the $F_{n}$ coefficients in (1.2) as $F_{n}=\partial_{\zeta}^{n}F|_{\zeta=0}$.

While similar approaches to study non-Gaussianity have been developed (for instance, see Refs. [44, 46], where the main focus is the analysis of non-Gaussianities arising from preheating), our purpose in this work is to provide a description of local non-Gaussianity in terms of the curvature perturbation field $\zeta$ in a way that explicitly reflects how the scales involved in the computation affect the result, and therefore, on how the primordial deviation from Gaussianity ${\cal F}_{\rm NG}$ can be related to predictions and observables in the CMB and LSS by adjusting the set of scales to those involved in the measurement at hand.

1.1 Main results

The object that determines the statistics of fluctuations, either in the temperature or in the density field, is the probability distribution function (PDF) from which one can compute moments, cumulative functions and several other objects such as mass functions, bias factors, etc. Since both the temperature and density fluctuations are sourced by curvature perturbations, the parametrisation (1.5) of $\zeta$ in terms of Gaussian random fields allows us to construct various classes of such functions characterising the statistics444There are several nonprimordial sources of non-Gaussianity, especially in late-time fields. Here, we always refer to its primordial component; see Sec. 2.5 for comments on the limitations posed by secondary effects. of CMB and LSS, hence, enabling us to use their distributions as a probe of the primordial universe. In this work we derive several classes of PDFs for the curvature perturbation $\zeta$, the CMB temperature fluctuation $\Theta\equiv\Delta T/T$ and the matter overdensity $\delta\equiv\delta\rho/\rho$, all following from the non-Gaussian ansatz (1.5) by using appropriate transfer functions.

Specifically, we achieve the following goals:

• $\star$

We derive the full probability distribution functional implied by the generalised ansatz (1.5). The result is shown in Eq. (2.10) and it explicitly depends on the function $F$ (which at the same time, is related to the landscape potential in the case where the primordial non-Gaussianity is produced during inflation).

• $\star$

We use this functional to derive various statistical tests for CMB and LSS, among them, 1- and 2-point PDFs, written in Eqs. (2.31) and (2.42). We also write the corresponding PDFs for the CMB temperature anisotropies in Eqs. (3.9) and (3.16).

• $\star$

We derive the NG halo mass function expressed in terms of the function $F$. This result is shown in Eqs. (4.8)-(4.10). We further comment on how statistical estimators probing the PDF can be used in LSS datasets in connection to primordial NG via Eq. (4.13).

• $\star$

We derive a generalised version of the linear halo bias resulting from the non-Gaussian deformation ${\cal F}_{\rm NG}$. The result is written in Eq. (4.27), where we show that the scale dependent bias is sensitive to the full function $F$ and not to any specific truncated version of it. In other words, all nonlinearity parameters appear in a degenerate manner.

To summarise, we wish to point out that since non-Gaussianity has not shown up in low $n$-point functions and given that its presence has solid foundations in our understanding of UV physics and quantum field theory/gravity, extending our study towards other estimators is part of the next step. The PDF is the object that fully encodes the statistical information of fluctuations enclosing all possible cases and as such it may serve as a central tool in the search for non-Gaussianity.

1.2 Structure of the paper

The article is organised as follows: In Sec. 2, we write down the probability functional describing the distribution of the curvature fluctuation $\zeta$, its Fourier dual —the partition function— as well as the 1-point and 2-point PDFs. We discuss several aspects regarding the use of window functions related to renormalisation, survey sky coverage and secondary NG effects. In Sec. 3, we derive the temperature distribution and focus on the 2-point PDF commenting on its use with CMB data. In Sec. 4, we derive the halo mass function and the halo bias for arbitrary local NG. We conclude in Sec. 5, while in the Appendices we present various details regarding the NG ansatz and the functionals emanating thereof presented in the main text.

2 $\zeta$ Statistics

In this section we derive the probability distribution functional (and other classes of distributions) for the primordial curvature perturbation $\zeta$. Our starting point is the generalised local ansatz shown in Eq. (1.5). The resulting PDF’s will set the stage for later sections, where we explore the consequences their non-Gaussian deviations may have on both CMB and LSS datasets.

2.1 Probability distribution functional: “bare” theory

The first task we can perform, and the most readily available to us thus far, is to attempt to use Eq. (1.5) to its full potential and derive directly the complete functional distribution that governs the $\zeta$ statistics: given that we know how a Gaussian random field is distributed, we may simply perform a change of variables to obtain the probability distribution functional ${\rm P}[\zeta]$ for $\zeta({\bm{x}})$ such that

[TABLE]

where $D\zeta$ is the functional measure of the PDF. Such an approach was used in [52] for the study of certain non-Gaussian initial conditions from inflation, in [53] in the context of nonlinear evolution of the correlation function of matter overdensities, while in [54, 55] it was shown how initial conditions with nonzero skewness/kurtosis can be described by the change of variables from a Gaussian PDF. Here, we extend this technique to the case of the generic NG field of Eq. (1.5).

The resulting distribution must account, to leading order in the perturbation ${\cal F}_{\rm NG}$, for every conceivable correlation function that may be constructed from the field $\zeta({\bm{x}})$ and for every expectation value of function(al)s of $\zeta({\bm{x}})$. To start with, the Gaussian random field $\zeta_{\rm G}$ is drawn from the following functional distribution:

[TABLE]

where $\mathcal{N}$ is an overall normalisation constant, while $\Sigma^{-1}({\bm{x}},{\bm{y}})$ and $P_{\zeta}(k)$ are the covariance matrix and the power spectrum respectively, related as

[TABLE]

In implementing the transformation $\zeta({\bm{x}})=\zeta_{\rm G}({\bm{x}})+{\cal F}_{\rm NG}[\zeta_{\rm G}]({\bm{x}})$, there are two contributions that emerge: one is composed by the terms that come out of the Taylor expansion of the Gaussian distribution by regarding the perturbation ${\cal F}_{\rm NG}$ as small, and the other arises from the functional determinant of the transformation. The latter is given by

[TABLE]

As long as the determinant exists (concretely, ${\rm d}F/{\rm d}\zeta>-1$) and is nonzero, we may, in principle, find an inverse to the relation $\zeta=\zeta_{\rm G}+{\cal F}_{\rm NG}[\zeta_{\rm G}]$ and denote it with $\zeta_{\rm G}({\bm{x}})=G[\zeta]({\bm{x}})$. Then we may change variables from the PDF associated to $\zeta_{\rm G}$, to find

[TABLE]

This is an exact result. However, if $F$ is small in comparison to the typical scales of the background theory on which the fluctuation field $\zeta$ lies, and so is ${\rm d}F/{\rm d}\zeta$, we can approximate the logarithm in the exponential with the first term in its power series expansion, and furthermore, we may approximate the inverse mapping by $G[\zeta]\approx\zeta-{\cal F}_{\rm NG}[\zeta]$. This yields a non-Gaussian exponential correction to the PDF, the exponent of which, to first order in ${\cal F}_{\rm NG}$, reads

[TABLE]

Using the definition of the covariance (2.3), we may finally write the result as

[TABLE]

This functional will serve as the guiding principle for all subsequent results. Note that even in a first order approximation, the probability distribution is always positive. However, to make computations tractable, we find it convenient to also perform a power series expansion of the non-Gaussian exponential factor. Then, to first order in ${\cal F}_{\rm NG}$, we find that

[TABLE]

The functional distribution (2.8) has support at all the scales where the underlying theory does, or at least, at the scales where the corresponding EFT is presumed to hold true. However, observable quantities do not typically involve all of the scales, and therefore it may be that the “bare” departure from Gaussianity ${\cal F}_{\rm NG}$ is not the most adequate quantity to describe them. We thus now turn to the discussion of using window functions and how to integrate scales out.

2.2 Probability distribution functional: running and renormalising

When making predictions, any EFT will force us to recognise certain scales at which the theory is no longer well-suited to describe physical observables. This typically implies a high-energy scale, where the theory has to be cut off. Thus, we will set $k_{\rm UV}$ as the maximum possible wavenumber the mode expansion of the curvature perturbation can have. Similarly, while it is not always introduced, one can make the same assertion with the very long wavelengths. As much as our theory may have predictions concerning phenomena happening beyond the present Hubble radius, they are currently unobservable. Therefore, in establishing predictions for observable quantities it seems natural to integrate out those scales, so that they are properly incorporated into the final, effective result. Because of this, we will take a conservative attitude and also define an infrared cutoff $k_{\rm IR}$, which can be thought of as the inverse of the current Hubble radius, thus bounding the domain of the theories we will be studying to $k\in(k_{\rm IR},k_{\rm UV})$.

Moreover, in realistic situations the experiment at hand may not allow us to access every value for the momentum scale $k$ evenly. In those cases we may wish to introduce a window function $W(k)$ to filter our results and give more weight to some scales. Accordingly, one would be interested in the statistics of the filtered field

[TABLE]

To derive the probability distribution functional of the field $\zeta_{W}({\bm{x}})$ it is enough to perform the change of variables $\zeta\to\zeta_{W}({{\bm{k}}})=W(k)\zeta({{\bm{k}}})$ in (2.8). One obtains

[TABLE]

where $P_{W}(k)\equiv W^{2}(k)P_{\zeta}(k)$. Leaving aside the argument of the deviation from Gaussianity ${\cal F}_{\rm NG}$ for a moment, this PDF has the same structure as the unfiltered PDF of Eq. (2.8).

This expression for ${\rm P}_{W}$ poses an interesting question: what if the window function of choice is defined (as usual) with hard cutoffs, just as if we were redefining the limits of our EFT? That is to say, how does ${\rm P}_{W}$ look if we have

[TABLE]

as the window function? It turns out that, for the functional integral to be well-defined, we have to integrate out of the theory all the modes which will not take part in our observable quantities. Upon doing so, we can avoid dealing with $\zeta_{W}({\bm{q}})/W({\bm{q}})$, which from the perspective of the theory with the window function would be ill-defined for the scales where $W=0$, but is perfectly finite (and equal to $\zeta({\bm{q}})$) from the perspective of the original theory. To this end, let us take the original functional distribution (2.8) (with $k$-space variables) and integrate out the modes outside the support of the window function $W$. We may write this as

[TABLE]

Upon integration over the prescribed range of modes, the purely Gaussian term in (2.8) gives a reduced Gaussian measure that considers only $k\in(k_{L},k_{S})$. To integrate the term containing ${\cal F}_{\rm NG}$ it is convenient to separate the integral over momentum space in (2.8) into two contributions:

[TABLE]

In this form, it is easy to see that upon integrating over $\zeta({\bm{k}})$, the term in the second line of (2.13) vanishes after performing a functional partial integration. Thus, we will only need to know how to deal with the quantity555Note that when considering the functional integration of the first line of Eq. (2.13), the differential operator $\left(\frac{\delta}{\delta\zeta(-{\bm{k}})}-\frac{\zeta({\bm{k}})}{P_{\zeta}(k)}\right)$ involves scales $|{\bm{k}}|\in(k_{L},k_{S})$, while the functional integration goes over modes with $|{\bm{k}}|\notin(k_{L},k_{S})$. We may thus pull it out of the integral.

[TABLE]

so as to see if, and how, the interaction is renormalised. It turns out that if we define

[TABLE]

then we may identify

[TABLE]

as the effective self-interaction, because we would have integrated out all the scales that are outside the range of interest and still leave the other scales within the measure ${\rm P}$ of the PDF, while maintaining its analytic structure.

The last equality of (2.15), the Weierstrass transform of $F(\zeta)$, can be obtained in numerous manners. If one were to follow standard diagrammatic perturbation theory, it arises from summing back every “loop” contraction performed by the Gaussian measure of ${\cal F}_{\rm NG}$ with itself. Since the momenta flowing through those loops are bounded by the range being integrated out, and there is no “external” momenta flowing through the diagrams, we have that their numerical value is the same for every loop and equal to the variance

[TABLE]

that was subtracted from the Gaussian field statistics when the modes $|{\bm{k}}|\not\in(k_{L},k_{S})$ were integrated out of the theory. Therefore, by using these results in Eq. (2.12), we obtain

[TABLE]

i.e., that ${\rm P}_{W}[\zeta]$ for the restricted variable (2.9) has the same functional form as the original PDF, with the only modification that now the departure from Gaussianity is given by a “filtered” interaction $\bar{\cal F}_{\rm NG}$ instead of the “bare” interaction ${\cal F}_{\rm NG}$.

In practice, there is more than one way of how to represent the running of $F$ depending on the scales one wants to include in the theory. Perhaps the most ethereal representation, but at the same time the most revealing of the theory’s structure, is through the differential expression of the Weierstrass transform, as in

[TABLE]

This expression makes clear how the theory runs by removing more or less scales, as well as the fact that the transformation rule between $F$ at different scales follows an adequate composition property: integrating out ranges of momenta $A$ and $B$ is implemented via $\sigma_{\rm out,A}^{2}$ and $\sigma_{\rm out,B}^{2}$, and doing so yields the same result independently of the order in which one subtracts the modes from the theory. Furthermore, this shows that the functional form of ${\cal F}_{\rm NG}$ and $\bar{\cal F}_{\rm NG}$ is the same, in the sense that the quantities that determine its concrete expression are exactly the same; the only thing that the window does is to restrict the range of modes entering in the observables.

Conversely, just as the PDF may be recast in an analogous manner to that of the original theory, the field with modes between $k_{L}$ and $k_{S}$ may also be written down as a local departure from Gaussianity

[TABLE]

merely because it follows statistics analogous to $\zeta$. Here we have to remind the reader that this is only so for $W(k)$ of the form (2.11). Other window functions still give rise to an explicit PDF, namely equation (2.10), but the deviations from Gaussianity may no longer be written as concisely as in (2.20). The difference lies in that a general window function does not render irrelevant some degrees of freedom of the theory; it only gives them dissimilar weights in the final result. However, in order to obtain the function ${\bar{\cal F}_{\rm NG}}$ it is crucial that we reduce the number of independent variables in our theory by integrating out their effects, as all of them will leave their signature, if small, in any given correlation function. Having made this point explicit, from now on, and for the rest of this paper, we will go back to considering general window functions.

Thus far, we have established how we may write the probability distribution functional of our theory depending on the scales under consideration and also how we may incorporate window functions into the distribution functional. Before passing to simpler statistical estimators stemming from this functional, we now discuss the partition function.

2.3 Partition function and $n$-point correlators

In practice, using the full probability distribution functional directly on cosmological data proves to be difficult, as there is only one realisation of our universe to probe and conduct measurements in, so it is not possible to take a frequentist approach to its statistics. While this suggests the use of Bayesian analysis to find the most probable ${\cal F}_{\rm NG}$ given the data by the means of the functional (2.8), it also reveals why one typically chooses to work with correlation functions to probe departures from Gaussianity: they can be computed from many Fourier modes on the sky, whose past history is presumably independent (at least if the nonlinearities are turned off), and therefore averages may be performed and compared with the theoretical predictions for the expectation values or correlations.

Fittingly, there is an object that encapsulates the information of all the correlation functions in a perhaps clearer way than the full probability density functional ${\rm P}$. This is the partition function $Z[J]$, which is the object that generates the $n$-point functions via functional differentiation:

[TABLE]

or equivalently, the functional Fourier transform of the PDF, which in the context of probability theory is called the characteristic function,

[TABLE]

Both expressions may be employed to obtain $Z[J]$: the first requires to know all the $n$-point functions beforehand and reconstruct the object that has them as its functional derivatives, while the second requires to know an explicit expression for the probability distribution functional. Since we have the latter, we may carry out this computation explicitly666The details of this derivation are presented in App. B. to first order in ${\cal F}_{\rm NG}$, obtaining

[TABLE]

Here we have defined $\sigma_{\zeta}^{2}\equiv\int_{k}P_{\zeta}(k)$ as the 1-point variance associated to the power spectrum for the relevant range of momenta. Window functions are easily incorporated by substituting $J({\bm{k}})$ with $J({\bm{k}})W(k)$, as this procedure will add a factor of $W(k)$ to every external leg in any given diagram.

Now that we have equation (2.23), we may compute the $n$-point functions directly, without having to resort to functional integration as we would with (2.8). Moreover, the structure that will emerge in these correlations is more closely related to (2.23), as is demonstrated by their explicit expressions in position space

[TABLE]

where the subscript $c$ indicates the result only considers the fully connected piece. For completeness, we write down their counterparts in momentum space

[TABLE]

where the coefficients777In a more standard notation, the first few terms would correspond to $f_{2}=f_{\rm NL}$, $f_{3}=g_{\rm NL}$, etc. $f_{n}$ are given by (probabilists’) Hermite moments of $F$:

[TABLE]

The coefficients $f_{n}$ are quantities of mass dimension $1-n$, which are invariant under the renormalisation procedure, discussed in Sec. 2.2, in a very fitting sense: because $\{f_{n}\}_{n=2}^{\infty}$ are the coefficients of a Hermite polynomial expansion888We omit $n=0,1$ in the Hermite expansion because we assume $\langle\zeta\rangle=0$ and $\langle\zeta\zeta\rangle$ to be set by the free theory matching the power spectrum of observations. To put it differently, with this definition of the local ansatz, due to orthogonality properties of the Hermite polynomials, the power spectrum of $\zeta$ is not modified to first order in the nonlinearity parameters [56]., we have

[TABLE]

where we have introduced the variance $\sigma_{\zeta}^{2}$ as an argument of $F$ in order to emphasise that the associated field has the corresponding amplitude for its fluctuations. This means that $\bar{F}$ in Eq. (2.15), or equivalently Eq. (2.19), takes the following form:

[TABLE]

where the coefficients $\{f_{n}\}_{n}$ remain unchanged; only the variance gets reduced to its new value after integrating some modes out.

Nonetheless, it is important to point out that each individual $n$-point function does not encapsulate all of the non-Gaussian information contained in $F$. Indeed, each correlation function only yields one term of an infinite series expansion of $F$, all of which are independent, at least in principle. Therefore, it is natural to try and find objects that are able to keep all of this information, without having to compute an infinite number of quantities. For this reason, we now turn to exploring 1- and 2-point probability density functions.

2.4 Fixed-point probability distribution functions

In the presence of a generic deviation from Gaussian statistics, involving both local and nonlocal terms, to assume that it is possible to capture all non-Gaussian information by looking at the single-point statistics of a field (i.e., correlations with all of the spatial coordinates at the same position) seems misguided, as the restriction to a single point is likely to mix local and nonlocal effects, making it difficult to disentangle them. However, if we restrict ourselves to local deviations from Gaussianity only, it is indeed possible to capture all such information.

In this subsection we write down explicit 1-point and 2-point PDFs for the curvature perturbation. In order to alleviate the discussion, we shall leave some of the details for App. C.

2.4.1 1-Point probability distribution function

Now we set ourselves to derive the simplest distribution that can be obtained within this framework: a density function for the 1-point statistics. It is defined as the distribution ${\rm P}(\zeta;{\bm{x}})$ that satisfies

[TABLE]

Given that we assume a homogeneous universe, ${\rm P}(\zeta;{\bm{x}})$ cannot depend on ${\bm{x}}$. Thus, we write ${\rm P}(\zeta;{\bm{x}})={\rm P}(\zeta)$. From a functional perspective, it is given by

[TABLE]

which may be evaluated in the same way as the partition function $Z[J]$ by writing down the Dirac delta as $\delta(\zeta({\bm{x}})-\bar{\zeta})=\int_{\gamma}e^{i\gamma(\zeta({\bm{x}})-\bar{\zeta})}$ and noticing that what will be left with in the functional integral is exactly $Z[J(-{\bm{k}})=\gamma e^{-i{\bm{k}}\cdot{\bm{x}}}]$. Then the remaining integral over $\gamma$ may be carried out by completing squares.

Thus, in the same spirit as the probability distribution functional we obtained earlier represents a first-order correction to Gaussian statistics, the 1-point probability density function may also be written as a slight departure from Gaussianity. Moreover, this density function resembles more closely the structure of $Z[J]$ than that of ${\rm P}[\zeta]$ because marginalising over all the other positions in presence of a finite range of wavelengths induces a filtering, which is manifest in $Z[J]$ but not so in ${\rm P}[\zeta]$. Given the various applications it may find in LSS or in primordial black hole formation [50], it is of interest to write it for an arbitrary window function.

Smoothing the field and its variance as in Eq. (2.9), the resulting expression is

[TABLE]

where

[TABLE]

Here we have written $W(x)=\int_{{\bm{k}}}e^{i{\bm{k}}\cdot{\bm{x}}}W(k)$, the position-space representation of the window function $W$, with $x=|{\bm{x}}|$, and we have also defined the (co)variances

[TABLE]

with a “bias” factor999This parameter is not the usual bias. If we think of the window as setting the scale of a tracer, we can define a bias as the ratio of field-tracer to the tracer-tracer correlation functions. ${\rm b}$ here is the ratio of the field-tracer to the 1-point tracer-tracer correlation functions. given by

[TABLE]

In words, this expression means that, given a range of modes defining the cutoffs of the theory, the perturbative correction to the PDF scans for the structure of the effective interaction at those scales through the action of $({\bar{\zeta}}-\sigma_{\zeta}^{2}\partial_{\bar{\zeta}})$, and then filters it, through the Weierstrass transform, according to the difference between the variance and the correlation implied by incorporating window functions. This observable can account for all the information contained within the local function $F$. This can be seen from the fact that all of the information concerning $F$ is stored within the $f_{n}$ coefficients, which can be retrieved by looking only at the 1-point statistics (c.f. Eq. (2.24)). Indeed, an analysis aiming to constrain the departure from Gaussianity ${\cal F}_{\rm NG}$ by the means of the 1-point CMB temperature distribution has already been performed in [49] .

Before passing to the 2-point functional, and in order to make contact with current literature, let us comment on the relation of this PDF to the Edgeworth representation. If we use the Hermite polynomial expansion of $F$, given in Eq. (2.27), the effect of the Gaussian filtering of Eq. (2.32) becomes transparent:

[TABLE]

that is, it replaces the field variable with the biased one. Using the fact that the Weierstrass transform commutes with derivatives, the latter now being evaluated at the biased field with its corresponding variance, and that $(x-\sigma^{2}\partial_{x}){\rm He}_{n}(x/\sigma)=\sigma{\rm He}_{n+1}(x/\sigma)$, we can rewrite the NG deviation of the 1-point PDF as

[TABLE]

where $\langle\zeta_{W}^{n}\rangle_{c}$ is given by Eq. (2.25) integrated over momenta. This is exactly the Edgeworth expansion of a non-Gaussian PDF truncated to first order in the couplings $f$, since we have restricted our derivation of the PDF (2.32) to first order101010In principle, there is nothing stopping us from computing the Edgeworth expansion to any order in the couplings; indeed, we can expand the exact PDF of Eq. (2.5) to any order in $F$. We choose, however, for simplicity to truncate the series to first order. in $F$.

Moreover, in this case it is possible to invert the Edgeworth expansion (2.36) in terms of the function $F$. Taking the $m$th Hermite moment of this expansion gives

[TABLE]

and by replacing the cumulant $\langle\zeta_{W}^{m}\rangle_{c}$ in terms of $f_{m}$ we get

[TABLE]

which means that we have obtained the coefficients of the Hermite polynomial expansion (2.27) of $F$. Therefore, all one needs to write an explicit expression for $F$ is to perform the correlation integrals $\int_{0}^{\infty}\!\!{\rm d}x\;4\pi x^{2}W(x){\rm s}^{2n}(x)$, and then sum back111111The precise scaling with $n$ of the correlation integrals depends on the window function of choice $W$, and thus it is not always possible to give the result for $F$ in a closed form. its Hermite expansion (2.27). This allows for a direct test of NG in a given dataset and provides a way of constraining the primordial quantities involved in its generation, in a manner that is identical to what has been previously done in [49].

However, this type of analysis integrates out information about the momentum dependence of the correlation functions, which can be used to place tighter constraints on the parameters of a given model. We now proceed to explore the two-point PDFs for the filtered primordial perturbations, which do contain such information.

2.4.2 2-Point probability distribution function

Now we would like to write down an observable able to account for the non-Gaussian statistics, without integrating out the information about correlations in the sky at different scales. Therefore it cannot be a single-point PDF. Thus, we try to do the next least complicated thing: a 2-point PDF ${\rm P}(\zeta_{1},\zeta_{2};{\bm{x}}_{1},{\bm{x}}_{2})$. This function satisfies

[TABLE]

where we have written ${\rm P}(\zeta_{1},\zeta_{2};{\bm{x}}_{1},{\bm{x}}_{2})={\rm P}(\zeta_{1},\zeta_{2};|{\bm{x}}_{1}-{\bm{x}}_{2}|)$ because we are assuming our universe to be statistically homogeneous. One way to obtain such an object is by conditioning in two points in a manner analogous to Eq. (2.30),

[TABLE]

where again, the final result can only depend on the spatial coordinates through the distance between the two positions ${\bm{x}}_{1}$ and ${\bm{x}}_{2}$. With this in mind, let us define the scalar variables

[TABLE]

Then, by inspecting the $n$-point functions of the theory (2.24) evaluated at the two points of interest, we obtain a similar expression as the 1-point PDF, but with two points defining the filtering instead of one:

[TABLE]

where ${\rm P}_{G,W}(\zeta_{1},\zeta_{2},r)$ is the bivariate Gaussian measure, with a covariance matrix given by the $2\times 2$ bottom right block of $\boldsymbol{\Sigma}$, defined below in Eq. (2.43). Let us go through this expression: the first thing to notice is the presence of two points, ${\bm{x}}_{1}$ and ${\bm{x}}_{2}$, defining a filtering through the same function as in the 1-point case. The second important aspect is that now the Gaussian that is convoluted with $F$ has a different mean and variance. However, they emerge in the same manner as $\sigma_{W}^{2}(x)$ and $\zeta_{W}(x)$ emerge in the 1-point case: $\sigma^{2}_{W}(r,r_{1},r_{2})$ and $\zeta_{W}(r,r_{1},r_{2})$ are the variance and mean of $\bar{\zeta}$ after conditioning on the values of $(\zeta_{1},\zeta_{2})$, starting from a joint Gaussian distribution for $({\bar{\zeta}},\zeta_{1},\zeta_{2})$ with covariance matrix

[TABLE]

where we have written the covariance between the two externally chosen points ${\bm{x}}_{1}$ and ${\bm{x}}_{2}$ as

[TABLE]

The functions $G_{ij}$ also appear in a similar way: $G_{i1}$ and $G_{i2}$ are “rotated” versions of $\sigma_{\zeta}^{2}$ and $\bar{\zeta}$, involving combinations of the free theory covariances that make the overall expression reduce to that of the 1-point PDF as ${\bm{x}}_{1}\to{\bm{x}}_{2}$. Their precise definitions are listed in App. C. In there, we delineate how to obtain the 2-point PDF in an alternative way: we start from correlators of the type $\langle\zeta_{W}^{n}({\bm{x}}_{1})\zeta_{W}^{m}({\bm{x}}_{2})\rangle$ and we deduce the function from which they emanate, corresponding to (2.39).

This PDF contains all the information of the free theory, as having two points allows to scan over all the range of distances in the sky, thus probing, among others, the 2-point correlation function completely, which is the defining object of a Gaussian theory. Even though Eq. (2.42) has its non-Gaussian features encoded in a perhaps more complicated fashion than its 1-point counterpart (2.31), both contain the same information about the underlying function ${\cal F}_{\rm NG}$. Indeed, one can obtain Eq. (2.31) by integrating over one of the field variables in Eq. (2.42). However, observationally, it might be more efficient to have information on the scale, since then we can, for example, disentangle different momentum shapes of correlation functions.

2.5 Linear transfer functions and their limitations

The above discussion sets a clear framework to study various cosmological observables in the late-time universe with linear perturbation theory. In particular, we can consider endowing the window function $W(k)$ with an additional argument, thus defining the transfer function $T({\bf a},{\bm{k}})$, to obtain an observable $A({\bf a})$, where $\bf a$ is a label for the new quantity/field that is linearly related to $\zeta$ as

[TABLE]

up to first order in cosmological perturbation theory. In the case of CMB and LSS, the temperature fluctuations map

[TABLE]

and the matter density contrast

[TABLE]

are well described by linear cosmological perturbation theory as long as the distance/momentum scales are sufficiently large/small so that late-time nonlinearities do not affect them. This is indeed the case for the CMB scales $10^{-4}\;{\rm Mpc}^{-1}\leq k\leq 0.3\;{\rm Mpc}^{-1}$ Planck observed [3, 5], and for the large-scale structure of our universe provided that $k\lesssim 0.1\;{\rm Mpc}^{-1}$ and the redshift $z$ satisfy $z\gtrsim 2$ [57] so that gravitational nonlinear effects remain suppressed. Sections 3 and 4 delve into statistical objects and probes that can be used to test for imprints of non-Gaussianity in the CMB and LSS, respectively, in this way.

However, even though for Planck and Wmap linear perturbation theory is indeed accurate enough, as the errorbars are expected to keep shrinking, this might not be the case for future surveys  [58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 30, 68, 69, 70, 71]. At this point, it is thus appropriate to clarify that our subsequent discussion, as it stands, is valid for perturbation theory with linear transfer functions, and as such, it will have some limitations inherited from this approach. For example, we will not take into account the fact that the gravitational potentials differ at second order, nor will we use the nonlinear transfer function corresponding to such higher order terms [72, 68, 69, 70, 73]. We also do not discuss nonlocal effects like the ISW to second order, which can be included in the window function via a convolution in time. That would complicate the formulas because it would force us to include additional integrals over the history of our universe to preserve the generality of the approach (and not only over a single spatial slice), which could in turn reduce the usefulness of defining transfer functions in our setup, as including integrals over time might often be equivalent to following the dynamical evolution of the perturbations up to some fixed order in time-dependent perturbation theory. The alternative would be to carry out a traditional perturbative computation, without resorting to a probability distribution function in the process. There remains, hence, the question of what the information extracted from measurements means for our understanding of the primordial universe: since the late-time nonlinear evolution is only treated perturbatively, at some order one might lose primordial information in the local ansatz because nonlinearities in the later evolution of the universe will modify the observable consequences of the primordial ${\cal F}_{\rm NG}$, and consequently, the primordial PDF we would infer for $\zeta$.

We note, however, that if the experimental precision demands to have control of the cosmological perturbations up to second order, one can intent to proceed further and include such effects by writing the full temperature field (or the matter density contrast $\delta$) as $\Theta=\Theta^{(1)}+\Theta^{(2)}+\cdots$. One could then induce a PDF for the fully nonlinear field taking into account both primordial and secondary NG to any desirable order121212By doing so, one would have to consider the field and its gradient as independent Gaussian variables, since the latter contributes to the second order curvature perturbation [72]., allowing for generic nonlocal mappings ${\cal F}_{\rm NG}$. In other words, one can parametrise the generalised ansatz as in Eq. (1.5), in the manner of an effective object that adequately describes the late-time observables, where the expansion coefficients account for nonlinearities in both the initial conditions and due to gravitational evolution [70]. In [67, 71] it was shown that second order effects provide for corrections to the skewness (albeit with distinguishable signals from the primordial component) around the sensitivity of future CMB surveys like CMB-S4, and one might, arguably, worry that this could also be the case for higher order correlators like the trispectrum; it is, however, conceivable that the secondary effects would keep diminishing as one considers higher order statistics. A way out could thus be a case-by-case evaluation of the secondary effects in higher order $n$-point functions [70], or a full characterisation of the secondary effects at the level of the PDF, via, e.g., simulations.

Even though, admittedly, we do not offer any crucial insight into the problem of disentangling primordial from late-time NG, our results set the stage for further study of NG signatures on PDF estimators, especially targeting models beyond single-field inflation (c.f. App. A.1), where primordial local NG is enhanced compared to other shapes of NG and secondary effects.

3 $\Theta$ Statistics

Having set the stage for studying late-time cosmological perturbations, we can now turn to discussing our specific scenarios of interest. To begin with, it is of particular relevance to write down testable quantities that we can obtain by looking at the primordial information that can be stored in spherical shells on the sky, such as the CMB. Namely, both the probability distribution functional ${\rm P}[\zeta]$ and its associated partition function $Z[J]$ may be projected onto the celestial sphere to yield distributions of the temperature fluctuations $\delta T(\hat{n})$.

Consider a generic linear transfer function $T({\bm{k}},\hat{n})$ from the primordial perturbations $\zeta$ to an observable defined on the sphere $\Theta(\hat{n})$ such that

[TABLE]

Of course, we are interested in the specific case where $\Theta(\hat{n})=\delta T(\hat{n})=(T(\hat{n})-T_{0})/T_{0}$ (where $T(\hat{n})$ is the CMB temperature measured in a given direction $\hat{n}$, and $T_{0}$ its average over all directions), however, the present discussion is also valid for more general observables that depend on the line of sight $\hat{n}$. Then, if we take $\Sigma(\hat{n},\hat{n}^{\prime})$ to be the covariance matrix defining the observable’s correlations between different directions in the sky $(\hat{n},\hat{n}^{\prime})$, i.e.,

[TABLE]

and $\Sigma^{-1}(\hat{n},\hat{n}^{\prime})$ as its inverse matrix, we find that the probability distribution functional for $\Theta$ is given by

[TABLE]

where now ${\rm Ker_{1}}({\bm{x}})$, ${\rm Ker_{2}}(\Theta;{\bm{x}})$, $\zeta_{\Theta}({\bm{x}})$, and $F_{\Theta}(\zeta;{\bm{x}})$ depend implicitly on the transfer function and, when noted, the variable of interest $\Theta(\hat{n})$. All of these quantities may be intuitively understood as the result of projecting a field defined on three spatial dimensions over a two-dimensional spherical surface. For instance, $\zeta_{\Theta}({\bm{x}})$ is given by

[TABLE]

which basically amounts to saying: take the original statistics of your theory, i.e., $e^{i{\bm{k}}\cdot{\bm{x}}}P_{\zeta}(k)$, project them onto the sphere by applying $T(\hat{n},{\bm{k}})$ and then correlate it with the field of interest $\Theta(\hat{n}^{\prime})$ by means of its inverse covariance matrix $\Sigma^{-1}(\hat{n},\hat{n}^{\prime})$. The integration kernels have similar definitions:

[TABLE]

However, the function $F_{\Theta}$ is a slightly different object than before. As a result of the projection, it acquires a spatial dependence, whose exact nature in terms of the primordial departure from Gaussianity $F$ is given by the Weierstrass transform:

[TABLE]

where $\sigma_{\zeta\Theta}^{2}({\bm{x}})$ may be understood as a position-dependent effective variance of $\zeta$, modified by projection effects:

[TABLE]

In all of the above, we have written $\Sigma(\hat{n},\hat{n}^{\prime})$ as a general function of the direction on the sphere. However, if we take into account that our universe is homogeneous, it must be possible to write it as a function of the angle between the two vectors, or equivalently, in terms of their scalar product $\hat{n}\cdot\hat{n}^{\prime}$. We will not overemphasize this in what follows, as the notation we deem natural to treat $\hat{n}$ and $\hat{n}^{\prime}$ is with them as separate directions, because they will have to be multiplied with another vector, the integration variable ${\bm{x}}$ that will appear in the NG kernel that modifies the PDF, which makes using $\hat{n}\cdot\hat{n}^{\prime}$ notationally heavier than just using $(\hat{n},\hat{n}^{\prime})$.

Now we have to find ways of using this. One option would be to compare how likely is our present-day CMB given a certain primordial deviation from Gaussianity of the local type $F$ by using Bayesian statistics. However, in order for this to be useful at its maximum capacity, it is likely that one would first have to establish a definitive imprint of primordial NG and be forced to introduce extra parameters into the effective description because more often than not a model comparison will favour the one with less parameters. Therefore, we will turn to the observable we outlined at the end of Sec. 2, which, to our knowledge, has been largely unexplored in this context131313The 2-point PDF has been used in LSS count-in-cells statistics [74]. and may offer valuable constraints on the nature of primordial NG: a 2-point PDF. We present this result after briefly discussing the 1-point statistics.

3.1 1-Point temperature distribution function

Given the discussion we have followed so far, it is now straightforward to reproduce the results of [49] and compute the 1-point PDF of the temperature map. The only difference is that the transfer functions map the fluctuations to the celestial sphere instead of the three-dimensional universe, but this does not prevent us from finding the distribution that generates the correlations. In analogy to (2.31), the result is given by

[TABLE]

where

[TABLE]

with the corresponding definitions of the window function

[TABLE]

and that of the variances

[TABLE]

with $b_{\Theta}$ given by

[TABLE]

The final result is independent of the angular direction $\hat{n}$ because of homogeneity. In practice this appears because the form of the transfer functions $T({\bm{k}},\hat{n})$ will imply that the only dependence on $\hat{n}$ will be through scalar products with ${\bm{x}}$ or ${\bm{k}}$, which as integration variables span the whole space, and thus they render all directions statistically homogeneous.

In Ref. [49], an analysis of the CMB statistics was carried out to search for non-Gaussian signals by estimating the PDF (3.9) directly from the 2015 Planck dataset [75], with an inconclusive result: if primordial non-Gaussianity of this type is present in the CMB, its statistical imprints are smaller than the intrinsic noise to which the 1-point statistics are subject to. With this and the forthcoming surveys in mind, we now proceed to outline another statistical object whose reconstruction from observations shows promise to have a lower intrinsic noise but still preserve all of the information contained within the 1-point PDF: a 2-point PDF.

3.2 2-Point temperature distribution function

It is interesting to write down expressions for quantities that are not typically used directly when characterising cosmological datasets. For instance, one rarely bothers to write down the full 2-point PDF for the CMB temperature map, as all of its information (in the Gaussian case) is already specified through the power spectrum. However, the scale dependence of this PDF may be a useful tool to search for non-Gaussianities. The result, analogously to what we had in the simpler case of a spatial window function, is

[TABLE]

As before, let us go through this expression: the first thing to notice is that the two filtering points ${\bm{x}}_{1}$ and ${\bm{x}}_{2}$ in (2.42) are now replaced by two directions in the sky, $\hat{n}_{1}$ and $\hat{n}_{2}$. Secondly, the Gaussian that is convoluted with $F$ now has a different mean and variance, obtained by conditioning the joint multivariate Gaussian PDF on the values of $\Theta_{1}$ and $\Theta_{2}$. The starting point to this is a joint Gaussian distribution for $({\bar{\zeta}},\Theta_{1},\Theta_{2})$ with covariance matrix

[TABLE]

where we have written the covariance between the two externally chosen points $\hat{n}_{1}$ and $\hat{n}_{2}$ as

[TABLE]

and $\sigma_{\Theta}^{2}=\sigma_{\Theta,{\rm ext}}^{2}(\hat{n},\hat{n})$, which is independent of the direction $\hat{n}$. The precise definitions of all the additional functions involved in this section are given in App. C. It is worth mentioning, as a reminder to the reader, that, as in the 2-point PDF for curvature fluctuations (2.42), both (3.16) (besides from the temperature variables) and $\Sigma(\hat{n}_{1},\hat{n}_{2})$ depend only on the angular distance between $\hat{n}_{1}$ and $\hat{n}_{2}$, or equivalently, on $\hat{n}_{1}\cdot\hat{n}_{2}$, and that this is a consequence of our universe’s homogeneity.

From this function, i.e., from the 2-point PDF (3.16), it is possible to obtain refined constraints on the local ansatz. Given a dataset, one can construct the 2-point PDF as follows: divide the angular distance into $N$ bins of width $\delta\vartheta$ and the temperature in $M\times M$ bins of size $\delta\Theta\times\delta\Theta$, as in the 1-point PDF but now with two axes for the temperature field. Now, for each bin associated to a given angular distance $\vartheta_{n}=n\cdot\delta\vartheta$, and for each value of $(i,j)$, count how many pairs of pixels separated by that angular distance $\vartheta_{n}$ have the values $(\Theta_{[i]},\Theta_{[j]})$ for the temperature in their respective positions. This process would generate $N$ two-dimensional histograms, with two temperature axes, which we label by $(\Theta_{1},\Theta_{2})$, whose value at coordinate $({\Theta_{1}}_{[i]},{\Theta_{2}}_{[j]})$ would give the number of pairs of pixels with temperatures in the $(i,j)$th bin, separated by angular distance in the $n$th bin. A $1/2$ symmetry factor must be included in the number counts for temperature bins with $i\neq j$, as the bin $({\Theta_{1}}_{[i]},{\Theta_{2}}_{[j]})$ is equivalent to $({\Theta_{1}}_{[j]},{\Theta_{2}}_{[i]})$.

How does this give refined constraints on the local ansatz? Let us appreciate that this set of PDFs contains information on the scale, or more concretely, on the temperature power spectrum and of its expansion in spherical harmonics (the standard $C_{\ell}$s) through the 2-point correlation $\Sigma(\hat{n}_{1},\hat{n}_{2})$. If NG is absent, then at each value of the angular distance the 2-point PDF will be a 2-variable Gaussian probability density with variances $\sigma_{\Theta}^{2}$ and covariance $\sigma_{\Theta,{\rm ext}}^{2}(\hat{n}_{1},\hat{n}_{2})=\Sigma(\hat{n}_{1},\hat{n}_{2})$. Then, in the presence of NG, each 2-point PDF will undergo a NG deviation induced by the same primordial mechanism $F$, but for each angular separation this deviation will be experienced differently because the covariance matrix implied by the Gaussian part is different. This means that for each angular distance, the kernel that acts upon $F$ in (3.16) gives a different deviation from Gaussianity, and therefore, each of the 2-point PDFs gives an independent estimator on the primordial NG field. For local NG, all of the $N$ 2-point PDFs at different angular separations in the sky should give consistent141414That is, within the experiment’s theoretical and systematic uncertainties. constraints/estimations of $F$. Conversely, if NG is measured and it does not adjust to the statistics implied by (3.16) at different angular scales, then purely local NG would be ruled out. Therefore, looking towards possible future directions to be explored, this type of object (a set of 2-point PDFs) shows promise to disentangle different shapes of NG, such as equilateral or orthogonal templates, in particular, from the local ansatz.

In order to search for non-Gaussianity within a 2-point PDF, many approaches are possible. Given a model, i.e., an explicit expression for $F$, and using it as a template with few adjustable parameters is usually the method that will give the best constraints. In the spirit of Eq. (2.27), however, another one is worth mentioning: one can use bivariate Hermite polynomials on the temperature variables $(\Theta_{1},\Theta_{2})$, so as to express the PDF in terms of a bivariate Edgeworth expansion [76, 77]. One then looks for any statistically significant nonzero coefficient in the expansion, in analogy to what was done for the 1-point case in [49]. This way, the existence of NG can be tested as a yes/no question, as any nonzero coefficient in an Edgeworth expansion implies a non-Gaussian distribution. This may be particularly useful when searching for NG in the next generation CMB surveys [78].

4 $\delta$ Statistics

Upcoming cosmological surveys will focus on the statistics of LSS, promising to bring precision cosmology to a new era. Indeed, the proliferation of observed modes due to the three-dimensional probe offered by LSS will enhance our statistics giving us invaluable information about the fundamental aspects of the early/late universe.

At small scales several sources of nonlinearity induce NG, like gravitational interactions and galaxy bias, obscuring the primordial contribution to the statistics. However, for relatively long modes, $k\lesssim 0.1\;$Mpc*-1* and higher redshifts [57], linear perturbation theory can be trusted, which makes it easier to identify primordial signatures. Perturbative techniques pushing our analytic control towards smaller, weakly nonlinear scales include several schemes like SPT [79] and more recently EFTofLSS [80] and TSPT [81, 82, 83], which are set within a hydrodynamics framework, while going even further requires full Boltzmann solvers via N-body simulations. In this work, we will focus on the purely linear regime, leaving weakly nonlinear evolution with NG initial conditions [83] for future study.

The main probes of non-Gaussianity are the bispectrum and/or trispectrum, number counts and bias. The spectra retain information about the shape of the 3- and 4-point functions in momentum space, which can be linked to the mechanism responsible for generating NG. Number counts probe directly the 1-point PDF, which even though loses the shape information, it serves as a complementary and equally powerful estimator of NG [54, 84, 85, 86, 87, 88, 89]. Finally, the halo bias serves as a third independent channel, which can give clear enhanced signatures of local NG at large scales.

In this section, we wish to track how primordial NG, in the form of the generic local ansatz (1.5), gets transmitted to the matter field in the linear regime. We extend our result, the non-Gaussian PDF of curvature fluctuations, in two directions: 1) we deduce a PDF for the matter density contrast $\delta$, and hence, a halo mass distribution; 2) we compute the effect of the generalised local ansatz on the halo bias. These are complementary probes of the non-Gaussian initial condition via cluster number counts and power spectra, respectively, which should be accessible by surveys such as Lsst, Euclid, Spherex and Ska.

4.1 Halo mass function

The matter overdensity $\delta({\bm{x}})=\delta\rho({\bm{x}})/\bar{\rho}$, with $\bar{\rho}=\Omega_{m}\rho_{\rm cr}$, is related to the primordial Newtonian potential, $\Phi=\frac{3}{5}\zeta$, as

[TABLE]

where $D(z)$ is the linear growth rate, $r_{H}$ the current Hubble radius and ${\mathcal{T}}(k)$ the transfer function [90]. We smooth the density field over a radius $R_{M}=\left(3M/4\pi\bar{\rho}\right)^{1/3}$ as in Eq. (2.9),

[TABLE]

using a top-hat filter $W_{M}(x)=V_{M}^{-1}{\rm H}(R_{M}-x)$, with H the Heaviside function and $V_{M}$ the volume enclosing a mass M. The probability distribution for the smoothed overdensity, ${\rm P}(\delta_{W})$, is then given by the 1-point PDF of Eq. (2.31) upon the replacements

[TABLE]

Now, in principle, having the matter distribution function one can compute the halo number density, that is, the number density of halos of mass between $M$ and $M+{\rm d}M$ at redshift $z$, and the observables derived from it. One way to do this is via the Press-Schechter (PS) scheme [91] extended to the NG case [92, 93]. Let

[TABLE]

(with $\nu\equiv\delta/\sigma_{W}$) be the tail distribution above some threshold value $\delta_{c}(z)$. Then, one can assert that the total fraction of mass collapsed into bound structures will be proportional to this cumulative PDF:

[TABLE]

where ${\rm d}n(M)$ is the number density of halos with masses in the range $M$ and $M+{\rm d}M$. The PS mass function then reads151515The fudge factor $2$ corrects for the cloud-in-cloud problem, that is, a collapsed object of mass $M_{1}$ tracing a volume $V_{1}$ can be missed by the PS function if it is part of an underdense region of larger volume $V_{2}$. The Gaussian value $2$ is a good approximation in case of a small NG deformation [54].

[TABLE]

The collapse threshold, through which the $z$ dependence arises, is taken to be given by the spherical model as $\delta_{c}(z)\simeq 1.686\;D(0)/D(z)$. The Gaussian PS function can be evaluated exactly by replacing ${\rm P}\to{\rm P}_{\rm G}$ in Eq. (4.4):

[TABLE]

For the NG case, the PS mass function is easy to compute from Eq. (2.31):

[TABLE]

where a prime stands for the derivative with respect to the mass, while $F\left({\rm b}\delta_{c};{\rm b}^{2}\sigma_{W}^{2}\right)$ is the Weierstrass transform of the local ansatz [see Eq. (2.35)] evaluated at the threshold $\delta_{c}(z)$ (note that here, $W(x)$ is the Fourier transform of the window function written in Eq. (4.3)). Upon using the Hermite expansion (2.27) of the function $F$, we may obtain a series representation of the mass function (4.8):

[TABLE]

where we have defined the reduced cumulants $\kappa_{n}(M)\equiv\frac{\langle\delta_{W}^{n}\rangle_{c}}{\sigma_{W}^{n}}$. Equation (4.8) offers a generalisation of the $f_{\rm NL}$, $g_{\rm NL}$ truncation161616For example, when truncated to $n=2$, the expansion (4.9) agrees with Eq. (4.19) of Ref. [84]. of the local ansatz (see e.g. Refs. [84, 94]) to arbitrary functions $F$. From its moment expansion (4.9), we can see that positive moments of the PDF ($\kappa_{n},\kappa^{\prime}_{n}>0$) lead to overabundance of collapsed objects at the high mass end, where $\nu,\nu^{\prime}\gg 1$, as long as $(\ln\kappa_{n+1})^{\prime}<\nu\nu^{\prime}$, which is satisfied since the cumulants $\kappa$ depend weakly on the mass [94].

However, due to the highly nonlinear character of the collapse, one cannot fully parametrise collapsed objects by a single threshold number $\delta_{c}$. Indeed, it has been shown that the PS prescription does not accurately estimate the halo abundance even in the Gaussian case (for example a small non-spherical perturbation can have a considerable impact) [95]. Hence, the extension to the NG case is guaranteed to also have errors with respect to simulations. What can be done though is to characterise the deviation from Gaussianity by comparing the ratio of G-to-NG densities to that of the PS scheme [96, 84], since the latter is expected to deviate equally in both cases:

[TABLE]

For the Gaussian mass function, \frac{{\rm d}n}{{\rm d}M}\big{|}_{\rm G}, we can adopt a Sheth-Tormen (ST) ansatz [95], which is better fitted to simulations than the Press-Schechter one, in which case the Gaussian ratio reads

[TABLE]

with the ST parameters $a=0.707,\;A=0.322184,\;p=0.3$. With the mass function (4.10) at hand, which is just Eq. (4.8) with the replacement $n_{\rm PS}|_{\rm G}\to n_{\rm ST}$, we may compute the number of clusters per redshift bin above some mass $M$ as [84]

[TABLE]

where $f_{\rm sky}$ is the fraction of the sky covered by each survey.

This formula can serve as a template for the number density of clusters. For example, in the context of light isocurvature axions one can compute the non-Gaussian PDF deformation $\Delta$ [48] and thus the mass function $\frac{{\rm d}n}{{\rm d}M}$ and create mock data via simulations. One can then pick an estimator and devise an overlap between the template and the data, in exactly the same manner that one uses the cosine estimator to measure the overlap of the local template with the actual bispectrum. The simplest thing to do is to take the number counts as an estimator. For example, given a model, one can count how many clusters of mass $M$ exist at redshift $z$ in the mock data and compare this number to the real data.

We may also try to reconstruct the PDF from data using other statistical estimators. This can be done by simply solving Eq. (4.10), together with (4.6), as an ODE for the cumulative PDF, $\mu^{>}$, to get

[TABLE]

with $r_{\rm G}$ given by Eq. (4.11). The left hand side now gives the tail distribution of $\delta$, that is, the probability of $\delta>\delta_{c}$ at redshift $z$. On the right hand side we have the total mass of collapsed objects with $m>M$ at redshift $z$, accounting for corrections to the spherical collapse model via assigning a mass value $m(z)=m_{\rm obs}/r_{\rm G}[\nu_{c}(z)]$ to an object of observed mass $m_{\rm obs}$ and multiplicity $\frac{{\rm d}n}{{\rm d}m}(m_{\rm obs},z)$; in principle, this can be deduced from number counts. One can now apply statistical estimators like Minkowski functionals to the dataset $\left\{\frac{m_{\rm obs}}{r_{\rm G}[\nu_{c}(z)]},\frac{{\rm d}n}{{\rm d}m}(m_{\rm obs},z)\right\}$, whose difference from their Gaussian estimate, as in the CMB [97, 98] and LSS [99], will be a direct probe the NG deformation of the primordial PDF. Let us note, however, that measuring accurately the mass and redshift of clusters is a hard task, which might complicate a direct connection between cluster counts and primordial NG in this manner.

To summarise, in the context of tomographic NG, one can 1) use Eqs. (4.9) and (4.10) as a template for the mass function and thus the number density up to arbitrary order in the Edgeworth expansion; 2) assume a microphysical model, compute the PDF (as in e.g. Ref. [48]) and use Eqs. (4.8) and (4.10) as a template without having to refer to moments; and finally, 3) consider statistical estimators on the LSS dataset probing the primordial PDF via the counting scheme implied by Eq. (4.13).

4.2 Scale dependent halo bias

Another powerful probe of primordial non-Gaussianity is the halo bias, which enters in the late-time power spectra. As already argued, non-Gaussian initial conditions alter the halo abundance in a nontrivial way by increasing the number of rare density peaks that collapse into halos. This is easier to visualise in the case of the quadratic local ansatz $\delta=\delta_{\rm G}+f_{\rm NL}\delta_{\rm G}^{2}$ [15]: a positive $f_{\rm NL}$ adds positive skewness to the density distribution; thus, the same probability now corresponds to higher values of $\delta$ with respect to the Gaussian field, leading to more probable enhanced peaks. In [15], it was shown how this is encoded in a scale dependent correction to the halo bias given by

[TABLE]

where $b_{\rm G}$ is the Eulerian Gaussian bias. The $k^{2}$ factor in $\alpha$ (see Eq. (4.1)) implies that the effects of non-Gaussianity will be accentuated in large scales (the transfer function goes to 1 for $k\to 0$), which brings surveys like Ska and Lsst to the frontline of NG searches, promising $\sigma(f_{\rm NL})\lesssim 1$ [100, 101]. This result was rederived and generalised in [102, 103] using different approaches and has been confirmed with N-body simulations (see e.g. Refs. [104, 105, 106, 107, 108, 109]), while in Refs. [56, 94], the scale dependent halo bias was computed for the case of cubic $g_{\rm NL},\tau_{\rm NL}$-type local NG —see Refs. [110, 111] for reviews. In what follows we extend it to the generalised local ansatz (1.5) and show how in the case of an isocurvature source, future surveys can probe the landscape potential via measurements of the bias factor.

To begin with, let us rewrite the NG ansatz (1.5) for the Newtonian potential $\Phi=\frac{3}{5}\zeta$ as

[TABLE]

where $\phi$ is a Gaussian random field with standard deviation $\sigma_{0}=\sqrt{\langle\phi^{2}\rangle}$. The $3/5$ factor is put explicitly to relate the primordial curvature perturbation $\zeta$ to the late-time (matter-dominated era) gravitational potential $\Phi$. Following the peak-background split method [112, 113], we now separate the Gaussian gravitational potential into long and short modes with respect to some characteristic halo scale $R_{\star}\sim R(M)$ as

[TABLE]

which induces a similar split in the variance, $\sigma_{0}^{2}=\sigma_{L}^{2}+\sigma_{S}^{2}=\langle\phi_{L}({\bm{x}})^{2}\rangle+\langle\phi_{S}({\bm{x}})^{2}\rangle$. With the help of the expansion in Hermite polynomials (2.27), and using known identities of these polynomials, we can write

[TABLE]

where

[TABLE]

with $f_{0}=f_{1}=0$. How is this expansion of $F$ useful? If one is mainly interested in the short-wavelength dynamics (for instance, to study the gravitational collapse of the matter distribution into galaxies), this expansion allows one to identify the functions $\beta_{m}$ of long-wavelength fluctuations as effective nonlinearity parameters for the short modes, as may be seen directly from (4.17): $m=1$ is a correction to the amplitude of the short modes, $m=2$ a nonlinearity of $f_{\rm NL}$-type, $m=3$ a nonlinearity of $g_{\rm NL}$-type, etc.

We now have to understand quantitatively how long wavelength density fluctuations affect the statistics of short modes and thus the halo number density per unit of halo mass, ${\rm d}n/{\rm d}M$, which from now on we denote by $n_{L}({{\bm{x}}})$, through each term in the expansion (4.17). Generically, the halo mass function is a function of the matter contrast and the amplitude of the short modes:

[TABLE]

Firstly, because of the dependence on $\rho$, irrespective of the presence of non-Gaussianity, a long wavelength perturbation $\delta_{L}$ will induce a linear background shift in $n$ as

[TABLE]

where $\bar{n}=n_{L}[\bar{\rho}]$. Moving to the non-Gaussian part, the purely long wavelenght contributions —$m=0$ coefficient in the expansion (4.17)— by definition will not affect the power spectra, and hence, the bias to first order in the nonlinearity parameters $f_{m}$, so we may disregard them. From a short-wavelength modes’ perspective, they are constant numbers that will only affect the background density through $\delta_{L}$ in the previous expression.

Next, we may observe that the coefficient of the term linear in $\phi_{S}$ ($m=1$) depends on $\phi_{L}$. Therefore, a long mode will induce a shift in $n({{\bm{x}}})$ through this term, since short modes feel a background perturbed by the local amplitude of the long wavelength perturbation as $\Delta_{\phi_{S}}\to\left(1+\beta_{1}\!\left(\phi_{L}\right)\right)\Delta_{\phi_{S}}$, leading to

[TABLE]

Finally, we need to take into account that the halo density depends on the local amplitude of the long modes through all the functions $\beta_{m}\left(\phi_{L}\right)$ for all $m\geq 2$, since each $\beta_{m}$ acts as a nonlinearity parameter assuming a local value set by the long wavelength fluctuation $\phi_{L}$ [103, 56]. That is, we should consider the halo density as a function of the form

[TABLE]

where $\beta_{1}$ controls the amplitude of the short modes, and the rest controls the amplitude of their nonlinearities. Thus, we can write

[TABLE]

From this expression we can now compute the linear bias, defined as

[TABLE]

where the matter-halo power spectrum is defined as

[TABLE]

and the matter-matter one as

[TABLE]

Thus, we will need to compute the matter-halo correlator and Fourier transform it ($\mathbb{F}$). The matter contrast $\delta$ is related —in subhorizon scales— to the Newtonian potential through the Poisson equation (4.1), so that $P_{\phi\delta}=P_{\delta\delta}/\alpha$. The first term in Eq. (4.23) is trivial since it just yields the $\delta$ propagator, which in Fourier space cancels the denominator in Eq. (4.24), resulting in the standard constant bias, while the second yields the scale dependent correction. Putting everything together, we get

[TABLE]

Evidently, the coefficient of the scale dependent correction contains a summation over all the nonlinearity parameters, or better put, the whole function $F$. In order to see this, we can write it down, equivalently, as

[TABLE]

where $\mathcal{N}\equiv\frac{e^{-\phi_{L}^{2}/2\sigma_{L}^{2}}}{\sqrt{2\pi}\sigma_{L}}\sum_{m=1}^{\infty}\frac{\partial\log n_{L}}{\partial\beta_{m}}\frac{1}{\sigma_{L}^{m}}{\rm He}_{m}\!\left(\frac{\phi_{L}}{\sigma_{L}}\right)$. Consequently, a scale-dependent halo bias can only signal the presence of some form of local NG but not of a specific parametrisation of it.

A downside of this is that the derivatives of $n_{L}$ w.r.t. $\beta_{m}$ in Eq. (4.28) would have to be computed from simulations if one is to obtain information about $F$, or alternatively, assume a model of collapse into haloes. In the first case, one may compute the derivatives of the halo mass function by varying the initial conditions of the simulations. However, the coefficients (4.18) of the scale dependent bias cannot be observed in real data in this manner, i.e, by repeating the collapse process as one would do in simulations, because they only come once and with the same initial condition for the curvature perturbation field $\zeta$. Therefore, it becomes necessary that the next step be a connection between the effective bias coefficient $\int\mathcal{N}\frac{{\rm d}F}{{\rm d}\phi}{\rm d}\phi$ and observable quantities by modelling the halo mass functional $n_{L}$ in some way.

It turns out that the previous section 4.1 provides enough tools to accomplish this. In particular, if one models local NG at short scales $\phi_{S}$ using the extended PS scheme, as in Eq. (4.10), the variation with respect to the effective nonlinearity parameter ${\rm d}\log n_{L}/{\rm d}\beta_{m}$ may be computed directly as ${\rm d}\log n_{L}/{\rm d}f_{m}^{S}$, where $f_{m}^{S}$ is the local nonlinearity parameter of the short-scale theory (and, to first order in the perturbation $F$, it also defines the corresponding nonlinearity parameter of the full theory because of Eq. (2.27)). This can be done systematically using Eqs. (2.36), (4.9) and (4.27).

The first scale dependent NG correction $(m=1)$, associated to cubic NG, was computed in Ref. [103] assuming a universal mass function, in agreement with the result of Ref. [15]:

[TABLE]

where $b_{\rm G}\equiv\frac{\partial\log n_{L}}{\partial\delta_{L}}$ is the Gaussian bias. The second term $(m=2)$, corresponding to a $g_{\rm NL}$-type of NG, was computed in Ref. [94] using an Edgeworth expansion of the halo mass function and was found to be given by

[TABLE]

where $\kappa_{n}(M)$, is the reduced cumulant defined below Eq. (4.9). It turns out that when (4.27) is expanded up to $f_{3}=g_{\rm NL}$, the resulting expression is in very good agreement with simulations [94]. Derivatives of the halo mass function with respect to the higher NG nonlinearity parameters $f_{4},f_{5},\ldots,$ may be easily computed from Eqs. (2.36) and (4.9). Concretely, they are given by

[TABLE]

where primes denote derivatives w.r.t. the mass of the halo $M$. These terms can (and should) also be tested with N-body simulations for $f_{m>3}^{S}$ NG initial conditions, but such a computation is beyond the scope of this paper.

From Eq. (4.27) it is clear that the bias alone cannot differentiate between $f_{\rm NL},\;g_{\rm NL}$ or any higher order NG but it can give an answer to the question of whether Gaussianity is present or not171717A proposal for partially disentangling the contributions of different cumulants is the density-in-spheres 1-point PDF [88], which does so by taking into account both overdense and underdense regions, thus breaking the degeneracy between odd/even moments.. However, given a specific model for primordial NG, this can be a powerful probe in the sense of matching a template against data. For example, in App. A, we show that one particular realisation of this situation is the presence of an isocurvature mode with a potential $\Delta V(\psi)$. In this case, the function $F$ of the generalised local ansatz (1.5), is related to the potential as $F\propto\Delta V^{\prime}$. Thus, a measurement of or a constraint on the bias would translate into a constraint on the parameters of the landscape potential. In such a context, one can choose a well-motivated potential $\Delta V$, which fixes the function $F$, depending on few parameters (two or three cover most physically motivated potentials). Then, with the help of simulations, one can use the PDF (2.10) with the replacements (4.3) to draw appropriate NG initial conditions for the density fluctuation field, and then look for a bias of the form (4.27). For example, within the axion parameter space, there are regions that lead to such cases [48] and may be probed in the near future with LSS surveys.

5 Concluding remarks

Non-Gaussianity, notwithstanding how small, is a robust prediction of inflation, and may be detectable if the inflaton had the chance to interact with other degrees of freedom during inflation. Up until now, NG does not show up in low $n$-point correlation functions of the temperature map [4], at least not in the form of a bispectrum nor a trispectrum, when their estimators are compared to the well-motivated local, equilateral, folded and orthogonal templates. However, primordial NG might be hidden in the data in the form of patterns that need to be revealed through different estimators.

In this work, we have focused on the object that contains the full information about the distribution of anisotropies of the temperature and density fields, that is, the probability density function. Our starting point has been the bottom-up parametrisation of the curvature perturbation as a Gaussian random field plus an arbitrary analytic function thereof. Instead of truncating the series expansion of this function to the first few terms, corresponding to the standard $f_{\rm NL}$, $g_{\rm NL}$ parametrisation, we have kept the entire series enabling us to derive a probability functional that encodes the full local ansatz.

The various kinds of distributions derived from this functional provide us with new estimators, like the 1-point and 2-point statistics, which constitute complementary channels for the search of NG. In particular, the 2-point PDF encoding the scale dependence of the temperature distribution can serve as a full probe of the power spectrum able to disentangle local NG from other shapes. Moreover, the community has been shifting the focus towards LSS statistics as the latter will probe much more modes enhancing the statistical power of the datasets. With a view towards the near-future LSS surveys, we have computed the distribution of halos of mass $M$ at redshift $z$ showing that it is directly related to the primordial NG deviation. Finally, we computed the NG scale dependent correction to the linear halo bias arguing that it is the full local ansatz that contributes to the scaling. Therefore, the bias offers a unique probe of local NG in full generality, which cannot be used to distinguish the set of nonvanishing nonlinearity parameters, since these contribute as a sum to the correction.

Our results can be used with the future CMB/LSS data in various ways. For example, one can obtain refined constraints on the local ansatz from an Edgeworth expansion of the 2-point CMB temperature PDF. In addition, upon assuming well-motivated microphysical models leading to NG, one can use the formulas for the halo mass function and the bias as templates to be compared with data. Finally, we may use other types of statistical estimators like Minkowski functionals on the CMB/LSS datasets to draw conclusions on the primordial PDF. Let us, nevertheless, stress once more that, as discussed in Sec. 2.5, our treatment is strictly valid only within linear perturbation theory and as such it is more accurate in cases where primordial NG dominates over secondary effects, especially when one considers higher $n$-point functions, which are well-motivated, for instance, from multi-field effects [48, 49] or reheating models [44, 46].

Let us finish by pointing out that even though this work is focused on local primordial non-Gaussianity, it should be possible to extend the discussion to other cases. An “equilateral” ansatz of the form discussed in the Introduction would lead to a different scale dependence in the 2-point PDF, which could in turn serve as an equilateral NG estimator. Note that constraining such an object via e.g. an appropriate Edgeworth expansion could be less costly than that of $n$-point functions using the standard templates. Finally, one may be able to further generalise this by including combinations of spatial derivatives acting on the curvature fluctuation that lead to enfolded correlators. In this more general case, by the same procedure of changing variables in the Gaussian, one would be able to fully describe all the four momentum shapes employed at the level of the PDF.

Acknowledgements

We wish to thank Yashar Akrami, Diego Blas, Jonathan Braden, Xingang Chen, Rolando Dunner, Nelson Padilla, Domenico Sapone, Eva Silverstein and Cora Uhlemann for useful discussions and comments. We are also grateful to the organisers of the workshop Inflation & Geometry at IAP, Paris, for creating a vibrant atmosphere, which helped us shape part of this work, as well as to all the participants for extensive discussions and feedback. GAP and BSH acknowledge support from the Fondecyt Regular project number 1171811 (CONICYT). BSH is supported by a CONICYT grant number CONICYT-PFCHA/MagísterNacional/2018-22181513. SS is supported by the CUniverse research promotion project (CUAASC) at Chulalongkorn University.

Appendix A Generalised ansatz from quantum fluctuations during inflation

During inflation, the primordial curvature perturbation $\zeta$ is sourced by quantum fluctuations of the inflaton field or possibly other degrees of freedom such as isocurvaton fields. As a general statement, one can write down the field operator $\zeta({\bm{x}})$ at the final time slice $t$ as

[TABLE]

where $U$ is the temporal evolution operator in the interaction picture of quantum mechanics, and $\zeta_{I}$ the interaction-picture field, which follows the dynamics of the free theory.

Naturally, the field $\zeta$ will generate a specific set of $n$-point functions when computing expectation values. One can then construct a PDF ${\rm P}[\zeta]$ that generates these statistics through functional integration

[TABLE]

over the field configurations $\zeta({\bm{x}})$. Quantum mechanics does provide the tools to determine ${\rm P}[\zeta]$ directly, at least in principle. The operation $\zeta=U^{\dagger}\zeta_{I}U$ can be reframed in terms of a functional expression

[TABLE]

that depends on the whole spacetime evolution of the interaction-picture fields, $\zeta_{I}$ and other degrees of freedom $\psi_{i}^{I}$, and that of their conjugate momenta, $\Pi_{I}$ and $\Pi_{i}^{I}$. In principle, one can compute correlations directly from this expression. However, if the dependence of $\mathcal{O}$ on the interaction-picture fields is known, one can determine the PDF of $\zeta$ by integrating over all possible configurations:

[TABLE]

where ${\rm P}_{\rm G}$ is a Gaussian measure, with appropriate prescriptions to take into account the ordering of operators in $\mathcal{O}$. The measure is guaranteed to be Gaussian because the free fields evolve linearly in time, and therefore the contraction of quantum field operators obeys Wick’s theorem, which is equivalent to saying that the statistics are Gaussian.

Once ${\rm P}[\zeta]$ is obtained, the statistics of $\zeta$ is fully determined. Computationally, however, it is useful to have a probability distribution from which one knows how to obtain expectation values. On the other hand, one knows that the observed statistics for $\zeta$ are consistent with Gaussianity, and that deviations, if any, must be small. This motivates finding a functional map $\zeta_{\rm G}({\bm{x}})=G[\zeta]({\bm{x}})$, with inverse $\zeta({\bm{x}})=\zeta_{\rm G}({\bm{x}})+\mathcal{F}[\zeta_{\rm G}]({\bm{x}})$, such that $\zeta_{\rm G}$ has Gaussian statistics, i.e., such that

[TABLE]

The difficulty, of course, lies in finding such a mapping. Afterward, one can include another mapping, that makes the power spectrum of $\zeta_{\rm G}$ to be consistent with current observations, if this is not the case already.

Thus, we have justified writing $\zeta$ in terms of a Gaussian field $\zeta_{\rm G}$,

[TABLE]

where the functional ${\cal F}_{\rm NG}$ is, in principle, arbitrary, and should be determined from the specifics of the inflationary model at hand. We now give an example of a model that motivates searching for non-Gaussianity wherein

[TABLE]

which is the main focus of our work.

A.1 A concrete example: multi-field inflation

Our current understanding of fundamental theories, such as string theory and supergravity, requires us to take into consideration the existence of many scalar fields. In such a framework, it does not make much sense to talk about an inflaton field; instead, one has an inflationary path meandering through a landscape potential. Consequently, curvature fluctuations may have been coupled to (many) other dynamically active degrees of freedom inevitably sourcing departures from Gaussianity. Even in the simplest extension —from single-field to multi-field inflation— one finds a rich phenomenology that can be tested by near-future surveys. The presence of multiple fields during inflation leaves unique imprints in the CMB, a subject that has gained a lot of attention in the last years with quasi-single field inflation [114, 115] and the Cosmological Collider program [116]. The “smoking gun” signature of massive degrees of freedom active at the inflationary energy scale is a bispectrum of the local shape that peaks for triplets of modes with one momentum going to zero.

In the case of two-field models of inflation, the Lagrangian describing the dynamics of the curvature fluctuation $\zeta$ interacting with an isocurvature field $\psi$ is found to have the following generic form:

[TABLE]

In the previous expression, $a$ is the scale factor, $\epsilon=-\dot{H}/H^{2}$ is the first slow-roll parameter (where $H=\dot{a}/a$ is the Hubble expansion rate during inflation), $\mu$ is the so-called entropy mass of the isocurvature field $\psi$, and $\alpha$ is a coupling that appears as a consequence of the shape of the background inflationary path in the multi-field space. It parametrises the bending of inflationary paths in field space [117, 118, 119], with $\alpha=0$ corresponding to the case of geodesic trajectories (a straight line in field space). The piece $\mathcal{L}_{\rm int}$ contains operators of cubic order (and higher) in terms of the fields $\zeta$ and $\psi$. A nonvanishing value of $\alpha$ makes $\psi$ act as a source for the amplitude of $\zeta$ via the quadratic mixing term

[TABLE]

Crucially, this mixing allows the transfer of non-Gaussian statistics from the isocurvature field to the curvature fluctuation (but not the other way around [120]). The generated non-Gaussianity will depend on the particular form of $\mathcal{L}_{\rm int}$.

A standard approach to deal with (A.8) consists in assuming that $\mathcal{L}_{\rm int}$ can be organised perturbatively in terms of powers of both fields $\zeta$ and $\psi$. Such an approach supposes that the nonlinear dynamics is dominated by cubic operators $\mathcal{L}_{\rm int}\supset\mathcal{L}_{\rm int}^{(3)}$, and implies that the main non-Gaussian departures will be parametrised by the bispectrum, the amplitude of the 3-point function in momentum space [26, 27, 28, 29, 30, 31, 32, 33]. Within this approach, one finds that the main operators leading to a nonvanishing bispectrum are given by a self-interaction term and a mixing term of the forms $\mathcal{L}_{\rm self}^{(3)}\propto a^{3}g\psi^{3}$ and $\mathcal{L}_{\rm mix}^{(3)}\propto a^{3}\alpha\dot{\zeta}\psi^{2}$, respectively. Any other term in $\mathcal{L}_{\rm int}^{(3)}$ turns out to be suppressed by powers of the slow-roll parameters. Both interaction terms lead to a non-analytic behaviour of the bispectrum: in the $\mu/H<3/2$ regime, one obtains an NG amplitude enhanced by $g$, with the bispectrum following an irrational power law in the squeezed limit leading to an intermediate shape between the equilateral and local ones [114, 115]; the other region $\mu/H>3/2$ leads to a bispectrum with a shape sensitive to the value of the entropy mass, displaying oscillatory patterns [121, 122, 123, 124], while in the $\mu/H\gg 3/2$ limit the heavy field may be fully integrated out, leading to equilateral shape NG [125].

However, one may conceive regimes characterised by interaction Lagrangians $\mathcal{L}_{\rm int}$ in which one is not allowed to disregard terms of higher powers in the fields. For instance, the field $\psi$ may have a potential $\Delta V(\psi)$ displaying a rich structure within a wide field range $\Delta\psi$ [126]. In this case, the Lagrangian (A.8) can be rewritten as

[TABLE]

where we have absorbed the mass term in $\Delta V(\psi)$ by means of the identification $\mu^{2}\equiv\Delta V^{\prime\prime}(0)$, where primes denote derivatives with respect to $\psi$. The elipses $\cdots$ denote terms that are either suppressed by slow-roll parameters or by $\alpha/H$, which is assumed to be small. This Lagrangian (A.10) describes the dynamics of perturbations in those cases where the landscape potential has a rich structure in directions orthogonal to the inflationary path. If the potential remains shallow, i.e. $\Delta V/H^{4}\ll 1$, the field $\psi$ will display characteristic fluctuations $\Delta\psi$ that could traverse many minima and maxima of the potential, therefore probing the structure of the landscape. In this case, one cannot just stick to computations involving the first two terms of the expansion $\Delta V=\frac{\mu^{2}}{3}\psi^{2}+\frac{g}{3!}\psi^{3}+\cdots$; instead, one has to perform perturbative computations where information about the entire potential $\Delta V$ is kept under control. As shown in [48, 49], as long as $\Delta V/H^{4}\ll 1$, it is indeed possible to compute every $n$-point correlation function for $\zeta$ and derive its probability distribution function, ${\rm P}(\zeta)$. Remarkably, this PDF turns out to be given by a Gaussian profile with small non-Gaussian corrections determined by the shape of the landscape potential $\Delta V$. Thanks to this property, it is in principle possible to reconstruct the shape of the landscape potential $\Delta V$ out of CMB observations [49], providing information about a section of $\Delta V$ during the period of inflation when the modes of wavelengths relevant to the CMB were exiting the horizon (hence the term, tomographic non-Gaussianity) and beyond.

The generation of tomographic non-Gaussianity may be traced back to the self-interactions of an isocurvature field $\psi$. These self-interactions are transferred to the curvature perturbations thanks to a linear interaction coupling both fields [see Eq. (A.10)]. To appreciate how self-interactions give rise to this form of non-Gaussianity, we may start by recalling that in the interaction picture, the evolution of the field $\psi({\bm{x}},\tau)$ is given by

[TABLE]

From a quantum-mechanical perspective in the Heisenberg picture, one can write

[TABLE]

Similarly, the interaction picture gives that to first order in the interaction the field is given by

[TABLE]

where the subscript $I$ informs us that the corresponding operator is in the interaction picture and evolves as a free field. If, for simplicity, we take de Sitter spacetime as a background, the interaction-picture Hamiltonian reads $H_{I}(\tau)=\int_{x}a^{4}(\tau)\Delta V(\psi_{I}({\bm{x}},\tau))$ with $a(\tau)=-1/H\tau$. With this, we may work on our previous equation to obtain

[TABLE]

With the help of canonical commutation relations for the appropriate field $v\equiv a\psi$ [120, 48], the commutator $[\psi_{I}({\bm{x}}^{\prime},\tau^{\prime}),\psi_{I}({\bm{x}},\tau)]$ is just a number and we can carry out the integral over $\tau^{\prime}$ explicitly provided that:

The quantum field $\psi_{I}({\bm{x}}^{\prime},\tau^{\prime})$ in the argument of $\partial\Delta V/\partial\psi$ may be treated as a constant over time $\tau^{\prime}$. 2. 2.

The range of modes under consideration involves superhorizon modes with $|k\tau^{\prime}|\lesssim 1$.

In fact, if the range of modes satisfies $|k\tau^{\prime}|\ll 1$ for all $k$, then the first condition is implied by the second.

Let us give some comments about these conditions: The first point seems natural in the sense that the statistics of $\psi_{I}$ do not evolve over time: it may be seen as a Gaussian random field with a definite covariance. The second point, although appealing, is both physically and mathematically suspect, since in principle the interaction-picture Hamiltonian involves every mode (i.e. every momentum scale). Nonetheless, from an EFT perspective this is perfectly acceptable, as long as the potential $\Delta V$ is responsible for describing the physics at those scales. Moreover, this is the appropriate course of action when studying CMB or LSS modes that spent a large number of e-folds outside the horizon, because they do satisfy $|k\tau^{\prime}|\ll 1$ throughout most of their history (practically for every time after horizon crossing this condition is fulfilled).

Using these considerations, one obtains

[TABLE]

where $\Delta N$ is the number of e-folds spent outside the horizon by the range of modes under consideration, which we take to satisfy $\Delta N\gg 1$. While the former may seem to be a heavy restriction, if we consider that the currently observable range of scales in the CMB satisfies $\ln(k_{S}/k_{L})\simeq 8$ [3] and $\Delta N\sim 60$, we see that approximating $\Delta N$ by a single value for all modes in the considered range is justified. Of course, as was earlier suggested, a rigorous proof of this requires to go through every $n$-point function, performing an adequate renormalisation procedure to make the computations consistent at every bounded range of momenta. This was thoroughly dealt with in previous works [48, 49], although stopping just short of writing down Eq. (A.15).

The curvature perturbation $\zeta({\bm{x}},\tau)$ may be obtained in a completely analogous manner, only that we now have to consider an extra commutator to account for the quadratic mixing term (A.9):

[TABLE]

where $H^{V}_{I}$ is the term of the interaction-picture Hamiltonian that contains the $\psi$ self-interactions and $H^{\alpha}_{I}$ contains terms associated to the quadratic mixing. The ellipsis $\cdots$ stand for higher order terms. It is of crucial importance to obtain the correct result to notice that the commutators in the last term of (A.16) only give a nonzero result when the pieces of the interaction Hamiltonian are written in that order. This, alongside the time ordering, yields an additional $1/2$ factor for the statistical transfer of the nonlinear perturbation $\Delta V$. After a calculation analogous to the one that led us to (A.15), with the same working assumptions, one obtains

[TABLE]

This equation can accommodate a variety of regimes. However, we choose to work in a situation wherein the linear transfer from $\psi$ to $\zeta$ dominates. Even though this last equation is a perturbative result, $\zeta=\frac{\alpha\Delta N}{H}\psi+\zeta_{0}$ is an exact solution of the equations of motion on superhorizon scales [120]. This allows us to neglect the first term. Furthermore, since conventionally $\zeta$ should satisfy $\langle\zeta\rangle=0$, from this point forward we will consider

[TABLE]

or equivalently, omitting the temporal coordinate and the projection integration $\int_{\bm{y}}\int_{{\bm{k}}}e^{i{\bm{k}}\cdot({\bm{x}}-{\bm{y}})}$,

[TABLE]

where we have written $\zeta_{\rm G}$ instead of $\psi_{I}$ to stress the nature of our result: $\zeta$ is made up from a Gaussian contribution plus a local non-Gaussian term. This result has the desired form $\zeta({\bm{x}})=\zeta_{\rm G}({\bm{x}})+{\cal F}_{\rm NG}[\zeta_{\rm G}]({\bm{x}})$. To identify ${\cal F}_{\rm NG}$ it is convenient to recognise that, because at linear order $\zeta\simeq\frac{\alpha\Delta N}{H}\psi$, the power spectrum of $\zeta$ satisfies $P_{\zeta}=\frac{\alpha^{2}\Delta N^{2}}{H^{2}}P_{\psi}$, where $k^{3}P_{\psi}/2\pi^{2}=H^{2}/4\pi^{2}$ (because in the free theory $\psi$ behaves as a massless field). This allows one to find [49] $\alpha^{2}\Delta N^{2}=4\pi^{2}A_{s}$ where $A_{s}=k^{3}P_{\zeta}/2\pi^{2}$ is the amplitude of the power spectrum $P_{\zeta}$ of $\zeta$. This finally leads to

[TABLE]

which is the result181818A similar expression with a generic NG deviation parametrised by some function was written in Refs. [44, 46] in the context of preheating. Here we see that another natural interpretation of such non-Gaussianity is an isocurvature potential, with the trigonometric case corresponding to an axion [48]. reported in the Introduction.

Appendix B The partition function

The defining property of a partition function is that upon functional differentiation as in Eq. (2.21), it should give the $n$-point functions:

[TABLE]

This is accomplished by the following functional:

[TABLE]

To evaluate this expression, let us use the expression for ${\rm P}[\zeta]$ as given in Eq. (2.8):

[TABLE]

We will now evaluate (B.2) by expanding $F$ in a power series and using Wick’s theorem. Since (B.2) may be read as computing the expectation value of

[TABLE]

over a Gaussian measure, we will do exactly that. Therefore, we will have to compute

[TABLE]

where the subscript $G$ instructs to take the expectation value over a Gaussian measure. Note that the second expectation value, per Wick’s theorem, can be written as the sum of two expectation values

[TABLE]

Now, using that $\Sigma({\bm{x}},{\bm{z}})=\int_{\bm{q}}P_{\zeta}(q)e^{i{\bm{q}}\cdot({\bm{x}}-{\bm{z}})}$, and replacing the last term into the corresponding term of (B.4), we get

[TABLE]

Then, integrating over ${\bm{x}}$ gives $|{\bm{q}}|=|{\bm{k}}|$, and thus yields

[TABLE]

which is equal but opposite in sign to the first term of (B.5). Therefore, those two cancel out, and we only have to compute

[TABLE]

If we define

[TABLE]

then the object of interest in the computation is

[TABLE]

There are three type of contractions in this expression: two self-contractions of fields originating from equivalent expressions and a mixed one. Performing the combinatorics gives

[TABLE]

which one can sum over $n,m$ to eliminate the constraints on the $m^{\prime},n^{\prime},\ell^{\prime}$ sums. This gives

[TABLE]

Identifying $\sigma_{\zeta}^{2}=\Sigma({\bm{z}},{\bm{z}})$ and using the Weierstrass transform, we get

[TABLE]

Finally, we use that

[TABLE]

to get, after integration by parts,

[TABLE]

which gives

[TABLE]

as shown in the main text.

Appendix C Details of the fixed-point PDFs

C.1 2-Point PDF: general case

In the main text we argued that the 2-point PDF can be deduced from the probability density functional upon conditioning in two points. Here we arrive at the same result starting from the corresponding 2-point correlators and deriving the PDF. To our knowledge, this has not been derived earlier and thus we outline the procedure in some detail.

As in the previous case, the PDF must include the non-fully connected contributions. This is a combinatorial mess, for we expect

[TABLE]

Let us calculate $\#_{n_{1},n_{2},m_{1},m_{2},m_{t}}$: we have an overall factor $n_{1}!n_{2}!$ from which we must divide out the overcounted terms. In this counting, we have $m_{t}!$ redundant permutations when connecting ${\bm{x}}_{1}$ and ${\bm{x}}_{2}$, plus $2^{m_{1}}m_{1}!2^{m_{2}}m_{2}!$ when pairing amongst themselves. Finally, the ones that are assigned to the fully connected contribution undergo no further permutation, thus we must also divide by $(n_{1}-2m_{1}-m_{t})!(n_{2}-2m_{2}-m_{t})!$. Thus,

[TABLE]

Careful inspection of this result reveals that (C.2) leads to a 2-point distribution analogous to what was obtained in [49], but with two points defining the filtering instead of one:

[TABLE]

where the coefficients $G_{ij}$, defined right below, depend on both the field variables and the spacetime positions ${\bm{x}},{\bm{x}}_{1},{\bm{x}}_{2}$ via the scalar variables

[TABLE]

In obtaining this expression, we have defined a number of functions that depend uniquely on the structure of the Gaussian theory. The variance $\sigma_{\zeta}^{2}(r,r_{1},r_{2})$ and $\zeta(r,r_{1},r_{2})$ are the regression coefficients obtained by conditioning a Gaussian distribution of $(\bar{\zeta},\zeta_{1},\zeta_{2})$ over $(\zeta_{1}^{W},\zeta_{2}^{W})$, with covariance matrix given by (2.43):

[TABLE]

Furthermore, the coefficients $G_{ij}$ are given by

[TABLE]

C.2 2-Point PDF: map to the sphere $S^{2}$

The only difference with the previous Appendix is that the window function should be a map onto an angular coordinate $\hat{n}$ instead of a three-dimensional (flat) space. In the text, we wrote

[TABLE]

Here we have

[TABLE]

while the regression coefficients are given by

[TABLE]

Finally, the functions $G_{ij}^{\Theta}$ read

[TABLE]