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

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.


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 [1][2][3][4], 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 [5]. 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 [6]. These features are strong hints that an inflationary mechanism indeed took place 1 . 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 (∼ 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 [9]. 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 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 [1][2][3]. There are countless good introductions to inflation and the production of primordial perturbations, including textbook material [31][32][33][34][35][36] and reviews [11,[37][38][39]. 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.

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 µν and slowly rolling down its potential. Its vacuum energy drives the accelerated expansion of the Universe. The corresponding action reads where M P is the reduced Planck mass, the first term is the Einstein-Hilbert action and V(φ) is the slow-roll potential. The background evolution is studied by assuming a flat Friedmann-Lemaître-Robertson-Walker (FLRW) metric where a(t) is the scale factor, for which the equations of motion are Here the prime on the potential V indicates derivation with respect to the background field φ 0 , while the dot denotes the time derivative with respect to cosmic time t. The expansion rate H =ȧ/a, known as the Hubble parameter, is determined by the first equation, while the evolution of the background field φ 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,Ḣ ≈ 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 N(φ) ≡ a e a dln a = t e t H dt (6) 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 [6] 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 ≡ −Ḣ H 2 1, η ≡˙ H < 1 (7) 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 Within the approximation of perfectly constant rate of expansion, the scale factor grows exponentially in time a ∝ e Ht and the space time is called De-Sitter (DS) space-time.

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 φ(x, t) = φ 0 (t) + δφ(x, t).
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 where we have defined conformal time τ as dt ≡ adτ. Note that τ ∈ (−∞, 0). The equations of motion in Fourier space for this action read δφ (k, τ) + 2Hδφ (k, τ) + k 2 δφ(k, τ) = 0, (11) where H is the Hubble parameter in conformal time. The generic solution to this differential equation is 3 To quantize the perturbations, we promote the field δφ and its conjugate momentum δπ ≡ ∂L/∂δφ to operators and write canonical commutation relations [δφ(x, τ), δπ(y, τ)] = iδ(x − y), (13) [δφ(x, τ), δφ(y, τ)] = 0, [δπ(x, τ), δπ(y, τ)] = 0.
It is useful to decompose the field defining time-independent creation and annihilation operators in Fourier space δφ(k, τ) = u(k, τ)a k + u * (−k, τ)a † −k , (15) δπ(k, τ) = a 2 u (k, τ)a k + a 2 u * (−k, τ)a † −k , (16) with commutation relations [a p , a † −q ] = (2π) 3 δ (3) (p + q), (17) [a p , which follow from the fact that a 2 u(k, τ)u * (k, τ) − 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 |c 1 | 2 − |c 2 | 2 = 1 (19) and the second from choosing the vacuum state. To do this, we notice that on sub-horizon scales, that is for kτ 1, the system changes on time scales much shorter than the typical time scale of 3 Here we used the approximation of perfectly DS spacetime, which implies τ = −1/Ha. Corrections are proportional to slow-roll parameters, which indeed we are neglecting at this stage.
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 [40]. With this choice of initial conditions, the solution reads While well within the horizon, kτ 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 δφ 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 where the approximate equality indicates again corrections suppressed by slow-roll parameters. An analysis extended to the coupled inflaton-metric fluctuations would show that ζ 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, ζ ≈ −5∆T/T. Correlation functions of this quantity therefore characterize perturbations produced during inflation. The power spectrum P ζ is defined as 0|ζ(k 1 , 0)ζ(k 2 , 0)|0 = (2π) 3 δ (3) (k 1 + k 2 )P ζ (k 1 ). (22) Combining Eqs. (21), (15) and (20) we get where it is customary to define a dimensionless power spectrum as P ζ (k) ≡ k 3 /2π 2 P ζ (k). Here * indicates that the quantities are computed at horizon crossing, kτ * = 1, and we reinstated slow-roll suppressed parameters in the last factor, giving the power spectrum a spectral index The amplitude A s = H 2 * /8π 2 M 2 P * and spectral index n s of the power spectrum have been measured by Planck [9] to be A s = (3.044 ± 0.014) × 10 −10 , n s = 0.9649 ± 0.0042. (26) 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 , 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.

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 ζ 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 ζ 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.

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 [41][42][43][44][45]. A generic form for an in-in correlator is where Q(τ) represents the operator of the product of n-curvature perturbation ζs so thatQ = ζ k 1 ζ k 2 ...ζ k n and Ω is the in state, the vacuum of the interacting theory. Note that, more in general, there could be combinations of both ζ and γ, the tensor perturbation, inQ. In this review we only focus on ζ, but tensor perturbations are equally important in constraining the inflationary scenario [46][47][48][49][50].
The time τ 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(τ) back to initial time τ 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 int . Eq. (27) is therefore evaluated formally as where T andT are time and reversed time-ordering symbols andQ I andĤ I int are evaluated using interaction picture operators. The i prescription allows to project the interacting vacuum state, Ω to the free one. The above expression can be expanded as a power series in H int . Interactions are then organized as usual using Feynman diagrams. As an example, let us take the three-point correlator of ζ expanding Eq. (28) to first order where H int = −L 3 + O(ζ 4 ) and L 3 is the perturbed Lagrangian up to cubic order in ζ and we have taken the superhorizon limit kτ → 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 It is customary to decompose B ζ as where we assumed statistical homogeneity and isotropy and we evaluate the dimensionless power spectrum, P ζ 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 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 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 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 NL and defined as where, as we can see, f NL can generically depend on K.

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 ζ. 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 ζ, 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.

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 and η, 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 [47] and we redirect directly to it for details on the computation. The dimensionless bispectrum that results from the computation reads 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 GF NL = 55/9 + 15/9 η. 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 Extensions and generalizations of this statement have been made in several subsequent efforts [47,49,[51][52][53][54]. 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 ζ 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 long /k short lim 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 φ and its derivative [55] where X ≡ −1/2g µν ∂ µ φ∂ ν φ. In the context of an effective field theory of inflation, the form of P(X, φ) is constrained by the following considerations: • Non-derivative operators, such as φ n /Λ n−4 , being Λ the largest energy scale where the effective description holds, contribute directly to the inflaton potential and are therefore strongly constrained by the background [56]. • Derivative operators of the form X n /Λ 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 /Λ 4 , gives Non-Gaussianities of order unity are therefore generated whenφ 0 ≈ Λ 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φ 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 [57,58], for which where f (φ) is the (squared) warp factor of the AdS-like throat related to these models 4 . The dimensionless bispectrum for a generic P(X, φ) model can be written in a rather model-independent way [60] S HD (k 1 , k 2 , where the subscript HD stands for "Higher-Derivative" terms and and we have defined the following quantities [55] 4 One should bear in mind that there is a degree of fine-tuning needed also for this model [59]. 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 As it is clear from Eq. 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 [61,62]. 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 [63][64][65][66][67] and they can also be studied in the context of the effective field theory of inflation [68]. Such features typically show up primarily in the power spectrum and constrains can be put with observations both of the CMB [9] and LSS [14], 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 µ ∼ ∆V/V and width σ, it can be easily shown on general grounds that observations of the power spectrum from the CMB constrain the ratio µ/ 1 [9]. The slow-roll parameter η on the other hand can change a lot if the change in occurs within a short time. If we take For most models with a sharp feature, the calculation of the bispectrum requires to use numerical solutions. Here we quote an approximated shape [61] S sharp (k 1 , 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.
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 [69,70]. In these models, the inflaton potential is characterized by an oscillatory term added to the usual slow-roll one where here Λ 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 ω > H, we might expect a resonance between couplings and modes which sources non-Gaussianities [71,72]. 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 being the slow roll parameter and f M P is predicted in string theory constructions with sub-Planckian decay constants [73,74]. The corresponding dimensionless bispectrum for this type of models reads where α = √ 2 / f M 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τ 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 [71] 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 ζ has therefore at least one contribution from them. Effectively, this translates into sending one of the three momenta from k → −k.

Solid inflation
Solids can be described in the context of field theory [75] by introducing three scalar fields φ I whose background values are identified with spatial coordinates φ I (x) = x I (52) and they are time-independent. Using this framework, [76] 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 [76].
Most interestingly for this review, the three-point function of scalar perturbations drastically violates [77] the standard single-field consistency relation of Eq. (34). The dimensionless bispectrum is calculated to be [76] where F and F Y are free parameters of the solid Lagrangian, 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, τ c is the conformal time at which the longest modes of observational relevance today exit the horizon, while τ 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 [76], respectively. It is important to notice that Eq. (53) is computed for not-too-squeezed momenta, that is, for k long /k short > √ . The squeezed limit was investigated in detail in [77], at leading order in slow-roll expansion it can be written in a much more compact form as which should be compared to the consistency relation of Eq. (34). Here θ is the angle between q and 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.

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 ζ.
Massive particles might be spontaneously created in the expanding space-time during inflation through non-perturbative effects [78][79][80]. Their production is particularly interesting as a probe to ultraviolet completions of inflation motivated by string theory constructions [12]. 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 H, decay exponentially fast after horizon crossing and their effect can be integrated out from the dynamics of ζ. 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 ζ 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 [81,82]. While the effect of massive scalar fields have been thoroughly investigated in the context of quasi-single field inflation [12,54,[83][84][85][86][87][88][89][90][91][92], higher-spin particles studies arose more recently [87,88,91,[93][94][95][96][97].

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 [98,99] for particles with spin s ≥ 1 principal series complementary series discrete series for t = 0, 1, 2, ..., s − 1. Similarly, scalar particles with mass m ≤ 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 ≥ 1 are required to have a minimal mass, m 2 ≥ s(s − 1)H 2 , unless they belong to the discrete series 5 . For these values, the system acquires an additional gauge invariance and the corresponding fields are called partially massless fields [101]. Partially massless fields produce sizeable tensor-scalar-scalar bispectra and the scalar trispectra [91,93,95]. 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 where ∆ = 3/2 ± iµ and 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 where η is conformal time and ∆ ± s = 3/2 ± iµ s is the conformal weight of the field and 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 µ s , for which the wavefunction oscillates logarithmically in conformal time, and imaginary µ 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 [88]. This extends earlier works focusing 5 Light particles with spin are allowed in inflationary set ups which breake DS isometries, as shown in [100]  only on scalar fields. At tree-level, [88] 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 ζζζ .
The corresponding dimensionless bispectrum for these three diagrams can be written in the following compact form where (κ) = (a), (b) and (c) represents the three diagrams, α s are functions of the angles between momenta (typically Legendre polynomials) and I (κ) s are complicated integrals whose expressions can be found in [88] 6 and c π 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: S 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 [87].
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 O( f NL ) 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. (33) and it does not contain information on the spin and mass of the particle at leading order, The non-analytic part for scalar particles has been derived in the context of quasi-single field inflation [83]. It has to be distinguished into two cases: for massive particles belonging to the principal series, µ ≥ 0 it reads while for particles belonging to the complementary series for which µ becomes imaginary, the scaling changes to 6 Note that there is a conversion factor of (k 1 k 2 k 3 ) 2 between I (κ) s and those defined in [88]. For instance, I where ν = −iµ is real. The case of higher-spin particles proceeds analogously: massive particles belonging to the principal series, µ s ≥ 0, have a squeezed limit of the form for even spins, where the L s is the Legendre polynomial of order s and φ s is a phase that is uniquely fixed in terms of µ s and c π [88]. 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 µ s becomes imaginary, the scaling changes to where ν s = −iµ 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 S ∝ (k 1 /k 3 ) 0 .

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 [39] for recent lectures on the topic). Observations of the CMB constrain the amount of isocurvature perturbations to be very small [9]. Nevertheless, there are several mechanisms with which isocurvature perturbations can be converted to curvature ones after inflation and evade CMB constraints [102][103][104][105][106][107][108][109]. 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 δN-formalism [110][111][112], 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 [102,[113][114][115] 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 x, The corresponding dimensionless bispectrum in Fourier space is of the form 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 δN-formalism
In the single-field model we introduced in §1, the comoving curvature perturbation ζ 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ζ . Each comoving superhorizon patch in the universe will have its own frozen value of ζ 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, δN, at different positions. This interpretation is dubbed δN-formalism.The generalization of this picture to the case of multi-field inflation is well summarized in [61], here we broadly follow its steps.
Let us consider a set of scalars φ i = φ 0i + δφ 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 (φ 0i ) as the number of e-folds between the initial and final slices for the unperturbed field φ 0i , while the one for the perturbed fields will be N(φ 0i + δφ i (x)). The curvature perturbation ζ can be expressed as a Taylor expansion of the variation of N around the initial values φ 0i as where the subscripts on N denote partial derivatives with respect to φ i and indices are summed using the Einstein convention, with i = 1, ..., N, being N the number of fields. Correlation functions of ζ are then related to correlation functions of δφ 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 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 and it can be easily shown that the corresponding dimensionless bispectrum S is the local one of Eq. (66) and 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 δN-formalism, such as the curvaton model, the modulated reheating and preheating scenarios (see [116] for a review).

Final remarks of this section
We have identified three main features that we can extract from the bispectrum of the curvature perturbation ζ: shape, running and amplitude. Among these three, the amplitude, parametrized via a single parameter f 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.

Figure 2.
A schematic representation of current and future limits on the primordial non-Gaussianity amplitude parametrized by f NL . Credits to [88].
The limit of f NL ∼ 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φ 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 ≈ k 2 ≈ k 3 ). Furthermore, a positive detection of an order unity f NL for squeezed bispectra (i.e. with k 1 k 2 ≈ k 3 ) would rule out single-field models of inflation in one shot. The ultimate goal is to measure f NL down to the level of slow-roll parameters, f NL 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 [6] f loc 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 [117,118]. 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 [119][120][121]. 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 [87] 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.

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 functions 7 .
The initial conditions for structure formation are set by connecting the primordial curvature perturbation ζ to the Newtonian potential Φ via a transfer function, T(k). At the linear level and using the Poisson equation, we can write the useful relation which defines the linear (L) matter overdensity. Here D(z) is the linear growth factor normalised to unity at present day and Ω m,0 and H are the matter density and hubble parameter today, respectively. The linear relation is reliable as long as the matter overdensity δ is much smaller than unity. In this regime, one can solve the (Newtonian) dynamical equations perturbatively (see [127] for a review and [128][129][130][131][132][133][134][135][136][137][138][139][140][141][142] 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 0.2 h/Mpc at redshift 0, 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 δ 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 and its variance To avoid clutter, we will condense the transfer function and the smoothing window into M R (k, z) ≡ M(k, z)W R (k) and drop the redshift dependence unless needed.

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 x, i.e. the value of the smoothed density field δ R at each 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 7 Weak lensing [122] and temperature anisotropies in the 21cm background from the pre-reionization epoch [123][124][125] can be used to trace DM directly, although the latter might also be biased [126] where δ N R c is the connected n-point moment. Note that, in the linear regime, σ R ∝ D(z) and therefore S N ∝ 1/D N−2 , which means that for N ≥ 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 where 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. (33) and gives σ R S 3,GF = A(R) + B(R)η, being A and B mild functions of the scale R and with amplitude A, B ∝ O(10 −4 ) for R ≈ 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 δ R , or more commonly the so-called peak height ν = δ R /σ R , to get where H N are Hermite polynomials. Note that, as a consequence of how the reduced cumulants scale with the linear growth factor, the combinations σS 3 and σ 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 [143][144][145][146][147].
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 ≈ 34/7 at lowest order [127], 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.

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, It is immediately clear from Eq. (81) that a bispectrum, or higher-order correlations, of ζ 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 [138].
In standard Eulerian perturbation theory (EPT) (see [127] for a review), one follows the evolution in time of the matter density contrast δ at position x. In the weakly non-linear regime, the Fourier mode of the density field reads where we simply denote δ as the unsmoothed, non-linear matter field and F 2 is the PT kernel and µ is the cosine between q 1 and q 2 . For the present analysis, it is not essential to smooth the matter field, so we work only with the unsmoothed δ. A contribution from a non-zero primordial bispectrum therefore arises already on the matter power spectrum at the 1-loop order and reads where notice that here we have used the statistical homogeneity and isotropy to write the momentum dependence of B ζ . The full 1-loop power spectrum is defined as where P 22 and P 13 are the 1-loop contributions to the Gaussian case [127]. 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 ∼ 0.1 h/Mpc are approached. Moreover, in [138] 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, (86) Extensions to higher loops have been computed in the context of perturbation theory (see [138,148] for a treatment in the effective field theory of LSS). It is customary to define a reduced bispectrum 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 [149][150][151]. 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 [152]. 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 [148].

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 [153][154][155], 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 1980s [156][157][158][159][160][161][162]. Almost as early, numerical methods have been used to test predictions and provide fits to data [163][164][165][166][167][168]. 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.

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 [169] to local-type non-Gaussianities in the initial conditions [170][171][172][173][174][175]. This extension was also studied for higher order primordial non-Gaussianities [172,176,177] and a range of other bispectrum shapes [178] ; the excursion set approach [179][180][181][182][183], which was introduced to solve problems suffered by the PS model, was thoroughtly also investigated with non-Gaussianity in the initial conditions [173][174][175][184][185][186][187][188][189][190][191][192][193][194][195]. Lastly, the peak model [196], which was recently combined with the excursion set approach in the Excursion Set Peaks (ESP) model [197,198], has been applied to non-Gaussian peak statistics [199][200][201][202][203][204]. Recently, methods based on spherically averaging cosmic densities have been shown to successfully disentangle primordial non-Gaussianities with late time ones [204].

The Press-Schechter mass function
The Press-Schechter (PS) model [169] 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 [205]. Let us first take the case of Gaussian initial conditions. The PDF for the density field δ R smoothed on a scale R ∝ M 1/3 , where M is the halo mass, is a simple Gaussian with zero mean and variance σ 2 R . In the PS approach the halo mass function reads whereρ is the mean matter density, F >δ sc is the level excursion probability that the density field smoothed on a scale M has overcome the threshold for spherical collapse, δ sc ≈ 1.687 8 , and ν(M, z) = δ sc /σ R (z) is the peak height 9 . 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 [179][180][181][182][183]. 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 ν, i.e. more massive halos are more rare. Despite the fact that the linear theory value of the threshold for spherical collapse is used, δ sc , the calculation is fully non-linear [205]. A feature of Eq. (88), and common to all universal mass functions 10 , is that it is entirely specified by a function of ν only, where f (ν) is called multiplicity function and we have dropped a subscript ν 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 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 [161] by performing a saddle-point approximation on an Edgeworth expansion on F >δ sc and computing the ratio of the non-Gaussian-to-Gaussian multiplicity functions where ν * = δ sc √ 1 − S 3 δ 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, [170] expanded directly the PDF, Eq. (80), to calculate F >δ sc and then the ratio where here we only considered the skewness. Combining the two predictions to better match numerical simulations, [172] got where ν * from Eq. (91) was expanded to first order in f 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 [208][209][210].
These formulae can be extended in case of primordial non-Gaussianity from trispectra of ζ. In this case, the kurtosis S 4 needs to be added [172,176,177]. 8 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. 9 The smoothing scale R is directly related to the mass of the halo M, we will use them interchangeably. 10 As 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 [206,207].

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 δ 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) = δ 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, With a similar meaning to the multiplicity function, the first crossing fraction f (s)ds is defined as where it is customary in this framework to use the variance s = σ 2 0 (R) as the reference variable for the random walk. As we said, the excursion set ansatz consists in requiring that δ < 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 δ < B(S − ∆s) for ∆s → 0. We can then Taylor expand in δ and B and get the following condition: with dδ/ds ≥ dB/ds at s = S, which is a condition on the 'velocity' of δ being greater than the tilt of the barrier [197,211,212]. This implies that now we have to deal with a bivariate distribution on δ and its velocity δ ≡ dδ/ds, P(δ, δ ), but with the simplification of having only one condition to impose, rather than an infinite number of them. The first crossing distribution reads If we now take the limit ∆s → 0 and change the integration variable δ → v = δ − B we get the general formula As an example, we can compute Eq. (99) for Gaussian initial conditions and a constant barrier, B(s) = δ sc . In this case, it is convenient to work with the following change of variables and we have defined We obtain therefore the following first crossing distribution, as a function of ν, where the subscript G indicates the Gaussian distribution and we have used the conditional probability There are several efforts that extend the excursion set to non-Gaussian initial conditions [173][174][175][184][185][186][187][188][189][191][192][193][194][195], 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 [187,188,213], which build up on results obtained for the Gaussian case [214][215][216] 1. They formulated the excursion set using top-hat filtering in real space, which is the preferred filter to compare with data and simulations [187]. 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 [217,218] for another approach). 2. They gave a physical explanation of the correction to the spherical threshold for collapse δ sc → √ qδ sc which is commonly used to improve fits to simulations (see [182,[219][220][221] for the earliest applications ) and showed how it is expected that the factor depends on the halo finder used [215]. 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 [187]. 4. They extended the formalism to a generic moving barrier B(σ) in [188].
These findings are summarized in the following prediction for the non-Gaussian multiplicity function where q = 1/(1 + D B ) being D B the diffusion coefficient of the barrier, which advocates 2. of the previous list, is the correction due to the use of top-hat filtering (addressing 1.), Γ(0, x) is the incomplete Gamma function and is the correction due to non-Gaussianity (addressing 3.), being U and 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 [188], where they also used an improved method based on the saddle-point approximation proposed by [191]. Another series of papers analyzes drifting diffusing barriers in the context of the excursion set approach [193], including also an extension to a class of primordial trispectra [195]. 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.

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 [196]. 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 ∝ (M/ρ) 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 [153][154][155]222]. Several authors have developed the technical tools to deal with local density maxima in the presence of non-Gaussian initial conditions [159,196,199,200,203,204]. 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 [197]. A full review and detailed description of the peak model and its extension to the ESP is nicely presented in [207]. 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 δ(x) as a point process, where δ D is the 3-dimensional delta function. The continuous field δ(x) can be Taylor expanded around a peak position x pk , where we imposed the peak condition on the first derivative, ∇ i δ(x pk ) = 0. In a similar manner, we can Taylor expand ∇ i δ(x) around the peak as Using Eq. (107) and Eq. (108), we can express the delta functions of Eq. (106) as provided that the hessian ∇ i ∇ j δ is definite negative, to ensure to be looking at maxima, and non singular in x 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 [223].
As a result, the number density of peaks of height ν at position x can be expressed in terms of the smoothed field δ s and its derivatives as [196] n pk (ν, R s , where R * = √ 3σ 1s /σ 2s is the characteristic radius of the peak and we have conveniently defined the normalized peak density and its derivatives as and we will adopt the following notation for the variance of the smoothed density field, (linearly extrapolated to present-day) and its derivatives, We have imposed explicitly that the stress tensor −ζ ij is definite negative through the positivity of λ 3 , its lowest eigenvalue.
The number density n 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 pk , ρ (N) 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 densitȳ n pk = n pk (ν, R s , x) = dν d 3 η d 6 ζ n pk (ν, R s , x)P 1 (ν, η, ζ A ) (116) being A = 1, ..., 6 since ζ is symmetric and the joint probability distribution P 1 is given by where M is the corresponding covariance matrix. Owing to rotational invariance, we have regrouped the set of variables into the following vector where J 1 = −tr(ζ ij ) is the peak curvature, J 2 = 3 2 tr(ζ 2 ij ) and J 3 = 9 2 tr(ζ 3 ij ), beingζ ij = ζ ij + δ ij J 1 /3 the components of the traceless part of the Hessian. We can then factorize the exponential into the following form and we define the correlation This decomposition allows to write Eq. (117) in a more compact way as a product of a bivariate Gaussian N in the variables ν and J 1 and 3-and 5-degrees of freedom χ-squared distributions in 3η 2 and 5J 2 respectively, and P(Ω) 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. [224,225]). The variable J 3 is uniformly distributed and constrained to satisfy J 2 3 ≤ J 3 2 by the fact that ζ 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, where we just adapted Eq. (96) to the variable R s , with the condition that δ − B ≤ 0, where now the prime indicates derivative with respect to R s . Hence, the corresponding ESP discrete number density can be written as n ESP (ω) = − µ σ 0s ν θ H (µ) n pk (y).
where we have defined the variable µ = B − δ and extended the vector of variables to ω = (y, µ). Combining Eq. (116) and Eq. (97) we get [201] n ESP (ν, R s )∆R s = n pk (y)P 1 (ω) (124) and using a similar approach to Eq. (A1), we can calculate the ESP multiplicity function to be where V = m/ρ, σ 0s = dσ 0s /dR s and G 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 [207] 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 [207] 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 ω = {ν, J 1 , µ, 3η 2 , 5J 2 , J 3 }. This expansion takes the name of Gram-Charlier expansion and it has been developed in a number of papers [199,203,226,227].
For generic non-Gaussian initial conditions we define P NG = P G [1 + δP NG ] and the Gram-Charlier 11 expansion up to third order in δ gives where we define the Hermite polynomials the generalized Laguerre polynomials and the correlations · GC as where we have assumed m ≤ 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 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 ω defines itself a bias parameter for the Gaussian case, for instance Hence we have found that each term in Eq. (131) defines a third order bias multiplied by the corresponding "skewness" of the type X 3 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 δn NG ESP , which is rather lengthy, but straightforward to calculate. Here we quote only the first three terms where we define the skewnesses as where g and χ vary for each variable X considered.

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 [208,209] focused on local-type non-Gaussianity with an amplitude of O( f NL ) ∼ 100s. In [209] they fitted the multiplicity function as where A, B, C and D each follow the functional form This fit showed agreement within ∼ 10% of the measurements. In [193,229], 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 δ sc , parametrizing the deviation with two quantities: the rate β at which the collapse threshold deviates on average from the spherical collapse δ 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 β 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 [230] 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 [231], 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 NL ≈ 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 [17]), the estimate of [231] is still preliminary and should be refined.
In [232], 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 [53, [233][234][235][236][237][238][239] (see discussion in the next Sec. §5.1.3), this type of searches should be investigated further.

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 [240,241] and x-ray detected clusters [242][243][244], 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 NL of the local type [245]. On the theoretical side, calculations within the context of the extended Press-Schechter mass function were perfomed in [175,246,247] and later using the excursion set formalism [190]. It is also possible to constrain directly the dark matter density field using weak lensing measurements [171,186,[248][249][250][251].
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.

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.

The breakthrough: scale-dependent bias from local PNG
In [208], 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 [252][253][254], multivariate bias expansions [255,256], and using a peak-background split (PBS) ansatz [257][258][259][260][261][262]. 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 [257].

Derivation with the peak-background split ansatz
The peak-background split approach [196,223] 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 δ. 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, 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 δ, but on Gaussian primordial fluctuations where indicates long modes and s short ones. Applying the split to Eq. (137) we get On sufficiently large scales, the relation holds, while the short-wavelengths modes of the density field are affected, at first order, as which can be interpreted as a local modulation of the amplitude of matter fluctuations by the long modes proportional to f NL , As a result, the local number density of halos of mass M is not a function of δ only, but also of this local modulationσ 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 δ , as From this expression we can define the first order bias in Fourier space in the presence of primordial non-Gaussianity the usual Gaussian Lagrangian bias and the amplitude of the so called scale-dependent bias. Here σ 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 ζ G .
The characteristic feature of this additional bias contribution at first order is that it scales as ∆b NG 1 ∝ 1/k 2 at large scales, since M(k) ∝ 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 loc NL is negative, with respect to the Gaussian case. In Fourier space, we therefore have where here we should remember that also the matter auto-power spectrum has non-Gaussian corrections, according to Eq. (84). As first argued in [257], 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 [257,263].

Universality of the mass function
The amplitude b NG 1 of Eq. (146) cannot be directly measured from galaxy surveys, as we do not have multiple realizations of the Universe with varying σ 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 σ 8 takes the simple form 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 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 [264][265][266]), 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.

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. (33), of order O( f NL ) ∼ O( , η) . 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, O( f NL ) ∝ n s − 1. This non-Gaussianity can be connected to a local-type one by an appropriate matching [47], 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 [53, [233][234][235][236][237][238][239]267,268], 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 [47], 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 [269,270]. An example in which this subtraction is done for the resonant running model of Eq. (51) is found in [271].

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.

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 we can expand the primordial bispectrum of ζ for k k s and get where L J are Legendre polynomials which appear only with even orders since in the squeezed limit k 1 ≈ −k 2 , removing the degree of freedom for odd number of exchanges of k s into −k s . The factor A J (k s , k ) 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. Table 1. List of values of A 2J from the squeezed limit expansion of Eq. (150) for interactions produced in single-field models of inflation.  Table 2. List of values of 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 and r s , where s indicates the spin of the particle, for principal and complementary series respectively. Full expressions for these amplitudes can be found in [88]. Notice that, as long as we are constraining non-Gaussianity using the halo power spectrum only, any information on the angle between k s and k is lost. One needs to look at the bispectrum [272][273][274][275][276][277][278][279][280][281][282][283][284][285][286][287] or, alternatively, to galaxy shapes [20,21] and intrinsic alignments [22] 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 for scale-invariant bispectra. We can then explicitly write the modulation of a specific long mode k on the local power spectrum around a small Lagrangian patch centered in q as where in the second term in the square brackets we have wrote explicitly the second order Legendre polynomial L 2 (l , 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 q, as compared to the Eulerian position at which we usually make observations of the (galaxy) density field, x [256,272,288]. Eq. (152) suggests a generalization of Eq. (142) by defining operators such as for J = 0 and for J = 2. The integration allows to account for the collective contribution from all long modes and has generically support for k < Λ, where Λ 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 k . These fields can be evaluated at the Eulerian position x by expanding around x(q, τ) = q + s(q, τ), being s the Lagrangian displacement, in perturbation theory 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 for the case of J = 0 at leading order. Here b ψ (∆) 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 where we have defined = a 0 ψ (∆) . Consequently, we get 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/5 f NL and ∆ = 0. For any ∆ > 0 the enhancement on large scales is reduced. For instance, non-Gaussianity from DBI inflation has ∆ = 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): ψ (2) ∝ ∇ 2 ζ is directly proportional to the matter density. Similarly, higher values of ∆ give contributions which are degenerate with bias parameters related to higher-derivative operators of the form ∇ ∆−2 δ. The intermediate range 0 ≤ ∆ ≤ 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 ψ ij (∆) needs to be multiplied by some operator O ij in order to appear in the bias expansion of Eq. (143), simply because the halo overdensity δ 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 σ 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 ∼ 0.02 h/Mpc. Given that observational data hardly constrain modes much larger than k ∼ 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 /k s and get 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 δ 2N to the bias expansion. For instance, for J = 0 we get for the halo auto-power spectrum 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 ψ (∆) with a bias generated by the operator ∇ 2 δ, we get a contribution approximately of the same order as b ψ (∆+2) . Moreover, the transfer function M(k) also contributes in this expansion when considering the full matter bispectrum in the squeezed limit [289] with terms which become relevant at the matter-radiation equality scale, k 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 [290]. 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 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 [68] 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 /k s ) 0 and (k /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 where we have also integrated over the angle between k s and k . We have now to distinguish two regimes: the ultra-squeezed limit, in which k k * , can be described by the scale-invariant ansatz of Eq. (151) with ∆ = 2. On the other hand, as k approaches k * , the product of k 2 s /k 2 * × k 2 /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 [271] 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.

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 [291].
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 ∼ R 3 * imposes that only second derivatives of the gravitational potential appear in the expansion. Therefore, we can write the expansion 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 δ and the tidal field 12 ∂ i ∂ j Φ. 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 [288,292]. 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 [225,[293][294][295][296][297]. Moreover, small deviations from locality can be accounted for perturbatively in terms of higher-derivative terms such as ∇ 2 δ 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 [298][299][300]. 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 [291], while original papers date back a few years earlier [301,302] 13 .
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 ζ, or equivalently Φ, itself are allowed in the expansion. The coupling in this case takes place at the Lagrangian position q, since it is imprinted in the initial conditions of structure formation. Following [291] we quote here the full bias expansion with PNG up to order J = 2 in the scale-invariant approximation of Eq. (151) 12 In the context of structure formation, it is conventional to work with the gravitational potential Φ, which in matter domination is simply related to the curvature perturbation by Φ = 3/5ζ. 13 See also [303] for a recent treatment in the context of quasi-single field inflation.
where we defined the tidal field The first line is leading order in δ and the non-Gaussian correction ψ, 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 ψ given that from current CMB observations we expect PNG corrections to be small [6]. We have only implicitly listed stochastic terms: at first order one needs to introduce an additional field, ψ with respect to the Gaussian case, which has only [288]. 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 δ 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 1/R * , higher-derivatives terms such as b ∇ 4 δ ∇ 4 δ can be neglected, as they typically scale as R 4 * k 4 δ [291].
Up to the order considered, Eq. (164) 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 [304][305][306][307], while analyses on universes with primordial non-Gaussianity are still missing.

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. (164) 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 70s. 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. (164). 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 [252,260,308], the excursion set approach [201,260,[308][309][310][311], the peak model [201] (see also [291] 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. (164).

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 [201,203,207,312,313]. The basic idea is to write the peak overdensity field δ L ESP as an effective perturbative expansion constructed from the rotational invariants of the system, namely the components of the vector ω defined in Eq.
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 δ L ESP (q)δ L ESP (q ) ⊃ c 2 ν that are generated when using the expansion to compute correlators of δ L ESP . The way to go is to write the expansion in terms of an entire set of orthogonal polynomials O n [ω(q)] where the indices n = {i, j, k, , m, n} correspond to the six variables of ω, σ [n] is an abbreviation for σ [n]  The polynomials of the form where L are Legendre Polynomials and x 3 = J 2 3 /J 3 2 .
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 [203]. The coefficients b L n are renormalized bias parameters which can be measured through 1-point ensemble averages [203] where the presence of the moments σ 0 , σ 1 , σ 2 and ς 0 ensures that the expansion is written in terms of the physical variables δ, η and ζ 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 δ where we have adopted the notation following previous literature. Let us make a few comments on this expansion. As compared to the Eulerian expansion of Eq. (164), 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 [275,[314][315][316][317]. This particular formulation of perturbation theory is specifically suited to embed the ESP framework. For the purpose of this review, working with Eq. (171) is enough to get the intuition of the physics at play. Secondly, Eq. (171) does not include any tidal shear, that is, scalar combinations of the field 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. (171).
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 ζ in the power spectrum, we need at least terms proportional to δ 3 R . 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 (175) where and 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. (177) allows to compute the leading order non-Gaussian contribution for any given primordial bispectrum.
It is instructive to take the squeezed limit k → 0 of Eq. (177) 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 where we defined 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 [227] it was proven that the two expressions are equivalent (180) 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 NG ESP = δ 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 ν c but also the spectral moments σ i . As argued in [201,202], the ESP mass function recovers the result for universal mass functions due to the fact that spectral moments appear only in ratios such as ν c /σ 0 or γ 1 = σ 2 1 /σ 0 σ 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 [154,155,193,194,210,222,[318][319][320][321][322][323][324] 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, being β 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 NG ESP = δ 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 [258]. In the context of the ESP model, in [202] 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 [206], 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 [220]. As remarked above, tidal shear is completely missing from the ESP model as we presented it up to now. As shown in [323], a significant part of the scatter in the collapse barrier is related to the local value of the tidal shear at halo formation. In [291], 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 σ 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 σ 0 . While [291] 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.

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 [119,121,206,[208][209][210]262,[325][326][327]. Extensions to other shapes than the local case have also been investigated [119,121,254,308,308]. The importance of this task is clear from Eq. (144): the non-Gaussian bias amplitude b NG 1 is degenerate in scale with f NL , so any theoretical or numerical uncertainty in its prediction affects the measurement of f 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 ∆ > 0. The degeneracy can be alleviated by measuring the 1/k ∆ scaling in the power spectrum of tracer populations with different mass and redshift using the multi-tracers technique originally proposed in [328], optimal weighting [326,329] and shot-noise suppression [330], but it is nonetheless important to have a handle on the uncertainty of b NG 1 . In this sense, the major feature of Lagrangian models presented in the previous section, that is, being able to predict the value of b NG 1 , 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.

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 ∆b NG 1 we are after, as first pointed out in [209,325]. At leading order, it reads where P NG and P 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 G hm is the linear Gaussian bias, which can be measured from simulations by measuring the ratio P G hm (k)/P G mm (k) at large scales, • ∆b NG I ( f NL ) is the scale independent correction which arises as a consequence of the change of the halo mass functionn 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 NL affects the high mass tale ofn h . For the same reason, the correction has opposite sign with respect to f NL , since the bias decreases (increases) whenever the halo mass function is enhanced (suppressed) with positive (negative) f NL , • β 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 [148] for a recent computation) or directly measured from simulations, • ∆b NG 1 (k) is the term we want to measure which is proportional to 1/k 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 NL [206].
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 ∆b NG 1 (k), because the non-Gaussian bias amplitude is proportional to the linear Gaussian bias b NG univ = δ sc b 1 . Since the earliest measurements [209,325], 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 [209], 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 [210], 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 [121,262,[325][326][327].
To settle this problem, [206] has compared the measurements of the non-Gaussian bias amplitude both with the universal prescription of b NG univ = δ sc b 1 and the peak-background split model-independent amplitude b NG PBS = ∂lnn h /∂lnσ 8 , see Fig. 3. The latter prescription was directly measured from simulations by running a set of Gaussian simulations with varying matter amplitude σ 8 and calculating the numerical derivative of the mass functionn 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 PBS NG , to the standard universal prediction δ sc b G 1 as a function of b 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  [331]) and right panel for a Friends-of-Friends algorithm with linking length λ = 0.2. Credits to [206]. 1 for each finder. Credits to [206].
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 δ sc < 1.687 at the high mass end of the halo mass function. This would possibly explain why b univ NG with δ sc = 1.687 overestimates b NG at high mass [210]. However, a change in δ sc cannot explain the fact that the amplitude b NG 1 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 σ 0 . Besides running simulations with different σ 8 as done in [206], one could also perform this measurement using separate universe simulations [327].

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 [262], 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 [232,301]. 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 [119,120], but also in CMB applications [117,118]. This is the reason why templates are often used in practical applications. In [121], 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 [333] 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 [334] generate a scale independent correction.
As a consequence of these considerations, [120,121] have introduced methods to implement a generic inflationary bispectrum in the initial conditions of N-body simulations. For instance, [121] has introduced the following ansatz 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 [120] is separable, but the power spectrum of Φ receives spurious corrections at large scales [121].

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 [16], 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.

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 [249][250][251]335] and galaxy intrinsic alignments [20][21][22]. • The Sunyaev-Zeldovich (SZ) effect [336] 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 [337]. The thermal SZ power spectrum [338] and kinetic SZ tomography [339] 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 [340][341][342]. • The scale-dependent bias outlined in Sec. §5 can be constrained using extrema counts as predictedy by peak theory [231] • The two-point statistics of voids, when combined with halos, give an order O(1) improvement in local-type non-Gaussianity constraints [343] • The covariance of galaxy number counts also helps in combination with number counts and variance [344] • Higher-redshift probes have been shown to be promising in constraining non-Gaussianity in the future: Lyman-alpha forest [345,346], 21-cm power spectra [23][24][25][26][27][28], CO and CII lines [29,30], cosmic reionization [347][348][349] and cross-correlations with CMB measurements [350,351]. • The response of the small-scale power spectrum to a squeezed bispectrum is particularly effective for the imprint of primordial non-Gaussianity [283,286,352].
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 [234,282,[353][354][355][356][357][358][359][360] and the importance of these relativistic corrections has been quantified for forthcoming surveys [361][362][363].
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 [272][273][274][275][276][277][278][279][280][281][282][283][284][285][286][287]364]. 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 [16].
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 [328,[365][366][367] 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 [368]. 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 [326,330,369].
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 loc NL = −0.9 ± 5.1, f equi NL = −26 ± 47 and f orth NL = −38 ± 24 at 65% confidence level, respectively. LSS searches have been also putting constraints since more than a decade [15,244,257,[370][371][372][373][374][375][376][377][378][379][380][381][382], 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 [15] and gives −51 < f loc NL < 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-Gaussianities 14 and O(10) for equilateral and orthogonal ones [16,19,21,30,280,281,285,343,. Several future surveys such as SPHEREx [16], Euclid [408], MeerKAT+DES [411] and SKA [19] might achieve these limits.

Conflicts of Interest:
The author declares no conflict of interest.