The Hunt for Primordial Interactions in the Large Scale Structures of the Universe

TL;DR

This paper reviews recent theoretical advances in understanding how primordial interactions during inflation influence the large-scale structure of the universe, especially the distribution of dark matter halos.

Contribution

It provides a comprehensive summary of models linking inflationary primordial interactions to observable large-scale structures, emphasizing the statistical distribution of dark matter halos.

Findings

Primordial interactions leave detectable imprints on large-scale structures.

A consistent framework connects inflationary correlations to halo distributions.

Insights could inform fundamental physics at high energies.

Abstract

The understanding of the primordial mechanism that seeded the cosmic structures we observe today in the sky is one of the major goals in cosmology. The leading paradigm for such a mechanism is provided by the inflationary scenario, a period of violent accelerated expansion in the very early stages of evolution of the Universe. While our current knowledge of the physics of inflation is limited to phenomenological models which fit observations, an exquisite understanding of the particle content and interactions taking place during inflation would provide breakthroughs in our understanding of fundamental physics at high energies. In this review, we summarize recent theoretical progress in the modelling of the imprint of primordial interactions in the large scale structures of the Universe. We focus specifically on the effects of such interactions on the statistical distribution of dark matter halos, providing a consistent treatment of the steps required to connect the correlations generated among fields during inflation all the way to the late-time correlations of halos.

Tables

Table 1. Table 1: List of values of A 2 ​ J subscript 𝐴 2 𝐽 A_{2J} from the squeezed limit expansion of Eq. ( 150 ) for interactions produced in single-field models of inflation.

Interactions in single-field models	$A_0$	$A_2$
Gravitational Floor (Eq. (2.2.1))	$2ϵ+\eta +\left(\frac{29}{6}ϵ+\frac{1}{4}\eta \right){\left(\frac{k_{\ell }}{k_s}\right)}^2$	$\left(\frac{1}{6}ϵ+\frac{5}{4}\eta \right){\left(\frac{k_{\ell }}{k_s}\right)}^2$
DBI (Eq. (39))	$\left[\frac{19}{12}\left(1-\frac{1}{c_s^2}\right)+\frac{3}{4}\left(\frac{1}{c_s^2}-1-\frac{2\lambda }{\Sigma }\right)\right]{\left(\frac{k_{\ell }}{k_s}\right)}^2$	$\frac{5}{6}\left(1-\frac{1}{c_s^2}\right){\left(\frac{k_{\ell }}{k_s}\right)}^2$
Solid Inflation (Eq. (54))	$0$	$-\frac{40}{9}\frac{F_Y}{F}\frac{1}{ϵ}\frac{1}{c_L^2}$

Table 2. Table 2: List of values of A 2 ​ J subscript 𝐴 2 𝐽 A_{2J} from the squeezed limit expansion of Eq. ( 150 ) for interactions produced in multi-field models of inflation. We have parametrized the amplitude of interactions involving the exchange of massive particles using α s subscript 𝛼 𝑠 \alpha_{s} and r s subscript 𝑟 𝑠 r_{s} , where s 𝑠 s indicates the spin of the particle, for principal and complementary series respectively. Full expressions for these amplitudes can be found in Lee et al. ( 2016 ) .

Interactions in multi-field models	$A_0$	$A_{2J}$
Local (Eq. (66))	$\frac{12}{5}{f}_{\mathrm{NL}}^{\mathrm{loc}}+\frac{3}{5}{f}_{\mathrm{NL}}^{\mathrm{loc}}{\left(\frac{k_{\ell }}{k_s}\right)}^2$	$3f_{\mathrm{NL}}^{\mathrm{loc}}{\left(\frac{k_{\ell }}{k_s}\right)}^2$
Massive particles
Even Spin (Principal Series)	${\alpha }_0{\left(\frac{k_{\ell }}{k_s}\right)}^{3/2}\cos \left[\mu \ln \left(\frac{k_{\ell }}{k_s}\right)\right]$	${\alpha }_{2J}{\left(\frac{k_{\ell }}{k_s}\right)}^{3/2}\cos \left[{\mu }_{2J}\ln \left(\frac{k_{\ell }}{k_s}\right)+{\varphi }_{2J}\right]$
Even Spin (Complementary Series)	$r_0{\left(\frac{k_{\ell }}{k_s}\right)}^{3/2-\nu }$	$r_{2J}{\left(\frac{k_{\ell }}{k_s}\right)}^{3/2-{\nu }_{2J}}$
Odd Spin (Principal Series)	${\alpha }_s{\left(\frac{k_{\ell }}{k_s}\right)}^2$	$𝒪\left[{\left(\frac{k_{\ell }}{k_s}\right)}^{5/2}\right]$
Odd Spin (Complementary Series)	${\alpha }_s{\left(\frac{k_{\ell }}{k_s}\right)}^2+r_s{\left(\frac{k_{\ell }}{k_s}\right)}^{5/2-{\nu }_s}$	$𝒪\left[{\left(\frac{k_{\ell }}{k_s}\right)}^{5/2}\right]$

Equations

$$S=\frac{M_P}{2}\int {\mathrm{d}}^{4}x\sqrt{-g}R+\int {\mathrm{d}}^{4}x\sqrt{-g}\left[-\frac{1}{2}{g}^{\mu \nu }{\partial }_{\mu }\varphi {\partial }_{\nu }\varphi -V(\varphi )\right],$$

$$\mathrm{d}{s}^2=-\mathrm{d}{t}^2+a^2(t)\mathrm{d}{\mathbf{x}}^2,$$

$$H^2=\frac{1}{3M_P^2}\left(\frac{1}{2}{\dot{\varphi }}_0^2+V({\varphi }_0)\right),$$

$$\dot{H}=-\frac{{\dot{\varphi }}_0^2}{2M_P^2},$$

$${\ddot{\varphi }}_0+3H{\dot{\varphi }}_0+V^{\prime }=0.$$

$$N(\varphi )\equiv {\int }_a^{a_e}\mathrm{d}\ln a={\int }_t^{t_e}H\mathrm{d}t$$

$$ϵ\equiv -\frac{\dot{H}}{H^2}\ll 1,\eta \equiv \frac{\dot{ϵ}}{ϵH}<1$$

$${ϵ}_V\equiv \frac{M_P^2}{2}{\left(\frac{V^{\prime }}{V}\right)}^2<1,{\eta }_V\equiv M_P^2\left|\frac{V^{\prime \prime }}{V}\right|<1.$$

$$\varphi (\mathbf{x},t)={\varphi }_0(t)+\delta \varphi (\mathbf{x},t).$$

$$S=\frac{1}{2}\int \mathrm{d}\tau {\mathrm{d}}^{3}x{a}^2\left[({\partial }_{\tau }\delta \varphi {)}^2-({\partial }_i\delta \varphi {)}^2\right],$$

$$\delta {\varphi }^{\prime \prime }(\mathbf{k},\tau )+2\mathcal{H}\delta {\varphi }^{\prime }(\mathbf{k},\tau )+k^2\delta \varphi (\mathbf{k},\tau )=0,$$

$$u(\mathbf{k},\tau )=c_1\frac{H}{\sqrt{2k^3}}(1+ik\tau )e^{-ik\tau }+c_2\frac{H}{\sqrt{2k^3}}(1-ik\tau )e^{ik\tau }.$$

$$[\delta \hat{\varphi }(\mathbf{x},\tau ),\delta \hat{\pi }(\mathbf{y},\tau )]$$

$$[\delta \hat{\varphi }(\mathbf{x},\tau ),\delta \hat{\varphi }(\mathbf{y},\tau )]$$

$$\delta \hat{\varphi }(\mathbf{k},\tau )$$

$$\delta \hat{\pi }(\mathbf{k},\tau )$$

$$[a_{\mathbf{p}},a_{-\mathbf{q}}^{†}]$$

$$[a_{\mathbf{p}},a_{-\mathbf{q}}]$$

$$|c_1{|}^2-|c_2{|}^2=1$$

$$u(\mathbf{k},\tau )=\frac{H}{\sqrt{2k^3}}(1+ik\tau )e^{-ik\tau }.$$

$$\zeta \approx H\delta t\approx H\frac{\delta \varphi }{{\dot{\varphi }}_0},$$

$$⟨0|\zeta ({\mathbf{k}}_1,0)\zeta ({\mathbf{k}}_2,0)|0⟩=(2\pi {)}^3{\delta }^{(3)}({\mathbf{k}}_1+{\mathbf{k}}_2)P_{\zeta }(k_1).$$

$${\mathcal{P}}_{\zeta }=\frac{H_{\ast }^2}{8{\pi }^2{M}_P^2{ϵ}_{\ast }}{\left(\frac{k}{k_{\ast }}\right)}^{n_s-1},$$

$$n_s-1\equiv \frac{\mathrm{d}\ln P_{\zeta }}{\mathrm{d}\ln k}=-2ϵ-\eta .$$

$${\mathcal{A}}_s$$

$$n_s$$

$$⟨\hat{Q}(\tau )⟩=⟨\Omega |\hat{Q}(\tau )|\Omega ⟩,$$

$$⟨\hat{Q}(\tau )⟩=⟨0|\bar{T}{e}^{i{\int }_{-\infty (1-iϵ)}^{\tau }{\hat{H}}_{\mathrm{int}}^I({\tau }^{\prime })\mathrm{d}{\tau }^{\prime }}{\hat{Q}}^I(\tau )T{e}^{i{\int }_{-\infty (1+iϵ)}^{\tau }{\hat{H}}_{\mathrm{int}}^I({\tau }^{\prime \prime })\mathrm{d}{\tau }^{\prime \prime }}|0⟩,$$

$$⟨{\zeta }_{{\mathbf{k}}_1}{\zeta }_{{\mathbf{k}}_2}{\zeta }_{{\mathbf{k}}_3}⟩\approx -i{\int }_{-\infty }^0\mathrm{d}\tau ⟨0|[{\hat{\zeta }}_{{\mathbf{k}}_1}{\hat{\zeta }}_{{\mathbf{k}}_2}{\hat{\zeta }}_{{\mathbf{k}}_3}(0),{\hat{H}}_{\mathrm{int}}(\tau )]|0⟩,$$

$$B_{\zeta }({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3)\equiv ⟨{\zeta }_{{\mathbf{k}}_1}{\zeta }_{{\mathbf{k}}_2}{\zeta }_{{\mathbf{k}}_3}⟩.$$

$$B_{\zeta }({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3)=(2\pi {)}^3{\delta }^{(3)}({\mathbf{k}}_1+{\mathbf{k}}_2+{\mathbf{k}}_3)\frac{(2{\pi }^2{)}^2}{(k_1{k}_2{k}_3{)}^2}{\mathcal{P}}_{\zeta }^2(k_{\ast })\mathcal{S}(k_1,k_2,k_3),$$

$$\mathcal{S}(k_1,k_2,k_3)\stackrel{k_1=k_2=k_3}{\to }\frac{18}{5}{f}_{\mathrm{NL}}(K),$$

$${\mathcal{S}}_{\mathrm{GF}}(k_1,k_2,k_3)=$$

Full text

\pubvolume

7 \issuenum3 \articlenumber5

\historyReceived: 31 May 2019; Accepted: 7 August 2019; Published: 8 August 2019 \TitleThe Hunt for Primordial Interactions in the Large Scale Structures of the Universe \AuthorMatteo Biagetti †\orcidA \AuthorNamesMatteo Biagetti

1. 1 Inflation and primordial perturbations

2. 1.1 Background evolution

3. 1.2 Quantum fluctuations of the inflaton

4. 2 Interactions from inflationary models

5. 2.1 Interactions as non-Gaussianities

6. 2.2 Interactions in models of inflation

7. 2.2.1 Interactions in single-field models

8. 2.2.2 Multi-field interactions during inflation

9. 2.2.3 Multi-field interactions after inflation

10. 2.3 Final remarks of this section

11. 3 From primordial interactions to matter overdensities

12. 3.1 The density probability distribution

13. 3.2 Dark matter correlation functions

14. 4 Imprints of primordial interactions on one-point halo statistics

15. 4.1 Analytic approaches

16. 4.1.1 The Press-Schechter mass function

17. 4.1.2 The excursion set approach

18. 4.1.3 The excursion set peaks model

19. 4.2 Numerical approaches

20. 4.3 Final remarks of this section

21. 5 Imprints of primordial interactions on two-point halo statistics

22. 5.1 The breakthrough: scale-dependent bias from local PNG

23. 5.1.1 Derivation with the peak-background split ansatz

24. 5.1.2 Universality of the mass function

25. 5.1.3 The single-field consistency relation

26. 5.2 Analytic approaches

27. 5.2.1 The squeezed limit of the primordial bispectrum

28. 5.2.2 The bias expansion

29. 5.2.3 Models of halo clustering

30. 5.3 Numerical approaches

31. 5.3.1 Local-type PNG

32. 5.3.2 N-body simulations with generic non-Gaussian initial conditions

33. 5.4 Final remarks of this section

34. 6 Observational prospects

35. A ESP model

36. A.1 The curvature function of density peaks

37. A.2 The non-Gaussian ESP mass function

Introduction

Cosmological observations reveal a Universe filled with structures over a wide range of scales. The last three decades of research in the field of cosmology have seen a huge development in our understanding of how these cosmic structures were formed throughout the history of the Universe. The current standard cosmological model is able to make a consistent timeline of the dynamical evolution of the Universe over the course of $13$ billion years. During the earliest known stage of cosmic evolution, an accelerated expansion phase known as cosmic inflation Guth (1981); Starobinsky (1980); Linde (1982); Albrecht and Steinhardt (1982), primordial perturbations are believed to have formed , providing the seed for the formation of all structures at later times. In the inflationary scenario, these perturbations come from quantum fluctuations of scalar fields in an expanding background Mukhanov and Chibisov (1981). As they are produced, they are stretched to very large scales, outside of the causal horizon, where they remain frozen. In subsequent stages of the cosmic evolution, these perturbations reenter the horizon, providing small inhomogeneities over the whole universe. The inhomogeneities grow due to gravitational instability and form the cosmic structures that we observe today.

The characteristics of the primordial perturbations have been best constrained by the statistical analysis of the temperature anisotropies in the Cosmic Microwave Background (CMB), the relic light that decoupled from all interactions in the moment in which electrons and protons combined to form neutral hydrogen atoms, 380,000 years after inflation. Observations of the CMB temperature anisotropies require that i) super-horizon, ii) nearly scale invariant, iii) very close to Gaussian and iv) adiabatic perturbations are produced in the early Universe Akrami et al. (2019). These features are strong hints that an inflationary mechanism indeed took place111We should not forget however that alternative models to inflation are also able to pass these tests (see for instance Ijjas and Steinhardt (2016) for a recent comparison of ekpyrotic models Khoury et al. (2001) with Planck data).. Despite the great success of observing such features, it is somewhat underwhelming to realize that the constraints from the CMB power spectrum are far from restrictive on inflationary models. Indeed, hundreds of models of inflation exist which are able to satisfy these constraints.

A deeper understanding of the physics of inflation would fill our knowledge of the Universe up to $10^{-30}$ seconds after the Big Bang. The theoretical and observational challenges pertaining such a search are somewhat hard to tackle. Other observational windows into the primordial Universe, such as primordial nucleosynthesis (happened around $3$ minutes after the Big Bang) and recombination ($\sim 380,000$ years), are governed by the laws of nuclear and atomic physics which have been established and extensively tested on Earth over the last century. On the other hand, the processes involved during inflation are still modeled using mostly simplified phenomenological mechanisms, because the typical energy scale at which inflation takes place is far above the TeV scale, thus inaccessible to experiments on Earth. Moreover, the inflationary environment most probably does not involve fields and interactions of our standard model of particle physics, the elementary particles we are made of being generated at a later stage of cosmic evolution. The bright side is that any new discovery is a clear window on new physics, being able to explore energies as high as $10^{14}$ GeV, just a couple of orders of magnitude away from the Planck energy scale, thus providing us the best hope of experimentally probing quantum gravity.

A strong prediction which is common to all inflationary constructions is the production of primordial gravitational waves. Unfortunately, we have not been able to observe them yet. The Planck data, combined with the BICEP/Keck array measurements, constrain primordial gravitational waves perturbations to have amplitude $6.4\%$ smaller than scalar ones Akrami et al. (2018). While a minimal amount of gravitational waves is necessarily produced by any model of inflation, lower bounds can be unobservably low. For instance, if the amplitude of tensor fluctuations scales as the fourth root of the energy scale of inflation Lyth (1997), as expected for a large class of models of inflation, primordial gravitational waves might be as far as $\sim 40$ orders of magnitude away from current limits222This limit is calculated by imposing the minimal reheating temperature which realizes successful big bang nucleosynthesis..

Beside the existence of primordial gravitational waves, a great deal of information is still hidden in the statistics of scalar perturbations. Indeed, the characteristics listed above concern only the power spectrum, which tells us about the underlying free-field theory of the inflaton. On the other hand, attempts at building inflationary models from string theory and even particle physics perspective are usually characterized by a very rich particle content and interactions (see Lyth and Riotto (1999); Baumann and McAllister (2015) for reviews). Current limits on higher-order statistics impose interactions among fields to be rather weak, at least on the inflationary direction in field space Akrami et al. (2019). Nevertheless, there is still plenty of room to explore possibilities. Moreover, even in the simplest model of single-field inflation a minimal coupling to gravity is present. This is usually called gravitational floor of interactions and it is a “must observe” feature of all inflationary models. Although the detection of such an imprint is still at least two orders of magnitudes away from the sensitivity of current experiments, it does represent a guaranteed discovery which should be a main focus in cosmological searches.

Given the importance of such a search, there has been a growing effort in finding new observables that can constrain inflation. A promising observational probe of inflation is the study of structure formation at large scales and late times. Primordial perturbations provide the initial conditions with which matter overdensities grew under the effect of gravitational instability and formed all the structures in the Universe. It is therefore natural to hope to extract information about inflation by studying how matter is distributed in the Universe. This task is complicated by the fact that gravitational instability is a non-linear process and most of the information on initial conditions is washed out once away from the linear regime. This is the reason why CMB searches for primordial non-Gaussianities have dominated the efforts of the last two decades: at the last scattering surface, perturbations are still mostly linear. However, recent developments in large-scale structure (LSS) theory and observations have demonstrated that LSS can provide better constraints than the CMB in the near future. A recent analysis showed the first such example: using Baryon Acoustic Oscillations data from the BOSS collaboration Beutler et al. (2017), the authors of Beutler et al. (2019) have put constraints on primordial features that are stronger than CMB ones Akrami et al. (2019). For primordial non-Gaussianity, results from eBOSS collaboration have given the most recent constraints on local type primordial non-Gaussianity Castorina et al. (2019) and are expected to reach CMB sensitivity including the full set of data in the next few months. In the close future, the recently funded SPHEREx mission is one of the most promising examples of this great potential Doré et al. (2014), along with Euclid Amendola et al. (2013), LSST Abate et al. (2012) and SKA Bacon et al. (2018). A strong feature of LSS observations is the fact that they explore a tridimensional volume, as opposed to the bidimensional photon last scattering surface probed by the CMB. The number of modes available to constrain the statistical distribution of perturbations is therefore greatly enhanced. Moreover, a variety of different probes can be exploited: the statistical distribution of galaxies can be mapped through spectroscopic and photometric surveys, while weak gravitational lensing can trace the dark matter distribution directly. In addition, galaxy intrinsic alignments can also provide constraints on primordial scalar and tensor perturbations Chisari and Dvorkin (2013); Chisari et al. (2016); Kogai et al. (2018). Higher-redshift probes include the 21-cm neutral hydrogen line Joudaki et al. (2011); Chongchitnan and Silk (2012); Chongchitnan (2013); Camera et al. (2013); Sekiguchi et al. (2019); Witzemann et al. (2019) and intensity mapping with other emission lines Moradinezhad Dizgah et al. (2019); Moradinezhad Dizgah and Keating (2019). While experimental efforts on these directions are still not as developed as for lower-redshift ones, their potential is huge given that they explore various ranges in red-shift and therefore observed volume.

Plan of the Review. In this review, we want to summarize recent theoretical progress in the modelling of the imprint of primordial interactions, taking place during inflation, in the clustering of dark matter halos. We follow primordial perturbations in chronological order from the very early Universe to the present day, considering their impact in the statistics of rare objects and of the lowest order clustering statistics. Particular attention is devoted to reviewing the latest theoretical developments in characterizing interactions taking place during, or right after, inflation as a source of primordial non-Gaussianity. We will therefore start by providing a quick overview of the inflationary mechanism with a focus on the generation of primordial perturbations in Section §1. The following Section §2 investigates interactions taking place during inflation and their observable signature in the three-point function of the primordial curvature perturbation. Section §3 is a bridge connecting the primordial curvature perturbation to perturbations in the dark matter density field. We then proceed with the two main sections of the review: the imprint of primordial interactions in one- and two-point statistics of dark matter halos, Sections §4 and §5 respectively. We conclude by providing an overview of observational prospects in Section §6.

Not in this Review. Before we begin, we find it useful to briefly stress what this review will not be about. On the inflationary side, this review will deal with interactions among fields, or self-interactions, which generate a non-zero three-point correlation function of primordial perturbations, neglecting all higher-order correlations which are also generated. This choice is partially justified by the fact that we expect most models to respect a perturbative expansion, and therefore higher-order correlations to be more and more suppressed. On the other hand, there are various examples in which the study of the four-point correlation function is interesting and can lead to observable imprints. Another limitation on the inflationary side is that, even restricting to interactions generating three-point correlations only, the list of models investigated here is not exhaustive. On the LSS side, the main limitation of ths review is that it does not deal with the real observable, i.e. galaxies and clusters of galaxies. While luminous galaxies are the main observable tracers of the dark matter statistical distribution, dark matter halos provide the building blocks for their formation. Their study is therefore a crucial step for connecting theoretical predictions from inflation to observable imprints in LSS. It is clearly not the last, nor the only step: indeed, this review will not deal with a number of interesting issues, such as, for instance, the way in which galaxies populate halos. Even within the modeling of dark matter halos, we will restrict mostly to the analytic treatment of the evolution of perturbations, from inflation all the way to the present Universe, although short paragraphs will be devoted to recent numerical progress in the study of structure formation. A further limitation is that we will not make any computation in red-shift space, which is where galaxies are observed. Finally, we will not investigate the three-point correlation function of halos (nor galaxies), which is a natural observable for primordial three-point functions.

1 Inflation and primordial perturbations

It is a remarkable achievement of the generic inflationary scenario to provide a mechanism to generate primordial perturbations, which source the formation of structure in the Universe, considering that it was not designed for it, but to solve the well-known hot big bang problems Guth (1981); Starobinsky (1980); Linde (1982). There are countless good introductions to inflation and the production of primordial perturbations, including textbook material Kolb and Turner (1990); Liddle and Lyth (2000); Dodelson (2003); Linde (1990); Mukhanov (2005); Weinberg (2008) and reviews Lyth and Riotto (1999); Malik and Wands (2009); Baumann (2011); Riotto (2018). The reader is encouraged to look at these references for a detailed analysis. Here we will give the main physical intuition without enter any technical detail.

1.1 Background evolution

A good simple model to start with is a single scalar field, called generically the inflaton, coupled to gravity through the metric $g_{\mu\nu}$ and slowly rolling down its potential. Its vacuum energy drives the accelerated expansion of the Universe. The corresponding action reads

[TABLE]

where $M_{P}$ is the reduced Planck mass, the first term is the Einstein-Hilbert action and $V(\phi)$ is the slow-roll potential. The background evolution is studied by assuming a flat Friedmann-Lemaître-Robertson-Walker (FLRW) metric

[TABLE]

where $a(t)$ is the scale factor, for which the equations of motion are

[TABLE]

Here the prime on the potential $V$ indicates derivation with respect to the background field $\phi_{0}$, while the dot denotes the time derivative with respect to cosmic time $t$. The expansion rate $H=\dot{a}/a$, known as the Hubble parameter, is determined by the first equation, while the evolution of the background field $\phi_{0}$ is determined by the third equation and the two are connected by the second equation. Inflation can be achieved when the potential $V$ dominates over the kinetic energy of the inflation. As a consequence, $\dot{H}\approx 0$ and the scale factor grows exponentially. Observations of the CMB require a minimum duration for inflation in order to solve the horizon and flatness problems. This is usually quantified in number e-folds $N$

[TABLE]

here measured from a time $t$ during inflation until the end of inflation, $t_{e}$. In this parametrization, inflation is required by CMB to run for $O(60)$ e-folds Akrami et al. (2019) and as a consequence the Hubble parameter has to stay almost constant within a typical Hubble time $H^{-1}$. The slow-roll conditions therefore are

[TABLE]

and have to be satisfied during inflation. Using Friedmann equations, we can impose the slow-roll conditions on the shape of the potential as well

[TABLE]

Within the approximation of perfectly constant rate of expansion, the scale factor grows exponentially in time $a\propto e^{Ht}$ and the space time is called De-Sitter (DS) space-time.

1.2 Quantum fluctuations of the inflaton

Let us now turn to perturbations. The inflaton field sets a “clock” for the amount of the expansion that the universe goes through during this phase and this amount is dictated, as a lower limit, by observations. A quantum-mechanical clock can not be infinitely precise though, as a consequence of Heisenberg’s uncertainty principle, instead it has a variance. The inflaton is therefore subject to spatially varying fluctuations, such that

[TABLE]

These spatial fluctuations determine small differences in the time at which inflation ends, so that the Universe inflates by different amounts in different regions. This physical phenomenon is at the basis of the generation of the primordial perturbations throughout the density distribution in the Universe.

Because of the quantum nature of these primordial perturbations, the full computation would involve the quantization of the coupled fluctuations of the inflaton and the metric. For the purpose of this review, it is sufficient to only show the case of a perturbed scalar field in DS without coupling to gravity. Counter-intuitively, this treatment captures most of the crucial points of the full coupled system and it is is somewhat technically simpler. We also ignore terms suppressed by the slow-roll parameters for now, so that, for example, the inflaton is effectively massless, being the second derivative of the potential constrained by Eq. (8). The action for the perturbed massless free field at second order reads

[TABLE]

where we have defined conformal time $\tau$ as ${\rm d}t\equiv a{\rm d}\tau$. Note that $\tau\in(-\infty,0)$. The equations of motion in Fourier space for this action read

[TABLE]

where $\mathcal{H}$ is the Hubble parameter in conformal time. The generic solution to this differential equation is333Here we used the approximation of perfectly DS spacetime, which implies $\tau=-1/Ha$. Corrections are proportional to slow-roll parameters, which indeed we are neglecting at this stage.

[TABLE]

To quantize the perturbations, we promote the field $\delta\phi$ and its conjugate momentum $\delta\pi\equiv\partial L/\partial\delta\phi^{\prime}$ to operators and write canonical commutation relations

[TABLE]

It is useful to decompose the field defining time-independent creation and annihilation operators in Fourier space

[TABLE]

with commutation relations

[TABLE]

which follow from the fact that $a^{2}u({\bf k},\tau)u^{\prime*}({\bf k},\tau)-\mbox{ c.c.}=t$ -independent constant. We would like now to fix the initial conditions $c_{1}$ and $c_{2}$. The first constraint comes from the normalization conditions

[TABLE]

and the second from choosing the vacuum state. To do this, we notice that on sub-horizon scales, that is for $k\tau\gg 1$, the system changes on time scales much shorter than the typical time scale of expansion, $H^{-1}$. In this limit, the first term of the solution Eq. (12) approaches the vacuum mode of Minkowski space-time and hence defines a preferable set of mode functions and a unique physical vacuum, usually referred to as Bunch-Davies state Bunch and Davies (1978). With this choice of initial conditions, the solution reads

[TABLE]

While well within the horizon, $k\tau\gg 1$, this solution is highly oscillatory, in the opposite limit, on super horizon scales, the amplitude asymptotes to a constant. This is the essential feature of models of inflation, because it states that patches of the size of the horizon or bigger evolve classically with spatially modulated amplitude given by the different values of $\delta\phi$ at horizon exit for that patch. As a consequence, the “inflation clock” stops at slightly different times in different patches of the Universe, or, in other words, the inflationary e-folding has space-dependent variations of size

[TABLE]

where the approximate equality indicates again corrections suppressed by slow-roll parameters. An analysis extended to the coupled inflaton-metric fluctuations would show that $\zeta$ can be defined as a Gauge-invariant quantity, called comoving curvature perturbation, which is conserved on superhorizon scales and it is directly related to the temperature anisotropies we observe in the CMB, $\zeta\approx-5\Delta T/T$. Correlation functions of this quantity therefore characterize perturbations produced during inflation. The power spectrum $P_{\zeta}$ is defined as

[TABLE]

Combining Eqs. (21), (15) and (20) we get

[TABLE]

where it is customary to define a dimensionless power spectrum as $\mathcal{P}_{\zeta}(k)\equiv k^{3}/2\pi^{2}P_{\zeta}(k)$. Here $*$ indicates that the quantities are computed at horizon crossing, $k\tau_{*}=1$, and we reinstated slow-roll suppressed parameters in the last factor, giving the power spectrum a spectral index

[TABLE]

The amplitude $\mathcal{A}_{s}=H^{2}_{*}/8\pi^{2}M^{2}_{P}\epsilon_{*}$ and spectral index $n_{s}$ of the power spectrum have been measured by Planck Akrami et al. (2018) to be

[TABLE]

It is remarkable to notice that the deviation of $n_{s}$ from $1$ is statistically significant, indicating the first direct measurement of time dependence in the inflationary evolution. On the other hand, since there is no direct constraint on the slow-roll parameter $\epsilon$, the value of the Hubble parameter during inflation, which is related to the energy scale of inflation, is still unknown. A measurement of the primordial gravitational waves power spectrum would break the degeneracy.

2 Interactions from inflationary models

In the previous section, we have shown how a simple phenomenological model of inflation with a single scalar field driving the accelerated expansion is able to produce primordial perturbations with characteristics that match CMB observations. Constructing an inflationary setting featuring these predictions is not as hard as it would seem: indeed a wide range of viable inflation models are available in the market. For this reason, it becomes crucial to extract information about higher-order primordial correlators, which would indicate that the statistics of $\zeta$ is not Gaussian and that nonlinear interactions are taking place during or after inflation. As we briefly show in this section, the nature of these nonlinear couplings is encoded in the higher-order correlators.

Plan of the section. We start by briefly introducing the basic points of the in-in formalism and connect it to the deviation from Gaussian statistics of the primordial perturbation $\zeta$ in Sec. §2.1, we then review different types of interactions as produced by various models, or classes of models, of inflation in Sec. §2.2 and we conclude with remarks about the current status and future prospects of detection of primordial non-Gaussianity, Sec. §2.3.

2.1 Interactions as non-Gaussianities

The computation of correlation functions in inflation is somewhat different than the usual quantum field theory methods applied to particle physics. In the latter case, scattering amplitudes are considered as non-interacting at some very early and very late times, far enough from the interaction region. In this way, the boundary conditions are taken to be vacuum states of the free theory on these limits, respectively called in and out states. On the other hand, during inflation the universe undergoes accelerating expansion and correlation functions are to be computed at fixed time. As we have seen in the previous section, only modes with wavelengths much smaller than the horizon can be approximated as living in a flat Minkowski space. This is the limit where the Bunch-Davies vacuum has been defined. Boundary conditions are therefore defined only at very early times, when most of the wavelengths are well within the horizon. The formalism to compute correlation functions in cosmological settings is called the in-in formalism Schwinger (1961); Bakshi and Mahanthappa (1963); Keldysh (1964); Calzetta and Hu (1987); Weinberg (2005). A generic form for an in-in correlator is

[TABLE]

where $Q(\tau)$ represents the operator of the product of $n$-curvature perturbation $\zeta$s so that $\hat{Q}=\zeta_{{\bf k}_{1}}\zeta_{{\bf k}_{2}}...\zeta_{{\bf k}_{n}}$ and $\Omega$ is the in state, the vacuum of the interacting theory. Note that, more in general, there could be combinations of both $\zeta$ and $\gamma$, the tensor perturbation, in $\hat{Q}$. In this review we only focus on $\zeta$, but tensor perturbations are equally important in constraining the inflationary scenario Arutyunov and Frolov (1999); Maldacena (2003); Maldacena and Pimentel (2011); Creminelli et al. (2012); Hinterbichler et al. (2014).

The time $\tau$ at which the correlator is computed is usually either at horizon crossing for the modes of interest or at the end of inflation. The strategy to compute correlators such as the one in Eq. (27) is to evolve $Q(\tau)$ back to initial time $\tau_{i}$ where the vacuum state is defined. In order to do that, the interaction picture is used, in which the background time dependence is determined by the quadratic Hamiltonian $H_{0}$, while interactions arise as corrections to $H_{0}$ through the interaction Hamiltonian $H_{\rm int}$. Eq. (27) is therefore evaluated formally as

[TABLE]

where $T$ and $\bar{T}$ are time and reversed time-ordering symbols and $\hat{Q}^{I}$ and $\hat{H}^{I}_{\rm int}$ are evaluated using interaction picture operators. The $i\epsilon$ prescription allows to project the interacting vacuum state, $\Omega$ to the free one. The above expression can be expanded as a power series in $H_{\rm int}$. Interactions are then organized as usual using Feynman diagrams. As an example, let us take the three-point correlator of $\zeta$ expanding Eq. (28) to first order

[TABLE]

where $H_{\rm int}=-L_{3}+\mathcal{O}(\zeta^{4})$ and $L_{3}$ is the perturbed Lagrangian up to cubic order in $\zeta$ and we have taken the superhorizon limit $k\tau\rightarrow 0$. This is how interactions connect to non-Gaussianities through the perturbed action: they generate higher-order correlators, such as the bispectrum in Eq. (29).

In this review, we exclusively consider non-Gaussian signatures which come from a non-zero three-point correlation function, or more frequently used, its equivalent in Fourier space, the bispectrum

[TABLE]

It is customary to decompose $B_{\zeta}$ as

[TABLE]

where we assumed statistical homogeneity and isotropy and we evaluate the dimensionless power spectrum, $\mathcal{P}_{\zeta}$ at pivot value, neglecting the small scale dependence $(k/k_{p})^{n_{s}-1}$ from now on. All the momentum dependence goes therefore in the function $\mathcal{S}$, which encodes all the crucial information about the bispectrum. This information can be categorized into three features:

• •

the shape of the bispectrum, which is usually expressed through the dependence of $\mathcal{S}$ on the ratio of the momenta, for instance $k_{2}/k_{1}$ and $k_{3}/k_{1}$.

• •

the running of the bispectrum, which refers to the dependence of $\mathcal{S}$ on the sum of the amplitude of the wave numbers $K=k_{1}+k_{2}+k_{3}$.

• •

the amplitude of the bispectrum, usually denoted as $f_{\rm NL}$ and defined as

[TABLE]

where, as we can see, $f_{\rm NL}$ can generically depend on $K$.

2.2 Interactions in models of inflation

Having made the connection between interactions and non-Gaussianities, we now classify models of inflation by the interactions, and consequently non-Gaussianities, they produce during inflation. All these models satisfy the minimal current constraints on inflation that we outlined in the previous section §1. Therefore, the only way to discriminate among them is to observe some degree of non-Gaussianity in the CMB temperature anisotropies or in the LSS.

In the phenomenological example we outlined in §1, inflation is run by a single, scalar field. Even in this simple case, non-Gaussianities through various types of interactions are produced, leading to a rich phenomenology in correlation functions of the curvature perturbation $\zeta$. On the other hand, it is reasonable to think that inflation was populated by many fields. The case of more than one field during inflation is clearly even richer, as it includes all the features of single-field models plus the possible interactions among the fields. These fields might have very diverse functions: i) contribute to the background, some/all of which ii) generate primordial perturbations, or some might simply iii) spectate, i.e. they do not give a substantial contribution to neither the accelerated expansion nor the perturbations. Here we are interested in the effects that a multi-field scenario has on the correlation functions of the curvature perturbation $\zeta$, and consequently on structure formation, thus we will restrict to the study of perturbations. We distinguish between two cases: first, the case of massive particles present during inflation, which decay right after horizon crossing and therefore their signature is left in the curvature perturbation during inflation. Second, the case in which perturbations from massless particles survive even on superhorizon scales and imprint non-Gaussianities after inflation.

2.2.1 Interactions in single-field models

In this section we review different ways to produce non-Gaussianities in single-field models. They range from very small (though never zero), i.e. of order of the slow-roll parameters $\epsilon$ and $\eta$, to possibly larger amplitudes when relaxing the assumption of minimal interaction of the inflaton with gravity. The list that follows is not complete, but rather it gives a schematic idea of the class of models that can generate non-Gaussianities in single-field inflation.

Gravitational floor

The example in Eq. (29) was not selected casually: it is the case of nonlinearities produced by the lowest order interaction of the inflaton with gravity. This is why it is called gravitational floor: every model of inflation is expected to have this minimal amount of interaction and non-Gaussianity produced. The integral in Eq. (29) was first performed in the seminal paper by Maldacena Maldacena (2003) and we redirect directly to it for details on the computation. The dimensionless bispectrum that results from the computation reads

[TABLE]

where $K=k_{1}+k_{2}+k_{3}$. By taking the limit where all the momenta are equal and using Eq. (32) we extract the amplitude $f_{\rm NL}^{\rm GF}=55/9\,\epsilon+15/9\,\eta$. The running of this shape is typically small, given that it is of order the running of the slow-roll parameters.

The squeezed limit of the bispectrum, i.e. the limit where one of the momenta is soft, provides the single-field consistency relation

[TABLE]

Extensions and generalizations of this statement have been made in several subsequent efforts Maldacena (2003); Creminelli and Zaldarriaga (2004); Cheung et al. (2008); Creminelli et al. (2011, 2012); Assassi et al. (2012). This relation is valid for any single-field model of inflation, not necessarily slow-roll, since it can be demonstrated on the sole assumption that the curvature perturbation $\zeta$ is constant on super-horizon scale, i.e. that it is adiabatic. It follows from the fact that a long-wavelength perturbation, if adiabatic, corresponds to a local rescaling of the background experienced by short-wavelength ones. The direct consequence of Eq. (34) is that any detection of non-Gaussianity for a bispectrum computed in the squeezed limit indicates that perturbations were generated by more than one field during or after inflation.

The consistency relation can be extended to finite long-wavelengths by an expansion in $k_{\rm long}/k_{\rm short}$

[TABLE]

where $b_{0}=1-n_{s}$ and $b_{1}$ is also fixed by symmetries, corresponding to a local constant gradient rescaling of short modes. Owing to the consistency relation, all powers of the expansion are integer. Interestingly, the presence of additional fields introduces non-analytic scalings, as we will show below for the case of quasi-single field inflation.

Higher-Derivative kinetic terms

A first extension of the model presented in the previous section, Eq. (1) is done by writing the most general Lorentz-invariant Lagrangian as a function of the inflaton $\phi$ and its derivative Garriga and Mukhanov (1999)

[TABLE]

where $X\equiv-1/2g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi$. In the context of an effective field theory of inflation, the form of $P(X,\phi)$ is constrained by the following considerations:

• •

Non-derivative operators, such as $\phi^{n}/\Lambda^{n-4}$, being $\Lambda$ the largest energy scale where the effective description holds, contribute directly to the inflaton potential and are therefore strongly constrained by the background Copeland et al. (1994).

• •

Derivative operators of the form $X^{n}/\Lambda^{4n-4}$ do not suffer this limitation. However, a simple estimation of the amplitude of the leading correction to the slow-roll Lagrangian, $X^{2}/\Lambda^{4}$, gives

[TABLE]

Non-Gaussianities of order unity are therefore generated when $\dot{\phi}_{0}\approx\Lambda^{2}$, which is the regime where the effective description brakes down. In other words, the effective Lagrangian becomes unstable to radiative corrections. Non-gaussianities of order unity of this type therefore represent a particularly well motivated target for observational searches, because $\dot{\phi}_{0}$ is already a relevant energy scale for the dynamics of the inflationary background, being the scale related to the breaking of exact DS background evolution.

Following these considerations, the only way to escape the break down of the effective description and to have sizeable non-Gaussianities is to write down a UV-complete model. One such example is provided by the Dirac-Born-Infeld (DBI) model Silverstein and Tong (2004); Alishahiha et al. (2004), for which

[TABLE]

where $f(\phi)$ is the (squared) warp factor of the AdS-like throat related to these models444 One should bear in mind that there is a degree of fine-tuning needed also for this model Chen (2008).. The dimensionless bispectrum for a generic $P(X,\phi)$ model can be written in a rather model-independent way Seery and Lidsey (2005)

[TABLE]

where the subscript HD stands for “Higher-Derivative” terms and

[TABLE]

and we have defined the following quantities Garriga and Mukhanov (1999)

[TABLE]

where $P_{X}$ denotes the derivative of $P$ with respect to $X$ and we stopped at the third derivative because we limit to the bispectrum and we quoted only the leading order in slow-roll parameters. In $P(X)$ models, $c_{s}$ is the speed of propagation of the scalar perturbations and can be typically different then unity, thus leading to sizeable bispectra. The bispectrum amplitude for these models indeed is given by

[TABLE]

As it is clear from Eq. (39), models with large non-Gaussianity, as the DBI model, need to have $c_{s}^{2}\ll 1$ and/or $\lambda/\Sigma\gg 1$. In the specific case of the DBI model, by combining Eq. (38) into Eq. (42) and (39), Eq. (45) vanishes so that only $S_{c}$ contributes with a sizable non-Gaussianity.

This bispectrum peaks in the equilateral triangle configuration, where all the momenta are similar. The physical intuition for this fact is straightforward: modes that are much longer than the others, once out of the horizon, cannot interact with those within the horizon. Large interactions can occur only when all momenta are similar and therefore exit the horizon at the same time. The running of this bispectrum is small.

Features during inflation

The presence of features in the primordial power spectrum is frequently linked to a sizeble bispectrum. From the observational point of view, a feature represents a breaking of scale-invariance of the primordial spectrum within a range of scales. From the point of view of building inflationary models based on quantum gravity, or string theory, constructions, scale invariance is often a result of various mechanisms, while features can appear rather naturally Chen (2010); Chluba et al. (2015). These features can be broadly classified into two categories: sharp features or periodic oscillations. Sharp features in the potential or in the internal field space can arise from a variety of models Adams et al. (2001); Chen et al. (2007); Bean et al. (2008); Achucarro et al. (2011); Miranda et al. (2012) and they can also be studied in the context of the effective field theory of inflation Bartolo et al. (2013). Such features typically show up primarily in the power spectrum and constrains can be put with observations both of the CMB Akrami et al. (2018) and LSS Beutler et al. (2019), so non-Gaussianities can be used as a cross-check.

*Sharp feature. *The feature causes the inflaton to momentarily exit the attractor phase and consequently the slow-roll parameters to vary over a few e-folds. If we characterize the feature by its relative height $\mu\sim\Delta V/V$ and width $\sigma$, it can be easily shown on general grounds that observations of the power spectrum from the CMB constrain the ratio $\mu/\epsilon\lesssim 1$ Akrami et al. (2018). The slow-roll parameter $\eta$ on the other hand can change a lot if the change in $\epsilon$ occurs within a short time. If we take $\Delta\epsilon\sim\Delta(\dot{\phi}_{0}^{2})/H^{2}\sim\mu$ and $\Delta t\sim\Delta\phi/\dot{\phi}_{0}\sim\sigma/\sqrt{V(c+\epsilon)}$ we get

[TABLE]

For most models with a sharp feature, the calculation of the bispectrum requires to use numerical solutions. Here we quote an approximated shape Chen (2010)

[TABLE]

where $1/k_{*}$ is the oscillatory frequency in Fourier space corresponding to the feature. The exponential cuts off long-wavelength modes which are much longer than $k_{*}$ and feature is smoothed with a power $n$ to be fit with numerical results, along with $m$. The most important property of this type of non-Gaussianity is the running of the bispectrum, which is explicit in the dependence on $K$.

Resonant running.

A different type of features might be generated by an oscillatory component in the background evolution, as predicted, for instance, by axion-monodromy inspired models of inflation McAllister et al. (2010); Flauger et al. (2010). In these models, the inflaton potential is characterized by an oscillatory term added to the usual slow-roll one

[TABLE]

where here $\Lambda$ is some high energy scale and $f$ is the axion decay constant. As we have already seen, each mode oscillates during inflation with decreasing frequency as it is stretched by the accelerated expansion, until it reaches $H$ where it becomes frozen. Therefore, for any oscillatory feature with frequency $\omega>H$, we might expect a resonance between couplings and modes which sources non-Gaussianities Chen et al. (2008); Flauger and Pajer (2011). This type of non-Gaussianity, differently from the other cases analyzed before, is generated on sub-horizon scales. It can be shown that the parameter space for this type of resonance can be large, since

[TABLE]

being $\epsilon$ the slow roll parameter and $f\ll M_{P}$ is predicted in string theory constructions with sub-Planckian decay constants Banks et al. (2003); Svrcek and Witten (2006). The corresponding dimensionless bispectrum for this type of models reads

[TABLE]

where $\alpha=\sqrt{2\epsilon}/fM_{P}$ is constrained to be large by Eq. (50) and $k_{*}$ is the pivot scale at which the amplitude of the dimensionless power spectrum is defined. Similarly to the sharp feature case, this type of non-Gaussianity has sizeble running.

Non Bunch-Davies vacua

The choice of the vacuum state during inflation is not unique. To address the ambiguity from basic principles one would need to know the full theory at the highest energies, where we expect the free-field approximation to break down, as well as the physics preceding inflation. Notwithstanding this, the issue can be addressed phenomenologically: any deviation from the attractor solution of the inflaton, such as the ones related to sharp features studied above, lead generically to a deviation from the standard Bunch-Davies vacuum. This is because the Bunch-Davies vacuum is chosen out the asymptotic limit of $k\tau\gg 1$ of the attractor solution. Non-Gaussianities produced by choosing different prescriptions for the vacua are rather model-dependent and have been explored in a number of papers (see Chen et al. (2008) for a summary and list of references). A common feature of non-Gaussianities resulting from non Bunch-Davies vacua is the fact that they are enhanced in the folded triangle limit, i.e. for $k_{1}+k_{2}-k_{3}=0$. The intuitive explanation can be understood looking at Eq. (12): in the Bunch-Davies case, only positive frequencies are considered and therefore the second term proportional to $c_{2}$ is neglected. Negative frequencies are instead produced when deviating from the attractor solution, and the leading order deviation in the in-in calculation of the bispectrum of $\zeta$ has therefore at least one contribution from them. Effectively, this translates into sending one of the three momenta from $k\rightarrow-k$.

Solid inflation

Solids can be described in the context of field theory Dubovsky et al. (2006) by introducing three scalar fields $\phi^{I}$ whose background values are identified with spatial coordinates

[TABLE]

and they are time-independent. Using this framework, Endlich et al. (2013) showed that inflation can be driven by a particular type of solid which has approximate dilation symmetry and exact rotational and translational internal symmetries. These symmetries allow to consistently build a solid that stretches during the accelerated expansion by many orders of magnitude. Even though the treatment of this model is done in the context of an effective field theory, the time-independence of the background fields implies that there is no breaking of time-translational invariance as in conventional effective field theory approaches of inflation. The computation of cosmological perturbations also shows peculiarities: for instance, adiabatic perturbations during inflation are absent Endlich et al. (2013).

Most interestingly for this review, the three-point function of scalar perturbations drastically violates Endlich et al. (2014) the standard single-field consistency relation of Eq. (34). The dimensionless bispectrum is calculated to be Endlich et al. (2013)

[TABLE]

where $F$ and $F_{Y}$ are free parameters of the solid Lagrangian, $\epsilon$ plays the role of the slow-roll parameter, $c_{L}$ and $c_{T}$ are the speeds of longitudinal and transverse phonons of the solid, respectively, $\tau_{c}$ is the conformal time at which the longest modes of observational relevance today exit the horizon, while $\tau_{e}$ is the conformal time at the end of inflation. The functions $Q$ and $U$ are given in Eq. (7.4) and (7.13) of Endlich et al. (2013), respectively. It is important to notice that Eq. (53) is computed for not-too-squeezed momenta, that is, for $k_{\rm long}/k_{\rm short}>\sqrt{\epsilon}$. The squeezed limit was investigated in detail in Endlich et al. (2014), at leading order in slow-roll expansion it can be written in a much more compact form as

[TABLE]

which should be compared to the consistency relation of Eq. (34). Here $\theta$ is the angle between ${{\bf q}}$ and ${{\bf k}}$. The fact that the angular dependence is that of a quadrupole and the overall amplitude not being constrained to be small are clear signals of the breaking of the consistency relation, despite the fact that solid inflation propagates only a single scalar mode.

2.2.2 Multi-field interactions during inflation

In this section, we provide a schematic summary of the interactions between massive particles and the inflaton taking place during inflation, which lead to non-Gaussian signatures in the correlation functions of $\zeta$.

Massive particles might be spontaneously created in the expanding space-time during inflation through non-perturbative effects Parker (1968, 1969, 1971). Their production is particularly interesting as a probe to ultraviolet completions of inflation motivated by string theory constructions Baumann and McAllister (2015). We must specify that they are massive since their mass cannot be arbitrarily light: to avoid back-reaction on the inflationary background, their typical mass must be of order $H$ or higher. On the other hand, the production rate is exponentially suppressed as a function of mass in De-Sitter space-time, roughly as $e^{-m/H}$. Therefore very massive particles, $m\gg H$, decay exponentially fast after horizon crossing and their effect can be integrated out from the dynamics of $\zeta$. In this case, inflation is effectively single-field, so that, for instance, the inflationary consistency relation of Eq. (34) has to be satisfied. Particles with mass of order $H$ instead produce characteristic non-local signatures in the correlation functions of $\zeta$ and will generically violate the single-field consistency relations. Here we want to stress that we are not only restricting to scalar fields: indeed a tower of high-spin states can arise in string theory constructions Rindani and Sivakumar (1985); Aragone et al. (1987). While the effect of massive scalar fields have been thoroughly investigated in the context of quasi-single field inflation Chen and Wang (2010); Baumann and Green (2012); Assassi et al. (2012); Noumi et al. (2013); Gong et al. (2013); Baumann and McAllister (2015); Arkani-Hamed and Maldacena (2015); Lee et al. (2016); Bordin et al. (2016); Flauger et al. (2017); Baumann et al. (2018); McAneny and Ridgway (2019), higher-spin particles studies arose more recently Arkani-Hamed and Maldacena (2015); Lee et al. (2016); Kehagias and Riotto (2017a, b); Baumann et al. (2018); Franciolini et al. (2018); Biagetti et al. (2017); Dimastrogiovanni et al. (2018).

Massive particles in De-Sitter

In the De-Sitter background space-time of inflation, particles can be classified as unitary irreducible representations of the De-Sitter group SO(1,4). Particles are characterized by their spin and mass, and the condition of unitarity imposes three allowed categories Wigner (1939); Bargmann and Wigner (1948) for particles with spin $s\geq 1$

[TABLE]

for $t=0,1,2,...,s-1$. Similarly, scalar particles with mass $m\leq 3/2H$ belong to the principal series, while lightest particles belong to the complementary series. Massless scalar particles are conformally invariant in DS. Notice that particles with spin $s\geq 1$ are required to have a minimal mass, $m^{2}\geq s(s-1)H^{2}$, unless they belong to the discrete series555Light particles with spin are allowed in inflationary set ups which breake DS isometries, as shown in Bordin et al. (2018). For these values, the system acquires an additional gauge invariance and the corresponding fields are called partially massless fields Deser and Waldron (2001). Partially massless fields produce sizeable tensor-scalar-scalar bispectra and the scalar trispectra Kehagias and Riotto (2017a); Baumann et al. (2018); Franciolini et al. (2018). In the case of study here, i.e. the scalar bispectrum, partially massless fields do not contribute, as their lack of a longitudinal degree of freedom kinematically prevents them to oscillate into a single scalar field.

As we already anticipated, massive particles decay at late times. The DS group acts as conformal group on the three-dimensional Euclidean space for the super-Hubble fluctuations, so that the asymptotic scaling at late times of a scalar particles is

[TABLE]

where $\Delta=3/2\pm i\mu$ and

[TABLE]

where the massless case $m=0$ corresponds to a conformally coupled field that does not decay at late times. We will investigate this case in the next section. Similarly, for a spin-$s$ field we get

[TABLE]

where $\eta$ is conformal time and $\Delta_{s}^{\pm}=3/2\pm i\mu_{s}$ is the conformal weight of the field and

[TABLE]

will be a crucial parameter in the non-Gaussian shapes, as we will show shortly. For finite mass fields, the decay scales as the conformal weight and can be distinguished in real values of $\mu_{s}$, for which the wavefunction oscillates logarithmically in conformal time, and imaginary $\mu_{s}$, for which particles belong to the complementary series and survive longer at late-times.

Effective approach to interactions of massive particles

The interactions of the massive particles described above in the context of inflation has been investigated using the effective field theory of inflation in Lee et al. (2016). This extends earlier works focusing only on scalar fields. At tree-level, Lee et al. (2016) distinguished three diagrams, illustrated in Fig. 1, which represent three different ways in which spin-$s$ fields can be exchanged in the three-point correlator $\langle\zeta\zeta\zeta\rangle$.

The corresponding dimensionless bispectrum for these three diagrams can be written in the following compact form

[TABLE]

where $(\kappa)=(a),(b)$ and $(c)$ represents the three diagrams, $\alpha^{(\kappa)}_{s}$ are dimensionless parameters, $L_{s}^{(\kappa)}$ are functions of the angles between momenta (typically Legendre polynomials) and $\mathcal{I}^{(\kappa)}_{s}$ are complicated integrals whose expressions can be found in Lee et al. (2016)666 Note that there is a conversion factor of $(k_{1}k_{2}k_{3})^{2}$ between $\mathcal{I}_{s}^{(\kappa)}$ and those defined in Lee et al. (2016). For instance, $\mathcal{I}^{(a)}_{s}=(k_{1}k_{2}k_{3})^{2}\mathcal{I}^{(s)}$. and $c_{\pi}$ is the sound speed of the Goldstone boson in the effective field theory of inflation. Although in a compact form, the importance of Eq. (59) is manifest: $\mathcal{S}_{\rm MP}$ depends on the mass and spin of the particle which mediates the exchange, so it provides a promising way to detect new particles at the high energies at which inflation takes place Arkani-Hamed and Maldacena (2015).

To get a more explicit idea of this type of bispectra, let us restrict to the squeezed limit of the single-exchange diagram $(a)$. Looking at the squeezed limit is relevant because of the single-field consistency relation, Eq. (34), which implies that in this limit non-Gaussianity of order $\mathcal{O}(f_{\rm NL})\gtrsim 0.1$ is only sourced by the presence of multi-field scenarios. In this limit, the bispectrum splits into an analytic part and a non-analytic part. The analytic part reflects local effects of massive particles with a scaling similar to the case of single-field models, Eq. (2.2.1) and it does not contain information on the spin and mass of the particle at leading order,

[TABLE]

The non-analytic part for scalar particles has been derived in the context of quasi-single field inflation Chen and Wang (2010). It has to be distinguished into two cases: for massive particles belonging to the principal series, $\mu\geq 0$ it reads

[TABLE]

while for particles belonging to the complementary series for which $\mu$ becomes imaginary, the scaling changes to

[TABLE]

where $\nu=-i\mu$ is real. The case of higher-spin particles proceeds analogously: massive particles belonging to the principal series, $\mu_{s}\geq 0$, have a squeezed limit of the form

[TABLE]

for even spins, where the $\mathcal{L}_{s}$ is the Legendre polynomial of order $s$ and $\phi_{s}$ is a phase that is uniquely fixed in terms of $\mu_{s}$ and $c_{\pi}$ Lee et al. (2016). This oscillatory scaling is the distinctive feature of this type of interactions and it tells information about the mass and spin of the particle involved in the exchange.

For particles belonging to the complementary series for which $\mu_{s}$ becomes imaginary, the scaling changes to

[TABLE]

where $\nu_{s}=-i\mu_{s}$ is real. In the case of odd spins, for symmetry reasons the squeezed limit of the non-analytic part is suppressed by at least an additional power of $k_{1}/k_{3}$, therefore the leading piece is always the analytic one of Eq. (60). Both in the case of scalar and higher-spin massive particles, the least suppressed scaling in powers of the squeezed ratio is the limiting case of $m=s(s-1)$ and corresponds to a constant shape $\mathcal{S}\propto(k_{1}/k_{3})^{0}$.

2.2.3 Multi-field interactions after inflation

Massive particles necessarily decay after horizon crossing during inflation, as seen in the previous section. This is not the case for massless particles: they survive at late times and typically have a non-linear evolution in multi-field space on superhorizon scales, which generate isocurvature perturbations (see Riotto (2018) for recent lectures on the topic). Observations of the CMB constrain the amount of isocurvature perturbations to be very small Akrami et al. (2018). Nevertheless, there are several mechanisms with which isocurvature perturbations can be converted to curvature ones after inflation and evade CMB constraints Salopek and Bond (1990); Bartolo et al. (2002); Bernardeau and Uzan (2002, 2003); Rigopoulos et al. (2006a, b); Seery and Lidsey (2005); Vernizzi and Wands (2006). We will not go into details of specific realizations of this conversion mechanism. Instead, we will summarize the results obtained in the framework of the $\delta N$-formalism Starobinsky (1985); Sasaki and Stewart (1996); Lyth and Rodriguez (2005), which is rather model-independent. It will quickly become clear that the shape of this type of non-Gaussianities is generic and goes under the name of local non-Gaussianity. Local-type non-Gaussianities are the most studied in the literature (see Salopek and Bond (1990); Gangui et al. (1994); Verde et al. (2000); Komatsu and Spergel (2001) for the earliest studies) and can be realized by a wide range of models. The name local comes from the fact that it can be expressed as a simple Taylor expansion around a Gaussian field at position ${{\bf x}}$,

[TABLE]

The corresponding dimensionless bispectrum in Fourier space is of the form

[TABLE]

and it has the nice feature of being already factorizable and simple to implement as optimal estimator in observed quantities. Moreover, it is a distinctive signature of multi-field models, since it peaks in the squeezed limit, where single-field models are necessarily producing very small non-Gaussianity as implied by the consistency relation Eq. (34).

The $\delta N$-formalism

In the single-field model we introduced in §1, the comoving curvature perturbation $\zeta$ is constant on superhorizon scales. Let us consider the gauge in which the inflaton perturbations are zero and only fluctuations on the metric are present. In this gauge, the spatial metric is simply given by $a^{2}e^{2\zeta}$. Each comoving superhorizon patch in the universe will have its own frozen value of $\zeta$ which determines the local evolution in space, independently of the other disconnected patches, through the scale-dependent scale factor. It can be shown that this is equivalent to say that each patch evolves as a separate universe each with slightly different number of e-folds, $\delta N$, at different positions. This interpretation is dubbed $\delta N$-formalism.The generalization of this picture to the case of multi-field inflation is well summarized in Chen (2010), here we broadly follow its steps.

Let us consider a set of scalars $\phi_{i}=\phi_{0i}+\delta\phi_{i}$ during inflation. We are interested in the behavior of the perturbations of these fields on superhorizon scales, looking at different causally disconnected patches. We choose an initial spatially flat slice where scalar metric fluctuations are zero. Note that modes we are interested in are already superhorizon on this slice, so we will evolve them classically. Moreover, we assume their statistics to be Gaussian. We now select uniform density final slices, that is, slices with the same energy density at each point of them. We define $N_{0}(\phi_{0i})$ as the number of e-folds between the initial and final slices for the unperturbed field $\phi_{0i}$, while the one for the perturbed fields will be $N(\phi_{0i}+\delta\phi_{i}({{\bf x}}))$. The curvature perturbation $\zeta$ can be expressed as a Taylor expansion of the variation of $N$ around the initial values $\phi_{0i}$ as

[TABLE]

where the subscripts on $N$ denote partial derivatives with respect to $\phi_{i}$ and indices are summed using the Einstein convention, with $i=1,...,N$, being $N$ the number of fields. Correlation functions of $\zeta$ are then related to correlation functions of $\delta\phi_{i}$, which we assumed to be Gaussian distributed massless fields. We have computed the mode functions of each of these fields in Eq. (20), their power spectrum is

[TABLE]

where the prime indicates that we are leaving the usual delta function implicit and $H_{*}$ is the Hubble parameter at horizon exit. Consequently, using Eq. (67) we can express the bispectrum of the curvature perturbation as

[TABLE]

and it can be easily shown that the corresponding dimensionless bispectrum $\mathcal{S}$ is the local one of Eq. (66) and

[TABLE]

The fact that this shape is local should not come as a surprise: this bispectrum is sourced by local interactions of fields on superhorizon scales. A number of multi-field models can be described using the $\delta N$-formalism, such as the curvaton model, the modulated reheating and preheating scenarios (see Bartolo et al. (2004) for a review).

2.3 Final remarks of this section

We have identified three main features that we can extract from the bispectrum of the curvature perturbation $\zeta$: shape, running and amplitude. Among these three, the amplitude, parametrized via a single parameter $f_{\rm NL}$, is the easiest to constrain from data and it tells us already a lot about what is the source of non-Gaussianity, as shown schematically in Fig. 2.

The limit of $f_{\rm NL}\sim\mathcal{O}(1)$ is particularly important for interactions generated by non-trivial higher-order kinetic terms in single-field inflation: it would signal that interactions come at an energy scale which coincides with $\dot{\phi}_{0}$, which is the scale at which the exact DS background evolution is broken. These are non-Gaussianities typically peaking on equilateral shapes of the bispectrum (i.e. with $k_{1}\approx k_{2}\approx k_{3}$). Furthermore, a positive detection of an order unity $f_{\rm NL}$ for squeezed bispectra (i.e. with $k_{1}\ll k_{2}\approx k_{3}$) would rule out single-field models of inflation in one shot. The ultimate goal is to measure $f_{\rm NL}$ down to the level of slow-roll parameters, $f_{\rm NL}\lesssim\mathcal{O}(0.01)$. At these limits, inflation predicts that non-Gaussianity has to be there no matter the inflationary model. Experimental sensitivity seems to be still far from this goal.

Current best limits on non-Gaussianity are set by observations of the bispectrum of CMB temperature anisotropy from the Planck satellite Akrami et al. (2019)

[TABLE]

at $68\%$ confidence level. While these limits refer directly to Eq. (66) for the local shape, the equilateral and orthogonal constraints come from templates, rather than one of the shapes presented above. This is because the analysis from data require optimal estimators of the bispectrum which need to sum over all modes available in the survey. This implies that only factorizable shapes, such as the local one of Eq. (66), are usable in practice Creminelli (2003); Fergusson et al. (2010). For this reason, templates need to be used for CMB constraints in place of the realistic predictions we have outlined above and the two can be related using so-called “fudge-factors”. This has been applied also to LSS Wagner et al. (2010); Fergusson et al. (2012); Wagner and Verde (2012). We will only briefly mention this issue in Sec. §5.3.2 in the context of N-body simulations.

In this section, we have shown that interactions during and after inflation necessarily take place even in the minimal scenario of single-field models through minimal coupling with gravity. This implies that there is a gravitational floor for detectable non-Gaussianities in the CMB and LSS and therefore a guaranteed signal. In non-minimal scenarios, consistency relations and constraints from current observations allow to classify interactions from realistic models of the early Universe and in some cases even open the possibility to use inflation as a cosmological particle collider Arkani-Hamed and Maldacena (2015) and explore new physics at energies as high as $10^{14}$ GeV. Such powerful predictions are rare in cosmology, and should be a primary target for experimental searches.

3 From primordial interactions to matter overdensities

In this section, we briefly review how non-Gaussianities are transferred from the primordial curvature perturbation to the distribution of dark matter and its correlation functions. We will show how non-Gaussian initial conditions affect directly the mass density probability function and also generate additional terms in the correlation functions of dark matter at different positions. We will argue that these effects are typically very small and their search is complicated by the fact that we do not have direct access to dark matter correlation functions777Weak lensing Bartelmann and Schneider (2001) and temperature anisotropies in the 21cm background from the pre-reionization epoch Loeb and Zaldarriaga (2004); Cooray (2005); Pillepich et al. (2007) can be used to trace DM directly, although the latter might also be biased Pritchard and Furlanetto (2007).

The initial conditions for structure formation are set by connecting the primordial curvature perturbation $\zeta$ to the Newtonian potential $\Phi$ via a transfer function, $T(k)$. At the linear level and using the Poisson equation, we can write the useful relation

[TABLE]

which defines the linear (L) matter overdensity. Here $D(z)$ is the linear growth factor normalised to unity at present day and $\Omega_{m,0}$ and $H$ are the matter density and hubble parameter today, respectively. The linear relation is reliable as long as the matter overdensity $\delta$ is much smaller than unity. In this regime, one can solve the (Newtonian) dynamical equations perturbatively (see Bernardeau et al. (2002) for a review and Baumann et al. (2012); Carrasco et al. (2012); Pajer and Zaldarriaga (2013); Carrasco et al. (2014); Mercolli and Pajer (2014); Carrasco et al. (2014); Carroll et al. (2014); Porto et al. (2014); Senatore and Zaldarriaga (2015); Baldauf et al. (2016); Assassi et al. (2015); Baldauf et al. (2015, 2016); Lewandowski and Senatore (2017); Senatore and Trevisan (2018) for an approach using effective field theory). Even pushing perturbation theories to their limits, theoretical predictions describing the evolution of matter are valid only in the weakly nonlinear regime, that is for $k\lesssim 0.2$ h/Mpc at redshift [math], even in the case of Gaussian initial conditions. Consequently, for all practical purposes, it is convenient to smooth small scale perturbations with a window function $W_{R}$, being $R$ the smoothing scale, and we indicate the smoothed field as $\delta_{R}$. Commonly used window functions in LSS are spherically symmetric functions, such as top-hat filters in real space or Gaussian windows. Consequently, the linear density power spectrum smoothed on a scale $R$ is given by

[TABLE]

and its variance

[TABLE]

To avoid clutter, we will condense the transfer function and the smoothing window into $\mathcal{M}_{R}(k,z)\equiv\mathcal{M}(k,z)W_{R}(k)$ and drop the redshift dependence unless needed.

3.1 The density probability distribution

The probability distribution function (PDF) of the smoothed density field at early times is usually assumed to be Gaussian at each fixed point in space ${\bf x}$,

[TABLE]

i.e. the value of the smoothed density field $\delta_{R}$ at each ${\bf x}$ is drawn from a Gaussian distribution. We know however from the previous section that small initial non-Gaussianities are always present, even in the minimal case in which only a single-field and gravity are present during inflation. These non-Gaussianities source higher-than-two reduced smoothed cumulants, defined as

[TABLE]

where $\langle\delta_{R}^{N}\rangle_{c}$ is the connected n-point moment. Note that, in the linear regime, $\sigma_{R}\propto D(z)$ and therefore $S_{N}\propto 1/D^{N-2}$, which means that for $N\geq 3$, higher order cumulants are increasingly suppressed by powers of the linear growth factor $D$. As a working example, let us compute the lowest order cumulant which is non-zero in the presence of non-Gaussianity, the skewness

[TABLE]

where

[TABLE]

It is easy to recognize the inflationary three-point function as a source of the skewness. The minimal amount of skewness produced in the initial conditions can be computed by making use of Eq. (2.2.1) and gives $\sigma_{R}S_{3,GF}=A(R)\epsilon+B(R)\eta$, being $A$ and $B$ mild functions of the scale $R$ and with amplitude $A,B\propto\mathcal{O}(10^{-4})$ for $R\approx 3Mpc/h$. Similar values apply also for the local type non-Gaussianities of Eq. (66). To get the physical intuition of how the PDF changes, it is enough to take the simplest approach and perform an Edgeworth expansion on the smoothed field $\delta_{R}$, or more commonly the so-called peak height $\nu=\delta_{R}/\sigma_{R}$, to get

[TABLE]

where $H_{N}$ are Hermite polynomials. Note that, as a consequence of how the reduced cumulants scale with the linear growth factor, the combinations $\sigma S_{3}$ and $\sigma^{2}S_{4}$ are redshift independent in the linear regime. More details on these expansion, and its refinements, can be found in several analyses dating back to the 90s Juszkiewicz et al. (1995); Fry and Scherrer (1994); Chodorowski and Bouchet (1996); Gaztanaga and Fosalba (1998); Durrer et al. (2000).

The generic effect of a positive (negative) skewness on the PDF is to produce more overdense (underdense) regions in the matter distribution. As perturbations grow from the initial time to later times, gravitational instabilities generate a positive skewness in the PDF which eventually dominates over primordial contributions. The effect of the skewness from gravitational evolution can be easily computed within perturbation theory to be $S_{3}\approx 34/7$ at lowest order Bernardeau et al. (2002), which is much larger than the one sourced by primordial bispectra respecting the current constraints from the CMB.

Beside studies in N-body simulations, the PDF of the matter density field is not a direct observable. Nevertheless, a PDF with non-Gaussian initial conditions has a direct impact on the abundance of clustered objects and on the clustering statistics. This is what we discuss in the following section §4.

3.2 Dark matter correlation functions

Primordial perturbations are transferred to the matter density field via Eq. (73) and correlation functions are necessarily affected. The matter $N$-point correlation functions reads, at the linear level,

[TABLE]

It is immediately clear from Eq. (81) that a bispectrum, or higher-order correlations, of $\zeta$ source at least the corresponding matter field correlator. As remarked above, gravitational evolution also sources secondary non-Gaussianities though, and these eventually dominate on the primordial ones. In the next section, we will show that the most promising observables for disentangling secondary non-Gaussianities from primary ones involve biased tracers of the matter field, rather than directly probing the dark matter field. Nevertheless, let us review a few essential points about dark matter correlation functions with non-Gaussian initial conditions, while we redirect the interested reader to the latest analyses in the context of the effective field theory for a wide range of non-Gaussianities Assassi et al. (2015).

In standard Eulerian perturbation theory (EPT) (see Bernardeau et al. (2002) for a review), one follows the evolution in time of the matter density contrast $\delta$ at position ${\bf x}$. In the weakly non-linear regime, the Fourier mode of the density field reads

[TABLE]

where we simply denote $\delta$ as the unsmoothed, non-linear matter field and $F_{2}$ is the PT kernel

[TABLE]

and $\mu$ is the cosine between ${{\bf q}}_{1}$ and ${{\bf q}}_{2}$. For the present analysis, it is not essential to smooth the matter field, so we work only with the unsmoothed $\delta$. A contribution from a non-zero primordial bispectrum therefore arises already on the matter power spectrum at the $1$-loop order and reads

[TABLE]

where notice that here we have used the statistical homogeneity and isotropy to write the momentum dependence of $B_{\zeta}$. The full 1-loop power spectrum is defined as

[TABLE]

where $P_{22}$ and $P_{13}$ are the $1$-loop contributions to the Gaussian case Bernardeau et al. (2002). It is interesting to notice that this term scales as $D(z)^{3}$, while $1$-loop terms from purely Gaussian initial conditions scale as $D(z)^{4}$. The non-Gaussian correction is suppressed at large scales, since the kernel $F_{2}$ vanishes in that limit. On the other hand, at small scales late-time non-linearities become quickly important as midly non-linear scales $k\sim 0.1$ h/Mpc are approached. Moreover, in Assassi et al. (2015) it was shown that the $2$-loop Gaussian terms are comparable to the non-Gaussian ones at $1$-loop even in the mildy non-linear regime. Therefore, deviations from the Gaussian case are hardly reaching the percent level for non-Gaussianities within current constraints.

The matter bispectrum is sensitive at tree-level to primordial corrections as well as to gravitational non-linearities,

[TABLE]

Extensions to higher loops have been computed in the context of perturbation theory (see Assassi et al. (2015); Welling et al. (2016) for a treatment in the effective field theory of LSS). It is customary to define a reduced bispectrum

[TABLE]

which is time- and scale-independent in the Gaussian case at tree-level in perturbation theory for equilateral configurations, so it is the appropriate quantity to study different types of non-Gaussianities. Several studies have tested theoretical predictions against N-body simulations Sefusatti et al. (2010); Figueroa et al. (2012); Lazanu et al. (2017). The prospects for observing the imprint of non-Gaussianities in the matter bispectrum are also rather weak for two main reasons: first, the dark matter field is not directly observable, with the exception of the case of weak leansing bispectrum measurements, for which, however, it was shown that the sensitivity to primordial non-Gaussianity is around two orders of magnitude away from current limits Schaefer et al. (2012). Secondly, the amplitude of late-time non-linearities largely surpasses the primordial contributions on all scales of interests for future observations. Moreover, similarly to the case of the one-loop power-spectrum, the contribution from the two-loop gaussian bispectrum is larger than the one-loop non-Gaussian one Welling et al. (2016).

4 Imprints of primordial interactions on one-point halo statistics

As anticipated in the previous section §3.1, any higher-order correlation function generated during or right after inflation sources non-Gaussian terms on the PDF of the dark matter field. At early times, this is the only source of non-Gaussianity on the matter distribution and it implies a deviation from an equally probable abundance of overdense and underdense regions. Even as gravitational evolution generates secondary non-Gaussianities, the trace of the initial conditions can be disentangled to a certain degree by looking at very massive halos, whose formation is highly sensitive to the tails of the initial PDF. Massive halos are also more likely to be associated to peaks of the early time density field as shown in multiple studies in N-body simulations Ludlow and Porciani (2011); Elia et al. (2011); Ludlow et al. (2014), confirming that they are sensitive to the initial PDF. The pioneers in the analytic treatment of non-Gaussianities in the mass density PDF and in the abundance of DM tracers date back to the $1980$s Adler (1981); Matarrese et al. (1986); Lucchin and Matarrese (1988); Catelan et al. (1988); Borgani and Bonometto (1989); Matarrese et al. (2000); Scoccimarro et al. (2004). Almost as early, numerical methods have been used to test predictions and provide fits to data Moscardini et al. (1991); Weinberg and Cole (1992); Matarrese et al. (1991); Park et al. (1991); Gooding et al. (1991); Borgani et al. (1994). Since then, the growing interest in primordial non-Gaussianity has pushed for multiple developments in both directions.

Plan of the section. This section has two main goals: i) broadly summarize the theoretical advances by providing a background with the earliest attempts and then focusing on a few of the most recent progresses (Section §4.1) and ii) presenting the most recent numerical attempts at testing current predictions and providing with semi-analytic or fully fitted phenomenological models in Section §4.2. We then conclude with a few remarks in Section §4.3. This plan is far from being exhaustive and will inevitably refer the reader to the available literature for many details.

4.1 Analytic approaches

In order to extract most efficiently information about the early PDF of the matter field, which would give constraints on primordial non-Gaussianity, a consistent theoretical framework of halo formation is needed. In practice, however, what we really measure is number counts of biased tracers, such as galaxies and clusters of galaxies. It is therefore sufficient to have a working model for the number density of halos of a certain mass per unite volume, known as the halo mass function. Since the earliest attempts of modeling the abundance of bound objects in the presence of non-Gaussian initial conditions, this search has progressed side to side with models of the simpler Gaussian case: for instance, several attempts tried to extend the Press and Schechter (PS) mass function Press and Schechter (1974) to local-type non-Gaussianities in the initial conditions LoVerde et al. (2008); Jimenez and Verde (2009); Desjacques and Seljak (2010); Lam and Sheth (2009a, b); Lam et al. (2009). This extension was also studied for higher order primordial non-Gaussianities Desjacques and Seljak (2010); Chongchitnan and Silk (2010); LoVerde and Smith (2011) and a range of other bispectrum shapes Barnaby and Shandera (2012) ; the excursion set approach Epstein (1983); Bond et al. (1991); Lacey and Cole (1993); Sheth (1998); Porciani et al. (1998), which was introduced to solve problems suffered by the PS model, was thoroughtly also investigated with non-Gaussianity in the initial conditions Inoue and Nagashima (2002); Avelino and Viana (2000); Lam and Sheth (2009a, b); Lam et al. (2009, 2010); Maggiore and Riotto (2010); De Simone et al. (2011); Musso and Paranjape (2012); D’Amico et al. (2011a, b); Paranjape et al. (2011); Achitouv and Corasaniti (2012); Achitouv et al. (2013); Achitouv and Corasaniti (2012). Lastly, the peak model Bardeen et al. (1986), which was recently combined with the excursion set approach in the Excursion Set Peaks (ESP) model Paranjape and Sheth (2012); Paranjape et al. (2013), has been applied to non-Gaussian peak statistics Gay et al. (2012); Codis et al. (2013); Desjacques et al. (2013); Biagetti and Desjacques (2015); Lazeyras et al. (2016); Uhlemann et al. (2018). Recently, methods based on spherically averaging cosmic densities have been shown to successfully disentangle primordial non-Gaussianities with late time ones Uhlemann et al. (2018).

4.1.1 The Press-Schechter mass function

The Press-Schechter (PS) model Press and Schechter (1974) is a good framework to make simple analytic computations, hence we use it here to show the main physical intuition on the problem and then briefly mention extensions which improve it. It is based on the fundamental assumption that collapse is spherical Gunn and Gott (1972). Let us first take the case of Gaussian initial conditions. The PDF for the density field $\delta_{R}$ smoothed on a scale $R\propto M^{1/3}$, where $M$ is the halo mass, is a simple Gaussian with zero mean and variance $\sigma^{2}_{R}$. In the PS approach the halo mass function reads

[TABLE]

where $\bar{\rho}$ is the mean matter density, $F_{>\delta_{sc}}$ is the level excursion probability that the density field smoothed on a scale $M$ has overcome the threshold for spherical collapse, $\delta_{sc}\approx 1.687$888 This value corresponds to a Universe completely dominated by matter. Small corrections need to be taken into account when including the late time domination of dark energy., and $\nu(M,z)=\delta_{sc}/\sigma_{R}(z)$ is the peak height999 The smoothing scale $R$ is directly related to the mass of the halo $M$, we will use them interchangeably.. Here we have corrected for a factor of $2$ which accounts for the cloud-in-cloud problem, that is, the fact that the PS mass function does not include the possibility that overdense regions can be contained in bigger, underdense, ones. Excursion set models solve this problem Epstein (1983); Bond et al. (1991); Lacey and Cole (1993); Sheth (1998); Porciani et al. (1998). The halo mass function of Eq. (88) represents the differential number density of halos per unit mass and volume and it is an exponentially decreasing function of $\nu$, i.e. more massive halos are more rare. Despite the fact that the linear theory value of the threshold for spherical collapse is used, $\delta_{sc}$, the calculation is fully non-linear Gunn and Gott (1972). A feature of Eq. (88), and common to all universal mass functions101010As we discuss in Sec. §5.1.2, the assumption of universality of the mass function, although well tested in N-body simulations with Gaussian initial conditions, is cause of concern in the presence of primordial non-Gaussianity Biagetti et al. (2017); Desjacques et al. (2018). , is that it is entirely specified by a function of $\nu$ only,

[TABLE]

where $f(\nu)$ is called multiplicity function and we have dropped a subscript $\nu_{R}$ to avoid clutter, but the smoothing on a scale $R$, corresponding to a halo of mass $M$, is understood. In the PS case we have

[TABLE]

The PS mass function of Eq. (88) provides a decent fit to data in the intermediate halo mass regime, but poorly predicts the high mass regime. Nevertheless, an estimation of the non-Gaussian mass function was first given in Matarrese et al. (2000) by performing a saddle-point approximation on an Edgeworth expansion on $F_{>\delta_{sc}}$ and computing the ratio of the non-Gaussian-to-Gaussian multiplicity functions

[TABLE]

where $\nu_{*}=\delta_{sc}\sqrt{1-S_{3}\delta_{sc}/3}$ is introduced to enforce the normalisation of the halo mass function. The non-Gaussian mass function was then obtained by multiplying Eq, (91) to the Gaussian PS mass function. Alternatively, LoVerde et al. (2008) expanded directly the PDF, Eq. (80), to calculate $F_{>\delta_{sc}}$ and then the ratio

[TABLE]

where here we only considered the skewness. Combining the two predictions to better match numerical simulations, Desjacques and Seljak (2010) got

[TABLE]

where $\nu_{*}$ from Eq. (91) was expanded to first order in $f_{\rm NL}$. Earliest checks of these predictions against N-body simulations showed that they tend to overestimate the effect of non-Gaussianity at increasing mass and redshift Dalal et al. (2008); Pillepich et al. (2010); Grossi et al. (2009).

These formulae can be extended in case of primordial non-Gaussianity from trispectra of $\zeta$. In this case, the kurtosis $S_{4}$ needs to be added Desjacques and Seljak (2010); Chongchitnan and Silk (2010); LoVerde and Smith (2011).

4.1.2 The excursion set approach

Press-Schechter models are known to suffer the cloud-in-cloud problem. The excursion set approach solves this issue and provides an elegant method to count regions above a certain threshold. The goal of this approach is precisely to find regions that are sufficiently overdense on a given smoothing scale, but not on larger ones. To perform this check, one needs to consider the density field $\delta$ at any given point as a function of the smoothing scale. This function looks similar to a random walk, the starting point being the limit of infinitely large smoothing scales, where the overdensity is zero. In this case, the critical density for collapse defines another curve, typically independent of the point of space, namely a“barrier” for the random walk. In the spherical collapse model, this barrier is constant and flat (in $R$), $B(R)=\delta_{sc}$, but one may consider more involved, and more realistic, models. The cloud-in-cloud problem is solved by identifying the largest smoothing scale on which the walk first crosses the barrier. The excursion set ansatz therefore relates the abundance of haloes of mass $M$ to the fraction of random walks that first cross a barrier $B(R)$ on the scale $R$, the halo mass being the one contained inside a sphere of radius $R$,

[TABLE]

With a similar meaning to the multiplicity function, the first crossing fraction $f(s){\rm d}s$ is defined as

[TABLE]

where it is customary in this framework to use the variance $s=\sigma^{2}_{0}(R)$ as the reference variable for the random walk. As we said, the excursion set ansatz consists in requiring that $\delta<B(s)$ for all $s<S$. This is an infinite set of conditions and, in a generic case, calculating the first crossing distribution can be very complicated. Let us assume that, when the walk is strongly correlated, we can instead impose the condition on the one preceding step, that is, that $\delta<B(S-\Delta s)$ for $\Delta s\rightarrow 0$. We can then Taylor expand in $\delta$ and $B$ and get the following condition:

[TABLE]

with ${\rm d}\delta/{\rm d}s\geq{\rm d}B/{\rm d}s$ at $s=S$, which is a condition on the ‘velocity’ of $\delta$ being greater than the tilt of the barrier Appel and Jones (1990); Paranjape et al. (2012); Paranjape and Sheth (2012). This implies that now we have to deal with a bivariate distribution on $\delta$ and its velocity $\delta^{\prime}\equiv{\rm d}\delta/{\rm d}s$, $P(\delta,\delta^{\prime})$, but with the simplification of having only one condition to impose, rather than an infinite number of them. The first crossing distribution reads

[TABLE]

If we now take the limit $\Delta s\rightarrow 0$ and change the integration variable $\delta^{\prime}\rightarrow v=\delta^{\prime}-B^{\prime}$ we get the general formula

[TABLE]

As an example, we can compute Eq. (99) for Gaussian initial conditions and a constant barrier, $B(s)=\delta_{sc}$. In this case, it is convenient to work with the following change of variables

[TABLE]

and we have defined

[TABLE]

We obtain therefore the following first crossing distribution, as a function of $\nu$,

[TABLE]

where the subscript ${\rm G}$ indicates the Gaussian distribution and we have used the conditional probability $p_{\rm G}(x,\delta_{sc})=p_{\rm G}(\delta_{sc})p_{\rm G}(x|\delta_{sc})$.

There are several efforts that extend the excursion set to non-Gaussian initial conditions Inoue and Nagashima (2002); Avelino and Viana (2000); Lam and Sheth (2009a, b); Lam et al. (2009, 2010); Maggiore and Riotto (2010); De Simone et al. (2011); Musso and Paranjape (2012); D’Amico et al. (2011b); Paranjape et al. (2011); Achitouv and Corasaniti (2012); Achitouv et al. (2013); Achitouv and Corasaniti (2012), some of which also relax the assumption of spherical collapse and consider, for instance, ellipsoidal barriers motivated by the fact that the tidal shear of the local density field contributes to the threshold for collapse. Here we will review the main results of Maggiore and Riotto (2010); De Simone et al. (2011, 2011), which build up on results obtained for the Gaussian case Maggiore and Riotto (2010a, b); Ma et al. (2011)

They formulated the excursion set using top-hat filtering in real space, which is the preferred filter to compare with data and simulations Maggiore and Riotto (2010). Such a choice of filter introduces non-Markovianity of the random walks, that is, walks are correlated not only with their direct predecessor, but with the whole preceding path (see also Musso and Sheth (2012, 2014) for another approach). 2. 2.

They gave a physical explanation of the correction to the spherical threshold for collapse $\delta_{sc}\rightarrow\sqrt{q}\delta_{sc}$ which is commonly used to improve fits to simulations (see Sheth (1998); Sheth and Tormen (1999); Sheth et al. (2001); Sheth and Tormen (2002) for the earliest applications ) and showed how it is expected that the factor depends on the halo finder used Maggiore and Riotto (2010b). 3. 3.

Non-Gaussianities also introduce non-Markovian corrections, they were able to include them in their framework using the results from the Gaussian case with non-Markovianity Maggiore and Riotto (2010). 4. 4.

They extended the formalism to a generic moving barrier $B(\sigma)$ in De Simone et al. (2011).

These findings are summarized in the following prediction for the non-Gaussian multiplicity function

[TABLE]

where $q=1/(1+D_{B})$ being $D_{B}$ the diffusion coefficient of the barrier, which advocates 2. of the previous list,

[TABLE]

is the correction due to the use of top-hat filtering (addressing 1.), $\Gamma(0,x)$ is the incomplete Gamma function and

[TABLE]

is the correction due to non-Gaussianity (addressing 3.), being $\mathcal{U}$ and $\mathcal{V}$ first and second derivatives of the skewness $S_{3}$, respectively.

Eq. (103) is not generalized to account for moving barriers ( improvement 4. in the list above ), such a generalization can be found in De Simone et al. (2011), where they also used an improved method based on the saddle-point approximation proposed by D’Amico et al. (2011b). Another series of papers analyzes drifting diffusing barriers in the context of the excursion set approach Achitouv and Corasaniti (2012), including also an extension to a class of primordial trispectra Achitouv and Corasaniti (2012). In general, the introduction of a consistent barrier for collapse is a non-trivial task because it requires a detailed modeling of halo formation, which at present day is not understood, being a fully non-linear process.

4.1.3 The excursion set peaks model

It is rather natural to think that, to some degree of approximation, virialized haloes at late time may have formed out of maxima of the initial density field Bardeen et al. (1986). This amounts to saying that, if we trace back halos of a certain mass $M$ to the initial conditions, they sit on a peak of the initial density field. The corresponding Lagrangian patch is called “proto-halo” and it defines the halo Lagrangian scale $R\propto(M/\bar{\rho})^{1/3}$. This assumption has been extensively analyzed in numerical N-body simulations and shown to work rather well for a wide range of masses Ludlow and Porciani (2011); Ludlow et al. (2014); Elia et al. (2011); Borzyszkowski et al. (2014). Several authors have developed the technical tools to deal with local density maxima in the presence of non-Gaussian initial conditions Bardeen et al. (1986); Catelan et al. (1988); Gay et al. (2012); Codis et al. (2013); Lazeyras et al. (2016); Uhlemann et al. (2018). Here we want to collect all their progresses and include them in the effort of combining the peak model and the excursion set approach. Only recently it has been realized that the two methods can be unified into one Excursion Set Peaks (ESP) model by imposing that peaks of the density field at a given smoothing scale satisfy a first crossing condition, in such a way as to have a first crossing distribution for peaks Paranjape and Sheth (2012). A full review and detailed description of the peak model and its extension to the ESP is nicely presented in Desjacques et al. (2018). Here we highlight the main results.

The Gaussian case

We can generally define the number density of points that are local maxima of the Gaussian field $\delta(\mathrm{\bf x})$ as a point process,

[TABLE]

where $\delta_{\rm D}$ is the 3-dimensional delta function. The continuous field $\delta(\mathrm{\bf x})$ can be Taylor expanded around a peak position $\mathrm{\bf x}_{\rm pk}$,

[TABLE]

where we imposed the peak condition on the first derivative, $\nabla_{i}\delta(\mathrm{\bf x}_{\rm pk})=0$. In a similar manner, we can Taylor expand $\nabla_{i}\delta(\mathrm{\bf x})$ around the peak as

[TABLE]

Using Eq. (107) and Eq. (108), we can express the delta functions of Eq. (106) as

[TABLE]

provided that the hessian $\nabla_{i}\nabla_{j}\delta$ is definite negative, to ensure to be looking at maxima, and non singular in $\mathrm{\bf x}_{\rm pk}$.

The identification of proto-halos inserts a scale into the problem, as we want to model the clustering of halo centres, and not their internal substructure. We therefore need to smooth the DM field on the Lagrangian size $R_{s}$. Moreover, we may impose a high enough threshold on the value of the density at the peak, that is, a peak height, to make sure that we are dealing with regions sufficiently dense as to favour virialization, in the spirit of Kaiser (1984).

As a result, the number density of peaks of height $\nu^{\prime}$ at position $x$ can be expressed in terms of the smoothed field $\delta_{s}$ and its derivatives as Bardeen et al. (1986)

[TABLE]

where $R_{*}=\sqrt{3}\sigma_{1s}/\sigma_{2s}$ is the characteristic radius of the peak and we have conveniently defined the normalized peak density and its derivatives as

[TABLE]

and we will adopt the following notation for the variance of the smoothed density field, (linearly extrapolated to present-day) and its derivatives,

[TABLE]

We have imposed explicitly that the stress tensor $-\zeta_{ij}$ is definite negative through the positivity of $\lambda_{3}$, its lowest eigenvalue.

The number density $n_{\rm pk}$ can be used, in principle, to compute any N-point correlation functions among the peaks of the density field by ensemble averaging products of $n_{\rm pk}$,

[TABLE]

and in the Gaussian initial conditions considered here, multivariate normal distributions are assumed to perform the ensemble average. The case of $N=1$ is the averaged peak number density

[TABLE]

being $A=1,...,6$ since $\zeta$ is symmetric and the joint probability distribution $P_{1}$ is given by

[TABLE]

where $M$ is the corresponding covariance matrix. Owing to rotational invariance, we have regrouped the set of variables into the following vector

[TABLE]

where $J_{1}=-{\rm tr}(\zeta_{ij})$ is the peak curvature, $J_{2}=\frac{3}{2}{\rm tr}(\bar{\zeta}^{2}_{ij})$ and $J_{3}=\frac{9}{2}{\rm tr}(\bar{\zeta}^{3}_{ij})$, being $\bar{\zeta}_{ij}=\zeta_{ij}+\delta_{ij}\,J_{1}/3$ the components of the traceless part of the Hessian. We can then factorize the exponential into the following form

[TABLE]

and we define the correlation

[TABLE]

This decomposition allows to write Eq. (117) in a more compact way as a product of a bivariate Gaussian $\mathcal{N}$ in the variables $\nu$ and $J_{1}$ and 3- and 5- degrees of freedom $\chi$-squared distributions in $3\eta^{2}$ and $5J_{2}$ respectively,

[TABLE]

and $P(\boldsymbol{\Omega})$ represents the probability distribution of the five remaining d.o.f.. Since they are all angular variables, they do not generate bias factors because the peak (and halo) overabundance can only depend on scalar quantities (see e.g. Catelan et al. (1998); McDonald and Roy (2009)). The variable $J_{3}$ is uniformly distributed and constrained to satisfy $J_{3}^{2}\leq J_{2}^{3}$ by the fact that $\zeta_{ij}$ is symmetric.

We now want to apply the excursion set ansatz to the statistics of peaks of the density field. This amounts to impose to select only those peaks which are found at a certain scale $R$ which have a smaller height on the next larger smoothing scale, that is,

[TABLE]

where we just adapted Eq. (96) to the variable $R_{s}$, with the condition that $\delta^{\prime}-B^{\prime}\leq 0$, where now the prime indicates derivative with respect to $R_{s}$. Hence, the corresponding ESP discrete number density can be written as

[TABLE]

where we have defined the variable $\mu=B^{\prime}-\delta^{\prime}$ and extended the vector of variables to $\boldsymbol{\omega}=(\mathbf{y},\mu)$. Combining Eq. (116) and Eq. (97) we get Desjacques et al. (2013)

[TABLE]

and using a similar approach to Eq. (184), we can calculate the ESP multiplicity function to be

[TABLE]

where $V=m/\bar{\rho}$, $\sigma^{\prime}_{0s}={\rm d}\sigma_{0s}/{\rm d}R_{s}$ and $\mathcal{G}^{(1)}_{1}(x)$ is defined in appendix A.1.

The multiplicity function of Eq. (125) is a theoretical prediction of the halo mass function which does not depend on any free parameter. Its limits of applicability are directly related to the assumptions it is based on, namely

• •

Halos form out of peaks of the Lagrangian density field,

• •

The barrier for collapse is described as a sole function of the density field.

• •

Walks are correlated enough on large scales such that we can impose the up-crossing condition on the one preceding step only, rather than the full walk.

The most important extension to be worked out on this model is the inclusion of the effect of tidal shear on the gravitational collapse of halos. The formalism has been wrote down in Desjacques et al. (2018) for a generic collapse barrier which depends on the tidal field shear field. The next step would be to understand the particular form which this barrier should have, in order for the model to provide a prediction on the halo mass function. As argued in Desjacques et al. (2018) and in Sec. §5.2.3, this extension is particularly relevant for primordial non-Gaussianity.

The non-Gaussian extension

The goal is now to embed a non-Gaussian initial PDF for peaks and therefore a non-Gaussian ESP mass function. The important point here is that in order to compute a non-Gaussian probability through an expansion similar to Eq. (80), we need to expand all the variables in $\boldsymbol{\omega}=\{\nu,J_{1},\mu,3\eta^{2},5J_{2},J_{3}\}$. This expansion takes the name of Gram-Charlier expansion and it has been developed in a number of papers Pogosyan et al. (2009); Gay et al. (2012); Gong and Yokoyama (2011); Lazeyras et al. (2016).

For generic non-Gaussian initial conditions we define $P_{\rm NG}=P_{\rm G}[1+\delta P_{\rm NG}]$ and the Gram-Charlier111111For cumulants up to fifth order it can be demonstrated that the Gram-Charlier expansion coincides with the Edgeworth expansion. It is worth noticing that, especially for large non-Gaussianities, the Gram-Charlier expansion might not converge, see for example Blinnikov and Moessner (1998). expansion up to third order in $\delta$ gives

[TABLE]

where we define the Hermite polynomials

[TABLE]

the generalized Laguerre polynomials

[TABLE]

and the correlations $\langle\cdot\rangle_{\rm GC}$ as

[TABLE]

where we have assumed $m\leq 1$ and the average is weighted by the ESP mass function to ensure that only peaks of the density field are selected. The non-Gaussian correction to the mass function is

[TABLE]

Here we notice an important feature of the expansion: each Hermite, or Laguerre, polynomial that multiplies the Gaussian mass function and PDF, when integrated over $\boldsymbol{\omega}$ defines itself a bias parameter for the Gaussian case, for instance

[TABLE]

Hence we have found that each term in Eq. (131) defines a third order bias multiplied by the corresponding “skewness” of the type $\langle X^{3}\rangle_{\rm GC}$. If we were to include terms involving the kurtosis, they would generate fourth order bias Gaussian parameters and so on. We refer to the Appendix A.2 for the full expression for $\delta n^{\rm NG}_{\rm ESP}$, which is rather lengthy, but straightforward to calculate. Here we quote only the first three terms

[TABLE]

where we define the skewnesses as

[TABLE]

where $g$ and $\chi$ vary for each variable $X$ considered.

4.2 Numerical approaches

In parallel with analytic approaches, N-body simulations have been exploited both to test theoretical predictions and to find semi-analytic or fully fitted phenomenological formulas to be used in real data. The earliest fits to the non-Gaussian mass function from simulations with non-Gaussian initial conditions Dalal et al. (2008); Pillepich et al. (2010) focused on local-type non-Gaussianity with an amplitude of $\mathcal{O}(f_{\rm NL})\sim 100$s. In Pillepich et al. (2010) they fitted the multiplicity function as

[TABLE]

where $A,B,C$ and $D$ each follow the functional form

[TABLE]

This fit showed agreement within $\sim 10$% of the measurements.

In Achitouv and Corasaniti (2012); Achitouv et al. (2014), the halo mass function from excursion set theory with a drifting diffusive barrier was tested against N-body simulations. Their prediction is an extension of Eq. (103) to include a stochastic barrier with linearly drifting average up to next-to-leading order. The idea is to allow for stochastic deviations from the usual spherical collapse threshold value $\delta_{sc}$, parametrizing the deviation with two quantities: the rate $\beta$ at which the collapse threshold deviates on average from the spherical collapse $\delta_{sc}$ and the amplitude $D_{B}$ of a constant scatter around the average collapse threshold at a given mass scale. They compare their model with measurements of the halo mass function in simulations with local-type non-Gaussianity in the initial conditions, first fitting $\beta$ and $D_{B}$ and then calibrating them on the Gaussian simulations. Using the first fit, their prediction agree within $5$% with simulations. When calibrating on the Gaussian simulations, they find dependence of the parameters on non-Gaussianity, indicating that the stochastic nature of the barrier is correlated with it.

An important study of the large-scale density field was performed in Mao et al. (2014) using N-body simulations in order to understand whether measurements of the moments of large-scale structure can yield constraints on primordial non-Gaussianity. They found that most of the information is included in the variance of the galaxy density field, because of the effect of the scale-dependent bias, which we will discuss in the next Sec. §5. The result on the variance was recently further refined in Nusser et al. (2018), where they use maxima and minima counts of the halo density field of N-body simulations to show that primordial non-Gaussianity of the local-type could be constrained at the level of $f_{\rm NL}\approx 10$ with the Euclid survey using this method. Although these figures are not competitive with direct constraints from the power spectrum of galaxies (see for instance Amendola et al. (2013)), the estimate of Nusser et al. (2018) is still preliminary and should be refined.

In Adhikari et al. (2014), N-body simulations with primordial non-Gaussianity generated by a range of two-field models of inflation were run. The measurements were compared to semi-analytic estimations of the non-Gaussian mass function based on the Edgeworth expansion. The type of primordial non-Gaussianity they consider is such that cumulants of higher order than the skewness are more important than in the single-field local case. Their comparison with data shows that different prescriptions for the scaling of higher moments give appreciably different results. The analysis of multi-field inflation in N-body simulations proves an important point on being careful when truncating the Edgeworth expansion. Given that recent theoretical progress has demonstrated that the local-type non-Gaussianity for single-field models is not observable Tanaka and Urakawa (2011); Baldauf et al. (2011); Creminelli et al. (2011); Pajer et al. (2013); Dai et al. (2013); Camera et al. (2015); Bartolo et al. (2016); de Putter et al. (2015) (see discussion in the next Sec. §5.1.3), this type of searches should be investigated further.

4.3 Final remarks of this section

In this section, we have focused on the effect of primordial non-Gaussianity on the probability distribution of dark matter halos. We have showed how the strongest effects are in the tails of the probability distribution. For this reason, studies on the very high mass tail of what observed in a galaxy survey, such as extreme-mass clusters Chongchitnan and Silk (2012); Chongchitnan (2015) and x-ray detected clusters Roncarelli et al. (2010); Sartoris et al. (2010); Shandera et al. (2013), have been investigated as well. Similarly, very large voids probe the low-density tail of the PDF. Early studies of N-body simulations showed that this could be a promising venue for large values of $f_{\rm NL}$ of the local type Grossi et al. (2008). On the theoretical side, calculations within the context of the extended Press-Schechter mass function were perfomed in Kamionkowski et al. (2009); Song and Lee (2009); Lam et al. (2009) and later using the excursion set formalism D’Amico et al. (2011a). It is also possible to constrain directly the dark matter density field using weak lensing measurements Jimenez and Verde (2009); Lam et al. (2010); Pace et al. (2011); Marian et al. (2011); Jeong et al. (2011); Hilbert et al. (2012).

Despite almost 40 years of research into the possibility of constraining primordial non-Gaussianity using the probability distribution of dark matter, and dark matter tracers, in the Universe, it remains a great challenge to obtain competitive figures as compared to other probes, such as the CMB and correlation functions of dark matter tracers (see the following section §5). The main reason for this is that most of the relevant information is in the high-mass tail of the halo (or possibly other tracers) distribution, where we lack the statistics to have high signal-to-noise ratios. Moreover, current limits on non-Gaussianity constrain the correction to the probability distribution to be very small: for instance, a local-type non-Gaussianity of order unity implies a correction at the subpercent level in the halo mass function.

5 Imprints of primordial interactions on two-point halo statistics

The interactions taking place during inflation that we introduced in Sec. §2.2 manifest themselves as non-trivial correlation functions of the primordial curvature perturbation, which provides the seed for structure formation. We have focused in particular on three-point correlation functions: the natural choice for looking for these primordial correlations in the large-scale structures of the late universe would be the three-point correlation function of galaxies, which is sourced by the primordial correlator at tree-level in a similar way as for the dark matter case of Eq. (86).

It came as a breakthrough the finding that the two-point function of halos on large scales contains a great deal of information about the class of primordial interactions whose bispectrum shape peak in the squeezed limit. The discovery of a characteristic scale dependence at large-scales in the linear bias of halos came as a surprise since it was a generic prediction of clustering models with Gaussian initial conditions that the bias at large scales was scale-independent. As we will show in this section, there are two main features that make this observable so promising to constrain primordial interactions: firstly, the presence of primordial non-Gaussianity manifests the strongest on the largest scales. In this regime, non-linear biasing, non-linear gravitational interaction and red-shift space distortions have negligible effects. Secondly, it manifests as a peculiar scale dependence, which is not produced by any other phenomenon on those scales. The particular scaling depends on the primordial bispectrum shape that generates it, hence a good constraint on the scale dependence could discriminate among different models of inflation.

This finding opened an entire new stream of research on the effects of PNG on the power spectrum of galaxies and it provides one of the most promising observables in LSS as of today. This section is devoted in reviewing recent efforts both on the theoretical side and on the numerical side.

Plan of the Section. We start by introducing the imprint of local-type non-Gaussianity on the halo power spectrum (almost) as it was discovered, with a derivation using the peak-background split ansatz in Sec. §5.1. We then proceed by generalizing to any model of primordial non-Gaussianity in Sec. §5.2, describing different methods for modeling structure formation and halo clustering in the presence of such an effect. We then overview the current status of numerical analyses in Sec§5.3, complementing the analytic part. We conclude with final remarks in Sec.§5.4.

5.1 The breakthrough: scale-dependent bias from local PNG

In Dalal et al. (2008), the effect of primordial non-Gaussianity in the halo power spectra was measured using N-body simulations with non-Gaussian initial conditions of the local type. The measurements exhibit a strong $1/k^{2}$ scaling at large scales. First physical interpretations and predictions were given in the context of the high peaks approximation Matarrese and Verde (2008); Valageas (2010); Shandera et al. (2011), multivariate bias expansions McDonald (2008); Giannantonio and Porciani (2010), and using a peak-background split (PBS) ansatz Slosar et al. (2008); Afshordi and Tolley (2008); Schmidt and Kamionkowski (2010); Desjacques et al. (2011); Smith et al. (2012); Scoccimarro et al. (2012). As a first step into the study of this topic, in this section we provide a pedagogical derivation for the case of local-type primordial non-Gaussianity in the context of the peak-background split ansatz, loosely following Slosar et al. (2008).

5.1.1 Derivation with the peak-background split ansatz

The peak-background split approach Kaiser (1984); Bardeen et al. (1986) gives a very clean physical interpretation of the effect and it provides a rather model independent prediction. The starting point of a peak-background split ansatz is separating long- and short-wavelengths modes of the density field $\delta$. A typical short mode can be thought of as the Lagrangian radius of a proto-halo, i.e. the scale relevant for the formation of virialized objects. The long modes act as a local rescaling of the amplitude of short ones. The basic assumption for being able to perform such a separation is that the two regimes are decoupled. This is the case for Gaussian distributed fields. Here, we want to consider the case of non-Gaussian initial conditions, which by definition introduce coupling between long and short modes already at the primordial level. This coupling is manifest when considering the local ansatz for the primordial perturbations,

[TABLE]

If the initial conditions are non-Gaussian, as in Eq. (137), the separation of scales has to be performed not at the level of the density field $\delta$, but on Gaussian primordial fluctuations

[TABLE]

where $\ell$ indicates long modes and $s$ short ones. Applying the split to Eq. (137) we get

[TABLE]

On sufficiently large scales, the relation

[TABLE]

holds, while the short-wavelengths modes of the density field are affected, at first order, as

[TABLE]

which can be interpreted as a local modulation of the amplitude of matter fluctuations by the long modes proportional to $f_{\rm NL}$,

[TABLE]

As a result, the local number density of halos of mass $M$ is not a function of $\delta_{\ell}$ only, but also of this local modulation $\hat{\sigma}_{s}$ and, more in general, on all cumulants of the density field. The expansion in Lagrangian space for the halo overdensity reads, at first order in $\delta_{\ell}$, as

[TABLE]

From this expression we can define the first order bias in Fourier space in the presence of primordial non-Gaussianity

[TABLE]

being

[TABLE]

the usual Gaussian Lagrangian bias and

[TABLE]

the amplitude of the so called scale-dependent bias. Here $\sigma_{8}$ is the variance of the density field for a smoothing at $R=8$ Mpc/h, the typical reference amplitude of matter fluctuations. This result is general, at first order, as we did not make assumptions other than the scale separation in $\zeta_{\rm G}$.

The characteristic feature of this additional bias contribution at first order is that it scales as $\Delta b_{1}^{\rm NG}\propto 1/k^{2}$ at large scales, since $\mathcal{M}(k)\propto k^{2}$ in this limit. It implies that the two-point correlation function of halos at large separations should show a clear enhancement, or suppression if $f_{\rm NL}^{\rm loc}$ is negative, with respect to the Gaussian case. In Fourier space, we therefore have

[TABLE]

where here we should remember that also the matter auto-power spectrum has non-Gaussian corrections, according to Eq. (84). As first argued in Slosar et al. (2008), the amplitude of this non-Gaussian correction depends on the halo formation history, suppressing the amplitude of recently formed halos and enhancing it if they are formed early. We will not discuss the effect here and defer the reader to the available literature Slosar et al. (2008); Reid et al. (2010).

5.1.2 Universality of the mass function

The amplitude $b_{1}^{\rm NG}$ of Eq. (146) cannot be directly measured from galaxy surveys, as we do not have multiple realizations of the Universe with varying $\sigma_{8}$. The simplest way around this problem is to assume universal mass functions of the form of Eq. (89), as in the case of the Press & Schechter mass function, Eq. (88). In this case, the logarithmic derivative of the mass function with respect to the normalization amplitude $\sigma_{8}$ takes the simple form

[TABLE]

The above expression also assumes spherical collapse. In this way, the linear bias in the presence of local-type PNG can be characterized via a single bias parameter $b^{\rm G}_{1}$.

The assumption of universality of the mass function has been thoroughly investigated for Gaussian initial conditions in the last decade (see for instance Tinker et al. (2008); Reed et al. (2013); Despali et al. (2015)), but not as precisely for non-Gaussian inital conditions. Indeed, studies of N-body simulations have shown that the assumption is not as reliable as needed by the increasing precision required by data. We expand on this studies in the next Sec. §5.3.1.

5.1.3 The single-field consistency relation

Single-field models of inflation generate a minimal amount of non-Gaussianity, as explained in §2.2.1. The minimal interaction of a single scalar field with gravity produces a minimal amount of non-Gaussianity, Eq. (2.2.1), of order $\mathcal{O}(f_{\rm NL})\sim\mathcal{O}(\epsilon,\eta)$ . For squeezed configurations of the bispectrum, the inflationary consistency relation Eq. (34) states that all single-field models of inflation are constrained to have small non-Gaussianity, $\mathcal{O}(f_{\rm NL})\propto n_{s}-1$. This non-Gaussianity can be connected to a local-type one by an appropriate matching Maldacena (2003), so that one would be tempted to state that we are expected to see a minimal amount of non-Gaussianity in the scale dependent bias. However, following an interesting debate Tanaka and Urakawa (2011); Baldauf et al. (2011); Creminelli et al. (2011); Pajer et al. (2013); Dai et al. (2013); Camera et al. (2015); Bartolo et al. (2016); de Putter et al. (2015); Bravo et al. (2018a, b), it has been shown that the single-field consistency relation is equivalent to the statement that long-wavelength perturbations are not physically coupled to small-scale ones and therefore there is no scale-dependent bias for single-field models of inflation. The argument follows from studying how to connect the result of Maldacena (2003), which was computed in comoving gauge at the epoch of inflation, with actual observations of galaxy number counts at the present day. The positive side of this finding is that any detection of a scale-dependent bias on large scales of the type of Eq. (144) would rule out the possibility that primordial perturbations were generated at early times by a single field.

In fact, the statement is not only restricted to local-type non-Gaussianity, but to all squeezed limits of inflationary bispectra: one way to show it is to remove the unphysical contributions from any bispectrum by calculating it in Conformal Fermi Coordinates Dai et al. (2015a, b). An example in which this subtraction is done for the resonant running model of Eq. (51) is found in Cabass et al. (2018).

5.2 Analytic approaches

After warming up with the effect of local PNG on the clustering of halos in the previous section, here we discuss how to extend it to any type of primordial interaction. We will begin by introducing a framework which allows for minimal assumptions on the physics of halo formation and collapse using an effective field theory approach to LSS. This treatment is extremely powerful in characterizing the most promising signatures from the early Universe. Subsequently, we proceed with applying a few well-motivated ansatz on the physics of structure formation at late times thereby increasing the predictive power of the theory, but losing control of theoretical errors. The trade-off between model-independence and predictivity plays a major role in constraining PNG using the power spectrum of tracers from galaxy surveys.

5.2.1 The squeezed limit of the primordial bispectrum

The finding that the presence of a coupling between long and short modes in the initial conditions modulates the local number density of halos and therefore their clustering at large scales has great importance: it implies that there is a specific kinematic regime of the initial conditions in which halo clustering at large scales is particularly sensitive to early Universe physics. This kinematic regime is the squeezed limit of the bispectrum, which is precisely the triangle shape by which the power spectrum on short modes is modulated by long ones.

By defining long and short modes for a bispectrum triangle as

[TABLE]

we can expand the primordial bispectrum of $\zeta$ for $k_{\ell}\ll k_{s}$ and get

[TABLE]

where $\mathcal{L}_{J}$ are Legendre polynomials which appear only with even orders since in the squeezed limit ${{\bf k}}_{1}\approx-{{\bf k}}_{2}$, removing the degree of freedom for odd number of exchanges of ${{\bf k}}_{s}$ into $-{{\bf k}}_{s}$. The factor $A_{J}(k_{s},k_{\ell})$ encodes information about the primordial interaction we want to describe. Tables 1 and 2 gathers the values of $A_{2J}$ for the models introduced in Sec. §2.2, excluding interactions involving features: these cases generically break scale-invariance and need to be handled with extra care. We expand on this in the next paragraph.

Notice that, as long as we are constraining non-Gaussianity using the halo power spectrum only, any information on the angle between ${{\bf k}}_{s}$ and ${{\bf k}}_{\ell}$ is lost. One needs to look at the bispectrum Baldauf et al. (2011); Pollack et al. (2012); Sefusatti et al. (2012); Yokoyama et al. (2014); Tasinato et al. (2014); Dizgah et al. (2015); Hashimoto et al. (2016); Tellarini et al. (2016); Hashimoto et al. (2016); Yamauchi et al. (2017); Di Dio et al. (2017); Chiang et al. (2017); An et al. (2018); Moradinezhad Dizgah et al. (2018); de Putter (2018); Moradinezhad Dizgah et al. (2018) or, alternatively, to galaxy shapes Chisari and Dvorkin (2013); Chisari et al. (2016) and intrinsic alignments Kogai et al. (2018) to be sensitive to it.

Scale-invariant primordial bispectra

The coefficients $A_{J}$ of the expansion contain information on the specific inflationary model from which the interaction giving rise to the bispectrum generates. It is instructive to start from the sub-class of models we defined in Section §2.2 which are scale-invariant. It is straightforward to show that

[TABLE]

for scale-invariant bispectra. We can then explicitly write the modulation of a specific long mode $k_{\ell}$ on the local power spectrum around a small Lagrangian patch centered in ${{\bf q}}$ as

[TABLE]

where in the second term in the square brackets we have wrote explicitly the second order Legendre polynomial $\mathcal{L}_{2}(l_{\ell},k_{s})$, factorizing out short modes, and the ellipses indicate higher order terms in $J$. Here it is important to stress that the modulation is imprinted in the initial conditions for structure formation, at the Lagrangian position ${{\bf q}}$, as compared to the Eulerian position at which we usually make observations of the (galaxy) density field, ${{\bf x}}$ Giannantonio and Porciani (2010); Baldauf et al. (2011); Assassi et al. (2015). Eq. (152) suggests a generalization of Eq. (142) by defining operators such as

[TABLE]

for $J=0$ and

[TABLE]

for $J=2$. The integration allows to account for the collective contribution from all long modes and has generically support for $k_{\ell}<\Lambda$, where $\Lambda$ is an arbitrary scale splitting long from short modes. Higher orders can be calculated by using $J-$th order trace-free projection operators with respect to ${{\bf k}}_{\ell}$. These fields can be evaluated at the Eulerian position ${{\bf x}}$ by expanding around ${{\bf x}}({{\bf q}},\tau)={{\bf q}}+{\bf s}({{\bf q}},\tau)$, being ${\bf s}$ the Lagrangian displacement, in perturbation theory

[TABLE]

Using a similar reasoning as in the local case, we then have a scale-dependent bias term in the halo auto-power spectrum of the form

[TABLE]

for the case of $J=0$ at leading order. Here $b_{\psi_{(\Delta)}}$ is the bias parameter corresponding to the operator of Eq. (153), which is a generalization of Eq. (146). Using a similar argument on the separation of scales leading to the expression of Eq. (146), the rescaling of the short-wavelengths modes by the long ones in this case reads

[TABLE]

where we have defined $\epsilon=a_{0}\psi_{(\Delta)}$. Consequently, we get

[TABLE]

The peak-background split is extremely powerful to derive expressions for bias parameters based on the response of the local number density of halos to the variation of a particular physical quantity. These relations are inherently renormalized, physical quantities that can be directly used to compute n-point correlation functions of halos and tracers in general. In order to use them as predicted paramaters, a specific modelling of halo formation needs to be implemented, as we discuss in the next section.

The local-type non-Gaussianity corresponds to $a_{0}=12/5f_{\rm NL}$ and $\Delta=0$. For any $\Delta>0$ the enhancement on large scales is reduced. For instance, non-Gaussianity from DBI inflation has $\Delta=2$, as indicated in Tab. 1. In this case, the non-Gaussian contribution to the halo power spectrum is scale-independent and thus degenerate with the linear bias $b_{1}$. This can be understood by looking at Eq. (153): $\psi_{(2)}\propto\nabla^{2}\zeta$ is directly proportional to the matter density. Similarly, higher values of $\Delta$ give contributions which are degenerate with bias parameters related to higher-derivative operators of the form $\nabla^{\Delta-2}\delta$. The intermediate range $0\leq\Delta\leq 2$ is explored by models with massive particles of various spins produced during inflation, as discussed in Sec. §2.2.2.

The second term in Eq. (156) is the leading order correction to the linear bias in the presence of generic scale-invariant primordial bispectra. For $J>0$, the correction comes at higher orders: for instance for $J=2$, the operator $\psi_{(\Delta)}^{ij}$ needs to be multiplied by some operator $O^{ij}$ in order to appear in the bias expansion of Eq. (5.1.1), simply because the halo overdensity $\delta_{h}$ is a scalar quantity. We make this statement more precise in a few paragraphs.

Let us first discuss how to extend the squeezed limit expansion of the bispectrum of Eq. (150). The extension is motivated by the fact that the limit is strictly valid until the modes $k$ we observe do not affect the local formation of halo. This can be diagnosed by checking how the variance $\sigma^{2}$ smoothed on the small-scale $R_{s}$ is affected by $k$, and it approximately coincides with modes at the peak of the matter power spectrum, $k\sim 0.02$ h/Mpc. Given that observational data hardly constrain modes much larger than $k\sim 10^{-3}$ h/Mpc, it is important to understand how much can we exploit from analytic methods before resorting to N-body simulations.

If we restrict to the scale-invariant case, we can perturbatively add terms in higher powers of $k_{\ell}/k_{s}$ and get

[TABLE]

where $N$ runs also only on even powers. Since the correction comes in the same form as the leading term at each order in $J$, we then just have to add operators of order $\delta^{2N}$ to the bias expansion. For instance, for $J=0$ we get for the halo auto-power spectrum

[TABLE]

As mentioned in the squeezed limit case, these higher order terms are potentially degenerate with other bias terms: in addition with bias operators from higher-derivative terms present even with Gaussian initial conditions, in this case there is degeneracy also when combining lower order scale-dependent bias from non-Gaussian initial condition with Gaussian ones. For instance, if we combine $\psi_{(\Delta)}$ with a bias generated by the operator $\nabla^{2}\delta$, we get a contribution approximately of the same order as $b_{\psi_{(\Delta+2)}}$. Moreover, the transfer function $\mathcal{M}(k)$ also contributes in this expansion when considering the full matter bispectrum in the squeezed limit Schmidt (2013) with terms which become relevant at the matter-radiation equality scale, $k_{\rm eq}$.

Beyond scale-invariance

In the previous paragraph we have made the approximation of a scale-invariant primordial bispectrum, which excludes non-Gaussianities generated by features during inflation, or any case in which a particular scale is involved in the primordial interaction which generates the bispectrum (or higher-order correlators). The effect of models with features on the scale-dependent bias was first investigated in Cyr-Racine and Schmidt (2011). The results above can be generalized to the case of features, here we consider a simple example to illustrate the calculation. Let us consider a simplified form of the sharp feature shape of Eq. (48) by defining

[TABLE]

being $A$ some overall amplitude and $k_{*}$ the scale of the feature. Here we are neglecting the constant phase of the oscillations, the power law and damping factor of Eq. (48) as they are not relevant for the argument. Let us also suppose that this bispectrum is generated by a single-field model of inflation, e.g. a sharp feature on the potential (see Bartolo et al. (2013) for a treatment in the context of the EFT of inflation). According to what discussed in §5.1.3, when expanding in the squeezed limit, terms of order $(k_{\ell}/k_{s})^{0}$ and $(k_{\ell}/k_{s})$ are not observable. Let us suppose therefore that we have made the calculation in CFC and removed such terms. The squeezed limit at leading order will then look like

[TABLE]

where we have also integrated over the angle between $k_{s}$ and $k_{\ell}$. We have now to distinguish two regimes: the ultra-squeezed limit, in which $k_{\ell}\ll k_{*}$, can be described by the scale-invariant ansatz of Eq. (151) with $\Delta=2$. On the other hand, as $k_{\ell}$ approaches $k_{*}$, the product of $k^{2}_{s}/k^{2}_{*}\times k^{2}_{\ell}/k^{2}_{s}$ from the second term of Eq. (162) becomes order unity. It is straightforward to show that going to higher orders does not help, and indeed each order becomes equally important. The method therefore fails unless we are able to resum the expansion. In Cabass et al. (2018) it was shown in the case of resonant running corresponding to the shape of Eq. (51) that this can be done under certain assumptions on the dependence of halo abundance on the statistics of short modes. We consider such methods in the Sec. §5.2.3.

5.2.2 The bias expansion

The extension of the scale-dependent bias signature from local-type non-Gaussianity to a wider range of models provides the basic information to write a generic bias expansion with which we can systematically compute correlators of the halo density field. Here we highlight the main steps required to include PNG to the bias expansion, while we defer details on the Gaussian terms to the recent review Desjacques et al. (2018).

A perturbative bias expansion connecting the distribution of halos (or tracers in general) to the dark matter field relies on the assumption that the formation of halos takes place over a very long period of cosmic time, but it is affected by a relatively short range of spatial scales. In other words, if we assume that the halo formation process develops within some scale $R_{*}$, which is typically the Lagrangian radius of the halo, and we study the halo distribution on much larger scales, then we can effectively write this process as local in space (but non-local in time). In this picture, the equivalence principle applied to the “free-falling” small region of size $\sim R_{*}^{3}$ imposes that only second derivatives of the gravitational potential appear in the expansion. Therefore, we can write the expansion

[TABLE]

where $b_{O}$ are bias parameters related to the operators $O$ which are all the scalar combinations that can be made, order by order, out of the matter density field $\delta$ and the tidal field121212In the context of structure formation, it is conventional to work with the gravitational potential $\Phi$, which in matter domination is simply related to the curvature perturbation by $\Phi=3/5\zeta$. $\partial_{i}\partial_{j}\Phi$. Because of non-locality in time, that is, the fact that the halo formation process happens over a long range of time, the terms of the expansion will in general include convective time derivatives Mirbabayi et al. (2015); Assassi et al. (2015). The square brackets indicate that the operators need to be renormalized, as typically, excepting the simple linear term, they involve products of fields evaluated at the same point, generating ultraviolate divergencies McDonald (2006); McDonald and Roy (2009); Assassi et al. (2014); Schmidt et al. (2013); Senatore (2015); Angulo et al. (2015). Moreover, small deviations from locality can be accounted for perturbatively in terms of higher-derivative terms such as $\nabla^{2}\delta$ to be added to the expansion.

The final ingredient to be added to this expansion is stochasticity, i.e. the impact of small-scale perturbations on the formation of galaxies Dekel and Lahav (1999); Taruya and Soda (1999); Matsubara (1999). We can include these effects using techniques from the effective field theory approach, introducing stochastic fields in the expansion which do not correlate with long-wavelength mode fields. Stochasticity in the presence of PNG becomes especially relevant in the case in which non-Gaussianity is sourced in the initial conditions by a superposition of two or more fields (see Secs. §2.2.2 and §2.2.3) . These models are characterized by large trispectra, which are not discussed here, therefore we will neglect stochasticity. A complete treatment, including the case of primordial non-Gaussian initial conditions, can be found in Desjacques et al. (2018), while original papers date back a few years earlier Smith and LoVerde (2011); Baumann et al. (2013)131313See also McAneny et al. (2018) for a recent treatment in the context of quasi-single field inflation..

The assumption of separation of long- and short-scales is no longer valid in the case of primordial non-Gaussianity: indeed, we have seen above that more operators derived from $\zeta$, or equivalently $\Phi$, itself are allowed in the expansion. The coupling in this case takes place at the Lagrangian position ${{\bf q}}$, since it is imprinted in the initial conditions of structure formation. Following Desjacques et al. (2018) we quote here the full bias expansion with PNG up to order $J=2$ in the scale-invariant approximation of Eq. (151)

[TABLE]

where we defined the tidal field

[TABLE]

The first line is leading order in $\delta$ and the non-Gaussian correction $\psi$, where we have considered the leading order in Eq. (155), while the second line is second order. We do not include terms which are second order in $\psi$ given that from current CMB observations we expect PNG corrections to be small Akrami et al. (2019). We have only implicitly listed stochastic terms: at first order one needs to introduce an additional field, $\epsilon_{\psi}$ with respect to the Gaussian case, which has only $\epsilon$ Assassi et al. (2015). We have truncated the expansion on three different levels:

• •

The squeezed limit expansion $a_{J}$ with $J>2$. This corresponds to contributions from higher-spin fields. As highlighted in Sec. §2.2.2 a particular feature of this interaction is that the primordial bispectrum has dependence on the angle between the long and short mode. This information though is integrated over in the halo power spectrum. Moreover, we have suppressed powers of $a_{J}$ as we expect $a_{J}$s to be small from current constraints.

• •

At third order in the density field. Since most of the effect of PNG is at large scales, where $\delta\ll 1$, the series expansion can be safely truncated at third order. As smaller and smaller scales are probed, higher orders need to be introduced.

• •

Higher-derivatives terms. At large enough scales $k\ll 1/R_{*}$, higher-derivatives terms such as $b_{\nabla^{4}\delta}\,\nabla^{4}\delta$ can be neglected, as they typically scale as $R_{*}^{4}k^{4}\delta$ Desjacques et al. (2018).

Up to the order considered, Eq. (5.2.2) is the most generic bias expansion for the halo density field. The unknowns in this expansion are the bias parameters $b_{O}$: similarly to other effective field theory expansions, these need to be measured from data, in this case either via direct observations of the galaxy statistical distribution or by analysing N-body simulations of the Universe. A series of recent efforts have focused on extracting these parameters from high-resolution simulations for Gaussian simulations Fujita et al. (2016); Perko et al. (2016); Abidi and Baldauf (2018); de la Bella et al. (2018), while analyses on universes with primordial non-Gaussianity are still missing.

5.2.3 Models of halo clustering

The results of the previous paragraphs do not rely on a specific modelling of halo formation and clustering. On the assumption that on large enough scales, the formation of structures is well described by the sole action of gravity, it is possible to describe the statistics of tracers with a finite set of operators, with related bias parameters. Given a set of initial conditions, the perturbative treatment of the clustering of halos (or tracers in general) allows to be rather agnostic about the small scale details of the late universe evolution, and predictions have a calculable theoretical error at each order in the expansion. The drawback of this approach is that it has to rely on observations: the parameters of the expansion are free to vary and Eq. (5.2.2) clearly shows that, especially in the case of non-Gaussianity in the initial conditions, there are many of them to fit to observations or N-body simulations.

For this reason, a certain degree of modelling can be combined to the perturbative expansion in order to provide reliable priors on the bias parameters. Chronologically, models of halo and galaxy formation and clustering date back even before the first perturbative methods were applied, with the models of spherical collapse first introduced in the $70$s. Since then, progress has developed into two main directions: the excursion-set approach and the peak model, which we briefly introduced in the previous section §4. The common underlying idea of these models is to build a correspondence between low redshift halos and their progenitors at early times based on a few simple, but well motivated, assumptions on the statistics and physical processes of the initial conditions of structure formation. They are therefore inherently Lagrangian models, as opposed to the Eulerian expansion of Eq. (5.2.2). The advantage of the Lagrangian description, in the context of primordial non-Gaussianity, is that no additional operators have to be introduced with respect to the Gaussian case. There are several efforts in the literature to model the two-point correlation function of tracers including primordial non-Gaussianity, using thresholded regions Matarrese and Verde (2008); Desjacques et al. (2011, 2011), the excursion set approach Desjacques et al. (2011, 2011); D’Aloisio et al. (2013); Adshead et al. (2012); Musso et al. (2012); Desjacques et al. (2013), the peak model Desjacques et al. (2013) (see also Desjacques et al. (2018) for an overview). As in the previous section §4, we choose the example of the excursion set peaks model to derive the non-Gaussian contribution to the halo power spectrum and conclude by making the connection with the generic bias expansion of Eq. (5.2.2).

Excursion Set Peaks

In order to compute the non-Gaussian correction to the power spectrum of peaks, we make use of the effective ESP bias expansion developed in a number of recent papers Desjacques (2013); Desjacques et al. (2013); Biagetti et al. (2014); Lazeyras et al. (2016); Desjacques et al. (2018). The basic idea is to write the peak overdensity field $\delta_{\rm ESP}^{\rm L}$ as an effective perturbative expansion constructed from the rotational invariants of the system, namely the components of the vector $\boldsymbol{\omega}$ defined in Eq. (118)

[TABLE]

where $c_{i}$ are bias parameters and we have not respected any particular order in the expansion, but rather just shown that any combination of the invariants enter the expansion. In order to refine the expansion and make it useful to calculate correlation functions, we have to i) find a way to predict bias parameters and ii) ensure to remove all zero-lag terms such as for example $\langle\delta^{\rm L}_{\rm ESP}({{\bf q}})\delta^{\rm L}_{\rm ESP}({{\bf q}}^{\prime})\rangle\supset c^{2}_{\nu}$ that are generated when using the expansion to compute correlators of $\delta^{\rm L}_{\rm ESP}$. The way to go is to write the expansion in terms of an entire set of orthogonal polynomials $O_{\bf n}[\boldsymbol{\omega}({{\bf q}})]$

[TABLE]

where the indices ${\bf n}=\{i,j,k,\ell,m,n\}$ correspond to the six variables of $\boldsymbol{\omega}$, $\sigma^{[{\bf n}]}$ is an abbreviation for $\sigma^{[{\bf n}]}=\sigma_{0}^{i}\sigma_{1}^{2\ell}\sigma_{2}^{j+2m+3n}\varsigma_{0}^{k}$, where we have defined

[TABLE]

as the variance of the up-crossing variable $\mu=-{\rm d}\delta/{\rm d}R$ and $O^{*}$ is the dual of $O$. The polynomials $O_{\bf n}[\boldsymbol{\omega}({{\bf q}})]$ are found by looking at the probability distributions of the rotational invariants, namely

• •

The trivariate Hermite polynomials $H_{ijk}(\nu,J_{1},\mu)$ are associated to the trivariate normal distribution $\mathcal{N}(\nu,J_{1},\mu)$

• •

The Laguerre polynomial $L_{i}^{(1/2)}\left(\frac{3}{2}\eta^{2}\right)$ is associated to the $\chi^{2}$-distribution with $3$ degrees of freedom of $\eta^{2}$ and $J_{2}$

• •

The polynomials of the form

[TABLE]

where $\mathcal{L}$ are Legendre Polynomials and $x_{3}=J_{3}^{2}/J_{2}^{3}$.

It is straightforward to verify that indeed these polynomials satisfy the orthogonality conditions and that they are a complete basis for the variables of the system Lazeyras et al. (2016). The coefficients $b^{\rm L}_{\bf n}$ are renormalized bias parameters which can be measured through 1-point ensemble averages Lazeyras et al. (2016)

[TABLE]

where the presence of the moments $\sigma_{0}$, $\sigma_{1}$, $\sigma_{2}$ and $\varsigma_{0}$ ensures that the expansion is written in terms of the physical variables $\delta,\eta$ and $\zeta$ rather than the normalized ones (cfr. Eqs. (111) - (113)). Putting everything together, we can write down the effective ESP expansion up to second order in $\delta$

[TABLE]

where we have adopted the notation

[TABLE]

following previous literature. Let us make a few comments on this expansion. As compared to the Eulerian expansion of Eq. (5.2.2), there are two crucial differences: first of all, this effective ESP expansion is done in the initial conditions in Lagrangian space, therefore it does not account for any non-linearity generated by gravitational evolution. The perturbative treatment in the context of Lagrangian space, as well as the connection to Eulerian space can be done using the formalism of Integrated Perturbation Theory Matsubara (2011, 2012, 2014); Yokoyama et al. (2014); Matsubara and Desjacques (2016). This particular formulation of perturbation theory is specifically suited to embed the ESP framework. For the purpose of this review, working with Eq. (5.2.3) is enough to get the intuition of the physics at play. Secondly, Eq. (5.2.3) does not include any tidal shear, that is, scalar combinations of the field $\mathcal{K}_{ij}$ defined in Eq. (165). In the presence of non-Gaussianity, including the tidal shear in the modeling of the ESP is tightly connected with the universality of the mass function and the way one models the scatter in the collapse barrier, so we dedicate the following paragraph to discuss it properly. For the time being, let us assume that the collapse is spherical and that we can neglect the tidal shear in the expansion of Eq. (5.2.3).

We want now to compute the non-Gaussian correction to the power spectrum of peaks as predicted by the ESP model. In order to get a contribution from the primordial bispectrum of $\zeta$ in the power spectrum, we need at least terms proportional to $\langle\delta_{R}^{3}\rangle$. At leading order, this is possible only by combining a first order bias term with a second order one, so that, in Fourier space, we get

[TABLE]

where

[TABLE]

and

[TABLE]

are the first and second order ESP bias parameters in Fourier space . Note that from now on we will neglect all the zero-lag correction terms from the calculation, leaving implicit that we subtract them whenever needed. The expression of Eq. (5.2.3) allows to compute the leading order non-Gaussian contribution for any given primordial bispectrum.

It is instructive to take the squeezed limit $k\rightarrow 0$ of Eq. (5.2.3) in order to make the connection with Eq. (156) of the previous section §5.2.1. For scale-invariant primordial bispectra of the type Eq. (151), we get

[TABLE]

where we defined

[TABLE]

where primes denote derivation with respect to $R$. At first sight, this result seems at odds with the model-independent prediction from the peak-background split ansatz of Eq. (158). Nevertheless, in Gong and Yokoyama (2011) it was proven that the two expressions are equivalent

[TABLE]

for any deterministic barrier for collapse. In the case of local PNG and in the approximation of spherical collapse, one recovers the well known results, $b^{\rm NG}_{\rm ESP}=\delta_{sc}\,b_{100}$, where now the linear bias is a direct prediction of the ESP model. Although this result holds in general for universal mass functions, in this case it is true also for the ESP mass function, Eq. (125), which is clearly not universal, since it depends non-trivially not only on $\nu_{c}$ but also the spectral moments $\sigma_{i}$. As argued in Desjacques et al. (2013); Biagetti and Desjacques (2015), the ESP mass function recovers the result for universal mass functions due to the fact that spectral moments appear only in ratios such as $\nu_{c}/\sigma_{0}$ or $\gamma_{1}=\sigma^{2}_{1}/\sigma_{0}\sigma_{2}$.

Tidal shear, the collapse barrier and primordial non-Gaussianity

The results above can be used to predict the halo power spectrum at late redshift under the (strong) approximation that all halos form around peaks of the early matter density field that overcome a flat spherical collapse barrier. Several analyses Robertson et al. (2009); Grossi et al. (2009); Tinker et al. (2010); Corasaniti and Achitouv (2011); Ludlow et al. (2014); Achitouv and Corasaniti (2012); Elia et al. (2011); Achitouv et al. (2013); Despali et al. (2013); Del Popolo et al. (2013); Borzyszkowski et al. (2014); Castorina et al. (2016); Reischke et al. (2018) of N-body simulations show that the collapse barrier is not really flat, but mass-dependent, and has significant scatter around the mean. The ESP model was therefore extended to allow for a non-spherical barrier with log-normally distributed scatter,

[TABLE]

being $\beta$ the log-normally distributed variable and the parameters of the distribution are fitted to what found in the N-body analyses. This phenomenological barrier was shown to fit rather well simulations with Gaussian initial conditions. In parallel, it was realized that non-Gaussian bias amplitude prediction for universal mass functions, $b^{\rm NG}_{\rm ESP}=\delta_{sc}\,b_{1}$, did not reproduce well some N-body simulation measurements, systematically underestimating the signal (see discussion in the next paragraph §5.3). Attempts at improving the agreement by considering ellipsoidal collapse barriers, so as to have moving barriers, but without scatter, have been proposed Afshordi and Tolley (2008). In the context of the ESP model, in Biagetti and Desjacques (2015) it was argued that the phenomenological barrier of Eq. (181) might improve the fitting to simulations. However, for this barrier, the equivalence of Eq. (180) does not hold. Since the PBS result was derived in full generality, the ESP prediction of the non-Gaussian bias using the barrier Eq. (181) is not consistent. The validity of the PBS prediction was also directly verified in N-body simulations Biagetti et al. (2017), as we will show in the next section.

These considerations bring us to the tidal shear. It has been long known that the local tidal shear has an impact on the collapse threshold Sheth et al. (2001). As remarked above, tidal shear is completely missing from the ESP model as we presented it up to now. As shown in Castorina et al. (2016), a significant part of the scatter in the collapse barrier is related to the local value of the tidal shear at halo formation. In Desjacques et al. (2018), it was demonstrated that if tidal shear is properly included in the ESP model, then the prediction of the non-Gaussian bias amplitude does agree with the PBS result. The crucial difference with respect to the phenomenological barrier of Eq. (181) is that a physically motivated barrier should be determined by the mean density and tidal shear within the collapse region, and therefore it should not explicitly depend on $\sigma_{0}$. The reason why this is particularly important in the non-Gaussian bias prediction is that, according to the PBS derivation, the effect precisely comes from the modulation of the halo mass function to local changes in $\sigma_{0}$. While Desjacques et al. (2018) clarifies this issue and lays down all the necessary ingredients to model tidal shear in the context of the ESP model, a full implementation is still to be done. If tidal shear constitutes most of the observed scatter of the collapse barrier, one should expect that the ESP model successfully predicts the non-Gaussian bias amplitude.

5.3 Numerical approaches

Having discussed in detail various predictions from analytical approaches, we turn now to numerical studies. Calibrating the amplitude of the non-Gaussian bias from local PNG has been one of the main goals in the last decade Dalal et al. (2008); Desjacques et al. (2009); Pillepich et al. (2010); Grossi et al. (2009); Wagner et al. (2010); Hamaus et al. (2011); Scoccimarro et al. (2012); Wagner and Verde (2012); Baldauf et al. (2016); Biagetti et al. (2017). Extensions to other shapes than the local case have also been investigated Shandera et al. (2011); Desjacques et al. (2011); Wagner et al. (2010); Wagner and Verde (2012); Desjacques et al. (2011). The importance of this task is clear from Eq. (144): the non-Gaussian bias amplitude $b_{1}^{\rm NG}$ is degenerate in scale with $f_{\rm NL}$, so any theoretical or numerical uncertainty in its prediction affects the measurement of $f_{\rm NL}$. Generalizations to different types of PNG, such as Eq. (156), show that we expect the degeneracy to remain, or even increase, since the non-Gaussian term becomes less and less distinguishable from the Gaussian one as the scaling in $k$ is weaker for $\Delta>0$. The degeneracy can be alleviated by measuring the $1/k^{\Delta}$ scaling in the power spectrum of tracer populations with different mass and redshift using the multi-tracers technique originally proposed in Seljak (2009), optimal weighting Hamaus et al. (2011); Mueller et al. (2017) and shot-noise suppression Seljak et al. (2009), but it is nonetheless important to have a handle on the uncertainty of $b_{1}^{\rm NG}$. In this sense, the major feature of Lagrangian models presented in the previous section, that is, being able to predict the value of $b_{1}^{\rm NG}$, loses part of its advantage since the theoretical error is unknown and therefore it cannot be reliably used in data analyses. In the next paragraphs, we briefly summarize the most recent efforts in understanding PNG using N-body simulations.

5.3.1 Local-type PNG

For N-body simulations with local-type primordial non-Gaussian initial conditions, measurements of the power spectrum at large scales should show the characteristic $1/k^{2}$ scaling which is absent in Gaussian ones. The measurement of the non-Gaussian power spectrum has several contributions, other than the non-Gaussian bias $\Delta b^{\rm NG}_{1}$ we are after, as first pointed out in Desjacques et al. (2009); Pillepich et al. (2010). At leading order, it reads

[TABLE]

where $P^{\rm NG}$ and $P^{\rm G}$ are the measured power spectra from simulations with non-Gaussian and Gaussian initial conditions respectively and the subscripts ’hm’ and ’mm’ indicate the halo-matter and matter-matter correlations, respectively. Let us go through each term separately:

• •

$b^{\rm G}_{\rm hm}$ is the linear Gaussian bias, which can be measured from simulations by measuring the ratio $P^{\rm G}_{\rm hm}(k)/P^{\rm G}_{\rm mm}(k)$ at large scales,

• •

$\Delta b^{\rm NG}_{I}(f_{\rm NL})$ is the scale independent correction which arises as a consequence of the change of the halo mass function $\bar{n}_{h}$ in the presence of PNG, which we extensively discussed in Sec. §4. This effect grows with increasing halo mass, since the presence of $f_{\rm NL}$ affects the high mass tale of $\bar{n}_{h}$. For the same reason, the correction has opposite sign with respect to $f_{\rm NL}$, since the bias decreases (increases) whenever the halo mass function is enhanced (suppressed) with positive (negative) $f_{\rm NL}$,

• •

$\beta^{\rm NG}_{m}(k)$ is the correction due to the change of the matter power spectrum in the presence of PNG, which can be computed using perturbation theory methods (see Welling et al. (2016) for a recent computation) or directly measured from simulations,

• •

$\Delta b^{\rm NG}_{1}(k)$ is the term we want to measure which is proportional to $1/k^{2}$,

• •

$\mathcal{O}(b^{\rm G}_{2},f_{\rm NL}^{2})$ indicates that the leading order has corrections from two different directions: non-linear biasing and higher-order primordial correlators. Next-to-leading order corrections can be canceled out combining non-Gaussian simulations with opposite sign $f_{\rm NL}$ Biagetti et al. (2017).

Most studies of the non-Gaussian bias from local-type PNG in N-body simulations have used the assumption that the mass function maintains universality even for mildly non-Gaussian initial conditions. Universality simplifies considerably the measurement of the components of $\Delta b^{\rm NG}_{1}(k)$, because the non-Gaussian bias amplitude is proportional to the linear Gaussian bias $b^{\rm NG}_{\rm univ}=\delta_{sc}b_{1}$. Since the earliest measurements Desjacques et al. (2009); Pillepich et al. (2010), it was noticed though that this prescription did not match the measurements with the expected precision, and moreover this discrepancy often changed depending on the halo finding algorithm used. In Pillepich et al. (2010), it was argued that fitted formulas based on the universality assumption deviate from the measured halo mass function from $10$% up to $30$% for low mass halos found using a FoF (Friends-of-Friends) algorithm. Moreover, the fit also degraded as a function of decreasing redshift, deviating $10$% even for more massive halos. In Grossi et al. (2009), the amplitude of the measured scale-dependent bias deviated about $25$% from the prediction of Eq. (148). Later analyses on both SO (Spherical Overdensity) and FoF halos confirmed similar discrepancies Desjacques et al. (2009); Hamaus et al. (2011); Scoccimarro et al. (2012); Wagner and Verde (2012); Baldauf et al. (2016).

To settle this problem, Biagetti et al. (2017) has compared the measurements of the non-Gaussian bias amplitude both with the universal prescription of $b^{\rm NG}_{\rm univ}=\delta_{sc}b_{1}$ and the peak-background split model-independent amplitude $b^{\rm NG}_{\rm PBS}=\partial{\rm ln}\bar{n}_{h}/\partial{\rm ln}\sigma_{8}$, see Fig. 3.

The latter prescription was directly measured from simulations by running a set of Gaussian simulations with varying matter amplitude $\sigma_{8}$ and calculating the numerical derivative of the mass function $\bar{n}_{h}$. The measurements clearly show that the PBS prediction works extremely well, while the universal mass function approximation systematically underestimates the amplitude. In order to clarify this mismatch, they also measured the ratio of the non-Gaussian bias amplitude as predicted by the PBS split ansatz, $b_{\rm NG}^{\rm PBS}$, to the standard universal prediction $\delta_{sc}b^{\rm G}_{1}$ as a function of $b^{\rm G}_{1}$ for a combination of three mass bins and at redshifts $z=0,1$ and $2$, see Fig. 4. The discrepancy between the two prescriptions is evident at the level of $10\%$ or more for all the halo finder algorithms used.

These results call again for a more accurate modeling of the collapse process, which we already discussed in the context of the ESP model in the previous section. Without invoking a specific prescription for the collapse barrier, several attempts in the literature have been made to partially explain these discrepancies. For instance, it is known that FoF halos with linking length $0.2$ tend to have an effective spherical collapse threshold $\delta_{sc}<1.687$ at the high mass end of the halo mass function. This would possibly explain why $b_{\rm NG}^{\rm univ}$ with $\delta_{sc}=1.687$ overestimates $b_{\rm NG}$ at high mass Grossi et al. (2009). However, a change in $\delta_{sc}$ cannot explain the fact that the amplitude $b_{1}^{\rm NG}$ changes sign, as a function of mass, at a different value for the two different prescriptions. The positive outcome of this analysis is that there is indeed a way to calibrate in a model independent way the non-Gaussian bias amplitude, by measuring the response of the halo mass function to a change in the local value of $\sigma_{0}$. Besides running simulations with different $\sigma_{8}$ as done in Biagetti et al. (2017), one could also perform this measurement using separate universe simulations Baldauf et al. (2016).

5.3.2 N-body simulations with generic non-Gaussian initial conditions

While the effects of local-type PNG in the power spectrum of halos have been extensively investigated in N-body simulations, we know that a plethora of other possible interactions may arise from inflation, see Sec. §2.2. In Scoccimarro et al. (2012), a code to generate initial conditions to run simulations with equilateral and orthogonal templates was developed and used to analyse power spectrum and bispectrum measurements. Multi-field inflation was also implemented in N-body simulations, allowing to make precise considerations about the stochasticity generated by the presence of more than one field in the initial conditions Smith and LoVerde (2011); Adhikari et al. (2014). The implementation of a generic primordial bispectrum in the initial conditions of an N-body simulation is not trivial. The main challenge is the fact that physical bispectrum shapes from inflation are not separable, i.e. they cannot be factorized as a product of functions of $k_{1}$, $k_{2}$ and $k_{3}$. Separability is a crucial feature for efficient computational algorithm for simulations not only in LSS Wagner et al. (2010); Fergusson et al. (2012), but also in CMB applications Creminelli (2003); Fergusson et al. (2010). This is the reason why templates are often used in practical applications. In Wagner and Verde (2012), it was argued though that, since the non-Gaussian contribution on the power spectrum is peaked in the squeezed limit, one should check that the templates reproduce correctly the physical shape in this limit, even in the case in which the physical shape does not peak in that limit. For instance, the orthogonal template Senatore et al. (2010) generates a non-Gaussian bias which scales as $1/k$ in the squeezed limit, while typical orthogonal physical shapes from the effective field theory of inflation Bartolo et al. (2010) generate a scale independent correction.

As a consequence of these considerations, Fergusson et al. (2012); Wagner and Verde (2012) have introduced methods to implement a generic inflationary bispectrum in the initial conditions of N-body simulations. For instance, Wagner and Verde (2012) has introduced the following ansatz

[TABLE]

which, however, is not separable and therefore is computationally expensive. Note that the integration is bounded from both above and below, since the simulation has a finite box size and resolution. The method of Fergusson et al. (2012) is separable, but the power spectrum of $\Phi$ receives spurious corrections at large scales Wagner and Verde (2012).

5.4 Final remarks of this section

The last decade has seen an intense theoretical and numerical work in analysing the effect of primordial non-Gaussianity in the statistics of biased tracers. The appearance of a scale-dependent feature at large scales, sourced by certain types of interactions during inflation, has motivated observational efforts and even pushed the funding of a tailored experiment, SPHEREx Doré et al. (2014), which is scheduled to be launched in 2023. We have argued the importance of two main aspects of this search: on the inflationary side, a selection of theoretically motivated models with observable imprints; on the structure formation modeling side, the need for combining effective approaches with model-dependent ansatz in order to be able to accurately handle the astrophysical uncertainties which hamper our ability to extract the primordial signature. This procedure becomes more and more important when refining the prediction to include halo occupation distributions algorithms for the modeling of galaxy statistics, red-shift space distortions and survey-related uncertainties and when considering higher-order statistics such as the three-point correlation function of halos and galaxies. Another important ingredient, as for all LSS studies, is the use of N-body simulations. Not only the theoretical modeling has to be meticulously tested against numerical results, but also an accurate matching of the N-body outputs against real data needs to be performed. In the case of primordial non-Gaussianity, these two tasks are complicated by the fact that each inflationary prediction effectively realizes a different cosmology and therefore requires a different set of N-body simulations to run and test. For this reason, numerical efforts have concentrated mainly on simple templates for known non-Gaussianities, such as the local type. More work is surely needed in order to properly model the large amount of possible interactions, most of which we have summarized in Sec. §2.2. Indeed, as remarked in the previous section, the risk of using simplified templates is to lose constraining power on more complicated, yet well motivated, inflationary signatures.

6 Observational prospects

While the focus of this review has been restricted to the theoretical study of the imprint of inflationary interactions on the statistics of dark matter halos, the discussion can, and has to, be extended towards several directions. First of all, the imprint of primordial interactions which we reviewed in Secs. §4 and 5 has effects on a number of observables other than halos and galaxies number counts, which we summarize in the following list

• •

Galaxy shapes contain this imprint in the shear and convergence field probed through weak-lensing Marian et al. (2011); Maturi et al. (2011); Jeong et al. (2011); Hilbert et al. (2012) and galaxy intrinsic alignments Chisari and Dvorkin (2013); Chisari et al. (2016); Kogai et al. (2018).

• •

The Sunyaev-Zeldovich (SZ) effect Sunyaev and Zeldovich (1970) can be exploited to observe dense clusters of galaxies which constrain the high tail of the density distribution function of galaxies which is sensitive to primordial non-Gaussianity Mak and Pierpaoli (2012). The thermal SZ power spectrum Hill and Pajer (2013) and kinetic SZ tomography Munchmeyer et al. (2018) have also been used to put constraints on local-type non-Gaussianity.

• •

The pairwise velocity distribution of galaxies is an additional probe to galaxy density statistics, as primordial non-Gaussianities induce a non-zero skewness and higher-order momenta in the distribution Catelan and Scherrer (1995); Schmidt (2010); Lam et al. (2011).

• •

The scale-dependent bias outlined in Sec. §5 can be constrained using extrema counts as predictedy by peak theory Nusser et al. (2018)

• •

The two-point statistics of voids, when combined with halos, give an order $O(1)$ improvement in local-type non-Gaussianity constraints Chan et al. (2018)

• •

The covariance of galaxy number counts also helps in combination with number counts and variance Cunha et al. (2010)

• •

Higher-redshift probes have been shown to be promising in constraining non-Gaussianity in the future: Lyman-alpha forest Seljak (2012); Chongchitnan (2014), 21-cm power spectra Joudaki et al. (2011); Chongchitnan and Silk (2012); Chongchitnan (2013); Camera et al. (2013); Sekiguchi et al. (2019); Witzemann et al. (2019), CO and CII lines Moradinezhad Dizgah et al. (2019); Moradinezhad Dizgah and Keating (2019), cosmic reionization D’Aloisio et al. (2013); Lidz et al. (2013); Mao et al. (2013) and cross-correlations with CMB measurements Raccanelli et al. (2015); Schmittfull and Seljak (2018).

• •

The response of the small-scale power spectrum to a squeezed bispectrum is particularly effective for the imprint of primordial non-Gaussianity Chiang et al. (2017); Yuan et al. (2017); de Putter (2018).

This review has also neglected any general relativistic effects. These become important when probing the largest scales of the galaxy power spectrum, where the imprint of sizable primordial bispectra in the squeezed configuration is strongest. Several analyses have computed these corrections Wands and Slosar (2009); Bruni et al. (2012); Baldauf et al. (2011); Jeong et al. (2012); Yoo et al. (2011); Lopez-Honorez et al. (2012); Yoo et al. (2012); Villa et al. (2014); Di Dio et al. (2017); Castiblanco et al. (2018) and the importance of these relativistic corrections has been quantified for forthcoming surveys Maartens et al. (2013); Raccanelli et al. (2014); Lorenz et al. (2018).

So far, we have been concerned only with one- and two-point statistics of LSS observables. Higher-order tracers statistics are able to directly trace primordial interactions. For instance, the bispectrum of galaxies is sourced at tree-level by all the non-Gaussian shapes in Sec. §2.2. Differently than for the power spectrum, where only the squeezed configurations contribute to the signature, the full shape of non-Gaussianity can be probed in the bispectrum. Theoretical and numerical work has been also making progress in this direction in the last decade Baldauf et al. (2011); Pollack et al. (2012); Sefusatti et al. (2012); Yokoyama et al. (2014); Tasinato et al. (2014); Dizgah et al. (2015); Hashimoto et al. (2016); Tellarini et al. (2016); Hashimoto et al. (2016); Yamauchi et al. (2017); Di Dio et al. (2017); Chiang et al. (2017); An et al. (2018); Moradinezhad Dizgah et al. (2018); de Putter (2018); Moradinezhad Dizgah et al. (2018); Eggemeier et al. (2018). The main challenges to overcome with the galaxy bispectrum is the fact that gravitational non-linearities developed at late times dominate over primordial ones. The accurate modeling of these non-linearities requires introducing a large number of new parameters, as compared to the case of the power spectrum. Nevertheless, combining power spectrum and bispectrum information can significantly improve constraints, see for instance the forecast for the SPHEREx mission Doré et al. (2014).

Besides improving the modeling and combining different statistics, several optimisation techniques have been proposed. One of the major problems to overcome, common to both CMB and LSS, is sampling variance. The multi-tracer technique Seljak (2009); Slosar (2009); McDonald and Seljak (2009); Gil-Marin et al. (2010) is based on the assumption that on large scales halos are biased, but not stochastic, tracers of the dark matter density field. Recently it has been argued that similar improvements on the sampling variance can be achieved by selecting tracers with no bias with respect to the dark matter density field Castorina et al. (2018). One can therefore eliminate the cosmic variance error by correlating a highly biased population of galaxies against an unbiased one. Another limitation is caused by the discrete nature of the galaxy distribution statistics. This is taken into account usually by adding a Poisson shot noise term to the galaxy power spectrum. This term is particularly relevant for populations of galaxies with a large mass, due to their low density. These are also the galaxies which are more sensitive to the imprint of primordial non-Gaussianity, as showed in Secs. §4 and 5. A mitigation of this problem is provided by optimally weighting populations of galaxies Seljak et al. (2009); Hamaus et al. (2010, 2011).

The current best constraints on primordial non-Gaussianities as set by CMB observations by the Planck satellite. Constraints are given in terms of the local, equilateral and orthogonal templates to be $f^{\rm loc}_{\rm NL}=-0.9\pm 5.1$, $f^{\rm equi}_{\rm NL}=-26\pm 47$ and $f^{\rm orth}_{\rm NL}=-38\pm 24$ at $65\%$ confidence level, respectively. LSS searches have been also putting constraints since more than a decade Slosar et al. (2008); Tseliakhovich et al. (2010); Xia et al. (2010, 2011); Benson et al. (2013); Ross et al. (2013); Mana et al. (2013); Giannantonio et al. (2014); Shandera et al. (2013); Ma et al. (2013); Karagiannis et al. (2014); Agarwal et al. (2014); Ho et al. (2015); Giannantonio and Percival (2014); Leistedt et al. (2014); Castorina et al. (2019), but they do not give competitive constraints with these figures for any of the observables listed above. The latest constraint was put by the eBOSS collaboration Castorina et al. (2019) and gives $-51<f_{\rm NL}^{\rm loc}<21$ at $95\%$ confidence level and the full data analysis of the experiment including quasars might reach the CMB sensitivity. A strong effort in forecasting the possible imprint of primordial non-Gaussianity in LSS has been pushing the limits to order $O(1)$ for local-type non-Gaussianities141414See Table 12 of Desjacques et al. (2018) for a summary table of future prospects on $f_{\rm NL}^{\rm loc}$. and $O(10)$ for equilateral and orthogonal ones Fedeli et al. (2011); Namikawa et al. (2011); Raccanelli et al. (2012); Giannantonio et al. (2012); Pillepich et al. (2012); Takeuchi et al. (2012); Sefusatti et al. (2012); Becker et al. (2012); Biagetti et al. (2013a, b); Yamauchi et al. (2014); Ferraro and Smith (2015); Byun and Bean (2015); Camera et al. (2015); Xu et al. (2015); de Putter and Dorè (2017); Doré et al. (2014); Amendola et al. (2017); Sartoris et al. (2016); Raccanelli et al. (2016); Alonso et al. (2015); Gariazzo et al. (2015); Alonso and Ferreira (2015); Fonseca et al. (2015); Raccanelli et al. (2017); Yamauchi and Takahashi (2016); Hashimoto et al. (2016); Amendola et al. (2018); Kovetz et al. (2017); Chisari et al. (2016); Adhikari et al. (2016); Fonseca et al. (2017); Yamauchi et al. (2017); de Putter et al. (2017); Gleyzes et al. (2017); Li and Ma (2017); Moradinezhad Dizgah and Dvorkin (2017); Moradinezhad Dizgah et al. (2018); Bellomo et al. (2018); Moradinezhad Dizgah and Keating (2019); Bacon et al. (2018); Chan et al. (2018). Several future surveys such as SPHEREx Doré et al. (2014), Euclid Amendola et al. (2018), MeerKAT+DES Fonseca et al. (2017) and SKA Bacon et al. (2018) might achieve these limits.

\funding

This research was funded by the Nederlands Organisation for Scientific Research (NWO) under the VENI grant 016.Veni.192.210.

Acknowledgements.

The author thanks Daniel Baumann, Guilherme Pimentel and Toni Riotto for useful discussions and Vincent Desjacques, Azadeh Moradinezhad Dizgah, Guilherme Pimentel and Toni Riotto for comments on a draft. \conflictsofinterestThe author declares no conflict of interest. \abbreviationsThe following abbreviations are used in this manuscript:

[TABLE]

\appendixtitlesno \appendixsectionsmultiple

Appendix A ESP model

A.1 The curvature function of density peaks

Without entering into the details of the full calculation (see Bardeen et al. (1986)), the integral Eq. (116) can be simplified into the only integration over $J_{1}$

[TABLE]

where we define

[TABLE]

and

[TABLE]

and we used Bayes’ theorem to write $\mathcal{N}(\nu,J_{1})=\mathcal{N}(\nu)\,\mathcal{N}(J_{1}|\nu)$ and consequently written the conditional gaussian distribution $P_{\rm G}(x-\mu;s)$ with shifted mean $\mu=\gamma_{1}\nu$ and variance $s=1-\gamma_{1}^{2}$. The integral can be computed analytically and gives the final result

[TABLE]

where we define the integrals

[TABLE]

and the [math]-th order is the curvature function of density peaks Bardeen et al. (1986)

[TABLE]

Note that Desjacques et al. (2010) introduced the extra variable $\alpha$ in order to get a closed form expression for their 2-point peak correlation, while Desjacques (2013) showed that $\alpha\neq 1$ can be interpreted as a long-wavelength perturbation in $J_{2}(\mathrm{\bf x})$.

A.2 The non-Gaussian ESP mass function

Performing the transformations Eqs. (128) for all the terms and defining

[TABLE]

we finally get

[TABLE]

\reftitle

References

\externalbibliography

yes