<figure><img src="image_1.png" /></figure>

# Hunting high and low: disentangling primordial and late-time non-Gaussianity with cosmic densities in spheres

C. Uhlemann, E. Pajer, C. Pichon, T. Nishimichi, S. Codis, F. Bernardeau

► To cite this version:

C. Uhlemann, E. Pajer, C. Pichon, T. Nishimichi, S. Codis, et al. “Hunting high and low: disentangling primordial and late-time non-Gaussianity with cosmic densities in spheres.” Monthly Notices of the Royal Astronomical Society, 2018, 474 (3), pp. 2853–2870. 10.1093/mnras/stx2623. cea-01677174

HAL Id: cea-01677174

https://cea.hal.science/cea-01677174

Submitted on 8 Jan 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Hunting high and low: disentangling primordial and late-time non-Gaussianity with cosmic densities in spheres

C. Uhlemann,<fnref n="1" />* E. Pajer,<fnref n="1" /> C. Pichon,<fnref n="2" />,<fnref n="3" /> T. Nishimichi,<fnref n="4" />,<fnref n="5" /> S. Codis<fnref n="2" />,<fnref n="6" /> and F. Bernardeau<fnref n="2" />,<fnref n="7" />

Accepted 2017 October 5. Received 2017 September 27; in original form 2017 August 8

## ABSTRACT

Non-Gaussianities of dynamical origin are disentangled from primordial ones using the formalism of large deviation statistics with spherical collapse dynamics. This is achieved by relying on accurate analytical predictions for the one-point probability distribution function and the two-point clustering of spherically averaged cosmic densities (sphere bias). Sphere bias extends the idea of halo bias to intermediate density environments and voids as underdense regions. In the presence of primordial non-Gaussianity, sphere bias displays a strong scale dependence relevant for both high- and low-density regions, which is predicted analytically. The statistics of densities in spheres are built to model primordial non-Gaussianity via an initial skewness with a scale dependence that depends on the bispectrum of the underlying model. The analytical formulas with the measured non-linear dark matter variance as input are successfully tested against numerical simulations. For local non-Gaussianity with a range from $f_{\text{NL}} = -100$ to $+100$, they are found to agree within 2 per cent or better for densities $\rho \in [0.5, 3]$ in spheres of radius $15 \text{ Mpc h}^{-1}$ down to $z = 0.35$. The validity of the large deviation statistics formalism is thereby established for all observationally relevant local-type departures from perfectly Gaussian initial conditions. The corresponding estimators for the amplitude of the non-linear variance $\sigma_8$ and primordial skewness $f_{\text{NL}}$ are validated using a fiducial joint maximum likelihood experiment. The influence of observational effects and the prospects for a future detection of primordial non-Gaussianity from joint one- and two-point densities-in-spheres statistics are discussed.

**Key words:** methods: analytical – methods: numerical – large-scale structure of Universe – cosmology: theory.

## 1 INTRODUCTION

Deviations from primordial Gaussianity are key to understand the physics of the early Universe, in particular inflation, and will soon be tested via a number of ambitious, deep, wide-field large-scale structure (LSS) surveys like Euclid (Laureijs et al. 2011), DES (Sánchez et al. 2014), DESI (Levi et al. 2013), LSST (LSST Science Collaboration 2009) and SPHEREx (Doré et al. 2014). To pin down the parameters that characterize primordial non-

Gaussianity (pNG) as accurately as possible, it is wise to consider the information contained within different probes of LSS in light of our ability to disentangle primordial from non-linear effects caused by gravitational evolution.

This paper revisits the question of how much information about pNG can be extracted from simple one- and two-point densities-in-spheres measurements within the LSS of matter, hence without shape information from the three-point correlation function or the bispectrum. More precisely, it will consider the full one-point probability distribution function (PDF) of dark matter densities in spheres thereby consistently including cumulants of all orders. It updates and generalizes the analyses based on perturbative predictions for

<footnote n="1">¹Center for Extreme Matter and Emergent Phenomena, Institute for Theoretical Physics, Utrecht University, Princetonplein 5, NL-3584CC Utrecht, the Netherlands</footnote>

<footnote n="2">²Institut d'Astrophysique de Paris, Sorbonne Universités, UPMC Univ Paris 6 et CNRS, UMR 7095, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="3">³School of Physics, Korea Institute for Advanced Study (KIAS), 85 Hoegiro, Dongdaemun-gu, Seoul 02455, Republic of Korea</footnote>

<footnote n="4">⁴Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa 277-85 83, Japan</footnote>

<footnote n="5">⁵CREST, JST, 4-1-8 Honcho, Kawaguchi, Saitama 332-0012, Japan</footnote>

<footnote n="6">⁶Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St. George Street, Toronto, ON M5S 3H8, Canada</footnote>

<footnote n="7">⁷CNRS and CEA, UMR 3681, Institut de Physique Théorique, F-91191 Gif-sur-Yvette, France</footnote>

<footnote>* E-mail: cu226@cam.ac.uk</footnote>
skewness in Fry & Scherrer (1994) and Durrer et al. (2000) and kurtosis in Chodorowski & Bouchet (1996), spherical-collapse inspired predictions for cumulants in Gaztanaga & Fosalba (1998) or the full PDF derived from a steepest descent approach in Valageas (2002b), as well as the purely simulation-based approach to a joint study of variance, skewness and kurtosis in Mao et al. (2014). It will relate the statistics of densities in spheres (a refined wording for densities-in-cells emphasizing the symmetry of the shape of the cells; see e.g. Bernardeau, Pichon & Codis 2014) to previous studies of PDFs of the dark matter density (Grossi et al. 2008) and related works on the impact of pNG on abundances of haloes (as well as galaxies and clusters residing therein) from Press–Schechter or excursion set based theories (Matarrese, Verde & Jimenez 2000; Scoccimarro, Sefusatti & Zaldarriaga 2004; LoVerde et al. 2008; Jimenez & Verde 2009; Lam & Sheth 2009; Lam, Sheth & Desjacques 2009; D’Amico et al. 2011a) and voids as well as galaxies residing therein (Kamionkowski, Verde & Jimenez 2009; Song & Lee 2009; D’Amico et al. 2011b).

Extending one-point statistics estimators, this paper also determines the two-point sphere bias describing the clustering of densities in spheres in the presence of pNG. Building upon previous works on scale-dependent halo bias in real and redshift space (in particular Dalal et al. 2008; Matarrese & Verde 2008; Desjacques, Seljak & Iliev 2009; Valageas 2010), it generalizes the underlying idea to the two-point bias of average densities in spheres that simultaneously probes intermediate and extreme density environments such as voids and haloes.

The outline of the paper is the following: Section 2 puts densities-in-spheres statistics in the context of current efforts to constrain pNG using LSS. Section 3 briefly reviews the pillars of pNG and describes how pNG affects the cumulants of the smoothed density field by the example of local non-Gaussianity. Section 4 recalls the large deviation principle that allows us to obtain one- and two-point statistics for the density in spheres. The relationship between the initial density PDF, the initial rate function and correlations at large separation is in particular extended to non-Gaussian initial conditions. The spherical collapse model and the saddle-point approximation applied to a log-transformed density are implemented in order to obtain the final non-Gaussian one-point PDF and two-point bias of densities in spheres. Section 5 presents the analytical predictions for the one-point density PDF and two-point sphere bias and compares them to measurements from N-body simulations performed with non-Gaussian initial conditions. Section 6 discusses how the analytical PDFs can be used to constrain pNG based on densities-in-spheres observables that are accessible from LSS and addresses promising perspectives. Section 7 wraps up. Appendix A gives some background for large deviation statistics while Appendix B presents generalizations and details of the treatment in the main text.

## 2 CONTEXT: pNG IN ALL SHAPES AND SIZES

For the purpose of constraining pNG, two aspects determine how advantageous a particular observable is. First, how much information about pNG the observable contains, and second how well theoretical errors and observational biases can be controlled in order to extract this information and disentangle primordial from evolution effects. The following provides a very short overview of the different probes used to constrain pNG using the cosmic microwave background (CMB) and LSS and motivates why densities-in-spheres statistics provides complementary information that can prove useful for constraining pNG.

For an in-depth discussion of inflationary models that generate pNG, relevant references and further details of the reviews in Bartolo et al. (2004), Chen (2010), Desjacques & Seljak (2010b), Verde (2010), Takahashi (2014), Alvarez et al. (2014) and Renaux-Petel (2015) are of interest.

### 2.1 CMB versus LSS

Currently, the best constraints on pNG come from the CMB, hence essentially from two-dimensional information at one redshift; it is however a rather direct and clean probe of density fluctuations, which only underwent linear evolution. Combining temperature and polarization in the latest Planck results (Planck Collaboration et al. 2016a) gives $1\sigma$ intervals for the amplitude of primordial skewness depending on the shape of the primordial bispectrum of $f_{\text{NL}}^{\text{loc}} = 0.8 \pm 5.0$, $f_{\text{NL}}^{\text{equil}} = -4 \pm 43$ and $f_{\text{NL}}^{\text{ortho}} = -26 \pm 21$. Following the standardized normalization method for $f_{\text{NL}}$ from Ferguson & Shellard (2009), observational limits and errors on $f_{\text{NL}}^{\text{loc}}$ can be consistently compared for different models, finding that despite appearances, upon this normalization current bounds from the CMB are of comparable order for the different shapes.

Reaching the theoretically interesting benchmark of $f_{\text{NL}} \lesssim 1$ will require to extract yet more information that is in principle available within LSS because it probes three-dimensional information for many modes over a wide range of redshifts (that can be obtained accurately from spectroscopic surveys). Unfortunately, this information is obscured by non-linear gravitational evolution, galaxy bias and redshift space distortions, and hence requires accurate theoretical modelling to disentangle dynamically induced from pNG. Note that most LSS-based probes are only sensitive to a particular range of bispectrum shapes, and hence their ability to constrain the different parameters varies. Eventually one should of course consider combined constraints from the CMB and LSS (see e.g. Giannantonio & Percival 2014; Giannantonio et al. 2014), especially to constrain scale-dependent non-Gaussianity; see e.g. LoVerde et al. (2008), Sefusatti et al. (2009), Schmidt & Kamionkowski (2010) and Becker, Huterer & Kadota (2012). For the next generation ‘Stage-4’ ground-based CMB experiment, CMB-S4 (Abazajian et al. 2016), the expected factor of improvement over Planck is slightly more than a factor of 2 which will be insufficient to reach the interesting theoretical benchmark for local $f_{\text{NL}}$.

### 2.2 Abundance and large-scale clustering of rare objects

Both the abundance and the clustering signal of rare objects, namely extrema of the density distribution such as haloes or voids, are suitable to probe pNG. The effect of pNG on the halo mass function and cluster number counts has been studied theoretically in Matarrese et al. (2000), LoVerde et al. (2008), Afshordi & Tolley (2008), Valageas (2010) and Musso & Paranjape (2012) within the framework of an extended Press–Schechter theory and compared to simulations in Grossi et al. (2008) and LoVerde & Smith (2011). Verde et al. (2013) considered multivariate PDFs that characterize the effect of $f_{\text{NL}}$ on different observables that can be obtained from CMB maps, but also the PDF of the linearly extrapolated and smoothed density field. In all those studies, it has been found that the presence of pNG modifies the abundances of haloes and voids increasingly with the rarity of the object. While in principle the presence of pNG might also alter the density profiles of rare objects, the impact on halo profiles has been observed to be small (Moradinezhad Dizgah, Dodelson & Riotto 2013) and a similar study for
voids that is underway shows no perceptible differences on void profiles.<fnref n="1" />

Furthermore, following an observation from Dalal et al. (2008), pNG leads to strongly scale-dependent bias that impacts the large-scale clustering of rare objects (peaks of the density distribution such as haloes) which hence can be used to constrain pNG from LSS (Slosar et al. 2008). Subsequently, pNG-induced halo bias and its dependence on scale, mass and redshift in the limit of high peaks and large separations has been theoretically investigated in Matarrese & Verde (2008), compared to simulations in Desjacques et al. (2009), Grossi et al. (2009), Shandera, Dalal & Huterer (2011) and refined in a treatment in real space in Valageas (2010).

Scale-dependent bias has already been used to put constraints on pNG based on galaxies and quasars within the SDSS in Slosar et al. (2008), Ross et al. (2013), Leistedt, Peiris & Roth (2014) finally reaching an individual constraint on $-49 < f_{NL} < 31$ (95 per cent confidence) for local pNG. The prospects of constraining pNG in future surveys with two-point statistics like the 3D galaxy power spectrum in redshift space, the angular galaxy power spectrum and the projected weak-lensing shear power spectrum have been investigated in Giannantonio et al. (2012) finding a forecasted marginalised error of $\sigma(f_{NL}) \simeq 3$ on local pNG for Euclid-like survey combining all probes. This is in line with current forecasts of $\sigma(f_{NL}) \simeq 4-5$ (competitive with Planck) for local pNG achievable with Euclid (Amendola et al. 2016) by combining galaxy clustering from spectroscopic and photometric data as well as weak lensing.

With recent advancements, it should be possible to improve constraints for local pNG from the two-point clustering of single populations of LSS tracers, which are already at the level of pre-Planck constraints from the CMB. To this end, Seljak (2009) demonstrated how the effect of cosmic variance for the estimation of pNG can be mitigated using two sets of tracers which, at large scales, are deterministically biased with respect to dark matter. This allows us to extract the scale-dependent relative bias due to pNG with an accuracy set by the noise due to the discrete (Poisson) sampling of the density field rather than cosmic variance. Hamaus, Seljak & Desjacques (2011) investigated constraints on local pNG from two-point statistics when using two optimization strategies – avoiding sampling variance by comparing multiple tracers of different bias and suppressing shot noise by optimally weighting haloes of different mass.

Byun & Bean (2015) investigated the dependence of the halo mass function and scale-dependent halo bias on the primordial bispectrum shape. They found that the scale-dependent bias on large scales probes general squeezed configurations from the primordial bispectrum, while the scale-dependent bias on smaller scales and the halo mass function are more sensitive to a broader range of shapes. Looking at the not-quite-large scales is also interesting because effects of general relativity (GR) become significant in the very large scale regime, precisely where the scale-dependent bias induced by pNG is strongest (see e.g. Bruni et al. 2012; Camera, Santos & Maartens 2015), although it has been shown that GR effects on scale-dependent bias are not degenerate with $f_{NL}$ (Yoo et al. 2012).

## 2.3 Bispectrum and trispectrum of the galaxy distribution

Since constraining pNG aims at probing the component of the initial bispectrum induced by non-Gaussianity, looking at its observable

late-time equivalent, the galaxy bispectrum, is a natural possibility that in particular provides means of distinguishing between different shapes. For local pNG with small amplitude $f_{NL}$, one can show that the primordial bispectrum estimator is optimal and equivalent to calculating the full likelihood of the data (Creminelli, Senatore & Zaldarriaga 2007). Unfortunately, the theoretical modelling of the late-time dark matter bispectrum induced by gravitational evolution already proves to be difficult. In addition to the modelling of the dark matter bispectrum, observational effects such as tracer bias (Kaiser 1984; Dekel & Rees 1987) and redshift space distortions (Kaiser 1987; Taruya, Nishimichi & Saito 2010) need to be included in the analysis.

Constraints on $f_{NL}$ from measurements of the galaxy bispectrum in redshift surveys have been considered as early as in Scoccimarro et al. (2004) concluding that planned galaxy surveys at high redshifts can in principle give pNG constraints comparable to, or even better than, those from CMB experiments. This was confirmed in Sefusatti & Komatsu (2007) by performing a Fisher matrix analysis for the galaxy bispectrum at high redshift with pNG and non-linear but local bias, including a study of how constraints of pNG improve with volume, redshift range, as well as the number density of galaxies. These works relying on tree-level perturbation theory for the matter bispectrum in the presence of pNG have been subsequently extended to include one-loop corrections (see e.g. Jeong & Komatsu 2009; Sefusatti 2009; Bernardeau, Crocce & Sefusatti 2010; Matsubara 2011) that improve the agreement with numerical simulations (Sefusatti, Crocce & Desjacques 2010). Following up on observations of scale-dependent bias in the two-point clustering of peaks, Jeong & Komatsu (2009) considered the bispectrum of galaxies as peaks finding that the effect of local $f_{NL}$ on the galaxy bispectrum cannot be obtained by replacing the linear bias in the galaxy bispectrum with the scale-dependent bias obtained for the power spectrum. Those perturbative analyses based on tree-level perturbation theory and a local bias model have been generalized in Baldauf, Seljak & Senatore (2011) to incorporate a multivariate bias expansion (Giannantonio & Porciani 2010) and the peak-background split method (Slosar et al. 2008). The validity of this complete tree-level approximation at large scales has been established by simulations with local pNG (Sefusatti, Crocce & Desjacques 2012) and used to get a rough Fisher-matrix-based estimate for an expected improvement of a few over the halo power spectrum for a combined power-spectrum and bispectrum analysis.

A more recent comparison of different perturbative and phenomenological models for the matter bispectrum with Gaussian and non-Gaussian initial conditions against numerical simulations has been performed in Lazanu et al. (2016) and Lazanu et al. (2017); see also references therein. It has been found that among the perturbative approaches, the Effective Field Theory of Large Scale Structure (EFTofLSS, as introduced in Baumann et al. 2012) extends the range of validity furthest on intermediate scales. However, this comes at the cost of introducing free extra parameters requiring calibration on simulations or marginalization. The EFTofLSS includes those free parameters to encode our ignorance of the small-scale physics that cannot be captured with perturbation theory (Bernardeau et al. 2002) and hence to trace theoretical errors.

Based on the matter bispectrum from EFTofLSS including pNG (Assisi et al. 2015), estimates for the amount of pNG while including theoretical errors in the modelling of the bispectrum were obtained in Baldauf et al. (2016) and Welling, van der Woude & Pajer (2016). They found that accounting for these theoretical errors can weaken constraints considerably, for example, degrading by a factor 5 from the idealized forecast for Euclid of $\sigma(f_{NL}) \simeq 0.45$.

<footnote n="1">¹ Private communication with N. Hamaus and K. Chan.</footnote>
(with floating bias) in Tellarini et al. (2016) down to $\sigma(f_{NL}) \simeq 1.8$ (with priors on EFT parameters) in Welling et al. (2016). In the latter, it has also been established that assuming the wrong shape for the theoretical error might lead to a false detection of pNG, highlighting how difficult the theoretical modelling proves to be. In addition, recent studies show that relativistic effects are relevant for the bispectrum at an order $\sigma(f_{NL}) \simeq 1$ and hint at a degeneracy with local pNG (Di Dio et al. 2017).

Regarding higher order galaxy statistics, Verde & Heavens (2001) argued, based on an idealized case with essentially linear evolution, that a measurement of the trispectrum, which has weaker dependence on non-linear clustering, may provide pNG constraints complementary to the bispectrum in the future.

## 2.4 Densities-in-spheres statistics

The statistics of densities in spheres (or counts-in-cells as their counterpart for discrete tracers) entail in particular the one-point PDF of finding a certain average density in a top-hat sphere of given radius and the two-point clustering of those spheres given their average density. While those statistics are related to the distribution of tracers and rare objects such as haloes and voids, they probe a different regime than the bispectrum of tracers or the abundance and clustering of rare objects and hence can improve our ability to probe pNG. One-point densities-in-spheres statistics already contains some of the information that is in all higher order correlation functions of the smoothed density field, namely their zero separation limit – characterizing the cumulants and thereby the shape of the PDF. The cumulants of the smoothed density field can be predicted accurately using perturbation theory or spherical collapse starting from Gaussian statistics (Bernardeau 1994a), but can also be extended to include pNG in an expansion around Gaussian initial conditions (Fry & Scherrer 1994; Chodorowski & Bouchet 1996; Gaztanaga & Fosalba 1998). An alternative approach to the statistics of densities in spheres is the excursion set model (Sheth 1998), which can be used to relate the halo mass function to the dark matter distribution in Eulerian space. Recently, also excursion set theory has been extended to non-Gaussian initial conditions (Musso & Sheth 2014). Since the PDF includes cumulants of all orders, it comprises more information than (a finite set of) single cumulants which is especially relevant in the tails that are generated through non-linear gravitational dynamics (Bernardeau 1994b), but also sensitive to pNG. The full density PDF and its cumulants have been considered for models with weak pNG arising from standard slow-roll inflation (Gaztanaga & Fosalba 1998; Valageas 2002b) and strong pNG arising from dimensional scaling models (Turok & Spergel 1991; Gaztanaga & Maehoenen 1996; White 1999; Scoccamarro 2000). Densities-in-spheres statistics include some higher order information beyond the bispectrum, but they only incorporate a very specific part of its rich shape-dependent information. This may seem to be a significant drawback when considering statistical information content for constraining pNG, but it is mitigated by our ability to build robust estimators through ab initio theoretical modelling. While cumulants of average densities in spheres can be predicted robustly from simple spherical collapse dynamics for both Gaussian (Bernardeau 1994b) and close-to-Gaussian initial conditions (Gaztanaga & Fosalba 1998; Valageas 2002b), modelling the bispectrum requires significantly more complex perturbative methods as discussed in the previous paragraph.

The one point PDF includes information about abundances of rare objects (haloes and voids) simultaneously in both tails, hence offering a unified treatment of all density environments. Similarly,

two-point densities-in-spheres statistics, which can be computed for large separations (Bernardeau 1996; Abbas & Sheth 2007; Codis, Bernardeau & Pichon 2016b; Uhlemann et al. 2017b) for Gaussian initial conditions, captures the large-scale clustering of rare objects (extreme density environments occurring in haloes or voids) in the tails but also encodes differences in clustering for more common density environments. Given that the most competitive pNG constraints from LSS to date come from the scale-dependence of halo bias, it is promising to generalize this idea to a more diverse range of density environments. The unified treatment of both high and low densities could in particular help to disentangle otherwise degenerate effects such as $f_{NL}$ and $g_{NL}$.

In spirit related to the PDF of densities in spheres is the pairwise velocity PDF whose sensitivity on local pNG has been investigated in Lam, Nishimichi & Yoshida (2011) where 5–10 per cent effects of pNG on the PDF were found for $f_{NL} = \pm 100$, but its accurate theoretical modelling proved to be difficult.

# 3 PRIMORDIAL NON-GAUSSIANITY AND CUMULANTS

Let us first express the leading-order cumulants in terms of local pNG parametrized by $f_{NL}$. We refer the readers to Appendix B1.1 for next-to-leading-order terms and non-local expansions.

## 3.1 Local primordial non-Gaussianity

For simplicity, the main text of the paper considers leading-order local (quadratic) pNG where the non-Gaussian field $\Phi_{NG}$ (the gravitational potential) is expressed in terms of a Gaussian field $\Phi_G$

$$ \Phi_{NG} = \Phi_G + f_{NL} (\Phi_G^2 - \langle \Phi_G^2 \rangle), \quad (1) $$

with constant parameter $f_{NL}$. For a generalization including the cubic order with $g_{NL}$, see Appendix B1.2. In principle, one can express the non-Gaussian density in terms of the Gaussian density and the gravitational potential using the identity

$$ \Delta\Phi_{NG} = \Delta\Phi_G + 2f_{NL}(\Phi_G\Delta\Phi_G + |\nabla\Phi_G|^2), \quad (2) $$

where $\Delta$ is the Laplacian and hence the left-hand side is proportional to the density. In practice, it is often more convenient to work in Fourier space<fnref n="2" /> where the Poisson equation becomes a simple multiplication and smoothing operates by multiplication rather than convolution.

Let us quantify how pNG influences the PDF of initial densities that are smoothed over a certain radius. Hence, one needs to relate the linearly evolved density field to the primordial curvature perturbation $\zeta$ encoding the information of non-linearities produced during and after inflation. We assume the dimensionless primordial power spectrum of the (comoving) curvature perturbation to be of the form

$$ \Delta_{\zeta}(k) = A_s(k_0)(k/k_0)^{(n_s-1)}, \quad (3) $$

with the amplitude of the scalar power spectrum $A_s \simeq 2.5 \times 10^{-9}$ measured at the pivot scale, $k_0 = 0.002/h$ Mpc$^{-1}$, and the primordial spectral index $n_s \simeq 0.96$ (consistent with recent limits from Planck Planck CollaborationXX 2016b), encoding initial conditions with a small departure from scale-invariance as predicted from slow-roll

<footnote n="2">² The Fourier convention $\delta(x) = (2\pi)^{-3} \int d^3k \exp(ik \cdot r) \delta(k)$ is used throughout. For the sake of simplicity, we will omit the hat on the Fourier transforms.</footnote>
models of inflation (Mukhanov & Chibisov 1981). The gravitational potential $\Phi^3$ is related to the curvature perturbation $\zeta$, which is preserved on super-Hubble scales, by a constant factor of 3/5 during the matter domination era such that its power spectrum is obtained as

$$P_{\Phi}(k) = \left(\frac{3}{5}\right)^2 \frac{2\pi^2}{k^3} \Delta_{\zeta}(k). \quad (4)$$

The gravitational potential $\Phi$ and the linear perturbation to the matter density $\delta_p$ at redshift $z$ are related through the Poisson equation on sub-Hubble scales

$$\delta_p(k, z) = g(k)T(k)D(z)\Phi_{\text{NG}}(k), \quad g(k) = \frac{-2}{3\Omega_m} \left(\frac{k}{H_0}\right)^2, \quad (5)$$

where $T(k)$ is the transfer function of perturbation (which goes to unity on very large scales $T(k) \xrightarrow{k \to 0} 1$) and encodes the suppression of power for modes that entered the horizon before matter-radiation equality, $D(z)$ is the linear growth factor which is related to the redshift as $D(z) = (1+z)^{-1}$ in matter domination, $\Omega_m$ is the present time fractional matter density and $H_0 = 100 \text{ h kms}^{-1} \text{ Mpc}^{-1}$ the Hubble constant. Finally, the spherical top-hat filtering on the density field can be applied to obtain

$$\delta_R(k) = W_{3D}(kR)\delta_p(k), \quad (6)$$

where $W_{3D}$ is the Fourier transform of the top-hat filter

$$W_{3D}(x) = \frac{3}{x^2}(\sin(x)/x - \cos(x)). \quad (7)$$

Now, pNG for the gravitational potential can be expressed conveniently as

$$\begin{aligned} \delta_{R,\text{NG}}(k, z) &= \alpha_R(k, z)\Phi_{\text{NG}}(k), \\ \alpha_R(k, z) &= W_{3D}(kR)g(k)T(k)D(z). \end{aligned} \quad (8)$$

The key assumption that allows for a general treatment of a mildly non-Gaussian field is that it can be expanded as a local functional of an underlying Gaussian field. On the other hand, any form of non-locality that can be expressed as a convolution in real space can also be handled in a similar way by suitably modifying our definition of the kernel functions $\alpha_R$. This means that our results can be easily extended to a primordial bispectrum of an arbitrary shape.

## 3.2 Leading-order cumulants of the smoothed density field

For one-point statistics, the skewness of the smoothed density field is the leading-order correction to Gaussian initial conditions. To this order, the effect of local pNG on the skewness is given by (see also Matarrese et al. 2000)

$$\begin{aligned} \kappa_3(R) &\equiv \langle \delta_{R,\text{NG}}^3 \rangle, \\ &\approx \frac{3 f_{\text{NL}}}{4\pi^4} \int dk_1 k_1^2 \alpha_R(k_1) P_{\Phi}(k_1) \int dk_2 k_2^2 \alpha_R(k_2) P_{\Phi}(k_2) \\ &\qquad \times \int_{-1}^1 d\mu_{12} \alpha_R \left( \sqrt{k_1^2 + k_2^2 + 2k_1 k_2 \mu_{12}} \right). \end{aligned} \quad (9)$$

For the non-local cases, the skewness (and also higher order cumulants) of the smoothed density field can be obtained from the primordial bispectrum $B_{\Phi}$ of arbitrary shape, as discussed in Appendix B1.1 based on equation (B4). Because of the smoothing

involved, the scale-dependence for different shapes is typically very similar, as demonstrated in Fig. B1, so the focus is here on local non-Gaussianity for simplicity.

It is often convenient to consider cumulants that are rescaled by certain powers of the variance, such as the reduced cumulants $\tilde{S}_n$ (which remain approximately constant through late-time gravitational evolution for Gaussian initial conditions; Colombi et al. 1997) or the rescaled cumulants $\epsilon_n$ (which are approximately independent of smoothing scale and robust against linear growth as it cancels in the ratio; D'Amico et al. 2011a). They enter any cumulant-based expansion – e.g. the Edgeworth expansion – of the cosmic PDFs around Gaussian kernels

$$\tilde{S}_n = \frac{\kappa_n}{\sigma^{2(n-1)}}, \quad \epsilon_n = \frac{\kappa_n}{\sigma^n} = \sigma^{n-2} \tilde{S}_n. \quad (11)$$

The tilde is used to distinguish the primordial reduced cumulants from the ones that are gravitationally evolved. The rescaled cumulants $\epsilon_n$ are to a very good approximation perturbatively ordered, as pointed out in (D'Amico et al. 2011a, where however a top-hat filter in k-space was used instead of a spherical top-hat filter in real space as considered here). Note that the smallness parameter is $\epsilon_3 \propto f_{\text{NL}}A^{1/2}$ where $A_s \simeq 2.510^{-9}$ is the amplitude of the primordial power-spectrum $P_{\Phi}$, similarly $\epsilon_4 \propto \epsilon_3^2 \propto f_{\text{NL}}^2 A$ or $\epsilon_4 \propto g_{\text{NL}} A$. Hence, for reasonably small $f_{\text{NL}} \simeq 100$ one can typically neglect next-to-leading-order contributions appearing in form of a kurtosis and a second-order correction to the variance. Appendix B1.2 presents the next-to-leading-order results (including the kurtosis and a second-order correction to the variance) for a local model including both $f_{\text{NL}}^2$ and $g_{\text{NL}}$.

For the two-point statistics, the leading-order mixed cumulant $\kappa_{12}(r)$ is also needed; it is obtained as

$$\begin{aligned} \kappa_{12}(r) &= \langle \delta_{R,\text{NG}}(\mathbf{x} + \mathbf{r}) \delta_{R,\text{NG}}^2(\mathbf{x}) \rangle \\ &= \frac{2 f_{\text{NL}}}{(2\pi)^6} \iint d^3 k_1 d^3 k_2 \alpha_R(k_1) \alpha_R(k_2) \\ &\qquad \times P_{\Phi}(k_1) [P_{\Phi}(k_2) + 2P_{\Phi}(k)] \exp[i\mathbf{k} \cdot \mathbf{r}], \end{aligned} \quad (12)$$

where $\mathbf{r}$ is the separation between the two points and $\mathbf{k} = \mathbf{k}_1 + \mathbf{k}_2$. Integrating over the angle between $\mathbf{k}$ and $\mathbf{r}$ yields

$$\begin{aligned} \kappa_{12}(r) = \frac{f_{\text{NL}}}{4\pi^4} &\int dk \, k^2 \alpha_R(k) \operatorname{sinc}(kr) \int dk_1 \, k_1^2 \alpha_R(k_1) P_{\Phi}(k_1) \\ &\times \int_{-1}^{1} d\mu_k \, \alpha_R(k_2) [P_{\Phi}(k_2) + 2P_{\Phi}(k)], \end{aligned} \quad (14)$$

where $k_2 = |\mathbf{k} - \mathbf{k}_1| = \sqrt{k_1^2 + k_2^2 - 2k_1 k_2}$ and $\mu_k$ is the cosine of the angle between $\mathbf{k}$ and $\mathbf{k}_1$. For a generalization of the formula for arbitrary bispectra see Appendix B1.1, in particular equation (B5). Naturally, when evaluating the leading-order mixed cumulant at zero separation, one correctly recovers the skewness from equation (9).

## 4 CONSTRUCTIONS OF THE DENSITY PDFS

Let us first briefly present the large deviation principle that allows us to obtain one-point statistics for the density in a sphere of fixed final radius, focusing on the extension to non-Gaussian initial conditions. Readers already familiar with the formalism for Gaussian initial conditions may skip to the extension to pNG given in Section 4.2.2 and then proceed with the phenomenological effects in Section 4.6.

<footnote n="3">³ Note that in conformal Newtonian gauge, the gravitational potential is related to the Bardeen potential as $\Phi_B = -\Phi$.</footnote>
## 4.1 Large deviation statistics in spheres

When considering a highly symmetric observable such as the density in spheres, one can argue that the most likely dynamics (amongst all possible mappings between the initial and final density field) is the one respecting the symmetry (Valageas 2002a).<fnref n="4" /> Spherical symmetry allows us to take advantage of the fact that non-linear solutions to the gravitational dynamics are known explicitly in terms of the spherical collapse model.

Let us denote $\rho_{SC}(\tau)$ is the non-linear transform of an initial fluctuation with linear density contrast, $\tau$, in a sphere of radius $R_{ini}$, to the final density $\rho$ (in units of the average density) in a sphere of radius $R$ according to the spherical collapse model

$$ \rho = \rho_{SC}(\tau), \quad \text{with} \quad \rho R^3 = R_{ini}^3, \qquad (15) $$

where the initial and final radii are connected through mass conservation. An explicit possible fit for $\rho_{SC}(\tau)$ is given by

$$ \rho_{SC}(\tau) = (1 - \tau/\nu)^{-\nu} \Leftrightarrow \tau_{SC}(\rho) = \nu(1 - \rho^{-1/\nu}), \qquad (16) $$

where $\nu$ can be adjusted to the actual values of the cosmological parameters ($\nu = 21/13$ provides a good description of the spherical dynamics for an Einstein-de Sitter background for the range of $\tau$ values of interest).

Thanks to this analytic spherical collapse model, the one-point PDF and bias functions of cosmic densities in concentric spheres, brought about by non-linear gravitational evolution, can be predicted explicitly from the given (close-to Gaussian) initial conditions.

## 4.2 Initial one-point PDF and decay-rate function

The large-deviation principle yields a formula for the PDF of finding a certain density in a sphere of a given radius given the statistics of the initial conditions. The initial decay-rate function $\Psi_{R_{ini}}^{\text{ini}}$ encodes the exponential decay of the initial PDF with increasing density contrast given the cumulants of the initial field. While in general the initial PDF is related to the cumulant generating function via an inverse Laplace transform, this relation is well approximated by a saddle-point approximation for close-to-Gaussian initial conditions. Then, the initial PDF can be written as

$$ \mathcal{P}_{R_{ini}}^{\text{ini}}(\tau) = \sqrt{\left| \frac{\partial^2 \Psi_{R_{ini}}^{\text{ini}}}{\partial \tau^2} \right| \frac{\exp[-\Psi_{R_{ini}}^{\text{ini}}(\tau)]}{(2\pi)^{1/2}}}. \qquad (17) $$

### 4.2.1 Gaussian initial conditions

For Gaussian initial conditions, the initial decay-rate function is simply given by a quadratic form in the initial density contrast $\tau$ in a sphere of radius $R_{ini}$ where the linear variance $\sigma_{\text{lin}}^2$ encodes all dependency with respect to the initial power spectrum for the Gaussian field

$$ \Psi_{R_{ini}}^{\text{ini}}(\tau) = \frac{1}{2} \frac{\tau^2}{\sigma^2(R_{ini})}, \qquad (18) $$

$$ \sigma_{\text{lin}}^2(R_{ini}) = \int \frac{dk}{2\pi^2} k^2 P_{\text{lin}}(k) W_{3D}^2(k R_{ini}). \qquad (19) $$

Note that in this case, equation (17) is merely an unusual rewriting of a Gaussian distribution, emphasizing the central role of the decay-rate function (18). The next step involves obtaining a suitable form of the initial decay-rate function for non-Gaussian initial conditions based on a cumulant expansion whose result is given in equation (21) below. For this analysis, the key ingredient of the decay-rate function is the scale-dependent skewness of the smoothed density field which can be computed for a local model using equation (9) and for any given bispectrum shape according to equation (B4).

### 4.2.2 Non-Gaussian initial conditions: cumulant expansions

The complete series of cumulants of a PDF is related to the cumulant generating function whose successive derivatives give back the cumulants. Taking the modifications in the cumulants of the smoothed density up to quadratic order in non-Gaussianity (including skewness, kurtosis and second-order variance) into account leads to the following truncated cumulant generating function

$$ \tilde{\phi}_{R_{ini}}(\lambda) = \frac{\lambda^2}{2} \left[ \kappa_{2, R_{ini}}^{(0)+(2)} + \frac{\lambda}{3} \kappa_{3, R_{ini}} + \frac{\lambda^2}{12} \kappa_{4, R_{ini}} \right], \qquad (20) $$

where (0) and (2) refer to the zeroth and next-to-leading order in $f_{\text{NL}}$, respectively. Note that this equation only holds for local quadratic pNG as is assumed throughout the main text. Performing a Legendre transformation of the cumulant generating function then gives the initial rate function

$$ \Psi_{R_{ini}}(\tau) \approx \frac{\tau^2}{2\sigma(R_{ini})^2} \left[ 1 - \frac{\varepsilon_3(R_{ini})}{3} \frac{\tau}{\sigma(R_{ini})} - \varepsilon_2^{(2)}(R_{ini}) + \left( \frac{\varepsilon_3(R_{ini})^2}{4} - \frac{\varepsilon_4(R_{ini})}{12} \right) \frac{\tau^2}{\sigma(R_{ini})^2} \right], \qquad (21) $$

using equation (9) to express the leading-order term $\varepsilon_3$, and equations (B7) and (B9) for the next-to-leading-order terms in terms of the underlying power spectrum and $f_{\text{NL}}$; see also equations 50 and 58 in Matarrese et al. (2000) and Valageas (2010), respectively. Since $\varepsilon$ are the dimensionless initial cumulants, the linear growth cancels out in this quantity but does enter in the variance $\sigma$. Since the focus is in the tails while not relying on perturbation theory, $\tau/\sigma_r$ need not be small. This is the main difference of our analysis which hence complements the perturbation-theory-based analysis for the expansions for the mass function in D'Amico et al. (2011a). Hence, the cumulant-expanded form of the rate function in equation (21) is to be preferred over an Edgeworth expansion around a Gaussian PDF that is discussed in Appendix B2. In order for the truncation in the series of primordial cumulants to be accurate, one must have that $1 \gg \varepsilon_3\tau/\sigma \gg \varepsilon_4(\tau/\sigma)^2$. This in principle limits how far into the tails one is allowed to probe when just considering primordial skewness. To obtain an estimate for the allowed range given an $f_{\text{NL}}$, one can use the spherical collapse solution (16) to obtain the initial overdensity $\tau \to \tau(\rho)$ and the radius $R_{ini} = R\rho^{1/3}$ connected to a final density $\rho$ in a sphere of radius $R$. For $f_{\text{NL}} \lesssim 100$, higher cumulants are still suppressed by one order of magnitude for typical densities $\rho \in [0.1, 10]$ and hence it will be sufficient to account for the primordial skewness for our analysis.

## 4.3 One-point PDF for an evolved non-Gaussian field

The final decay-rate function is obtained from re-expressing the initial decay-rate (21) in terms of the final densities and radii

$$ \tau_{R_{ini}} = \tau_{\text{SC}}(\rho_R), \quad R_{ini} = R \rho^{1/3}, \qquad (22) $$

<footnote n="4">⁴ This is a result of the so-called contraction principle in the context of large deviation theory as explained in Bernardeau & Reimberg (2016), which formalizes the idea that amongst all unlikely fates (in the tail of the PDF) the least unlikely one (i.e. the spherical collapse solution) dominates.</footnote>
using the spherical collapse mapping equation (16), which corresponds to a saddle-point approximation<fnref n="5" /> and rescaling the linear variance at scale R to its correct non-linear value

$$ \sigma_{\text{lin}}(R_{\text{ini}}) \rightarrow \sigma(R_{\text{ini}}) \approx \frac{\sigma(R)}{\sigma_{\text{lin}}(R)} \sigma_{\text{lin}}(R \rho^{1/3}). \quad (23) $$

Altogether the decay-rate function reads

$$ \Psi_R(\rho) = \Psi_{R_{\text{ini}}=R\rho^{1/3}}^{\text{ini}} \left( \frac{\sigma_{\text{lin}}(R)}{\sigma(R)} \tau_{\text{SC}}(\rho) \right), \quad (24) $$

where $\Psi_{R_{\text{ini}}}$ is given by equation (21). From that decay-rate function, one can obtain the cumulant generating function and from there the PDF. Appendix A provides a concise explanation of the cumulant generating function and its relation to the decay-rate function and the PDF. As discussed in Uhlemann et al. (2016), one can use a saddle point approximation to accurately evaluate the corresponding expressions and obtain a direct relation between the final decay-rate function and the PDF. This requires a suitable choice of variable such that the final decay-rate function is convex, i.e. $\Psi''_{R} > 0$ for all densities of interest. One particularly advantageous choice is the logarithmic density $\mu = \log \rho$ for which one can apply the saddle-point approximation to predict the PDF of the logarithmically mapped density and translate this to the PDF of the density field via a simple change of variables

$$ \begin{aligned} \mathcal{P}_R(\rho) &= \mathcal{P}_{\mu,R}[\mu(\rho)] \left|\frac{\partial \mu}{\partial \rho}\right| = \sqrt{\frac{\partial^2 \Psi_R}{\partial \mu^2}} \left|\frac{\partial \mu}{\partial \rho}\right| \frac{\exp[-\Psi_R]}{(2\pi)^{1/2}} \\ &= \sqrt{\frac{\Psi''_R(\rho) + \Psi'_R(\rho)/\rho}{2\pi}} \exp[-\Psi_R(\rho)]. \end{aligned} \quad (25) $$

Note that since the decay-rate function $\Psi_R$ appearing here has been written down for the log-density, the variance $\sigma(R) = \sigma_\mu(R)$ that is entering equation (24) is the variance of the log-density $\mu = \log \rho$. To ensure a unit mean density after the mapping, one has to shift the log-density appropriately which is equivalent to considering the renormalized density $\tilde{\rho} = \rho/\langle\rho\rangle$ with the shorthand notation $\langle f(\rho) \rangle = \int_0^\infty d\rho f(\rho) \mathcal{P}_R(\rho)$. Furthermore, since the saddle-point method yields only an approximation to the exact PDF, the PDF obtained from equation (25) is not necessarily perfectly normalized (although this effect is usually subpercent). In practice, one can account for both effects by considering the normalized PDF

$$ \hat{\mathcal{P}}_R(\rho) = \mathcal{P}_R \left( \rho \frac{\langle \rho \rangle}{\langle 1 \rangle} \right) \frac{\langle \rho \rangle}{\langle 1 \rangle^2}. \quad (26) $$

Equations (20)–(26) yield the general evolved non-Gaussian PDF of the density field.

## 4.4 Initial two-point PDF and sphere bias

The initial sphere bias used in previous works on densities in spheres (see e.g. Bernardeau 1996; Codis et al. 2016b; Uhlemann et al. 2017b) is defined as the ratio between the conditional mean density induced by an initial density contrast $\tau$ in a sphere of radius $R_{\text{ini}}$ at separation $r$ and the average correlation

$$ b_{\text{ini}}(\tau, r) = \frac{\langle \tau'(\mathbf{x} + \mathbf{r}) | \tau(\mathbf{x}) \rangle}{\langle \tilde{\tau}'(\mathbf{x} + \mathbf{r}) \tilde{\tau}(\mathbf{x}) \rangle} = \frac{\int d\tau' \tau' \mathcal{P}_{R_{\text{ini}}} (\tau, \tau'; r)}{\mathcal{P}_{R_{\text{ini}}} (\tau) \xi_{\text{lin}} (R_{\text{ini}}, r)}, \quad (27) $$

and can be expressed in terms of the joint PDF of densities in spheres of radius $R_{\text{ini}}$ at separation $r$ and their correlation function. The initial correlation function is obtained from the linear power spectrum as

$$ \xi_{\text{lin}}(R_{\text{ini}}, r) = \int \frac{dk}{2\pi^2} k^2 P_{\text{lin}}(k) W_{3D}^2(k R_{\text{ini}}) \text{sinc}(kr). \quad (28) $$

Note that this expression only depends on the density $\tau$, the radius of the spheres $R_{\text{ini}}$ and their separation $r$ while the densities $\{\tau', \tilde{\tau}', \tilde{\tau}\}$ are dummy variables for the evaluation of the spatial expectation values. The sphere bias can be used to express the joint PDF of densities in spheres of radius $R$ at large separation $r \gg R$

$$ \frac{\mathcal{P}_{R_{\text{ini}}} (\tau, \tau'; r)}{\mathcal{P}_{R_{\text{ini}}} (\tau) \mathcal{P}_{R_{\text{ini}}} (\tau')} \approx 1 + \xi_{\text{lin}}(R_{\text{ini}}, r) b_{\text{ini}}(\tau, r) b_{\text{ini}}(\tau', r). \quad (29) $$

### 4.4.1 Gaussian initial conditions

To derive the initial bias, let us follow the procedure described in appendix C of Uhlemann et al. (2017b). The covariance matrix of $(\tau, \tau')$ is given by

$$ \Sigma_{\text{lin}} = \begin{pmatrix} \sigma_{\text{lin}}^2(R_{\text{ini}}) & \xi_{\text{lin}}(R_{\text{ini}}, r) \\ \xi_{\text{lin}}(R_{\text{ini}}, r) & \sigma_{\text{lin}}^2(R_{\text{ini}}) \end{pmatrix}, \quad (30) $$

where the linear variance is given by (19) and the correlation function by (28). This covariance matrix can be diagonalized by transforming $(\tau, \tau')$ to a set of independent variables $(\nu, \zeta)$

$$ \nu = \frac{\tau}{\sigma}, \quad \zeta = \frac{\sigma^2 \tau' - \xi \tau}{\sigma \sqrt{\sigma^4 - \xi^2}}, \quad (31) $$

which are built to be decorrelated $\langle \nu \zeta \rangle = 0$ and have unit variance $\langle \nu^2 \rangle = \langle \zeta^2 \rangle = 1$. For Gaussian initial conditions, thanks to diagonalization, $(\nu, \zeta)$ now follow a standard normal distribution, such that it is easy to check that the sphere bias reads

$$ b_{\text{ini}}^{G}(\tau) = \frac{\langle \tau'(\zeta, \nu) | \nu = \tau/\sigma \rangle}{\xi(r)} = \frac{\tau}{\sigma_{\text{lin}}^2(R_{\text{ini}})}, \quad (32) $$

which is proportional to the initial overdensity $\tau$ as expected from Kaiser (1984). Furthermore, it is independent of $r$ such that the separation and density dependence in equation (29) can be factorized.

### 4.4.2 Non-Gaussian initial conditions: cumulant expansions

For initially non-Gaussian fields, one has to redo the computation for the conditional mean that defines the initial bias function, equation (27), which now might acquire a dependence of the separation that was not present for Gaussian initial conditions. To obtain an expression for the two-point correlation for initially non-Gaussian fields, the joint cumulant generating function was expanded following equation (20). The leading-order correction to a bivariate Gaussian PDF will be the first term of a bivariate Edgeworth expansion (Chambers 1967; McCullagh 1984; Kendall, Stuart & Ord 1987) given in terms of third-order cumulants, the skewness $\kappa_3 = \langle \delta^3 \rangle$ given in equation (9) and the mixed cumulant $\kappa_{12} = \langle \delta' \delta^2 \rangle = \kappa_{21}$ given in equation (12).<fnref n="6" /> In the decorrelated variables $\nu$ and $\zeta$ introduced in equation (31), this expansion is given in terms of products of Hermite polynomials (see e.g.

<footnote n="5">⁵ Note that according to Valageas (2010) the effect of a realistic $|f_{\text{NL}}| \lesssim 100$ on the density profile of the saddle point is negligible such that also the onset of shell-crossing remains practically unchanged compared to the Gaussian case presented in Valageas (2009).</footnote>

<footnote n="6">⁶ Note that this treatment resembles what is performed in Chongchitnan & Silk (2011), where the scale-dependence of the bias is computed based on a bivariate Edgeworth expansion for models with pure $g_{\text{NL}}$.</footnote>
Barndorff-Nielsen & Cox 1979), which for the first non-Gaussian correction gives

$$ \mathcal{P}(\nu, \zeta) \approx \frac{\exp\left(-\frac{\nu^2}{2} - \frac{\zeta^2}{2}\right)}{2\pi} \left[ 1 + \frac{1}{6} \left( \langle \nu^3 \rangle H_3(\nu) + \langle \zeta^3 \rangle H_3(\zeta) \right) + \frac{1}{2} \left( \langle \nu \zeta^2 \rangle H_1(\nu) H_2(\zeta) + \langle \nu^2 \zeta \rangle H_2(\nu) H_1(\zeta) \right) \right], \quad (33) $$

with $H_1(x) = 1$, $H_2(x) = x^2 - 1$ and $H_3 = x^3 - 3x$. For an Edgeworth expansion up next-to-leading order that includes the joint kurtosis see Appendix B2.2 and Appendix B1.2 for the associated cumulants induced by $f_{\text{NL}}^2$ and $g_{\text{NL}}$. The conditional mean that determines the sphere bias is

$$ \langle \tau' | \tau \rangle = \frac{\int d\nu \int d\zeta \, \tau'(\nu, \zeta) \mathcal{P}(\nu, \zeta) \delta_D(\tau(\nu, \zeta) - \tau)}{\int d\nu \int d\zeta \, \mathcal{P}(\nu, \zeta) \delta_D(\tau(\nu, \zeta) - \tau)}, $$

where inverting equation (31) gives

$$ \tau = \nu\sigma, \quad \tau' = \frac{\nu\xi + \zeta\sqrt{\sigma^4 - \xi^2}}{\sigma}. \qquad (34) $$

After some algebra and using the properties of Hermite polynomials, and expressing the moments of $\nu$, $\zeta$ back via the moments of $\tau$, $\tau'$, the sphere bias at leading order in $f_{\text{NL}}$ becomes

$$ b_{\text{ini}}^{\text{NG}}(\tau, r) = \frac{\tau}{\sigma^2} \left[ 1 + \frac{1}{2} \frac{\tau}{\sigma^2} \left( \frac{\kappa_{12}(r)}{\xi(r)} - \frac{\kappa_3}{\sigma^2} \right) \right] - \frac{1}{2\sigma^2} \left( \frac{\kappa_{12}(r)}{\xi(r)} - \frac{\kappa_3}{\sigma^2} \right), \quad (35) $$

where the cumulants $\kappa_3$, $\kappa_{12}$, $\sigma$ and $\xi$ from equations (9), (14), (19) and (28), respectively, are computed with smoothing radius $R_{\text{ini}}$ and the term in the second row ensures a zero mean bias $\langle b(\tau) \rangle = 0$. For an analogous result for a pure $g_{\text{NL}}$ model, see equation (B16) in Appendix B2.2 and for the generalization of cumulants to arbitrary primordial bispectra, see Appendix B1.1. In contrast to the result for Gaussian initial conditions, the result is now separation dependent and the average two-point correlation function explicitly appears. Let us observe that the non-Gaussian sphere bias is a sum of the Gaussian sphere bias and a separation-dependent correction term

$$ b_{\text{ini}}^{\text{NG}}(\tau, r) = b_{\text{ini}}^{\text{G}}(\tau) + \Delta b_{\text{ini}}^{\text{NG}}(\tau, r). \quad (36) $$

This expression shows that in analogy to the strong scale-dependent bias for overdense regions (haloes), one expects an analogous scale-dependent bias for underdense regions (voids). This could be a promising idea for constraining pNG; while it has been mentioned in Sekiguchi & Yokoyama (2012), no exhaustive study of this effect has been performed yet, but is underway.<fnref n="7" />

## 4.5 Two-point sphere bias for an evolved non-Gaussian field

The two-point sphere bias for the time-evolved density field is defined in complete analogy to equation (27) as

$$ b_R(\rho, r) = \frac{\langle \rho' | \rho; r \rangle - 1}{\langle \bar{\rho} \rho'; r \rangle - 1}, \quad (37) $$

where the numerator is the conditional mean sphere density given a sphere density $\rho$ at separation $r$ and the denominator is the average sphere correlation function. As was done for the one-point PDF, a

replacement rule for the initial density and radius is used following the spherical collapse (22), and the variance is rescaled according to equation (23). For Gaussian initial conditions and at large separations, this recovers the separation-independent expression given in Uhlemann et al. (2017b)

$$ b_R^G(\rho) = \frac{\sigma_{\text{lin}}^2(R)\tau_{\text{SC}}(\rho)}{\sigma^2(R)\sigma_{\text{lin}}^2(R\rho^{1/3})}, \quad (38) $$

where the redshift dependence is encoded in the non-linear variance $\sigma^2(\zeta, R)$ and hence bias grows like $D^{-2}(z)$. The pNG-induced extra terms of equation (35) can be rewritten in terms of reduced cumulants giving the combination $\tilde{C}_{12} - \tilde{S}_3$, where $\tilde{C}_{12} = \kappa_{12}/(\xi\sigma^2)$. Those reduced cumulants are robust against non-linearities such that no further rescaling is needed. Due to the power of (linear) densities appearing in the reduced cumulants, the redshift dependence of this term is $\propto D^{-1}(z)$. For Gaussian initial conditions, the sphere bias is scale independent, meaning that there is no residual dependence on pair separation. For non-Gaussian initial conditions, the situation is different and one needs to relate separation-dependent quantities at early and late times, more precisely the initial (Lagrangian) and final (Eulerian) separation of the two spheres (of identical density). Valageas (2010) suggest to treat each sphere as a test particle that falls into the potential well caused by another sphere and hence to use the mass conservation mapping for both positions to obtain

$$ r_{\text{ini}} \approx r \left( 1 + \frac{2}{3} \tau_{\text{profile}}(\rho, R_{\text{ini}}, r) \right), \quad (39) $$

where $\tau_{\text{profile}}$ is the linear density contrast within radius $r$ around the sphere of radius $R_{\text{ini}}$ with density $\rho$. It probes the radial density profile of the spherical saddle point which has been obtained in Bernardeau (1994b) and Valageas (2009),

$$ \tau_{\text{profile}}(\rho, R_{\text{ini}}, r) = \tau_{\text{SC}}(\rho) \frac{\sigma_{\text{lin}}(R_{\text{ini}}, r)^2}{\sigma_{\text{lin}}(R_{\text{ini}})^2}. \quad (40) $$

For large separations $r \gtrsim 60 \text{ Mpc h}^{-1}$ and the range of radii $R \simeq 15 \text{ Mpc h}^{-1}$ and densities of interest here $\rho \in [0.5, 3]$, the density contrast according to this profile is at the 5–10 per cent level. While there are in principle modifications of the density profile in the presence of pNG (as calculated in Valageas 2010), their effect is negligible to leading order. Altogether the non-Gaussian bias is now given as

$$ b_R^{\text{NG}}(\rho, r) = b_R^G(\rho) \left[ 1 + \frac{\tau_{\text{SC}}(\rho)}{2} (\tilde{C}_{12} - \tilde{S}_3)(R_{\text{ini}}, r_{\text{ini}}) \right], \quad (41) $$

where the result is normalized to ensure $\langle b_R^{\text{NG}}(\rho, r) \rangle = 0$ and $\langle \rho b_R^{\text{NG}}(\rho, r) \rangle = 1$ as follows:

$$ \hat{b}_R^{\text{NG}}(\rho, r) = \frac{b_R^{\text{NG}}(\rho, r) - \langle b_R^{\text{NG}}(\rho, r) \rangle}{\langle (\rho - 1) b_R^{\text{NG}}(\rho, r) \rangle}. \quad (42) $$

### 4.5.1 Scale-dependent halo bias in Fourier space

In Fourier space, scale-dependent halo bias has been derived using several approaches, as summarized in Grossi et al. (2009); this includes peak theory (Dalal et al. 2008), the bispectrum (Matarrese & Verde 2008), peak-background split (Slosar et al. 2008), ellipsoidal collapse (Afshordi & Tolley 2008), perturbation theory (McDonald 2008; Taruya, Koyama & Matsubara 2008) and excursion sets (Adshead et al. 2012). All of them lead to consistent results for the scale-dependent correction to the linear bias factor induced by leading-order local non-Gaussianity

$$ b_{\text{h,NG}}(k) = b_{\text{h,G}} + \Delta b_{\text{h,NG}}(k), \quad (43) $$

<footnote n="7">⁷ Private communication with Hamaus N. and Chan K.</footnote>
$$ \Delta b_{h,NG}(k) = 2(b_{h,G} - 1)f_{NL}\delta_c \frac{3\Omega_m}{2D(z)r_H^2k^2}, \quad (44) $$

where $r_H$ is the Hubble radius, $\delta_c$ is the critical threshold for collapse and $D(z)$ is the growth rate. A benchmark of theoretical predictions against numerical simulations has been performed in Desjacques et al. (2009), Pillepich, Porciani & Hahn (2010), Giannantonio & Porciani (2010) and Smith, Desjacques & Marian (2011). The presence of an extra term, corresponding to the term proportional to the initial skewness $\epsilon_3$ in equation (41), has been noted in Slosar et al. (2008) and Afshordi & Tolley (2008). While it was argued in Slosar et al. (2008) that this scale-independent contribution does not cause a problem when fitting the bias to LSS data, Desjacques et al. (2009) found that it leads to substantial improvement when comparing theoretical prediction to simulations.

In this context, the generalization to other types of pNG can be done by promoting the constant $f_{NL}$ to a scale-dependent function $f_{NL}(k) \propto (kR_{L,h})^\Delta$ where $\Delta = 0$ gives back local pNG, $\Delta = 2$ corresponds to equilateral pNG and $\Delta \in [0, 2]$ includes quasi-single field inflation and some of its modifications. For a discussion of scale-dependent bias for bispectrum templates beyond local pNG, see e.g. Scoccolmarro et al. (2012) and Gleyzes et al. (2017).

### 4.5.2 Scale-dependent halo bias in real space

In real space, which is more natural for the treatment of densities in spheres, the scale dependence of halo bias has been derived in Matarrese & Verde (2008) based on an expansion of the two-point correlation for peaks as regions above a certain (high) threshold using a path-integral approach in Grinstein & Wise (1986), Matarrese, Lucchin & Bonometto (1986) and in Valageas (2010) with a steepest decent approach for leading-order local non-Gaussianity $f_{NL}$.

According to Matarrese & Verde (2008), the equal density two-point correlation for initially non-Gaussian fields in the limit of high peaks and large separations is given by the result for Gaussian initial conditions, $\xi_{h,G} = \xi_h(r)b_{ini}^{h,G}(\tau)^2$, and a correction term, also called scale-dependent bias, is given by

$$ \Delta\xi_h^{NG}(r_{ini}, \tau) = \kappa_{12}(r_{ini})b_{ini}^{G}(\tau)^3, \quad (45) $$

with $r = |\mathbf{x}_2 - \mathbf{x}_1|$ and with the leading-order mixed cumulant $\kappa_{12}(r)$ given by equation (14). Hence, this yields to a bias function with only the first two terms in equation (35). However, following the arguments in Valageas (2010), one has to include another (separation-independent) term proportional to $\epsilon_3\sigma$ that stems from the skewness of the initial PDF and leads to the first line of (35), but is missing the constant terms.

As explained in Desjacques et al. (2009), a contribution from the non-Gaussian modification to the (linear) power spectrum is also expected; it is however much smaller than the terms presented before and will therefore be neglected, as was done in Matarrese & Verde (2008) and Valageas (2010).

## 4.6 Phenomenological effect of primordial non-Gaussianities

Fig. 1 shows a comparison of the theoretical predictions for a fixed value of $f_{NL} = 100$ when varying the underlying non-linear variance in a way that corresponds to time evolution $\sigma(z) = \sigma_{15}D(z)$ with $\sigma_{15} \simeq 0.5$ and $D(z)$ the linear growth function. Note that although the non-linear variance does not agree with the linearly evolved variance at late times, the late-time non-linear variance on the scales of interest grows approximately with the linear growth function.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** (upper panel) The PDF of the log density $\mu = \log \rho$ in a sphere of radius $R = 15$ Mpc h$^{-1}$ comparing a Gaussian model (solid lines) with one with local pNG according to equation (1) for the gravitational potential $\Phi_\rho$ with $f_{NL}^\Phi = \pm 100$ (dashed and dotted lines) for variances from $\sigma_\mu = 0.05$ (blue inner lines) to 0.5 (red outer lines) in steps of 0.05 (lower panel). The corresponding ratio of the PDFs with and without pNG as a function of the density.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** The ratio of the one-point PDFs with $f_{NL} = 0$ and $f_{NL} = +100$ for variances from $\sigma_\mu = 0.05$ (blue) to 0.5 (red) in steps of 0.05 as a function of the rarity of the event for overdensities (dashed upper lines) and underdensities (dotted lower lines).</figcaption></figure>

One can see that the overall effect of pNG in the PDF is very weak except for the tails of the distribution. Note also that the apparent asymmetry between over and underdensities hints at a deviation from lognormality for high variances. Fig. 2 shows the relative effect of pNG on the PDF at early (low variances, blue curves) and late times (high variances, red curves) when fixing the probability of the to-be-found density value and looking at overdensities (dashed upper lines) and underdensities (dotted lower lines). Note the three properties as follows.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Density dependence of the scale-dependent two-point sphere bias $b_R(\rho, r)$ at redshift $z = 1$ with radius $R = 15$ Mpc h$^{-1}$ (where $\sigma_\mu \simeq 0.31$) and separations $r = 32, 64, 96$ Mpc h$^{-1}$ (solid, dotted, dashed lines) induced by pNG with $f_{\text{NL}} = \pm 100$ (lower blue and upper red lines, respectively).</figcaption></figure>

(i) For low variances corresponding to early times, the effect of pNG is significantly stronger for overdensities.

(ii) For high variances corresponding to late times, the effect of pNG is slightly stronger for underdensities.

(iii) The time (or variance) dependence is very mild for under-densities but substantial for overdensities. This is in line with the usual claim that underdensities are more pristine objects and better preserve the initial conditions.

While those observations are interesting from a theory point of view, they should not be naively translated to practical applications for which several complications (in particular possible inaccuracies of the PDF in the tails, the limited available volume and a discrete sampling of the density field in observations) need to be considered.

Fig. 3 displays a schematic picture of the scale dependence of the sphere bias $b_R(\rho, r)$ that is induced by pNG and describes density-dependent clustering. The scale-dependent sphere bias correction given in equation (41) scales like $\Delta b_R(\rho, r, f_{\text{NL}}) \propto f_{\text{NL}} b_{R,G}(\rho) r^2$ such that it increases quadratically with separation $r$ (see Fig. B2), the bias of the density $\rho$ in the sphere (hence the rarity of the event) and the amplitude of primordial skewness $f_{\text{NL}}$. As discussed before, the redshift dependence of scale-dependent spheres bias correction is roughly $\Delta b_R \propto 1/D(z)$.

# 5 VALIDATION ON SIMULATIONS

Let us now compare the prediction of equation (25) to density-in-spheres measurements in dark matter simulation encoding primordial non-Gaussianities.

## 5.1 Simulation

The simulations contain 2048³ particles in a box of length 4096 Mpc h$^{-1}$ and a given realization has been run for models with $f_{\text{NL}} = +100, 0, -100$. The relevant cosmological parameters of the simulations (based on WMAP5) are summarized in Table 1 and snapshots are created at redshifts $z = 3, 2, 1, 0.35$. The non-Gaussian initial condition generator was developed in Nishimichi (2012) based on a parallel code by Valageas & Nishimichi (2011),

<table><caption>**Table 1.** Cosmological parameters of the simulation.</caption><thead><tr><th>Ω<sub>m</sub></th><th>Ω<sub>Λ</sub></th><th>h</th><th>A<sub>s</sub>(k<sub>0</sub> = 0.002/h Mpc<sup>-1</sup>)</th><th>n<sub>s</sub></th><th>σ<sub>8</sub></th></tr></thead><tbody><tr><td>0.279</td><td>0.721</td><td>0.701</td><td>2.486 × 10<sup>-9</sup></td><td>0.96</td><td>0.8157</td></tr></tbody></table>

<table><caption>**Table 2.** Parameters of the simulation that characterize primordial skewness in terms of $f_{\text{NL}}$ and the non-linear variance of the log-density $\mu = \log \rho$ for $R = 15$ Mpc h$^{-1}$ at redshifts $z = 0.35, 1$ and $2$ as measured from the simulation and as fitted using the PDF template from equation (25). Note that differences in the non-linear variance between different non-Gaussian models are subpercent.</caption><thead><tr><th></th><th colspan="3">0.35</th><th colspan="3">1.00</th><th colspan="3">2.00</th></tr></thead><tbody><tr><td>f<sub>NL</sub></td><td>-100</td><td>0</td><td>+100</td><td>-100</td><td>0</td><td>+100</td><td>-100</td><td>0</td><td>+100</td></tr><tr><td>σ̂<sup>PDF</sup><sub>μ</sub></td><td>0.419</td><td>0.418</td><td>0.417</td><td>0.315</td><td>0.314</td><td>0.313</td><td></td><td>0.218</td><td></td></tr><tr><td>σ̂<sub>μ</sub></td><td>0.415</td><td>0.414</td><td>0.413</td><td>0.312</td><td>0.312</td><td>0.311</td><td></td><td>0.218</td><td></td></tr><tr><td>Ŝ<sup>μ</sup><sub>3</sub></td><td>0.231</td><td>0.294</td><td>0.358</td><td>0.198</td><td>0.286</td><td>0.374</td><td>0.148</td><td>0.277</td><td>0.406</td></tr><tr><td>-Ŝ<sup>μ</sup><sub>4</sub></td><td>0.421</td><td>0.367</td><td>0.293</td><td>0.529</td><td>0.461</td><td>0.354</td><td>0.631</td><td>0.547</td><td>0.376</td></tr><tr><td>σ̂<sub>ρ</sub></td><td>0.441</td><td>0.442</td><td>0.444</td><td>0.323</td><td>0.323</td><td>0.324</td><td></td><td>0.222</td><td></td></tr><tr><td>Ŝ<sup>ρ</sup><sub>3</sub></td><td>3.31</td><td>3.40</td><td>3.5</td><td>3.24</td><td>3.34</td><td>3.45</td><td>3.16</td><td>3.30</td><td>3.45</td></tr><tr><td>Ŝ<sup>ρ</sup><sub>4</sub></td><td>19.6</td><td>20.9</td><td>22.3</td><td>18.4</td><td>19.8</td><td>21.4</td><td>17.3</td><td>19.2</td><td>21.2</td></tr></tbody></table>

which computes the particle displacements using the second-order Lagrangian perturbation theory (e.g. Scoccimarro 1998; Crocce, Pueblas & Scoccimarro 2006). The gravitational evolution is followed using the Tree-PM code GADGET2 (Springel, Yoshida & White 2001; Springel 2005).

Table 2 provides the values for the variance, the driving parameter of the theory, measured in the simulation for both the log-density $\mu = \log \rho$ and the density $\rho$. The differences in the variances measured at $R = 15$ Mpc h$^{-1}$ are at the subpercent level which is in qualitative agreement with the finding in Mao et al. (2014).

The PDF was measured with a top-hat filter of radius $R = 15$ Mpc h$^{-1}$ using an FFT-based method with a clouds-in-cells mass assignment. The convergence against the mass assignment schemes [either nearest grid point (NGP) or clouds-in-cells (CIC) with different number of grid points] has been tested yielding stable results for intermediate and large densities $\rho$ where the difference between the two schemes is less than 1 per cent and potential problems in the highly underdense regions with $\rho < 0.4$, where there is a difference between the results obtained with NGP and CIC. The results shown are based on a CIC mass assignment with $1280^3$ grid points, with both window and aliasing corrections implemented before multiplying the top-hat function in k-space.

To measure the sphere bias, encoding the excess correlation between densities in spheres according to equation (37), a given separation $r$ is chosen, giving a grid of non-overlapping spheres. In that grid, the densities of the neighbouring spheres with right separation $r$ are then collected in bins of width $\Delta \rho = 0.1$; precise formulas are given by equations 19 and 20 in Uhlemann et al. (2017b).

## 5.2 Comparison of analytical predictions against simulations

### 5.2.1 One-point PDF

Fig. 4 shows a comparison of the theoretical predictions for the one-point PDF $P_R(\rho)$ against the simulations when plugging in the correct values of $f_{\text{NL}}$ and the measured non-linear log-variance $\sigma_\mu$ that is given in Table 2. It shows that the prediction is in excellent agreement with the measurements, in particular it is accurate at the 1 per cent level between densities of 0.5 and 2.5 and therefore can access a mildly non-linear regime that is inaccessible to perturbative
<figure><img src="image_6.png" /><figcaption>**Figure 4.** (Upper panel) Comparison between measurement in the simulation (data points) for $f_{\text{NL}} = -100, 0, +100$ (blue, green, red) and the saddle point prediction (lines) at redshifts $z = 0.35, 1, 2$ (solid, dotted, dashed) for the PDF with variances $\sigma_{\mu}^{\text{PDF}}$ as given in Table 2. (Lower panel) Corresponding residuals between the prediction and the measurements. For better visibility the $\rho$-values for $f_{\text{NL}} = \pm 100$ have been displaced from the ones of $f_{\text{NL}} = 0$ by $\pm 0.02$.</figcaption></figure>

methods. Note that, while the values of $f_{\text{NL}}$ probed here are rather large, the per cent-level accuracy is expected to extend to all values between $f_{\text{NL}} = -100$ and $f_{\text{NL}} = 100$, in particular to values around $f_{\text{NL}} = 0$ (previously also tested in Uhlemann et al. 2016, 2017b). Fig. 5 shows the ratio of the Gaussian and non-Gaussian PDFs as predicted and measured in the simulation; an exquisite agreement is found that shows that the inaccuracy of our PDF template acts the same way on the Gaussian and non-Gaussian version and hence comparing their ratio is in even better agreement with the measurements than what one would naively expect from the accuracy of the PDF which degrades for densities beyond what is shown in Fig. 4. Based on this observation, the constraining power of the one-point

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Residuals between the PDFs with pNG $f_{\text{NL}} = \pm 100$ (red and blue) and the fiducial PDF with $f_{\text{NL}} = 0$ as measured in the simulation (data points with error bars) and predicted from the saddle point approximation (lines) at redshifts $z = 0.35, 1, 2$ (solid, dotted, dashed). The grey horizontal lines indicate pNG-induced deviations from the fiducial PDF at the 1 per cent (solid) and 5 per cent (dashed) level.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Two-point sphere bias for local pNG $f_{\text{NL}} = \pm 100$ (upper red and lower blue) and the fiducial PDF with $f_{\text{NL}} = 0$ (middle green) as measured in the simulation (data points with error bars) and predicted from the saddle point approximation (lines) at redshift $z = 1$ for spheres of radii $R = 15 \text{ Mpc h}^{-1}$ at separation $r = 64 \text{ Mpc h}^{-1}$.</figcaption></figure>

density PDF for $f_{\text{NL}}$ is determined in Section 6 based on a maximum likelihood (ML) approach.

### 5.2.2 Two-point sphere bias

Fig. 6 shows a comparison of the theoretical prediction for the density- and separation-dependent two-point sphere bias $b_R(\rho, r)$. It appears clearly that the effect of pNG is strongest in the low- and high-density tails, in good agreement with the theoretical prediction from equations (41) and 42. The predictions shown are evaluated using the approximation $r_{\text{ini}} \simeq r$ giving very similar results to the approximation from equation (39) and only shows differences in the extreme density regions which are however within the error bars. The differences between the measurements with $f_{\text{NL}} = \pm 100$
and $f_{\text{NL}} = 0$ indicating the impact of pNG on the clustering of certain density environments are in excellent agreement with the theoretical expectation. Note that, since the simulations have been run from the same realization, error bars for the differences between non-Gaussian and Gaussian models are displayed as the individual error bars rather than their sum.

# 6 NON-GAUSSIANITY PARAMETER ESTIMATORS

## 6.1 Maximum likelihood with one-point statistics

Having validated that the functional form of the one-point PDF (equation 25) against simulations, let us now use it to do an ML estimate of the amount of pNG. The density PDF $\mathcal{P}(\rho|\sigma, f_{\text{NL}})$ at a given redshift depends on two parameters, the (gravitationally induced) non-linear variance, $\sigma$, and the amplitude of pNG encoded in $f_{\text{NL}}$.

Let us conduct the following fiducial experiment, following the assumptions of Codis et al. (2016a) where the approximately linear growth of the non-linear variance was used to constrain dark energy that changes the growth of fluctuations. Here, one can also take advantage of the fact that the non-linear variance, when measured for spheres of radii $R > 10$ Mpc $h^{-1}$, evolves with the linear growth to a good approximation but keep the fiducial cosmology fixed except for allowing for pNG and changes in the overall amount of clustering. In order to mimic a Euclid-like survey, let us consider redshifts between $z_{\text{min}} = 0.6$ and $z_{\text{max}} = 2.0^8$ binned so that the comoving distance of one bin is $d = 40$ Mpc $h^{-1}$. For every redshift bin, a number of spheres that correspond to regularly drawn spheres of radius $R = 15$ Mpc $h^{-1}$ separated by $d = 40$ Mpc $h^{-1}$ are considered.<fnref n="9" /> For a 15 000 square degree survey, it yields 55 bins of redshift ($z_i$) with a number of spheres ranging from about $N_1 \simeq 7000$ (at $z_1 = 0.6$) to $N_{55} \simeq 40000$ (at $z_{55} = 2.0$) for a total of about 1.2 million supposedly independent spheres. In this experiment, note that samples are drawn directly from the PDF, hence assuming that the model for the density PDF is exact, meaning that the local amplitude $f_{\text{NL}}$ well parametrizes the amount of pNG in terms of skewness and non-linear variance $\sigma$ well describes the outcome of non-linear gravitational evolution from close to Gaussian initial conditions.

Note that the present analysis is conservative with respect to the number of spheres one can get given a certain volume, but optimistic with respect to the modelling assuming no theoretical error. On the one hand, shrinking the radius $R$ of the spheres increases the available information and hence tightens constraints. On the other hand, the accuracy of the theoretical prediction decreases with the radius such that the theoretical error becomes relevant. Hence, a good balance of the two must be struck. One possible improvement on the following analysis would be to decrease the radius of the spheres with increasing redshift (for example using linear

<table><caption>**Table 3.** Collection of mean ML results determined from 10 samples for $f_{\text{NL}}$ parametrizing the amount of primordial skewness $f_{\text{NL}}$ when keeping the non-linear variance $\sigma_{15}$ fixed at its fiducial value and vice versa.</caption><thead><tr><th></th><th>fid</th><th>ML</th><th>1σ</th><th>2σ</th><th>3σ</th></tr></thead><tbody><tr><td>$f_{\text{NL}}$</td><td>0.0</td><td>-2.5</td><td>±9.5</td><td>±16.5</td><td>±24.5</td></tr><tr><td>$\sigma_{15}$</td><td>0.514 50</td><td>0.514 45</td><td>±0.0004</td><td>±0.0007</td><td>±0.001</td></tr></tbody></table>

growth) to keep the variance, which controls the theoretical accuracy, constant over the redshift slices. In practice, the chosen treatment to scale the number of drawn spheres with the available volume (and hence not probing the deep PDF tails, especially at low redshifts) ensures that the PDFs entering the ML estimation are within their regime of validity (as determined from the comparison with simulations).

### 6.1.1 Estimate of primordial skewness

In order to get constraints on the parameter $f_{\text{NL}}$ of pNG, let us compute the log-likelihood of 10 randomly samples for the 1.2 million measured densities $\{\rho_{i,j}\}_{1\le i \le N_j, 1\le j \le 55}$ given models for which $f_{\text{NL}}$ varies

$$ \mathcal{L}(\{\rho_{i,j}\}|f_{\text{NL}}) = \sum_{j=1}^{55} \sum_{i=1}^{N_j} \log \mathcal{P}(\rho_{i,j}|z_j, f_{\text{NL}}), \quad (46) $$

where $\mathcal{P}(\rho|z_j, f_{\text{NL}})$ is the theoretical density PDF at redshift $z_j$ for a primordial non-Gaussianity model parametrized by $f_{\text{NL}}$. Optimizing the probability of observing densities $\{\rho_{i,j}\}_{1\le i \le N_j, 1\le j \le 55}$ at redshifts $\{z_j\}_{1\le j \le 55}$ with respect to $f_{\text{NL}}$, yields a ML estimate for the primordial non-Gaussian parameter.

The resulting mean ML values (averaged over 10 samples) and the corresponding $\alpha = 1\sigma, 2\sigma, 3\sigma$ confidence intervals given in Table 3 correspond to the models for which $\mathcal{L}(\{\rho_{i,j}\}|f_{\text{NL}}) = \max_{f_{\text{NL}}} \mathcal{L}(\{\rho_{i,j}\}|f_{\text{NL}}) + \log(1 - \text{Erf}(\alpha/\sqrt{2}))$. Modulo our assumptions, this ML method allows us to detect non-Gaussianity with $\sigma(f_{\text{NL}}) \simeq 10$.

### 6.1.2 Joint estimate of primordial skewness and variance

In order to get joint constraints on primordial skewness parametrized through $f_{\text{NL}}$ and the non-linear log-variance at given radius at present time $\sigma_{15} = \sigma_{\mu}(R = 15 \text{ Mpc } h^{-1}, z = 0)$ (hence equivalent to $\sigma_8$), let us re-compute the log-likelihood of the 1.2 million measured densities $\{\rho_{i,j}\}_{1\le i \le N_j, 1\le j \le 50}$ given models for which both $f_{\text{NL}}$ and $\sigma_{15}$ vary

$$ \mathcal{L}(\{\rho_{i,j}\}|\sigma_{15}, f_{\text{NL}}) = \sum_{j=1}^{50} \sum_{i=1}^{N_j} \log \mathcal{P}(\rho_{i,j}|z_j, \sigma_{15}, f_{\text{NL}}). \quad (47) $$

Fig. 7 depicts the result for the joint ML estimate computed for 10 different samples. It shows that the joint estimation of $f_{\text{NL}}$ and $\sigma_{15}$ shows no degeneracy, hence the signature of pNG is qualitatively different from a change in the normalization of the density fluctuation field, associated with $\sigma_8$ (or in our case $\sigma_{15}$). Having access to both (overdense and underdense) wings of the PDF is what helps distinguishing these two effects that essentially have different symmetry signatures (variance versus skewness). Indeed, the marginalization over $\sigma_{15}$ decreases the accuracy that can be obtained for $f_{\text{NL}}$ only marginally compared to the run where $\sigma_{15}$ was kept fixed at its fiducial value. Our analysis probes a range

<footnote n="8">⁸ Limiting the redshift range to $z_{\text{max}} = 2$ also ensures that the effect of shot noise, which essentially cuts off the signal at high $z$ (see e.g. fig. 3 in Welling et al. 2016) but is neglected here, is not too extreme.</footnote>

<footnote n="9">⁹ At high redshifts, the flat-sky approximation is not used; instead the volume of spherical shells of thickness $d$ is computed for all redshifts. To estimate the number of non-overlapping spheres that can be fit in the shells, the average cut surface of the cosmic volume is divided by the cut surface of cubes with side length $d$. This estimate is close to the one obtained from assuming a hexagonal close packing with a density of $\eta = \pi/(3\sqrt{2}) \simeq 0.74$ which is the densest possible packing of equal spheres.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Result of joint ML for non-linear variance $\sigma_{15} = \sigma(R = 15, z = 0)$ and pNG $f_{\text{NL}}$ with the mean of the maximum (point) and the corresponding $1\sigma$, $2\sigma$ and $3\sigma$ contours (dark blue to cyan) averaged over 10 samples drawn from the PDF at different redshift with fiducial values ($\sigma_{15}, f_{\text{NL}}$) = (0.5145, 0) (dashed cross).</figcaption></figure>

of matter densities that (a) is symmetric around mean density in the probability of finding them and (b) increases with redshift due to the larger cosmic volume available for the sampling. While the dependence on the choice of observable tracers in surveys has not been examined yet, the inclusion of tracer bias (proceeding along Uhlemann et al. 2017a) is not expected to qualitatively change the most relevant region of the density PDF.

## 6.2 Observational effects

**Tracer bias.** It has been demonstrated recently in Uhlemann et al. (2017a) that tracer bias can be incorporated in densities-in-spheres statistics by a local mapping between average densities in spheres. Based on mass-weighted halo densities in both real and redshift space, a local one-to-one mapping has been shown to provide accurate predictions for one- and two-point densities-in-spheres statistics. The relationship between dark matter and mass-weighted halo densities has been observed to admit a quadratic model in the log-densities as efficient parametrization. The main effect of tracer bias is to cause a strong degeneracy between the dark matter variance and bias parameters for the one-point PDF which can only be lifted using two-point statistics. Since, however, the dark matter variance has been shown not to be degenerate with a local $f_{\text{NL}}$, bias is not expected to diminish the constraints except for finite sampling effects. This is in line with the finding in Se fusatti & Komatsu (2007) that for a bispectrum analysis, a local $f_{\text{NL}}$ is not degenerate with (local) galaxy bias. Even discreteness should be reasonably controlled, given that for typical tracer number densities of order $n \simeq 10^{-4} (h \, \text{Mpc}^{-1})^3$ (as expected for Euclid) spheres of radius $R = 15 \, \text{Mpc} \, h^{-1}$ should contain a sufficient number of objects. In Mao et al. (2014), the signal of pNG in the first three one-point cumulants of the smoothed non-linear density field has been studied on numerical simulations and mock galaxy cata logues resembling LRGs. For dark matter densities, it was shown that the pNG signal in the skewness and kurtosis of the dark matter densities is much larger than that induced by the non-linear variance. Note that in the present formalism, given that the one-point PDF includes all cumulants, the skewness and kurtosis are both taken into account when jointly constraining the dark matter variance and $f_{\text{NL}}$ via the ML approach, which does improve the accuracy of the parameter estimation compared to sample variance (as shown in Codis et al. 2016a). However, for number-weighted galaxy densi ties from mocks, the finding in Mao et al. (2014) is in stark contrast, and the variance almost completely dominates the pNG signal. Re cently, Gleyzes et al. (2017) investigated the scale-dependent halo bias for equilateral and quasi-single field inflation generated pNG, and looked at how they can or cannot be distinguished from a gen-

eral biasing model. For local (and some quasi-single field) shapes, they found that significant improvements over Planck are achievable with scale-dependent bias for a large volume LSS surveys, while the sensitivity to equilateral non-Gaussianity is heavily suppressed by marginalizing over a general halo biasing model.

**Redshift space distortions.** The effect of redshift space distortions has been assessed in Mao et al. (2014) where it has been found that they affect the variance and skewness very similarly in the Gaussian and non-Gaussian case and therefore do not affect the detectability of non-Gaussianity from these measurements, see fig. 4 in Mao et al. (2014). Recently, Uhlemann et al. (2017a) showed that redshift space distortions can be incorporated naturally in a mean local bias model for densities-in-spheres statistics such as one-point PDF and the two-point sphere bias.

**GR effects.** There are also GR effects that generate non-Gaussianity. According to Bruni, Hidalgo & Wands (2014), the non-linear relation between the spatial curvature and the metric perturbation translates into a specific non-Gaussian contribution to the initial comoving matter density. Note also that there are subtleties in the definition of bias depending on rest frames and relativistic effects can enter through projection or past light-cone effects (see section 9 of Desjacques, Jeong & Schmidt 2016, for a review). Relativistic effects are most important on the very large scales that are relevant for the scale-dependent bias in Fourier space (Camera et al. 2015), while the intermediate separations probed here should not be affected as much. However, in Yoo et al. (2012) it has been shown that the GR effect is not degenerate with the pNG signature in galaxy bias, and the ability to detect pNG is little compromised.

## 6.3 Outlook

For local pNG, one should eventually constrain pNG using jointly the one-point PDF and two-point sphere bias statistics. This could in particular be interesting for distinguishing between $f_{\text{NL}}$ and $g_{\text{NL}}$ which lead to a very similar signature in scale-dependent halo bias (Desjacques & Seljak 2010a; Giannantonio & Porciani 2010; Desjacques, Jeong & Schmidt 2011a,b; Smith, Fer raro & LoVerde 2012), see the extensions discussed in Appendix B. Indeed, joint analyses of abundances and clustering have been performed for galaxy clusters in Sartoris et al. (2010) and Mana et al. (2013) and shown to improve constraints based on individ ual probes. The prospects of constraining $f_{\text{NL}}$ from sphere bias and disentangling $f_{\text{NL}}$ from $g_{\text{NL}}$ will be explored in a forthcoming work. To render densities-in-spheres statistics applicable to counts-in-cells measurements in galaxy surveys, one should furthermore include tracer bias and redshift space distortions. While this is be yond the scope of this paper, the present procedure to include pNG in densities-in-spheres statistics can be combined with the biased tracer formalism described in Uhlemann et al. (2017a) which also demonstrates a joint parameter estimation using the one-point PDF and the two-point sphere bias function.

Recently, Chiang et al. (2015) introduced the position-dependent correlation function (or equivalently power spectrum) as the correlation between two-point functions of galaxy pairs within different (large) subvolumes with a given (small) mean density contrast at their location. Based on this, it has been demonstrated that for local $f_{\text{NL}}$ the position-dependent correlation function can yield comparable constraints to the full bispectrum, despite its inability to distinguish between linear and quadratic bias. This is encouraging because the spirit of the position-dependent correlation function is to capture the density dependence of two-point clustering and hence
similar to sphere bias that relies on different (small) spheres with a potentially large average density contrast.

For non-local pNG such as the equilateral and orthogonal types, the expected signal in sphere bias is much smaller than in the local case, see Fig. B2. Hence, one has to rely on the one-point PDF for which achievable constraints on $f_{\text{NL}}$ should be of the same order than for the local case because the main qualitative change is in the scale dependence of the skewness, see Fig. B1.

## 7 CONCLUSION

Relying on analytical predictions for the one-point PDF and the two-point bias of spherically averaged cosmic densities, we disentangled dynamically generated non-Gaussianity from primordial ones. The pNG is encoded in the form of a small initial skewness (controlled by $f_{\text{NL}}$) with a scale dependence that in general depends on the bispectrum of the underlying model and is very small for local pNG.

We successfully benchmarked the analytical predictions for the one-point PDF from local non-Gaussianity in a range from $f_{\text{NL}} = -100$ to +100 against numerical simulations finding excellent agreement, achieving 1 per cent accuracy within the central density region $\rho \in [0.5, 2.5]$ and about 5 per cent for the adjacent high-/low-density tails. Similarly, we have tested the predicted impact of pNG on two-point statistics encoded in the scale dependence of the sphere bias against the simulation finding very good agreement within the error bars of the simulation. Most notably, we have observed a scale-dependent sphere bias in extreme density environments that includes overdensities – a well-known result for halo bias – but also underdensities suggesting that void bias could be used to constrain pNG as well.

Using a simple joint ML estimator for the amplitude of the nonlinear variance $\sigma_8$ and local primordial skewness $f_{\text{NL}}$, we obtained an estimate for the constraining power of the one-point density-in-spheres PDFs for a Euclid-like survey finding $\sigma(f_{\text{NL}}) \simeq 10$. We then discussed the influence of tracer bias, redshift space distortions and relativistic effects and provided an outlook for constraining pNG jointly using the one-point PDF and two-point bias while including bias. Given the clear scale-dependent bias signal, we observed for both high- and low-density spheres, one can hope to improve upon pNG constraints obtained from halo bias alone in the near future and better disentangle between $f_{\text{NL}}$ and $g_{\text{NL}}$.

## ACKNOWLEDGEMENTS

This work is partially supported by the grants ANR-12-BS05-0002 and ANR-13-BS05-0005 of the French Agence Nationale de la Recherche. CU and EP are supported by the Delta-ITP consortium, a programme of the Netherlands organization for scientific research (NWO) that is funded by the Dutch Ministry of Education, Culture and Science (OCW). CU thanks IAP for hospitality when this project was initiated. TN acknowledges financial support from JSPS KAKENHI Grant Number 17K14273 and JST CREST Grant Number JPMJCR1414. Numerical simulations presented in this paper were carried out on Cray XC30 at Center for Computational Astrophysics, National Astronomical Observatory of Japan.

CU thanks Guido D'Amico, Matteo Biagetti, Kwan Chuen Chan, Vincent Desjacques, Nico Hamaus, Eiichiro Komatsu, Mikhail Ivanov, Marcello Musso, Marilena LoVerde, Mark Neyrinck, Emiliano Sefusatti, Fabian Schmidt, Ravi Sheth, Sergey Sibiryakov, Yvette Welling, Drian van der Woude for discussions and Toyokazu Sekiguchi together with Shuichiro Yokoyama for information.

## REFERENCES

* Abazajian K. N. et al., 2016, preprint (arXiv:1610.02743)

* Abbas U., Sheth R. K., 2007, MNRAS, 378, 641

* Adshead P., Baxter E. J., Dodelson S., Lidz A., 2012, Phys. Rev. D, 86, 063526

* Afshordi N., Tolley A. J., 2008, Phys. Rev. D, 78, 123507

* Alvarez M. et al., 2014, preprint (arXiv:1412.4671)

* Amendola L. et al., 2016, preprint (arXiv:1606.00180)

* Assassi V., Baumann D., Pajer E., Welling Y., van der Woude D., 2015, J. Cosmol. Astropart. Phys., 11, 024

* Baldauf T., Seljak U., Senatore L., 2011, J. Cosmol. Astropart. Phys., 4, 006

* Baldauf T., Mirbabayi M., Simonović M., Zaldarriaga M., 2016, preprint (arXiv:1602.00674)

* Barndorff-Nielsen O., Cox D. R., 1979, J. R. Stat. Soc. B, 41, 279

* Bartolo N., Komatsu E., Matarrese S., Riotto A., 2004, Phys. Rep., 402, 103

* Baumann D., Nicolis A., Senatore L., Zaldarriaga M., 2012, J. Cosmol. Astropart. Phys., 7, 051

* Becker A., Huterer D., Kadota K., 2012, J. Cosmol. Astropart. Phys., 1212, 034

* Bernardeau F., 1994a, A&A, 291, 697

* Bernardeau F., 1994b, ApJ, 427, 51

* Bernardeau F., 1996, A&A, 312, 11

* Bernardeau F., Kofman L., 1995, ApJ, 443, 479

* Bernardeau F., Reimberg P., 2016, Phys. Rev. D, 94, 063520

* Bernardeau F., Uzan J.-P., 2003, Phys. Rev. D, 67, 121301

* Bernardeau F., Colombi S., Gaztañaga E., Scoccimarro R., 2002, Phys. Rep., 367, 1

* Bernardeau F., Crocce M., Sefusatti E., 2010, Phys. Rev. D, 82, 083507

* Bernardeau F., Pichon C., Codis S., 2014, Phys. Rev. D, 90, 103519

* Bernardeau F., Codis S., Pichon C., 2015, MNRAS, 449, L105

* Blinnikov S., Moessner R., 1998, A&AS, 130, 193

* Bruni M., Crittenden R., Koyama K., Maartens R., Pitrou C., Wands D., 2012, Phys. Rev. D, 85, 041301

* Bruni M., Hidalgo J. C., Wands D., 2014, ApJ, 794, L11

* Byun J., Bean R., 2015, J. Cosmol. Astropart. Phys., 3, 019

* Camera S., Santos M. G., Maartens R., 2015, MNRAS, 448, 1035

* Chambers J. M., 1967, Biometrika, 54, 367

* Chen X., 2010, Adv. Astron., 2010, 638979

* Chiang C.-T., Wagner C., Sánchez A. G., Schmidt F., Komatsu E., 2015, J. Cosmol. Astropart. Phys., 9, 028

* Chodorowski M. J., Bouchet F. R., 1996, MNRAS, 279, 557

* Chongchitnan S., Silk J., 2011, Phys. Rev. D, 83, 083504

* Codis S., Pichon C., Bernardeau F., Uhlemann C., Prunet S., 2016a, MNRAS, 460, 1549

* Codis S., Bernardeau F., Pichon C., 2016b, MNRAS, 460, 1598

* Colombi S., Bernardeau F., Bouchet F. R., Hernquist L., 1997, MNRAS, 287, 241

* Creminelli P., Senatore L., Zaldarriaga M., 2007, J. Cosmol. Astropart. Phys., 3, 019

* Crocce M., Pueblas S., Scoccimarro R., 2006, MNRAS, 373, 369

* D'Amico G., Musso M., Noreña J., Paranjape A., 2011a, J. Cosmol. Astropart. Phys., 2, 001

* D'Amico G., Musso M., Noreña J., Paranjape A., 2011b, Phys. Rev. D, 83, 023521

* Dalal N., Doré O., Huterer D., Shirokov A., 2008, Phys. Rev. D, 77, 123514

* Dekel A., Rees M. J., 1987, Nature, 326, 455

* Desjacques V., Seljak U., 2010a, Phys. Rev. D, 81, 023006

* Desjacques V., Seljak U., 2010b, Adv. Astron., 2010, 908640

* Desjacques V., Seljak U., Iliev I. T., 2009, MNRAS, 396, 85

* Desjacques V., Jeong D., Schmidt F., 2011a, Phys. Rev. D, 84, 061301

* Desjacques V., Jeong D., Schmidt F., 2011b, Phys. Rev. D, 84, 063512

* Desjacques V., Jeong D., Schmidt F., 2016, preprint (arXiv:1611.09787)

* Di Dio E., Perrier H., Durrer R., Marozzi G., Moradinezhad Dizgah A., Noreña J., Riotto A., 2017, J. Cosmol. Astropart. Phys., 3, 006

* Doré O. et al., 2014, preprint (arXiv:1412.4872)
Durrer R., Juszkiewicz R., Kunz M., Uzan J.-P., 2000, Phys. Rev. D, 62, 021301

Fergusson J. R., Shellard E. P. S., 2009, Phys. Rev. D, 80, 043510

Fry J. N., Scherrer R. J., 1994, ApJ, 429, 36

Gaztanaga E., Fosalba P., 1998, MNRAS, 301, 524

Gaztanaga E., Maehoenen P., 1996, ApJ, 462, L1

Giannantonio T., Percival W. J., 2014, MNRAS, 441, L16

Giannantonio T., Porciani C., 2010, Phys. Rev. D, 81, 063530

Giannantonio T., Porciani C., Carron J., Amara A., Pillepich A., 2012, MNRAS, 422, 2854

Giannantonio T., Ross A. J., Percival W. J., Crittenden R., Bacher D., Kilbinger M., Nichol R., Weller J., 2014, Phys. Rev. D, 89, 023511

Gleyzes J., de Putter R., Green D., Doré O., 2017, J. Cosmol. Astropart. Phys., 4, 002

Grinstein B., Wise M. B., 1986, ApJ, 310, 19

Grossi M., Branchini E., Dolag K., Matarrese S., Moscardini L., 2008, MNRAS, 390, 438

Grossi M., Verde L., Carbone C., Dolag K., Branchini E., Iannuzzi F., Matarrese S., Moscardini L., 2009, MNRAS, 398, 321

Hamaus N., Seljak U., Desjacques V., 2011, Phys. Rev. D, 84, 083509

Jeong D., Komatsu E., 2009, ApJ, 703, 1230

Jimenez R., Verde L., 2009, Phys. Rev. D, 80, 127302

Juszkiewicz R., Weinberg D. H., Amsterdamski P., Chodorowski M., Bouchet F., 1995, ApJ, 442, 39

Kaiser N., 1984, ApJ, 284, L9

Kaiser N., 1987, MNRAS, 227, 1

Kamionkowski M., Verde L., Jimenez R., 2009, J. Cosmol. Astropart. Phys., 1, 010

Kendall M. G., Stuart A., Ord J. K. eds, 1987, Kendall's Advanced Theory of Statistics. Oxford Univ. Press, Oxford

Lam T. Y., Sheth R. K., 2009, MNRAS, 398, 2143

Lam T. Y., Sheth R. K., Desjacques V., 2009, MNRAS, 399, 1482

Lam T. Y., Nishimichi T., Yoshida N., 2011, MNRAS, 414, 289

Laureijs R. et al., 2011, preprint (arXiv:1110.3193)

Lazanu A., Giannantonio T., Schmittfull M., Shellard E. P. S., 2016, Phys. Rev. D, 93, 083517

Lazanu A., Giannantonio T., Schmittfull M., Shellard E. P. S., 2017, Phys. Rev. D, 95, 083511

Leistedt B., Peiris H. V., Roth N., 2014, Phys. Rev. Lett., 113, 221301

Levi M. et al., 2013, preprint (arXiv:1308.0847)

LoVerde M., Smith K. M., 2011, J. Cosmol. Astropart. Phys., 8, 003

LoVerde M., Miller A., Shandera S., Verde L., 2008, J. Cosmol. Astropart. Phys., 4, 014

LSST Science Collaboration, 2009, preprint (arXiv:0912.0201)

McCullagh P., 1984, Biometrika, 71, 461

McDonald P., 2008, Phys. Rev. D, 78, 123519

Mana A., Giannantonio T., Weller J., Hoyle B., Hütsi G., Sartoris B., 2013, MNRAS, 434, 684

Mao Q., Berlind A. A., McBride C. K., Scherrer R. J., Scoccimarro R., Manera M., 2014, MNRAS, 443, 1402

Matarrese S., Verde L., 2008, ApJ, 677, L77

Matarrese S., Lucchin F., Bonometto S. A., 1986, ApJ, 310, L21

Matarrese S., Verde L., Jimenez R., 2000, ApJ, 541, 10

Matsubara T., 2011, Phys. Rev. D, 83, 083518

Moradinezhad Dizgah A., Dodelson S., Riotto A., 2013, Phys. Rev. D, 88, 063513

Mukhanov V. F., Chibisov G. V., 1981, ZhETF Pisma Redaktsiiu, 33, 549

Musso M., Paranjape A., 2012, MNRAS, 420, 369

Musso M., Sheth R. K., 2014, MNRAS, 439, 3051

Nishimichi T., 2012, J. Cosmol. Astropart. Phys., 8, 037

Pillepich A., Porciani C., Hahn O., 2010, MNRAS, 402, 191

Planck Collaboration XVII, 2016, A&A, 594, A17

Planck Collaboration XX, 2016, A&A, 594, A20

Renaux-Petel S., 2015, C. R. Phys., 16, 969

Ross A. J. et al., 2013, MNRAS, 428, 1116

Roth N., Porciani C., 2012, MNRAS, 425, L81

Sánchez C., Carrasco Kind M., Lin H., Miquel R., Abdalla F. B., Amara A., 2014, MNRAS, 445, 1482

Sartoris B., Borgani S., Fedeli C., Matarrese S., Moscardini L., Rosati P., Weller J., 2010, MNRAS, 407, 2339

Schmidt F., Kamionkowski M., 2010, Phys. Rev. D, 82, 103002

Scoccimarro R., 1998, MNRAS, 299, 1097

Scoccimarro R., 2000, ApJ, 542, 1

Scoccimarro R., Sefusatti E., Zaldarriaga M., 2004, Phys. Rev. D, 69, 103513

Scoccimarro R., Hui L., Manera M., Chan K. C., 2012, Phys. Rev. D, 85, 083002

Sefusatti E., 2009, Phys. Rev. D, 80, 123002

Sefusatti E., Komatsu E., 2007, Phys. Rev. D, 76, 083004

Sefusatti E., Liguori M., Yadav A. P. S., Jackson M. G., Pajer E., 2009, J. Cosmol. Astropart. Phys., 12, 022

Sefusatti E., Crocce M., Desjacques V., 2010, MNRAS, 406, 1014

Sefusatti E., Crocce M., Desjacques V., 2012, MNRAS, 425, 2903

Sekiguchi T., Yokoyama S., 2012, preprint (arXiv:1204.2726)

Seljak U., 2009, Phys. Rev. Lett., 102, 021302

Shandera S., Dalal N., Huterer D., 2011, J. Cosmol. Astropart. Phys., 3, 017

Sheth R. K., 1998, MNRAS, 300, 1057

Slosar A., Hirata C., Seljak U., Ho S., Padmanabhan N., 2008, J. Cosmol. Astropart. Phys., 8, 031

Smith R. E., Desjacques V., Marian L., 2011, Phys. Rev. D, 83, 043526

Smith K. M., Ferraro S., LoVerde M., 2012, J. Cosmol. Astropart. Phys., 3, 032

Song H., Lee J., 2009, ApJ, 701, L25

Springel V., 2005, MNRAS, 364, 1105

Springel V., Yoshida N., White S. D. M., 2001, New Astron., 6, 79

Takahashi T., 2014, Prog. Theor. Exp. Phys., 2014, 06B105

Taruya A., Koyama K., Matsubara T., 2008, Phys. Rev. D, 78, 123534

Taruya A., Nishimichi T., Saito S., 2010, Phys. Rev. D, 82, 063522

Tellarini C., Ross A. J., Tasinato G., Wands D., 2016, J. Cosmol. Astropart. Phys., 6, 014

Touchette H., 2009, Phys. Rep., 478, 1

Touchette H., 2011, preprint (arXiv:1106.4146)

Turok N., Spergel D. N., 1991, Phys. Rev. Lett., 66, 3093

Uhlemann C., Codis S., Pichon C., Bernardeau F., Reimberg P., 2016, MNRAS, 460, 1529

Uhlemann C. et al., 2017a, MNRAS, preprint (arXiv:1705.08901)

Uhlemann C., Codis S., Kim J., Pichon C., Bernardeau F., Pogosyan D., Park C., L'Huillier B., 2017b, MNRAS, 466, 2067

Valageas P., 2002a, A&A, 382, 412

Valageas P., 2002b, A&A, 382, 431

Valageas P., 2009, A&A, 508, 93

Valageas P., 2010, A&A, 514, A46

Valageas P., Nishimichi T., 2011, A&A, 527, A87

Verde L., 2010, Adv. Astron., 2010, 768675

Verde L., Heavens A. F., 2001, ApJ, 553, 14

Verde L., Jimenez R., Alvarez-Gaume L., Heavens A. F., Matarrese S., 2013, J. Cosmol. Astropart. Phys., 6, 023

Welling Y., van der Woude D., Pajer E., 2016, J. Cosmol. Astropart. Phys., 8, 044

White M., 1999, MNRAS, 310, 511

Yoo J., Hamaus N., Seljak U., Zaldarriaga M., 2012, Phys. Rev., D86, 063514

## APPENDIX A: OBTAINING THE DENSITY PDF FROM LARGE DEVIATION STATISTICS

Considering the one-point PDF of densities in spheres, $P(\rho)$, a highly symmetric configuration allows us to take advantage of the spherical collapse model for gravitational dynamics. To obtain the PDF, we use the cumulant generating function of densities in spherical cells, $\varphi(\lambda)$, defined via a Laplace transform of the density PDF

$$ \varphi(\lambda) = \log \left[ \int d\rho \exp(\lambda\rho) P(\rho) \right]. \quad (A1) $$
This relationship is useful because, in the limit of zero variance, the cumulant generating function is obtained analytically from the decay-rate function $\Psi(\rho)$ via a Legendre transformation

$$ \varphi(\lambda) = \lambda\rho - \Psi(\rho), \quad \lambda = \Psi'(\rho), \qquad (A2) $$

where the conjugate variable $\lambda$ is a function of the density $\rho$ via the stationary condition on the decay-rate function. This is a key result in large deviation statistics, which can heuristically be understood in terms of a steepest decent method and more formally derived from the Gärtner-Ellis theorem. The final decay-rate function can then be obtained from the initial decay-rate function by a simple remapping according to spherical collapse as described in equation (24), as a result of the so-called contraction principle. For a pedagogical introduction to large deviation statistics we refer to Touchette (2011) and Touchette (2009); Bernardeau & Reimberg (2016) for a more detailed and thorough discussion of its underlying principles and theorems. The PDF of the density is then given as an inverse Laplace transform of the cumulant generating function $\varphi(\lambda)$

$$ \mathcal{P}(\rho) = \int \frac{d\lambda}{2\pi i} \exp[-\lambda\rho + \varphi(\lambda)]. \qquad (A3) $$

Hence, the PDF can be obtained from a numerical integration in the complex plane as done in Bernardeau et al. (2014) and Bernardeau, Codis & Pichon (2015) or analytically evaluated using a saddle point approximation for the log-density (which has a close-to-optimal range of validity) as described in Uhlemann et al. (2016).

## APPENDIX B: BEYOND LEADING-ORDER AND LOCAL

This appendix generalizes the formalism presented for leading-order local pNG to general primordial bispectra and next-to-leading-order local pNG. As prime example for arbitrary primordial bispectra, it will consider equilateral and orthogonal pNG and show how the model ingredients vary compared to local pNG. For next-to-leading-order local pNG, the clear hierarchy of primordial rescaled cumulants that is relevant for the validity of the Edgeworth-expansions for the primordial statistics used in the main text will be exhibited.

### B1 Scale dependence of primordial cumulants

#### B1.1 Non-Gaussianity beyond local shape

Beyond the local non-Gaussianity, whose bispectrum is a simple product of two power spectra

$$ B_{\Phi}^{\text{loc}}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = 2 f_{\text{NL}}^{\text{loc}} [P_{\Phi}(\mathbf{k}_1) P_{\Phi}(\mathbf{k}_2) + 2 \text{perm.}], \qquad (B1) $$

there are many other shapes of the primordial bispectrum, for example of equilateral and orthogonal type

$$ B_{\Phi}^{\text{eq}}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = 6 f_{\text{NL}}^{\text{eq}} \left[ -[P_{\Phi}(\mathbf{k}_1)P_{\Phi}(\mathbf{k}_2) + 2 \text{perm.}] \right. \\ \left. -2P_{\Phi}^{2/3}(\mathbf{k}_1)P_{\Phi}^{2/3}(\mathbf{k}_2)P_{\Phi}^{2/3}(\mathbf{k}_3) \right. \\ \left. + P_{\Phi}(\mathbf{k}_1)P_{\Phi}^{2/3}(\mathbf{k}_2)P_{\Phi}^{1/3}(\mathbf{k}_3) + 5 \text{perm.} \right], \qquad (B2) $$

$$ B_{\Phi}^{\text{orth}}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = 6 f_{\text{NL}}^{\text{orth}} \left[ \begin{aligned} & -3[P_{\Phi}(\mathbf{k}_1)P_{\Phi}(\mathbf{k}_2) + 2 \text{perm.}] \\ & -8P_{\Phi}^{2/3}(\mathbf{k}_1)P_{\Phi}^{2/3}(\mathbf{k}_2)P_{\Phi}^{2/3}(\mathbf{k}_3) \\ & +3[P_{\Phi}(\mathbf{k}_1)P_{\Phi}^{2/3}(\mathbf{k}_2)P_{\Phi}^{1/3}(\mathbf{k}_3) + 5 \text{perm.}] \end{aligned} \right]. \qquad (B3) $$

Note that, in contrast to the strictly local case, for the equilateral bispectrum there is no a priori reason for $f_{\text{NL}}^{\text{eq}}$ to be scale-independent (see e.g. LoVerde et al. 2008). While single-field inflationary models generically predict the local-type, single-field models generate predominantly other forms. Local non-Gaussianity peaks for squeezed configurations $k_1 \simeq k_2 \ll k_3$, while equilateral peaks for $k_1 \simeq k_2 \simeq k_3$, note that a fair comparison of the amplitudes $f_{\text{NL}}^{\text{loc}}$ and $f_{\text{NL}}^{\text{equi}}$ requires careful normalization that takes the different shape-dependence into account, as argued in Fergusson & Shellard (2009).

Extending the leading-order mixed cumulant given in equations (9) and (12) to a general form for the primordial bispectrum, $B_{\Phi}$ yields

$$ \kappa_{3,\text{NG}} = \int \frac{d^3 k_1}{(2\pi)^3} \alpha_R(k_1) \int \frac{d^3 k_2}{(2\pi)^3} \alpha_R(k_2) \\ \times \alpha_R(|\mathbf{k}_1 + \mathbf{k}_2|) B_{\Phi}(\mathbf{k}_1, \mathbf{k}_2, -(\mathbf{k}_1 + \mathbf{k}_2)), \qquad (B4) $$

and

$$ \kappa_{12,\text{NG}}(r) = \int \frac{d^3 k_1}{(2\pi)^3} \alpha_R(k_1) \int \frac{d^3 k}{(2\pi)^3} \alpha_R(k) \\ \times \alpha_R(|\mathbf{k} - \mathbf{k}_1|) B_{\Phi}(\mathbf{k}_1, -\mathbf{k}_1 - \mathbf{k}, \mathbf{k}) \exp[i\mathbf{k} \cdot \mathbf{r}]. \qquad (B5) $$

Since the present formalism does not make any assumption on the scale-dependence of the primordial cumulants, it can be applied to any known bispectrum shape and naturally incorporates models of fixed shape with running non-Gaussianity parameter as discussed in LoVerde et al. (2008). Note that for identical scale-independent pre-factors $f_{\text{NL}}^{\text{loc}} = f_{\text{NL}}^{\text{equi}}$, the relative amplitude of the skewness differs $\kappa_3^{\text{loc}} \simeq 3\kappa_3^{\text{equi}} \simeq -5\kappa_3^{\text{orth}}$ for radii around 10–15 Mpc h$^{-1}$, while the density scale-dependence is shown in Fig. B1. The separation scale-dependence of the leading-order mixed cumulant is compared in Fig. B2, which shows that the orthogonal model has a sign opposite to the local and equilateral models. This figure shows that the signal of scale-dependent bias is largest for the local model, suppressed

rescaled skewness $\frac{\epsilon_3(R * \rho^{1/3})}{\epsilon_3(R)}$ for R=10 Mpc/h

<figure><img src="image_10.png" /><figcaption>Figure B1. Scale-dependence of the rescaled skewness for different primordial bispectra: local (flat purple line), equilateral (ascending blue line) and orthogonal (descending red line). All models are close to having a scale-invariant skewness which is best achieved for a local model.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure B2.** Scale-dependence of the ratio of the leading-order mixed cumulant $\kappa_{12,R}(r)$ and the linear correlation function $\xi_R(r)$ that enters the two-point sphere bias $b(\rho)$ evaluated at redshift $z=1$ with radius $R=15$ Mpc h$^{-1}$ as a function of separation $r$ [Mpc h$^{-1}$] normalized by the amplitude pNG given by $f_{\text{NL}}$ for different primordial bispectra: local (upper purple line), orthogonal (middle red line) and equilateral (lower blue line).</figcaption></figure>

roughly by a factor 3 for the orthogonal model, but by one to two orders of magnitude for the equilateral model.

## B1.2 Next-to-leading-order local non-Gaussianity

One can also consider cubic local pNG that can be written as

$$ \Phi_{\text{NG}} = \Phi_G + f_{\text{NL}}^{\Phi} (\Phi_G^2 - \langle \Phi_G^2 \rangle) + g_{\text{NL}}^{\Phi} (\Phi_G^3 - 3\Phi_G \langle \Phi_G^2 \rangle), \quad (\text{B6}) $$

with an extra constant $g_{\text{NL}}^{\Phi}$ compared to (1). The corresponding kurtosis is

$$ \begin{aligned} \kappa_4(R) &\equiv \langle \delta_{R,\text{NG}}^4 \rangle_c \\ &\approx \frac{48 f_{\text{NL}}^2}{(2\pi)^9} \int d^3 k_1 \alpha_R(k_1) P_\Phi(k_1) \int d^3 k_2 \alpha_R(k_2) P_\Phi(k_2) \quad (\text{B7}) \\ &\quad \times \int d^3 k_3 \alpha_R(k_3) P_\Phi(k_3) \alpha_R(k_{123}), \\ &\approx \frac{3 f_{\text{NL}}^2}{2\pi^6} \int dk_1 k_1^2 \alpha_R(k_1) P_\Phi(k_1) \int dk_2 k_2^2 \alpha_R(k_2) P_\Phi(k_2) \\ &\quad \times \int dk_3 k_3^2 \alpha_R(k_3) \int_{-1}^1 d\mu_{13} P_\Phi(k_{13}) \int_{-1}^1 d\mu_{132} \alpha_R(k_{123}), \quad (\text{B8}) \end{aligned} $$

where $k_{13} = |\mathbf{k}_1 + \mathbf{k}_3| = \sqrt{k_1^2 + k_3^2 + 2k_1 k_3 \mu_{13}}$ and $k_{123} = |\mathbf{k}_1 + \mathbf{k}_2 + \mathbf{k}_3| = k_{13}^2 + k_2^2 + 2k_{13} k_2 \mu_{123}$. Note that at next-to-leading order also the variance receives a contribution from pNG

$$ \begin{aligned} \sigma_{\text{NG}}^2(R) &\equiv \kappa_2(R) \equiv \langle \delta_{R,\text{NG}}^2 \rangle \\ &\approx \sigma_G^2(R) + \frac{f_{\text{NL}}^2}{2\pi^4} \int dk_1 k_1^2 P_\Phi(k_1) \int dk_2 k_2^2 P_\Phi(k_2) \\ &\quad \times \int_{-1}^1 d\mu \, \alpha_R \left( \sqrt{k_1^2 + k_2^2 + 2k_1 k_2 \mu} \right)^2. \quad (\text{B9}) \end{aligned} $$

Roth & Porciani (2012) remark that non-zero $f_{\text{NL}}$ and $g_{\text{NL}}$ might cancel for the scale-dependent bias case and that the assumption of a one-parameter model (setting either $f_{\text{NL}}$ or $g_{\text{NL}}$ to 0) can significantly bias the estimation of the pNG parameters when both $f_{\text{NL}}$ and $g_{\text{NL}}$ do not vanish and the typically expected ordering $g_{\text{NL}} \simeq f_{\text{NL}}^2$ does not hold.

Let us now consider the effect of the cubic term in the local transformation (1) that is parametrized by $g_{\text{NL}}$ which is usually assumed to be of second-order $g_{\text{NL}} \simeq f_{\text{NL}}^2$ but can be generated independently of $f_{\text{NL}}$ (see e.g. Bernardeau & Uzan 2003). Again, such a cubic term can be accounted for by the modifications it induces in the lowest order cumulants, see Matarrese et al. (2000).

The subleading term of the variance would in fact be modified by the addition of a $g_{\text{NL}} \Phi^3$ term but is cancelled by the extra term included in the bracket in equation (B6). The subleading skewness would get an extra contribution that is however of third order and hence will be ignored. The kurtosis is already modified to leading order, gaining the extra piece

$$ \kappa_4^{f_{\text{NL}},g_{\text{NL}}} = \kappa_4^{f_{\text{NL}}} + \frac{1}{2} \frac{g_{\text{NL}}}{f_{\text{NL}}^2} \kappa_4^{f_{\text{NL}}} \quad (\text{B10}) $$

which still has the nice feature of appearing as a renormalization of the previously calculated leading-order kurtosis. For two-point statistics, the two joint kurtosis terms $\kappa_{22}$ and $\kappa_{13}$ are needed, given by the formulas 53) and 54) in Chongchitnan & Silk (2011) and reducing to A5 and A6 in Desjacques & Seljak (2010a) for pure $g_{\text{NL}}$. The full expressions read

$$ \begin{aligned} \kappa_{31}(r) ={}& \langle \delta_{R,\text{NG}}(\mathbf{x})^3 \delta_{R,\text{NG}}(\mathbf{x} + \mathbf{r}) \rangle_c \\ &\approx \frac{6}{(2\pi)^9} \prod_i \left( \int dk_i k_i^2 \alpha(k_i) P_\Phi(k_i) \int_{-1}^1 d\mu_i \int_0^{2\pi} d\phi_i \right) \quad (\text{B11}) \\ &\quad \times \alpha(k_{123}) \exp[ir(k_1\mu_1 + k_2\mu_2 + k_3\mu_3)] \\ &\quad \times \left[ g_{\text{NL}} \left( 1 + \frac{P_\Phi(k_{123})}{P_\Phi(k_3)} \right) \right. \\ &\quad \left. + 4 f_{\text{NL}}^2 \frac{P_\Phi(k_{12})}{P_\Phi(k_2)} \left( 1 + \frac{P_\Phi(k_2)P_\Phi(k_{123})}{P_\Phi(k_1)P_\Phi(k_3)} \right) \right] \quad (\text{B12}) \end{aligned} $$

and

$$ \begin{aligned} \kappa_{22}(r) ={}& \langle \delta_{R,\text{NG}}(\mathbf{x})^2 \delta_{R,\text{NG}}(\mathbf{x} + \mathbf{r})^2 \rangle_c \\ &\approx \frac{2}{(2\pi)^9} \prod_i \left( \int dk_i k_i^2 \alpha(k_i) P_\Phi(k_i) \int_{-1}^1 d\mu_i \int_0^{2\pi} d\phi_i \right) \\ &\quad \times \alpha(k_{123}) \exp[ir(k_1\mu_1 + k_2\mu_2)] \\ &\quad \left\{ \begin{aligned} &3g_{\text{NL}} \left( 1 + 2 \frac{P_\Phi(k_{123})}{P_\Phi(k_1)} + \frac{P_\Phi(k_{123})}{P_\Phi(k_3)} \right) \\ &+ 4 f_{\text{NL}}^2 \left[ \frac{P_\Phi(k_{13})}{P_\Phi(k_3)} \left( 1 + \frac{P_\Phi(k_3)P_\Phi(k_{123})}{P_\Phi(k_1)P_\Phi(k_2)} \right) \right. \\ &\quad \left. + \frac{P_\Phi(k_{12}) + P_\Phi(k_{23})}{P_\Phi(k_2)} \left( 1 + \frac{P_\Phi(k_2)P_\Phi(k_{123})}{P_\Phi(k_1)P_\Phi(k_3)} \right) \right] \end{aligned} \right\} \quad (\text{B13}) \end{aligned} $$

with $\mu_i = \cos \angle(\mathbf{r}, \mathbf{k}_i)$ and $k_{123} = |\mathbf{k}_1 + \mathbf{k}_2 + \mathbf{k}_3| = (k_1^2 + k_2^2 + k_3^2 + 2k_1 k_2 \mu_{12} + 2k_1 k_3 \mu_{13} + 2k_2 k_3 \mu_{23})^{1/2}$, where $\mu_{i,j} = \cos \angle(\mathbf{k}_i, \mathbf{k}_j) = [(1 - \mu_i^2)(1 - \mu_j^2) \cos(\phi_i - \phi_j) + \mu_i \mu_j]^{1/2}$.

## B2 Edgeworth expansion

### B2.1 Univariate Edgeworth expansion for one-point PDF

Since pNG is small, one typically systematically expands the initial PDF of the smoothed density field in an Edgeworth expansion (Bernardeau & Kofman 1995; Juszkiewicz et al. 1995) around a Gaussian distribution. The Edgeworth series $E_n$ is an asymptotic expansion to approximate a probability distribution using its cumulants $\kappa_n$. With the Gaussian distribution as reference function it can be written as, see equation (43) in Blinnikov & Moessner (1998),
<figure><img src="image_12.png" /><figcaption>**Figure B3.** The scale dependence of the rescaled cumulants $\epsilon_n = \kappa_n/\sigma^n$ divided by the appropriate power of the amplitude of pNG $f_{NL}$ in the relevant region for different radii that are relevant for typical densities in spheres of radii $R \simeq 10-15$ Mpc $h^{-1}$.</figcaption></figure>

$$E_n(x) = \frac{1}{\sqrt{2\pi\kappa_2}} \exp\left(-\frac{(x-\kappa_1)^2}{2\kappa_2}\right) \times \left[ 1 + \sum_{s=1}^{n} \sum_{r=1}^{s} \frac{B_{s,r}(\lambda_3, \dots, \lambda_{s-r+3})}{s!} H_{s+2r}\left(\frac{x-\kappa_1}{\sqrt{\kappa_2}}\right) \right], \quad (B14a)$$

where $\lambda_n$ are the normalized and rescaled cumulants

$$\lambda_n \equiv \frac{\epsilon_n}{n(n-1)}, \qquad (B14b)$$

$B_{s,r}$ are the Bell polynomials and $H_n$ are the probabilists' Hermite polynomials. Up to second order, the relevant Bell polynomials are

$$B_{1,1}(\lambda_3) = \lambda_3, \quad B_{2,1}(\lambda_3, \lambda_4) = \lambda_4, \quad B_{2,2}(\lambda_3) = \lambda_3^2, \qquad (B14c)$$

and the higher order Hermite polynomials read

$$H_4(x) = x^4 - 6x^2 + 3, \\ H_6(x) = x^6 - 15x^4 + 45x^2 - 15. \qquad (B14d)$$

When the Edgeworth expansion is evaluated up to next-to-leading order and expressed through the reduced primordial cumulants $\tilde{S}_n$ the PDF reads

$$\mathcal{P}(\tau) = \frac{1}{\sqrt{2\pi\sigma(r)}} \exp\left(-\frac{\tau^2}{2\sigma^2(r)}\right) \left[ 1 + \frac{\epsilon_3(r)}{3!} H_3\left(\frac{\tau}{\sigma(r)}\right) + \frac{\epsilon_4(r)}{4!} H_4\left(\frac{\tau}{\sigma(r)}\right) + \frac{\epsilon_3(r)^2}{72} H_6\left(\frac{\tau}{\sigma(r)}\right) \right]. \quad (B15)$$

While the Edgeworth expansion is usually accurate in the region that is close to the peak, it leads to bigger deviations in the tails. Indeed,

at least to next-to-leading order the Edgeworth expansion can be obtained from the saddle-point approximation by expanding the exponential with the non-quadratic part of the rate function given in equation (21). Note that the modification of the variance at second order should in principle be included here by further expanding $1/\sigma^2(r) = (1 - \epsilon_2(r))/\sigma_G^2(r)$ in the Gaussian term when the next to leading order in $f_{NL}$, or the leading order $g_{NL}$, is considered, but in practice turns out to be further suppressed; see Fig. B3.

## B2.2 Bivariate Edgeworth expansion for two-point sphere bias

The extension of the bivariate Edgeworth expansion for the two-point PDF from the leading order given in equation (33) to next-to-leading order is

$$\mathcal{P}(v, \zeta) \approx \frac{\exp\left(-\frac{v^2}{2} - \frac{\zeta^2}{2}\right)}{2\pi} \left[ 1 + \frac{1}{3!} (\epsilon_{30} H_3(v) + \epsilon_{03} H_3(\zeta)) \right. \\ + \frac{1}{2} (\epsilon_{12} H_1(v) H_2(\zeta) + \epsilon_{21} H_2(v) H_1(\zeta)) \\ + \frac{1}{4!} (\epsilon_{40} H_4(v) + \epsilon_{04} H_4(\zeta)) + \frac{1}{4} \epsilon_{22} H_2(v) H_2(\zeta) \\ + \frac{1}{3!} (\epsilon_{31} H_3(v) H_1(\zeta) + \epsilon_{13} H_1(v) H_3(\zeta)) \\ + \frac{1}{72} (\epsilon_{30}^2 H_6(v) + \epsilon_{03}^2 H_6(\zeta)) + \frac{1}{4} \epsilon_{12} \epsilon_{21} H_3(v) H_3(\zeta) \\ + \frac{1}{12} (\epsilon_{30} \epsilon_{21} H_5(v) H_1(\zeta) + \epsilon_{12} \epsilon_{03} H_1(v) H_5(\zeta) \\ + \epsilon_{30} \epsilon_{12} H_4(v) H_2(\zeta) + \epsilon_{21} \epsilon_{03} H_2(v) H_4(\zeta)) \\ \left. + \frac{1}{8} (\epsilon_{21}^2 H_4(v) H_2(\zeta) + \epsilon_{12} H_2(v) H_4(\zeta)) \right],$$

using the notation $\epsilon_{ij} = \langle v^i \zeta^j \rangle_c$ for the joint cumulants.

Note that the initial sphere bias for a pure $g_{NL}$ model where $\kappa_3 = \kappa_{12} = 0$, in complement to equation (35) for pure $f_{NL}$, is given as

$$b_{\text{ini}}^{\text{NG,gNL}}(\tau, r) = \frac{\tau}{\sigma^2} \left\{ 1 + \frac{1}{6} \left[ \left(\frac{\tau}{\sigma}\right)^2 - 3 \right] \left( \frac{\kappa_{13}(r)}{\sigma^2 \xi(r)} - \frac{\kappa_4}{\sigma^4} \right) \right\}. \quad (B16)$$

Interestingly, $g_{NL}$ predicts an antisymmetric signature of the sphere bias for high and low densities, meaning that the sign of the sphere bias changes when comparing overdensities and underdensities. This is in contrast to the result for $f_{NL}$ from equation (35) which has the same sign for both high and low densities. Hence, sphere bias for both high and low densities is promising to disentangle $f_{NL}$ and $g_{NL}$, which appear degenerate when looking at halo bias alone.

This paper has been typeset from a TeX/LA TeX file prepared by the author.